Estimation of exciton diffusion lengths of organic semiconductors in random domains

Exciton diffusion length plays a vital role in the function of opto-electronic devices. Oftentimes, the domain occupied by an organic semiconductor is subject to surface measurement error. In many experiments, photoluminescence over the domain is measured and used as the observation data to estimate this length parameter in an inverse manner based on the least square method. However, the result is sometimes found to be sensitive to the surface geometry of the domain. In this paper, we employ a random function representation for the uncertain surface of the domain. After non-dimensionalization, the forward model becomes a diffusion-type equation over the domain whose geometric boundary is subject to small random perturbations. We propose an asymptotic-based method as an approximate forward solver whose accuracy is justified both theoretically and numerically. It only requires solving several deterministic problems over a fixed domain. Therefore, for the same accuracy requirements we tested here, the running time of our approach is more than one order of magnitude smaller than that of directly solving the original stochastic boundary-value problem by the stochastic collocation method. In addition, from numerical results, we find that the correlation length of randomness is important to determine whether a 1D reduced model is a good surrogate for the 2D model.


Introduction
From a practical perspective, measurement error or insufficient data in many problems inevitably introduces uncertainty, which however has been overlooked for a long time. In materials science, recent adventure in manufacturing has reduced the device dimension from macroscropic/mesoscropic scales to nanoscale, in which the uncertainty becomes important [4]. In the field of organic opto-electronics, such as organic light-emitting diodes (LEDs) and organic photovoltaics, a surge of interest has occurred over the past few decades, due to major advancements in material design, which led to a significant boost in the materials performance [28,24,31]. These materials are carbon-based compounds with other elements like N, O, H, S, and P, and can be classified into small molecules, oligomers, and polymers with atomic mass units ranging from several hundreds to at least several thousands and conjugation length ranging from a few nanometers to hundreds of nanometers [13,24].
At the electronic level, exciton, a bound electron-hole pair, is the elementary energy carrier, which does not carry net electric charge. The characteristic distance that an exciton travels during its lifetime is defined as the exciton diffusion length, which plays a critical role in the function of opto-electronical devices. A small diffusion length in organic photovoltaics limits the dissociation of excitons into free charges [33,22], while a large diffusion length in organic LEDs may limit luminous efficiency if excitons diffuse to non-radiative quenching sites [1]. Generally, there are two types of experimental methods to measure exciton diffusion length: photoluminescence quenching measurement, including steady-state and time-resolved photoluminescence surface quenching, time-resolved photoluminescence bulk quenching, and exciton-exciton annihilation [20], and photocurrent spectrum measurement [27]. Exciton generation, diffusion, dissociation, recombination, exciton-exciton annihilation, and exciton-environment interaction, are the typical underlying processes. Accordingly, two types of models are used to describe exciton diffusion, either differential equation based or stochastic process based. The connections between these models are systematically discussed in [9].
We focus on the differential equation model in this paper. Accordingly, the device used in the experiment includes two layers of organic materials. One layer of material is called donor and the other is called acceptor or quencher due to the difference of their chemical properties. A typical bilayer structure is illustrated in Figure 1. These materials are thin films with thicknesses ranging from tens of nanometers to hundreds of nanometers along the x direction and in-plane dimensions up to the macroscopic scale. Under the illumination of solar lights, excitons are generated in the donor layer, and then diffuse. Due to the exciton-environment interaction, some excitons die out and emit photons which contribute to the photoluminescence. The donor-acceptor interface serves as the absorbing boundary while other boundaries serve as reflecting boundaries due to the tailored properties of the donor and the acceptor. As derived in [9], such a problem can be modeled by a diffusion-type equation with appropriate boundary conditions, which will be introduced in §2. Since the donor-acceptor interface is not exposed to the air/vacuum and the resolution of the surface morphology is limited by the resolution of atomic force microscopy, this interface is subject to an uncertainty with amplitude around 1 nm. At a first glance, this uncertainty does not seem to affect the observation very much since its amplitude is much smaller than the film thickness. However, in some scenarios [20], the fitted exciton diffusion lengths are sensitive to the uncertainty, which may affect a chemist to determine which material should be used for a specific device. Therefore, it is desirable to understand the quantitative effect of such an uncertainty on the exciton diffusion length and provide a reliable estimation method to select appropriate models for organic materials with different crystalline orders.
Uncertainty quantification is an emerging research field that addressing these issues [35,19,30]. Due to the complex nature of the problems considered here, finding analytical solutions is almost impossible, so numerical methods are very important to study these solutions. Here we give a briefly introduction of existing numerical methods, which can be classified into nonintrusive sampling methods and intrusive methods.
Monte Carlo (MC) method is the most popular non-intrusive method [16]. For the randomness in the partial differential equations (PDEs), one first generates N random samples, and then solves the corresponding deterministic problem to obtain solution samples. Finally, one estimates the statistical information by ensemble averaging. The MC method is easy to implement, but the convergence rate is merely O( 1 √ N ). Later on, quasi-Monte Carlo methods [7] and multilevel Monte Carlo methods [15] have been developed to speed up the MC method. Stochastic collocation (SC) methods explore the smoothness of PDE solutions with respect to random variables and use certain quadrature points and weights to compute solution realizations [36,2,25]. Exponential convergence can be achieved for smooth solutions, but the quadrature points increase exponentially fast as the number of random variables increases, known as the curse of dimensionality. Sparse grids were introduced to reduce the quadrature points to some extent [6]. For high-dimensional PDEs with randomness, however, the sparse grid method is still very expensive.
In intrusive methods, solutions of the random PDEs are represented by certain basis functions, e.g., orthogonal polynomials. Typical examples are the Wiener chaos expansion (WCE) and polynomial chaos expansion (PCE) method. Then, Galerkin method is used to derive a coupled deterministic PDE system to compute the expansion coefficients. The WCE was introduced by Wiener in [34]. However, it did not receive much attention until Cameron provided the convergence analysis in [8]. In the past two decades, many efficient methods have been developed based on WCE or PCE; see [14,37,38,3,18] and references therein.
When dealing with relatively small input variability and outputs that do not express high nonlinearity, perturbation type methods are most frequently used, where the random solutions are expanded via Taylor series around their mean and truncated at a certain order [21,11]. Typically, at most second-order expansion is used because the resulting system of equations are typically complicated beyond the second order. An intrinsic limitation of the perturbation methods is that the magnitude of the uncertainties should be small. Similarly, one also chooses the operator expansion method to solve random PDEs. In the Neumann expansion method, we expand inverse of the stochastic operator in a Neumann series and truncate it at a certain order. This type of method often strongly depends on the underlying operator and is typically limited to static problems [39,35].
In this paper, we employ a diffusion-type equation with appropriate boundary conditions as the forward model and the exciton diffusion length is extracted in an inverse manner. Surface roughness is treated as a random function. After nondimensionalization, the forward model becomes a diffusion-type equation on the domain whose geometric boundary is subject to small perturbations. Therefore, we propose an asymptotic-based method as the forward solver with its accuracy justified both analytically and numerically. It only requires solving several deterministic problems over the regular domain without randomness. The efficiency of our approach is demonstrated by comparing with the SC method as the forward solver. Of experimental interest, we find that the correlation length of randomness is the key parameter to determine whether a 1D surrogate is sufficient for the forward modeling. Precisely, the larger the correlation length, the more accurate the 1D surrogate. This explains why the 1D surrogate works well for organic semiconductors with high crystalline order.
The rest of the paper is organized as follows. In §2, a diffusion-type equation is introduced as the forward model and the exciton diffusion length is extracted by solving an inverse problem. Domain mapping method and the asymptotic-based method are introduced in §3 with simulation results presented in §4. Conclusion is drawn in §5.

Model
In this section, we introduce a diffusion-type equation over the random domain as the forward model and the extraction of exciton diffusion length is done by solving an inverse problem.
Refer to Figure 1. The donnor-acceptor interface, Γ, is described by x = h(z, ω), a random field with period L: where {θ k } are i.i.d. random variables, φ k (z) = sin(2kπ z L ), and λ k > 0 are eigenvalues that control the decay speed of physical mode φ k (z). In principle, one could also add the cosine modes in the basis functions {φ k }. We here only use the sine modes for simplicity. In the experiment,h ∼ 1 nm due to the surface roughness limited by the resolution of atomic force microscopy. The thickness d varys between 10 ∼ 100 nm in a series of devices. Therefore, the dimensionless parameter characterizing the ratio between measurement uncertainty and film thickness =h/d, ranges around [0.01, 0.1]. So, it is assume that the amplitudeh d in our models. The in-plane dimensions of the donor layer are of centimeters in the experiment, but we choose L ∼ 100 nm and set up the periodic boundary condition along the z direction based on the following two reasons. First, the current work treats exciton diffusion length as a homogeneous macroscopic quantity, which is a good approximation for ordered structures. For example, small molecules are the simplest and can form crystal structures under careful fabrication conditions [12,29]. Second, the light intensity and hence the exciton generation density is a single variable function depending on x only.
Define the domain D : Here σ is the exciton diffusion length which is an unknown parameter, and the σ 2 term in (2a) describes the exciton diffusion. Exciton-environment interaction makes some excitons emit phonons and die out, which is described by the term −u in (2a). The normalized exciton generation function G is R + -valued, and is smooth on R + ∪{0}. By solving the Maxwell equation over the layered device, one can find that G(x) is a combination of exponential functions which decay away from 0 [5]. x = d is served as the reflexive boundary and homogeneous Neumann boundary condition is thus used there, while x = h(z, ω) is served as the absorbing boundary and homogeneous Dirichlet boundary condition is used in (2b). Periodic boundary condition is imposed along the z direction in (2c). It is not difficult to see that the solution u to (2) is strictly positive in D by the maximum principle. The (normalized) photoluminescence is computed by the formula If the interface Γ is random but entirely flat, i.e., h(x, ω) = ξ(ω) for some random variable ξ, then the domain is a rectangle (ξ(ω), d) × (0, L). Notice that in (2), G is a function of x only. Then, (2) actually reduces to the following 1D problem For the 1D model (4), when L → 0, the photoluminescence defined by (3) reduces to This is why the normalized factor 1/L is used in (3). Due to the simple analytical formula, the 1D model given by (4) and (5) has been widely used to fit experimental data for photoluminescence measurement [20] and photocurrent measurement [17].
Since the roughness of the interface is taken into account, problem (2) with the random interface Γ is viewed as a generalized and more realistic model. The 1D model (4) still has the uncertainty of the boundary but fails to include the spatial variety of the donor-interface interfacial layer. We are interested in identifying under which condition the 1D model can be viewed as a good surrogate for the 2D model and how this condition can be related to the property of organic semiconductors.

2.2.
Inverse problem: Extraction of exciton diffusion length. In the experiment, photoluminescence data { I i } N i=1 are measured for a series of bilayer devices with different thicknesses Here i denotes the i-th observation in the experiment with d i the thickness of the donor layer. σ is the unknown parameter, and the optimal σ is expected to reproduce the experimental in a proper sense. To achieve this, we propose the following minimization problem in the sense of mean square error We use the Newton's method to solve (6) for σ. Given σ (0) , for n = 1, 2, . . . , until convergence, we have Here α n ∈ (0, 1] is given by a line search [26]. Details are given in Appendix A.

Methods for solving the forward model
In the photoluminescence experiment, the surface roughness is very small compared to the film thickness, i.e.,h ∼ 1 nm and 10 ≤ d ≤ 100 nm. Based on this observation, we propose an asymptotic-based method for solving the diffusion-type equation over the random domain. For comparison, we first describe the domain mapping approach [38].

Domain mapping method.
To handle the random domain D , we introduce the following transformationỹ so that D becomes the unit square D s = (0, 1) × (0, 1). Under this change of variables, Eq.
(2) becomes the following PDE with random coefficients (still use y and z to represent y and z, respectively) where the spatial differentiation operator is defined for a random element ω in the probability space and g(y, z, ω) The boundary condition is The photoluminescence defined in (3) is then transformed into Remark 3.1. In 1D, changing of variable y = x−ξ d−ξ also transforms (4) to a differential equation with random coefficients over the unit interval.
with (14) and the boundary condition Accordingly, the photoluminescence can be written as Remark 3.2. The generation term in (10) depends on both y and z after changing of variables. We expect some dimensional effect on the estimation of σ, which will be carefully examined in §4.

3.2.
Finite difference method for the model problem. We use finite difference method to discretize the forward model (8) .., N y + 1 and j = 1, ..., N z + 1, respectively. For the discretization in space, we use a second-order, centereddifference scheme [23]. We introduce the difference operators The operators D z 0 , D z − , and D z + are defined similarly. For each ω ∈ Ω and each interior mesh point (i, j) with 2 i N y , 2 j N z , we discretize the forward model (8) as where h, h , and h are evaluated at (y i , z j ). We then discretize the boundary conditions (11) on ∂D s . The Dirichlet boundary condition on y = 0 gives u 1,j = 0, 1 j N z + 1. For the Neumann boundary condition on y = 1, we introduce ghost nodes at (y −1 , z j ) and obtain a second order accurate finite difference approximation = 0. Then, the values of the u −1,j at the ghosts nodes are eliminated by combining with Eq. (17). Finally, the periodic boundary condition along the z direction gives u i,Nz+1 = u i,1 . We solve a system of N y (N z + 1) linear equations for {u i,j } with 2 i N y + 1 and 1 j N z + 1.
The equations have a regular structure, each equation involving at most nine unknowns. Thus the corresponding matrix of the system is sparse and can be solved efficiently using existing numerical solvers. After obtaining {u i,j }, we use the 2D trapezoidal quadrature rule to compute the photoluminescence I(σ, d) defined in (12).
In this paper, we choose the sparse-grid based SC method [6,25] to discretize the stochastic dimension in Eq. (8). As such the expectation of u(y, z, ω) is computed by where s q are sparse-grid quadrature points, w q are the corresponding weights, and Q is the number of sparse-grid points. Other functionals of u(y, z, ω) can be computed in the same way. When the solution u(y, z, ω) is smooth in the stochastic dimension, the SC method provides very accurate results.

3.
3. An asymptotic-based method. If we rewrite Eq. (2) in the nondimensionalized form with the change of variablesx = x/d andz = z/L, the domain D becomes where K is the mode number in the interface modeling. As discussed in §2, ∼ 0.01 − 0.1. Therefore, it is meaningful to derive the asymptotic equations when → 0. For ease of description, we list the main results below. The main idea is: (1) we rewrite Eq. (2) over D s, ; (2) with appropriate extension/restriction of solutions on the fixed domain D s , we obtain a Taylor series with each term satisfying a PDE of the same type with the boundary condition involving lower order terms; (3) we apply the inverse transform for each term and change the domain D s back to D 0 = (0, d) × (0, L). Detailed derivation can be found in Appendix B for self-consistency.
The interested readers can find the systematic study on asymptotic expansions for more general problems in [10]. The asymptotic expansion over the fixed domain D 0 is of the form The equation for each w n can be derived in a sequential manner. Only the first three terms are listed here. More details are included in Appendix B.
The leading term w 0 (x, z) is the solution to the boundary value problem and w 2 (x, z, ω) is the solution to the following boundary value problem Remark 3.3. As demonstrated in Eqs. (20), (21), and (22), the asymptotic expansion in (19) requires a sequential construction from lower order terms to high order terms and the partial derivatives of lower terms appear in the boundary condition for high terms. Numerically, we use the second-order finite difference scheme for (20), (21), and (22). For boundary conditions, we use the one-sided beam warming scheme to discretize ∂ x w 0 (0, z) and ∂ x w 1 (0, z, ω) so the overall numerical schemes are still of second order accuracy.
Define v [n] = n k=0 k w k . Note that w 0 is a function of (x, z) only. The zeroth order approximation of the photoluminescence is and so For k = 1, 2, . . . , K, let w 1,k (x, z) be the solution to (21) with φ k (z) in place ofh(z, ω) Then by linearity, the solution w 1 to (21) withh given by (18) can be expressed as Hence the first order approximation of the photoluminescence becomes and so Next, we consider the second order approximation of the photoluminescence. Sinceh and w 1 are given by (18) and (26), the boundary condition for w 2 at x = 0 can be written as Introduce w 2,j,k (x, z) as the solution to the boundary value problem (22) with the boundary condition at x = 0 replaced by Then and consequently, the second order approximation of the photoluminescence is and we have In general, w n can be written as the sum of K n functions, each of which solves a deterministic problem.
The approximation accuracy of a finite series in (19) is given by the following theorem. Proof can be found in [10].
where u is the solution to (2) and v [n] = n k=0 k w k .
To proceed, let us recall the definition of Bochner spaces.

Proof. From Theorem 3.4, for
Since w n , n 0 satisfies the same elliptic equation (20) with a boundary condition depending on w k , k n − 1. By the Lax-Milgram's theorem, we have w n ∈ L 2 P (Ω, H 1 (D 0 )), n 0. Therefore, v [n] ∈ L 2 P (Ω, H 1 (D 0 )), n 0 and the desired result is obtained. A direct consequence of Proposition 3.6 is Based on the above assertions, we have Corollary 3.7. For (24), (28), and (30), we have the following approximation errors In summary, by using the asymptotic expansion solution, we circumvent the difficulty of sampling the random function and solving PDEs on irregular domains for each sample. In our approach, there is no statistical error or errors from numerical quadratures as in MC method, SC method, and PCE method. However, our method is applicable only for small perturbation of the random interface, where a small n is sufficient in practice. The computational cost depends on the approximation order n and the number of modes K used to represent the random interface, and increases proportionally to K n .

Numerical Results
In this section, we numerically investigate the accuracy and efficiency of the asymptotic-based method in computing photoluminescence and the efficiency in estimating the exciton diffusion length. In addition, we study of the validation of the diffusion-type model, i.e., under which condition the 1D model can be viewed as a good surrogate for the 2D model.  , ω), d), z ∈ (0, L)}. Recall that the random interface h(z, ω) between the donor and the acceptor is parameterized by h(z, ω) = h K k=1 λ k θ k (ω) sin(2kπ z L ), where θ k (ω) are i.i.d. uniform random variables and K is the number of random variables in the model.
We first solve (8) over the fixed domain D s = (0, 1) × (0, 1) in the domain mapping method using the SC method. Note that the spatial differentiation operator in (9) depends on the random variables in a highly nonlinear fashion, which makes the WCE method and PCE method extremely difficult. In the asymptotic-based method, we solve deterministic boundary value problems (20), (21), and (22) over the fixed domain D 0 = (0, d) × (0, L), respectively. Recall that in the asymptotic-based method, =h/d and the random interface becomesh(z, ω) = K k=1 λ k θ k (ω) sin(2kπz). In our simulation, the random interface h(z, ω) is parameterized by K = 5 random variables. The accuracy of the asymptotic-based method is verified by two numerical tests. In the first test, θ k ∼ U (0, 1), while in the second one θ k ∼ U (−1, 1).
To compute the reference solution, we employ the finite difference method to discretize the spatial dimension of Eq. (8) with a mesh size H = 1 128 , and use the sparse-grid based SC method to discretize the stochastic dimension. We choose level six sparse grids with 903 quadrature points. After obtaining solutions at all quadrature points, we compute the expectation of the photoluminescence, which provides a very accurate reference solution. In the asymptotic-based method, we use the finite difference method to discretize the spatial dimension of boundary value problems (20), (25), and (22) for w 2,j,k with a mesh size H = 1 64 . Expectations E[θ k ] in (28) and E[θ j θ k ] in (30) can be easily computed beforehand. Therefore, given the approximate solutions w 0 , w 1,k , and w 2,j,k , we immediately obtain different order approximations of the expectation of the photoluminescence. This provides the significant computational saving over the SC method. For = 2 −i , i = 2, ..., 7, Figure 2 shows the approximation accuracy of the asymptotic-based method. In Figure 2(a), θ k ∼ U (0, 1). The approximated expectation of the photoluminescence obtained by using the zeroth, first and second order approximations are shown in the lines with circle, star, and triangle, with convergence rates 1.21, 1.99, and 3.81, respectively. In Figure  2 We conclude this subsection with a discussion on the computational time of our method. In these two tests, on average it takes 164.5 seconds to compute one reference expectation of the photoluminescence. If we choose a low level SC method to compute the expectation of the photoluminescence, it takes 27.3 seconds to compute one reference expectation of the photoluminescence that gives a comparable approximation result to our asymptotic-based method. However, our method with the second order approximation only takes 1.56 seconds to obtain one result. We achieve a 18X speedup over the SC method. Generally, the ratio of the speedup is problem-dependent. It is expected that higher ratio of speedup can be achieved it one solves a problem where the random interface is parameterized by high-dimensional random variables.

4.2.
Estimation of the exciton diffusion length. In this section, we estimate the exciton diffusion length in an inverse manner with the asymptotic-based method as the forward solver. Since only limited photoluminescence data from experiments are available, we solve the forward model (2) to generate data in our numerical tests. Specifically, given the exciton diffusion length σ, the exciton generation function G, the in-plane dimension L, and the parametrization of the random interface h(z, ω), we solve Eq. (2) for a series of thicknesses {d i }, and calculate the corresponding expectations of the photoluminescence data p{Ĩ i } according to Eq. (3). Therefore, {d i ,Ĩ i } serves as the "experimental" data. We then solve the minimization problem (6) based on our numerically generated data {d i ,Ĩ i } to estimate the "exact" exciton diffusion length σ in the presence of randomness, denoted by σ exact and will be used for comparison later.
We fix L = 4 in all our numerical tests since it is found that this minimizer is not sensitive to the in-plane dimension L. We show the convergence history of exciton diffusion lengths for various in Figure 3, where the photoluminescence data are generated with σ = 5, σ = 10, and σ = 20, respectively. Here the relative error is defined as E n, = | σexact−σ n, σexact |, where n is the iteration number, σ exact is the "exact" exciton diffusion length, and σ n, is the numerical result defined in Eq. (7). To show more details about the accuracy of our asymptotic-based method, in Tables 1, 2, and 3, we list the relative errors of our method for plotting Figures 3(a), 3(b), and 3(c). In all numerical tests, we choose the same termination criteria |σ (n) − σ (n−1) | < 10 −4 in the Newton's method. Our asymptotic-based method performs well in estimating the exciton diffusion length. In general, the smaller amplitudes the random interface, the more accurate the exciton diffusion length and the smaller the iteration number. Additionally, for larger exciton diffusion lengths σ exact , a faster convergence in the optimization approach is observed.   Table 3. Relative errors E n, = | σexact−σ n, σexact | for iteration number n = 1, 2, 3, ..., and various . The prescribed σ is 20.
Again, only limited photoluminescence data from experiments are available and we have to solve the forward model to generate data in our numerical tests. Specifically, given the exciton diffusion length σ, the exciton generation function G, and the parametrization of the random interface h(ω), we solve Eq. (4) for a series of thicknesses {d i }, and calculate the corresponding expectations of the photoluminescence data {Ĩ i } according to Eq. (5). Therefore, {d i ,Ĩ i } serves as the "experimental" data generated by the 1D model. We then solve the minimization problem (6) based on our numerically generated data {d i ,Ĩ i } to estimate the "exact" exciton diffusion length σ in the presence of randomness, denoted by σ exact and will be used for comparison.
In our numerical tests, we use the 1D model (4) with σ = 5 and σ = 10 to generate photoluminescence data. d i = 10i, i = 1, ..., 10,h = 1, and =h/d i . We use K = 10 random variables to parameterize the random interface. We set λ k = k β , where β 0 controls the decay rate of λ k . The random interface therefore takes the form 1]. Figure 4 plots the covariance function of the random interface defined by Eq. (36) for β = 0 and β = −2. It is clear that the smaller the β, the larger the correlation length.  The convergence history of the exciton diffusion length for various β is plotted in Figure 5, where the photoluminescence data is generated by the 1D model (Eqs. (4) and (5)) with σ = 5 and σ = 10. Again, the relative error is defined as E n,β = | σexact−σ n,β σexact |, where n is the iteration number, σ exact is the "exact" exciton diffusion length, and σ n,β is the numerical result defined in Eq. (7). Note that σ n,β depends also on implicitly but we omit its dependence for convenience. Tables 4 and 5 list the relative errors of our method for plotting Figure 5. The same criteria |σ (n) − σ (n−1) | < 10 −4 is used here. The numerical exciton diffusion length obtained by our method converges to the reference one with the relative error less than 1% when β −1.
Our numerical results show that a faster decay of the eigenvalues λ k leads to a better agreement between the results of the 1D model and the 2D model. The smaller the β, the better the agreement. On the other hand, the smaller the β, the larger the correlation length. Therefore, the larger the correlation length, the better the agreement. Our study sheds some light on how to select a model as simple as possible without loss of accuracy for describing exciton diffusion in organic materials. In the chemistry community, it is known that under careful fabrication conditions [12,29], organic semiconductors, including small molecules and polymers, can form crystal structures, which have large correlation lengths. As a consequence, exciton diffusion in these materials can be well described by the 1D model [27,20,32]. For organic materials with  Table 4. Relative errors E n,β = | σexact−σ n,β σexact | for iteration number n = 1, 2, 3, ..., and various β. The prescribed σ is 5. low crystalline order, i.e., small correlation length, however, our result suggests that the 1D model is not a good surrogate of the high dimensional models.

Conclusion
In this paper, we model the exciton diffusion by a diffusion-type equation with appropriate boundary conditions over a random domain. The exciton diffusion length is extracted via minimizing the mean square error between the experimental data and the model-generated data. Since the measurement uncertainty for the domain boundary is much smaller compared to the  Table 5. Relative errors E n, = | σexact−σ n,β σexact | for iteration number n = 1, 2, 3, ..., and various β. The prescribed σ is 10. device thickness, we propose an asymptotic-based method as the forward solver. Its accuracy is justified both analytically and numerically and its efficiency is demonstrated by comparing with the SC method as the forward solver. Moreover, we find that the correlation length of randomness is the key parameter to determine whether a 1D surrogate is sufficient for the forward modeling.
The discussion here focuses on the photoluminescence experiment. For the photocurrent experiment, from the modeling perspective, the forward model is the same but the objective function is different. An exciton either contributes to the photoluminescence or the photocurrent, so the photocurrent is defined as the difference between a constant (total exciton contribution) and the photoluminescence [9]. Therefore, the proposed method can be applied straightforwardly with very little modification.
, (37) where α k ∈ (0, 1] is given by the line search technique. For example, we take (6) as the minimization problem and the domain mapping formulation in §3.1 as the forward problem. Other combinations can be worked out similarly. In 2D, for the first derivatives, we have and Denote the derivatives of u(y, z) with respective to the parameter σ by u 1 (y, z) := ∂u ∂σ (y, z), and u 2 (y, z) := ∂ 2 u ∂σ 2 (y, z). Differentiating (8) with respect to σ directly, we have and u 1 shares the same boundary condition as u.
For the second derivatives, we have and and u 2 satisfies Again, the same boundary condition applies for u 2 .
Note that Eq. (42) is in fact a Cauchy problem of the time evolution equation not a boundaryvalue problem of the elliptic PDE. The velocity is specified on the interface Γ by ∂ n u and the wave travels along the normal n. So, the solution of (42) exists for 0 x h (z, ω) [10].
Particularly, we have the existence of the value of u lay at x = 0.
B.1. The solution on regular domain and its asymptotic expansion. Now the solutions u and u lay are both well-defined on D s, and L by (41) and (42), respectively. In the next, we introduce a function piecewisely defined by these two functions on the regular domain D and want to find the correct equation for this function on D in order to carry our asymptotic method. Let w be defined on D = D s, ∪ L as follows This definition is justified by (42) and immediately implies the following obvious but important fact which arises from the boundary condition on the interface Γ of u : w (x, z) = 0 on Γ .
It is easy to see that w is the unique solution to the following problem where u lay at x = 0 is given a prior : w (0, z) = u lay (0, z), for 0 z 1, ∂ x w (1, z) = 0, for 0 z 1, w (x, z + 1) = w (x, z), for (x, z) ∈ D. (45) We start with the following ansätz for w , w (x, z) = ∞ n=0 n w n (x, z) for (x, z) ∈ D.
Plug this ansätz into the equation (45), and match the terms at the same order of , then we obtain the following equations for w n in D: L w 0 + g(x) = 0, Lw n = 0, n 1.
Next, we discuss the boundary conditions for these PDEs. The two of the boundary conditions in (45), ∂ x w (1, z) = 0 and w (x, z + 1) = w (x, z), do not depend on u lay . Thus, the ansätz (46) simply gives us the same boundary conditions for each w n : ∂ x w n (1, z) = 0, and w n (x, z + 1) = w n (x, z).
The boundary condition of (45) at x = 0, i.e., on Γ 0 , depends on the data u lay on this boundary. If one works on this boundary condition, it is possible to solve the Cauchy problem (42) for small analytically so that u lay (x = 0, z) can be obtained in terms of u (i.e., w ), and eventually certain connections for w n can be built. But the use of the very original boundary condition (44) on Γ ⊂ ∂D s, , not on ∂D, actually significantly simplifies the calculations and finally offers more friendly results. The details follow below.