A Fractional Drift Diffusion Model for Organic Semiconductor Devices

: Because charge carriers of many organic semiconductors (OSCs) exhibit fractional drift diffusion (Fr-DD) transport properties, the need to develop a Fr-DD model solver becomes more apparent. However, the current research on solving the governing equations of the Fr-DD model is practically nonexistent. In this paper, an iterative solver with high precision is developed to solve both the transient and steady-state Fr-DD model for organic semiconductor devices. The Fr-DD model is composed of two fractional-order carriers (i.e., electrons and holes) continuity equations coupled with Poisson’s equation. By treating the current density as constants within each pair of consecutive grid nodes, a linear Caputo’s fractional-order ordinary differential equation (FrODE) can be produced, and its analytic solution gives an approximation to the carrier concentration. The convergence of the solver is guaranteed by implementing a successive over-relaxation (SOR) mechanism on each loop of Gummel’s iteration. Based on our derivations, it can be shown that the Scharfetter–Gummel discretization method is essentially a special case of our scheme. In addition, the consistency and convergence of the two core algorithms are proved, with three numerical examples designed to demonstrate the accuracy and computational performance of this solver. Finally, we validate the Fr-DD model for a steady-state organic field effect transistor (OFET) by fitting the simulated transconductance and output curves to the experimental data. This work aimed to develop a Fr-DD model solver for simulating the anomalous dynamics of OSC devices. Two algorithms based on a novel discretization scheme and successively over-relaxed Gummel’s iteration are proposed here to solve the transient and steady-state Fr-DD model equations. This study has identified the consistency of the two algorithms by showing that the truncation error from the discretized divergence of current density functions will vanish with the spatial step size to a positive fractional order of 1 − β . The convergence analysis reveals that the Gummel mapping is a contraction mapping if we consider the successive over-relaxation mechanism in the Gummel’s iteration, which thus completes the proof of convergence. Three numerical examples, including one benchmark example and two others constructed from the perspective of engineering applications are employed to demonstrate the algorithms’ accuracy and computational performance. It is found in the first example that altering α and β can impact the spatial convergence order but only varying α will affect the convergence order in time. The increase rate of CPU time is less than the shrinking rate of temporal step size and lower than the growth rate of spatial grid points. These findings suggest that our solver has high precision and fast computational speed as it limits the computational error to a predefined satisfactory level (from ∼ 10 − 4 to ∼ 10 − 6 ) at a relatively small expense of CPU time (from ∼ 20 s to ∼ 100 s). The results reported in two numerical examples reveal the prediction and characterization of the transient-state and steady-state dynamics for any type of organic semiconductor device. Finally, we provide experimental verification for the fractional drift diffusion model of a DNTT-OFET. We stipulate that this is the first to date exploration of the Fr-DD model solver laying the groundwork for future research into fractional drift diffusion modeling of flexible organic electronics.

in various forms. In the Van Roosbroeck representation of the DD model, the current density equation can be augmented by Einstein's relation, which gives a fixed proportional relationship between the diffusion coefficients D p , D n and the drift mobilities μ p , μ n [7]. In this paper, the Van Roosbroeck representation of the DD model can be expressed in a closed-form as Eqs. (1)- (3), where current densities are given by J n = −qμ n n∇ϕ + qD n ∇n and J p = −qμ p p∇ϕ − qD p ∇p, T is the thermal temperature, q is the charge of an electron, ε is semiconductor's absolute dielectric permittivity, N + D and N − A are ionized donor and acceptor concentrations. G n and G p are the net electron and hole generation rates, respectively. Previous research data collected from Silicon/Germanium test experiments have revealed the effectiveness of the DD model for modeling the charge carrier transports in ISCs [8]. In the past several decades, plentiful numerical algorithms have been developed for solving Eqs. (1)-(3), including the finite element method [9], finite difference fractional step method [10], mixed finite volume and modified upwind fractional difference method [11], and monotone iterative method based on the adaptive finite element discretization [12,13], etc. All of those numerical methods have one thing in common: an efficient iterative method, e.g., Newton's iteration, Gauss-Seidel iteration, or Gummel's iteration was utilized to decouple Eqs. (1)- (3). Among these iteration methods, Gummel's approach is generally more effective than other methods due to its flexibility in finding its initial guess and customizing the update formulas to improve the convergence speed and computational performance. Moreover, the effectiveness, stability and convergence of Gummel's decoupling method and iteration for its application to DD simulations were also thoroughly and rigorously proved by mathematicians [14][15][16][17]. Recent research revealed that the conventional (integer-order) DD model may not be able to characterize the charge carrier transports in organic semiconductors (OSCs), evident from the long-tail behavior of the photocurrent curve observed in OSCs [18]. Based on the DD model, Reference [19] showed that the mean squared displacement (MSD) of the carrier trajectory should be proportional to its diffusion time, i.e., E x 2 (t) ∝ t. However, the long-tail behavior of the photocurrent curve observed in OSCs implies that the MSD in this scenario is given by E x 2 (t) ∝ t α , for α termed as the dispersive parameter of the OSC, 0 < α < 1, depending on the temperature and band structure disorders [20][21][22][23]. This long-tail photocurrent phenomenon was first observed by using time-of-flight measurements [24], and the mechanism that underpins the dispersive carrier transports can be precisely explained by the "multiple trapping model" [23,25], the "single trapping model" [26] and the "hopping model" [27][28][29], respectively. Based on the "multiple trapping model," the charge carriers in OSCs can be classified as free (delocalized) charge carriers p f , n f and trap (localized) charge carriers p t , n t . The free charge carrier is the carrier that can hop freely between two trap centers and the trap charge carrier is the carrier that is permanently captured by a localized trap center. References [18,30] proved that the free hole density and the trap hole density in the p-type OSCs have a relationship as given in Eq. (4), where RL 0 D α t is Riemann-Liouville (RL) fractional derivative of order 0 < α < 1, τ 0 is the mean time of delocalization (mean free time for a free charge carrier moving between two entrapments), c is the charge carrier capture coefficient defined as c = ω 0 [sin (πα) /π α] α , ω 0 is the capture rate of the trap charge carriers, and α = kT/E 0 is the dispersive parameter which depends on the temperature T and the expected value of the exponential density of trap states E 0 . The 1D continuity equation for free charge carriers in p-type OSCs was also derived by References [18,30] as given in Eq. (5), where E (x, t) = − ∂ϕ(x, t) ∂x is the intensity of the electric field in the 1D domain, and p (x, 0) δ(t) is the initial charge carriers agitated by impacting of photon beams. Consider that p = p f + p t and p t p f in OSCs. Substituting Eq. (4) into Eq. (5) can derive the continuity equation for total charge carrier density, as given in Eq. (6), where F α (x, t) = τ 0 c α μ p E (x, t) is the anomalous advection coefficient and D α = τ 0 c α D p is the anomalous diffusion coefficient. The hole mobility μ p and hole diffusion coefficient D p satisfies the Einstein relation. Eq. (6) coupled with the 1D Poisson equation forms the 1D fractional drift diffusion (Fr-DD) model. 1 A discretization scheme, which discretizes the time-fractional derivative with backward finite difference method and the integer-order spatial derivatives with finite center difference method, was proposed to solve the 1D Fr-DD model [20][21][22][23][24][25][26][27][28][29][30][31]. Reference [20] showed that the photocurrent curves obtained from the 1D Fr-DD model are in good agreement with the recorded transient photocurrents from regio-random OSCs poly(3-hexylthiophene) (RRa-P3HT) and regio-regular poly(3-hexylthiophene) (RR-P3HT). In addition, Reference [32] introduced the fractional reduced differential transform method to solve the 1D Fr-DD model and also suggested the existence of a more general Fr-DD model with both time derivative and spatial derivatives fractionalized. As the Fr-DD model emerges as a useful tool for understanding the dispersive transport behavior of the charge carriers in OSCs, investigating how to solve it is instrumental for predicting the steady-state and transient electrical responses of OSC devices. Up to now, far too little attention has been paid to the development of a general Fr-DD model solver. Although a certain number of researches have been carried out on developing the solvers for the conventional DD model, the resulting solvers often have low accuracy and high computational complexity. The goal of this research is to develop a solver with high precision and computational performance for the Fr-DD model. The Fr-DD model is described by a group of coupled fractional-order PDEs as presented in Eqs. (7)- (9), where C 0 D α t is the Caputo's time-fractional derivative of order 0 < α ≤ 1, the fractional-order electron current density I and hole current density J are given by I = −qμ n n∇ϕ + qD n C ∇ β r n and is the Caputo's fractional gradient operator of order 0 < β ≤ 1 in 3D. Since the OSCs are typically treated as intrinsic materials without dopants, the ionized donors N + D and acceptors N − A can be omitted in Poisson's equation [33]. If the fractional-order gradient operator is defined with Riemann-Liouville's (RL) fractional derivatives, the divergence of current density function ∇ · J in Eq. (9) can be expanded [34,35]. As that defining the fractional derivatives using RL's definition may not be convenient for us to specify and assign physical meanings to the initial values and boundary conditions, we adopt Caputo's (C) fractional derivatives to replace RL's in Eqs. (7)- (9). Since composing Caputo's fractional derivative with an integer-order gradient operator from its left side may not result in the same expansion as defined by RL's derivative (See Lemma 2.2), it requires us to show that the truncation error from the discretized divergence term ∇ ·J will vanish as spatial step size reduces to zero to guarantee the consistency of the discretized divergence terms under two different definitions of fractional derivatives (RL's and Caputo's). In this paper, we propose Theorem 4.1 to illustrate the consistency of the discretized divergence term ∇ ·J under RL's and Caputo's fractional definitions.
Here, we set up a general-form Fr-DD model to simulate the anomalous transport behavior of charge carriers in OSCs. Equipped with proper initial values and boundary conditions, the Fr-DD model can handle the transient or steady-state dynamics of any-type OSC device. In addition, we develop an iterative solver for the Fr-DD model based on two novel algorithms and propose Theorem 4.2 to show the convergence of the model solver. It can be shown that the discretized DD model equation via our discretization scheme coincides with the discrete-form Fr-DD equations derived from the Scharfetter-Gummel (SG) discretization method [9], which implies that our Fr-DD model solver has wider applicability than the DD model solver based on SG method. Finally, we design three numerical examples to demonstrate the high accuracy and computational performance of the Fr-DD model solver, and experimentally validate the Fr-DD model for a steady-state OFET.
The remainder of the paper is organized as follows. Section 2 presents preliminaries in fractional calculus. Section 3 develops the solver in detail. Section 4 discusses the consistency and convergence analysis of the algorithms. Three numerical examples are provided in Section 5 to support our theoretical analysis and demonstrate the computational performance of our method. In Section 6, we adjust the parameters in the Fr-DD model to fit the experimental characteristic curves measured from a DNTT-based OFET [36]. In Section 7, we show the conclusions of this work.
Both RL and Caputo's fractional derivatives can be considered as interpolation to integerorder derivatives, which means RL The relationship between RL can Caputo's fractional derivatives can be expressed as the following lemma. ([a, t]) and n − 1 < γ ≤ n, then the following equality holds Proof. See [35].
By directly employing the definitions, the composition rules for fractional derivatives can be given as the following lemma.
Proof. Eqs. (13) and (14) can be inferred from RL and Caputo's definitions, and the proof for Eq. (15) is given in [35]. From Lemma 2.1, we have which completes the proof for Eq. (16).
It can be observed that both RL's and Caputo's fractional derivatives can be composed with an integer-order derivative from two sides, but the composition is not commutative. Next, let us give the Laplace transformation on RL and Caputo's fractional derivatives as the following lemma.
One important formula relating the Laplace transform and two-parameter Mittag-Leffler function is given in Eq. (19), and its proof can be found in [37].
Subsequently, we will present the analytic solution for Caputo's fractional linear time-invariant (LTI) state equation.

Analytic Solution of Caputo's Linear Fractional-Order ODE
If we let 0 < γ ≤ 1, the analytic solution of Caputo's linear fractional-order ODE is given in the following theorem.

Theorem 2.4 Consider the Caputo's linear fractional-order ODE defined in a discrete 1D space
Then, its solution is given by where (x) = E γ (Ax γ ) is the generalized state transition function, E γ (t) is the one-parameter Mittag-Leffler function, and the fictitious input functionv(y) is obtained byv (x) = L −1 V (s) s 1−γ .
Proof. Apply Laplace transform on both sides of Eq. (20), it gives Rearrange both sides of Eq. (22) and take the inverse Laplace transform, we have where Id = 1 or Id is the identity matrix in case u, v are vectorized variables. The last step comes from Eq. (19), i.e., the inverse Laplace transform of the Mittag-Leffler function, by letting α = γ and β = 1. We also apply the convolution theorem in the last step derivation, and this completes our proof.
The theorem we presented above establishes the precise relationship between states on two consecutive grid points in a 1D discrete space Let us assume that the input function v (t) = 1, the fictitious input function is then evaluated , and considering the commutative property of the convolution integral, then Eq. (23) in this special case can be reformulated as We notice that the second term on the RHS of Eq. (24) involves a fractional integral of order γ . The fractional integral (or the left Riemann-Liouville integral) of order γ is defined as Eq. (25). [38] In the next section, we present discrete approximation formula for the left Riemann-Liouville integral and fractional derivatives.

Discrete Approximation of Fractional Integrals and Fractional Derivatives
The fractional integral of the generalized state transition function cannot be evaluated through an analytic formula. The following lemma gives the composite Simpson's rule for approximating a left Riemann-Liouville integral.
x]) then the following (3 + γ )-th order approximation of the left Riemann-Liouville integral can be obtained where x j = j x with a positive integer j and step size x. Since we let x = x 2n , the coefficients b Proof. See [39].
For the transient-state Fr-DD model, the discretization of the time-fractional derivative is necessary, the following lemma gives a first-order approximation for Caputo's fractional time derivative of order 0 < γ ≤ 1.

]) then the following first-order approximation of the Caputo's time-fractional derivative can be obtained
where t j = j t with a positive integer j and step size t. Since we let T = t k+1 , the coefficients b Proof. The discrete approximation in Eq. (29) can be constructed by applying the piecewise quadrature, the error estimates can be derived as follows, where the mean value theorem is applied for f (t) and f (t) with ξ 1 , ξ 2 ∈ (t m , t m+1 ), and the continuous function f (2) on a compact domain [0, T] assumes its maximum value, we let M = max t∈[0, T] f (2) (t).

Derivation of the Computational Scheme
Without loss of the generality, we implement the discretization schemes in the twodimensional spatial domain, the equations and algorithms derived afterward can be extended to one-dimensional and three-dimensional scenarios. Let the spatial step size in the x direction be x = L x /(N x + 1) and in the y direction be y = L y /(N y + 1). The two-dimensional grids are given by x i = i x, i = 0, 1, . . . , N x + 1 and y j = j y, j = 0, 1, . . . , N y + 1. In the temporal domain, we let the temporal step size be t = T/N, and the temporal girds are then determined by t n = n t, n = 0, 1, 2, . . . , N. For the sake of a clear demonstration, a set of notations for the charge carrier concentrations and current densities in the 2D discrete domain are presented in Fig. 1. As shown in Fig. 1a, the electron concentration n(t k , x i , y j ) and the two components of the electron current density vector IX t k , x i+1/2 , y j , IY (t k , x i , y j+1/2 ) in the discretized electron continuity equation are denoted by n k i, j and IX k i+1/2, j , IY k i, j+1/2 , respectively. Also, Fig. 1b records a similar set of notations employed in the discretized hole continuity equation for the hole concentration p and two components of the hole current density vector (JX , JY ). In contrast to the charge carrier concentration, it is required in the scheme that the current density should be specified only in the semi-grids (the midpoints of two horizontally or vertically consecutive normal grids). In the discretized Poisson equation, the electrostatic potentials ϕ k i, j are specified in the normal grids. However, in case the underpinned material is inhomogeneous, the dielectric permittivity is not a constant within the material domain and its spatially distributed values ε i+1/2, j+1/2 should be assigned in the semi-grids (the midpoints of two diagonally consecutive normal grids) in our discretization scheme. The notations for the discrete electrostatic potentials and the spatially distributed dielectric permittivity values are not depicted in Fig. 1 due to the space limitation for demonstration. In the following sections, we derive the computational schemes for the transient and steady-state Fr-DD model separately.

Discretization of Fr-DD Model in Transient State
For Poisson equation, we directly apply the second-order finite central difference on the Laplace operator, then Eq. (7) becomes where the generalized dielectric coefficients are given byε at each time step t k , k = 1, . . . , N, rearranging Eq. (31) gives a matrix equation is the unknown electrostatic potentials inside the boundary and b (k) ϕ ∈ R N x N y is the vector composed of electric charges and known electrostatic potentials on the boundary. Depending on the types of boundary conditions, the electrostatic potentials on the boundary should be either updated repeatedly at each time step or updated just once at the initial time step. The Dirichlet boundary conditions are given as follows: where f ϕ (∂ D ) is a predefined function on the Dirichlet boundary ∂ D . The Neumann or Robin boundary conditions, as given in Eq. (33), should be updated repeatedly at each time step.
where g ϕ (ϕ k−1 ) is a function defined on the interior discrete points, and the form of g ϕ is given by discretizing the continuous Neumann or Robin boundary conditions.
For the electron continuity equation, the diffusion coefficient D p , D n and carrier mobility μ p , μ n can be spatially distributed, so they can be parameterized as similar as the dielectric permittivity and form four independent parameter matrices. Nevertheless, in the ensuing derivations, we treat D p , D n , μ p and μ n as constants within the whole domain in order to reduce the indicial complexity.
For Caputo's space-fractional gradient operator in the electron continuity equation, firstly we will treat the current density flowing through the interval of two consecutive normal grids as a constant, which can result in two Caputo's linear fractional-order ODEs (assume that the time step is at k + 1): Referring to Eq. (24) and Theorem 2.4, the solutions to these two fractional-order ODEs can be obtained as follows: where the generalized state transition functions are defined by In Eq. (8), the gradient of current density is approximated by (8) and approximating the time-fractional derivative using Lemma 2.6, we can get the discrete form of Eq. (8) as follows: i+1, j and C k+1 i, j+1 are given in Eqs. (41)-(45), where the generalized state transition functions 1 (t) and 2 (t) are given in Eqs. (38) and (39), while 3 (t) and 4 (t) are given by For the left Riemann-Liouville integrals appearing in Eqs. (41)-(45), the three-point Simpson's rule inferred by Lemma 2.5 can be applied to obtain their approximated values to order of 3 + β.
For instance, J β 0 + 1 ( y) is approximated by The other Riemann-Liouville integrals J , and J β 0 + 4 ( y) can thus be approximated accordingly. As a result, Eq. (40) can then be arranged in a matrix form is a pentadiagonal matrix, n (k+1) is the unknown interior electron concentrations at the (k + 1)-th time step, b (k+1) n is a column vector composed of the known boundary electron concentration at the (k + 1)-th time step and the known interior electron concentrations solved from the previous k time steps.
The hole continuity equation in Eq. (9) can also be discretized with a similar procedure, and its discrete form is obtained in Eq. (49), i+1, j and D k+1 i, j+1 are given in Eqs. (50)-(54).
where we denote the generalized reversed state transition functions byˆ 1 (t),ˆ 2 (t),ˆ 3 (t) and 4 (t), and their definitions are given by Eqs. (55)-(58). Similarly, the left Riemann-Liouville integrals of the generalized state transition functions can also be approximated by three-point Simpson's rule.
As remarkably similar in format to Eqs. (40), (49) can also be transformed into a compact . The boundary conditions for these linear matrix equations are then specified as follows: where f n , f p are the predefined functions in Dirichlet boundary conditions, and g n , g p are the functions derived from the Neumann or Robin boundary conditions. Furthermore, the initial value conditions are specified in Eq.
Step-1.3 Update (g+1) n k+1 i, j by damping results to improve the convergence rate.

Special Case when α = 1 and β = 1
When α = 1 and β = 1, the Fr-DD model becomes the conventional DD model which is universally employed in the modeling of crystalline semiconductors (e.g., Si, Ge, etc.). In this case, through simple calculations and substitutions, it can be verified that Eqs. (40) and (49) where the new coefficients are defined as B k+1 (n, m), j = B is the Bernoulli function. Here the discretized system of equations formed by Eqs. (31), (64) and (65) is identical to the discretized system of equations derived from the well-known Scharfetter-Gummel method [9]. Likewise, Algorithm 1 and 2 can be easily modified to solve the numerical solution of the conventional DD model.

Consistency and Convergence Analysis
The proposed discretization scheme is consistent if the truncation error terms can be made to vanish as the mesh and time step size is reduced to zero. First of all, the consistency of the finite center difference scheme applied to the Poisson equation can be easily proved [40]. Furthermore, it can be inferred from Lemma 2.5 and Lemma 2.6 that the truncation error of the discretized carrier continuity equations in Eqs. (40) and (49) will vanish as the spatial and time step sizes shrink to zero. Nevertheless, Eq. (16) hints that an additional truncation error can be generated by composing Caputo's fractional derivative terms in the current density with an integer-order gradient operator on the left side of the equation. To test the influence of this truncation error on the consistency of Eqs. (40) and (49), we propose Theorem 4.1, which gives the shrinking order of this truncation error with the spatial step sizes. Theorem 4.1 Consider the two-dimensional divergence terms ∇ · I and ∇ · J in Eqs. (8) and (9) with I = −qμ n n∇ϕ + qD n C ∇ β r n and J = −qμ p p∇ϕ − qD p C ∇ β r p, then the following equations hold for 0 < β < 1.
where i x , i y are unit vector in the x and y direction.
Proof. Observing that both equations are similar in structure and symmetric in x and y directions, it is sufficient to prove Eq. (66) in only the x direction. From Eq. (12), we obtain If we take the first derivative of both sides of Eq. (68), we can obtain According to Eq. (12), RL 0 D β+1 x (n(x, y)) can be expanded as Substituting Eq. (70) into Eq. (69) and observing that (1 − β) = −β (−β), we can obtain Then we can see Eq. (66) as a corollary to Eq. (71). This completes the proof.
In the derivation of Eqs. (40) and (49), we treat current I and J as constants and solve Caputo's linear fractional-order ODE within two consecutive grid points. Therefore, from Theorem 4.1, we can get where I = IX i+1/2, j , IY i, j+1/2 is the augmented electron current density vector. Eq. (72) shows that the discretized current density terms can be composed with an integer-order gradient operator from the left with an inclusive truncation error. By forcing x, y → 0, these truncation error terms will decay with the spatial step size to a fractional order 1 − β. This then consolidates our claims on the consistency of the proposed discretization scheme.
For convenience, the convergence analysis is only performed on Algorithm 2, but the conclusions of the analysis also apply to Algorithm 1 due to its structural similarity to Algorithm 2 within each step of time advancement. Let us begin our analysis by setting up a finite dimensional vector space V := ϕ ∈ R N x N y : ϕ ∞ < ∞ and a product vector The Gummel map is a mapping A : Y −→ V ×Y that relates a pair [n 1 , p 1 ] to a triplet (ϕ, n 2 , p 2 ). Therefore, the Gummel mapping in Algorithm 2 can be represented by a series of linear matrix computations in Eqs. (73)-(75), We need to show that the Gummel mapping is a contraction mapping with contraction constant L < 1, and the fixed point theorem guarantees the convergence of the algorithms.

Theorem 4.2 The Gummel mapping A is a contraction mapping if we consider the successively over relaxation mechanism in the Gummel iteration.
Proof. substitute (g) ϕ into Eqs. (74) and (75), we have By damping the intermediate results, we can get Taking quotient on both sides and applying triangular inequality yield Since the relative sizes of (g+1)n and (g+1)p to (g) n and (g) p are indeterminable, we need to discuss the following two cases. For (g+1)n ≥ (g) n and (g+1)p ≥ (g) p , Steps 1.4-1.6 in Algorithm 2 imply that ω n −→ 0 and ω p −→ 0 if divergence count d exceeds 1000. Therefore, A ≤ 1 and iteration converge in large g. For (g+1)n < (g) n or (g+1)p < (g) p , it can be inferred from Eqs. (80) and (81) that A < 1. This completes the proof.

Numerical Examples
In this section, we consider three numerical examples to evaluate the accuracy and demonstrate the computational performance of our Fr-DD model solver. All the numerical computations below are based on a MATLAB (R2019b) subroutine and performed on a laptop (MacBook Pro 2019) with Intel Core i9 CPU and 16 GB of RAM.
Example 5.1 Consider the following single-carrier transport problem with fractional derivatives in both time and space for (x, y) ∈ and t > 0.
where = (0, 1)× (0, 1), 0 < α ≤ 1, and 0 < β ≤ 1. The exact solution to this problem is prescribed as ϕ = − 1 2π 2 exp −2π 2 t sin (π x) sin(π y) and n = exp(−2π 2 t) sin(π x) sin(π y), where the nonlinear term on the RHS of Eq. (83) is then given by + sin 2 (π x) cos (2π y) − exp(−2π 2 t) π β+1 sin (π x + (β + 1) π 2 sin(π y) + π β+1 sin π y + (β + 1) π 2 sin(π x) Example 5.1 is a benchmark problem constructed by the method of manufactured solutions [41]. The ground truth is known with its solutions at t = 0.02 s sketched in Fig. 2. The ground truth is compared to the numerical solutions to evaluate the convergence order of our algorithms. The error in this example is calculated by a variant form of the Frobenius norm acting on the error matrix, i.e., e (τ , Here, the spatial step sizes in x and y dimensions are both given by x = 1/(N + 1), where N is the number of internal grid points in one dimension.
To verify the convergence order in time, we make the spatial step size x small enough, such as x = 0.01 in this case, to ensure that the spatial discretization error is much smaller than the time discretization error. With different temporal step sizes τ , the numerical errors and the CPU times can be obtained and are shown in Tabs. 1 and 2, respectively. In Tab. 1, we compare three different combinations of α and β for fixed β = 1, and it is observed that most of the numerical convergence orders in time lie within (1,2). This observation gives an estimate for the convergence order in time and shows its dependence on the time-derivative order α. In Tab. 2, we compare three different combinations of α and β for fixed α = 1. With space-derivative order β varying, the convergence orders in time approach zero, implying the negligent effect of β on the convergence order in time. The overall error of our discretization scheme is dominated by the approximation error of the Riemann-Liouville integral (See Eq. (48)), as the approximation error of the Riemann-Liouville integral is only determined by the β value and the spatial step size x. Therefore, as the temporal step size τ shrinks the error in Tab. 2 remains relatively unchanged. The overall error can be further decreased by reducing the spatial step size x, which is consistent with the trend observed in Tab. 4. However, it should be mentioned that the approximation of the Riemann-Liouville integral will not limit the overall accuracy of the scheme in the case of β = 1 (See Tab. 1) since we can analytically calculate the Riemann-Liouville integral when β = 1. In addition, it can also be observed in Tabs. 1 and 2 that the speed at which the CPU time increases is less than the speed at which the temporal step size shrinks suggesting that the computational error of the solver can be reduced to any pre-set magnitude at the cost of a relatively small increase in CPU time.   To check the spatial convergence order, we take a sufficiently small temporal step size τ = 0.00001 to guarantee that the temporal discretization errors can be neglected compared with the spatial errors. Similar to Tabs. 1 and 2, we record the errors and CPU times calculated under different spatial step sizes in Tabs. 3 and 4. In Tab. 3, the space-derivative order β is fixed to 1 and then forms three different combinations with α. It is observed that the distribution of the spatial convergence order is not uniform in α, and the order decreases dramatically as the spatial step size decreases. In Tab. 4, we set three different combinations of α and β for fixed α = 1. It can be noted that the spatial convergence order is very close to 1 regardless of the change in β, which reveals the linear dependency of the scheme error on the spatial step size when β < 1. The CPU time under different spatial step sizes does not exhibit a specific growth trend as what we observed in Tabs. 1 and 2. However, considering that the growth rate of the number of discrete spatial grids is the square of the reduction rate of the spatial step size, the CPU time still grows at a slower rate relative to the growth of the number of discrete spatial grids. According to these observations, we can infer that the solver precision can indeed be raised to a certain level at the expense of a relatively small increase in computational complexity (CPU time). Table 3: The errors, numerical spatial convergence orders and CPU times for different spatial step sizes x with fixed temporal step size τ = 1e − 5 and fixed space-derivative order β = 1  Example 5.2 Consider the following steady-state single-carrier transport problem in a 2D ptype organic field effect transistor (OFET).
where the effective hole mobility μ p and diffusion coefficient D p are constants for homogeneous materials. The net generation-recombination rate G p ≈ 0 since the generation and recombination activities are relatively unimportant in OFETs as a majority carrier device [33,37]. The solution domain is defined in Fig. 3, where the size of the organic semiconductor (OSC) layer is 500 μm × 30 nm and the size of the dielectric layer is 500 μm× 64 nm. The parameters for OFET simulation are presented in Tab. 5, and the diffusion coefficient for OSC is determined by Einstein's relation D p = V T μ p . The geometric sizes and the material types of the OFET domain are the same as the OFET fabricated in [36] and the material parameters are taken from [42,43]. The encapsulating layer (Parylene) in [36] is not considered in this numerical example in order to simplify boundary conditions.  Similarly, the boundary potential on AH is specified by ϕ AH = −3.0 V. These electric potential boundary values are all reasonably selected, and the OFET has been proven to work normally under these boundary potentials [36,42]. If we assume the dielectric layer is a perfect isolator, the hole concentration on region ABGH should be 0. Applying the above boundary conditions and Algorithm 2, we can obtain the simulated steady-state electric potentials and hole concentrations within the solution domain for β = 1 and 0.8, as shown in Fig. 4. However, it should be noted that in Example 5.1 and 5.2 only cases where β ≥ 0.7 can be simulated in our program since four generalized (reversed) state transition functionŝ i ( x) (i = 1, 2, 3, 4) may blow up to infinity as x is small and β < 0.7. For the cases of β = 1 and β = 0.8, we impose the specified boundary conditions and obtain the simulated surface plots of electric potentials and hole concentrations in Fig. 4. It is observed that the distribution of electric potentials is not significantly affected by the selection of different β. Nevertheless, the profiles of hole concentration under different β are obviously shifted along the thickness direction, which implies that different β values are related to the intensity of charge carriers' diffusion motions.

Figure 4:
The simulated steady-state electric potentials and hole concentrations within the thinner solution domain for an OFET when space-derivative order β = 1 and β = 0.8, respectively more within each step of time advancement (i.e., p ≈ Cτ α ). For the second group of numerical experiments, we fix α = 0.9 and set β = 1, 0.8 and 0.6 to solve for solutions at t = 1e − 5 s. It is found in Fig. 5b that the decay rate for hole concentration reduces as the decrease of β, and this phenomenon can be ascribed to the more inactive diffusive motions of charge carriers under smaller β. In Example 5.2, we simplify boundary conditions to better discuss the influence of β values on the charge carriers' diffusive motions in the steady-state OFET. As an extension to Example 5.2, this section provides the experimental validation for the Fr-DD model. The fabrication and the experimental characterization of the OFET that we model was discussed in [36]. As shown in Fig. 6, compared to the simplified structure of the OFET in Example 5.2, the complete structure of the OFET is encapsulated in a polymer layer made of Parylene and all the metallic electrodes have a thickness of 30 nm. The material parameters are specified in Tab. 5. In addition to the boundary conditions given in Example 5.2, we should also treat the encapsulating layer as a perfect insulator, where no charge carrier is transmitted, and the dielectric displacement should be continuous on its borders. Consider that the length of the source electrode L S and the drain electrode L D are both 200 μm and the width of the OFET (out-of-plane dimension) W is 1000 μm, the net current flowing through the drain electrode is evaluated by where J y is the y-component of the continuous current density at the Au-DNTT interface, i D is the discrete grid index in the x-direction and j D is the discrete grid y index at the Au-DNTT interface. In our case, the spatial step sizes x = 5 μm and y = 1 nm, so we can fix j D = 184 and calculate the summation over i D = 1, . . . , 40. The fractional Riemann-Liouville integral J β 0 +ˆ 1 ( y) is then evaluated using Eq. (48).
The β value in the fractional drift diffusion OFET model depends on the spatial coordinates and electrode potentials, i.e., β = B(x, y, V gs , V ds ). This inhomogeneity of β for different spatial coordinates and boundary conditions is caused by the irregular crystalline structure of OSCs and the electronic polarization under different boundary conditions [44,45]. If we ignore the dependence of β on spatial coordinates and only consider its dependence on V gs and V ds , we can obtain the relationship curves for β = B(V gs , V ds ) in Figs. 7b and 8b by fitting the experimental data.
When β is adjusted for a fixed V ds and varying V gs according to Fig. 7b, we can utilize Eq. (88) to calculate the drain current I ds under different V gs and obtain the theoretical transconductance. In Fig. 7a, it is found that the theoretical transconductance curve (solid black line) is in good agreement with the experimentally measured transconductance (red circles). Similarly, if we adjust β for varying V ds and four fixed V gs according to dependence curves in Fig. 8b, we can notice that the theoretical output curves (black solid lines) can well fit with the experimentally characterized outputs (Symbol, i.e., red circles, blue squares, etc.) as shown in Fig. 8a. These results not only confirm the validity of the Fr-DD model for predicting the OFET characteristics, but also suggest the highly nonlinear dependence of β value on the electrode potentials V gs and V ds . The relationship curves between β and electrode potentials can be constructed very quickly by writing a simple optimization subroutine to minimize the error between the experimental data and the theoretical predictions. Due to the flexibility of adjusting β, the Fr-DD model avoids considering the involuted trap states in OSCs, thus greatly improving the modeling efficiency of OSC devices compared with the conventional OFET analytic models involving the trap or impurity states.

Conclusion
This work aimed to develop a Fr-DD model solver for simulating the anomalous dynamics of OSC devices. Two algorithms based on a novel discretization scheme and successively overrelaxed Gummel's iteration are proposed here to solve the transient and steady-state Fr-DD model equations. This study has identified the consistency of the two algorithms by showing that the truncation error from the discretized divergence of current density functions will vanish with the spatial step size to a positive fractional order of 1 − β. The convergence analysis reveals that the Gummel mapping is a contraction mapping if we consider the successive over-relaxation mechanism in the Gummel's iteration, which thus completes the proof of convergence. Three numerical examples, including one benchmark example and two others constructed from the perspective of engineering applications are employed to demonstrate the algorithms' accuracy and computational performance. It is found in the first example that altering α and β can impact the spatial convergence order but only varying α will affect the convergence order in time. The increase rate of CPU time is less than the shrinking rate of temporal step size and lower than the growth rate of spatial grid points. These findings suggest that our solver has high precision and fast computational speed as it limits the computational error to a predefined satisfactory level (from ∼10 −4 to ∼10 −6 ) at a relatively small expense of CPU time (from ∼20 s to ∼100 s). The results reported in two numerical examples reveal the prediction and characterization of the transient-state and steady-state dynamics for any type of organic semiconductor device. Finally, we provide experimental verification for the fractional drift diffusion model of a DNTT-OFET. We stipulate that this is the first to date exploration of the Fr-DD model solver laying the groundwork for future research into fractional drift diffusion modeling of flexible organic electronics.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.