Computers, Materials & Continua DOI:10.32604/cmc.2021.017439
Article

A Fractional Drift Diffusion Model for Organic Semiconductor Devices

Yi Yang*, Robert A. Nawrocki, Richard M. Voyles and Haiyan H. Zhang

School of Engineering Technology, Purdue University, West Lafayette, 47906, USA
*Corresponding Author: Yi Yang. Email: yang1087@purdue.edu
Received: 30 January 2021; Accepted: 8 March 2021

Abstract: 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.

Keywords: Fractional drift diffusion model; Gummel’s iteration; Caputo’s fractional-order ordinary differential equation; organic field effect transistor

1  Introduction

The mathematical modeling of the electrons and holes transports in an inorganic semiconductor (ISC) is established by a system of coupled partial differential equations (PDEs), which are formulated by Gauss’ law applied to the electrical potential φ, and the continuity of the electron and hole current densities, Jn and Jp, respectively [1,2]. Besides modeling of ISCs, this system of coupled PDEs, namely, the drift diffusion (DD) model, has also found extensive applications in modeling other diffusion-reaction processes, such as ion exchanges in the electrochemical solvents [3,4], and the transports of positive/negative particles within cell membranes [5,6]. Depending on different application scenarios, the DD model can be represented 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 Dp, Dn 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),

−Δ(εφ)=q(p−n+ND+−NA−)(1)

∂n∂t=1q∇⋅Jn+Gn(2)

∂p∂t=−1q∇⋅Jp+Gp(3)

where current densities are given by Jn=−qμnn∇φ+qDn∇n and Jp=−qμpp∇φ−qDp∇p, Einstein’s relations are Dnμn=Dpμp=VT=kTq, VT is the thermal voltage, k is Boltzmann constant, T is the thermal temperature, q is the charge of an electron, ε is semiconductor’s absolute dielectric permittivity, ND+ and NA− are ionized donor and acceptor concentrations. Gn and Gp 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–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(x2(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(x2(t))∝tα, for α termed as the dispersive parameter of the OSC, 0<α<1, depending on the temperature and band structure disorders [20–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–29], respectively. Based on the “multiple trapping model,” the charge carriers in OSCs can be classified as free (delocalized) charge carriers pf,nf and trap (localized) charge carriers pt, nt. 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),

∂pt(x,t)∂t=1τ0cα0RLDtα(pf(x,t))=1τ0cα1Γ(1−α)∂∂t∫0tpf(x,s)(t−s)αds(4)

where 0RLDtα 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/E0 is the dispersive parameter which depends on the temperature T and the expected value of the exponential density of trap states E0. The 1D continuity equation for free charge carriers in p-type OSCs was also derived by References [18,30] as given in Eq. (5),

∂pf(x,t)∂t+1τ0cα0RLDtα(pf(x,t))+∂∂x[μpE(x,t)pf(x,t)]−Dp∂2pf(x,t)∂x2=p(x,0)δ(t)(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 = pf + pt and pt≫pf in OSCs. Substituting Eq. (4) into Eq. (5) can derive the continuity equation for total charge carrier density, as given in Eq. (6),

0RLDtα(p(x,t))+∂∂x[Fα(x,t)p(x,t)]−Dα∂2p(x,t)∂x2=p(x,0)t−αΓ(1−α)(6)

where Fα(x,t)=τ0cαμpE(x,t) is the anomalous advection coefficient and Dα=τ0cαDp is the anomalous diffusion coefficient. The hole mobility μp and hole diffusion coefficient Dp 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–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),

−Δ(εφ)=q(p−n)(7)

0CDtα(n(x,y,z,t))=1q∇⋅I+Gn(8)

0CDtα(p(x,y,z,t))=−1q∇⋅J+Gp(9)

where 0CDtα 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μnn∇φ+qDnC∇rβn and J=-qμpp∇φ-qDpC∇rβp, and C∇rβ=(0CDxβ0CDyβ0CDzβ) 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 ND+ and acceptors NA- 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 as ∇⋅J=-qμp(pΔφ+∇p⋅∇φ)-qDpRL∇rβ+1p, similarly for ∇⋅I in Eq. (8) [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.

2  Preliminaries

2.1 Definition of Fractional Operators

The Riemann–Liouville (RL) fractional derivative with order γ>0 is defined in Eq. (10) [34],

0RLDtγ(f(t))=1Γ(n−γ)dndtn∫0tf(τ)(t−τ)γ+1−ndτ,n−1<γ<n(10)

where Γ(n-γ) denotes the gamma function. Similarly, Caputo’s fractional derivative with order γ>0 is defined as Eq. (11) [34].

0CDtγ(f(t))=1Γ(n−γ)∫0tf(n)(τ)(t−τ)γ+1−ndτ,n−1<γ<n(11)

Both RL and Caputo’s fractional derivatives can be considered as interpolation to integer-order derivatives, which means 0RLDtn(f(t))=0CDtn(f(t))=f(n)(t). The relationship between RL can Caputo’s fractional derivatives can be expressed as the following lemma.

Lemma 2.1 Assumef∈Cn-1([a,t])andn-1<γ≤n, then the following equality holds

aCDtγ(f(t))=aRLDtγ(f(t))−∑k=0n−1f(k)(a)Γ(k−γ+1)(t−a)k−γ(12)

Proof. See [35].

By directly employing the definitions, the composition rules for fractional derivatives can be given as the following lemma.

Lemma 2.2 Assumef∈Cn+m-1([a,t]), n-1<γ≤n, and m > 0 is an integer, then the following relations hold

dmdtm[aRLDtγ(f(t))]=aRLDtγ+m(f(t))(13)

aCDtγ(dmdtmf(t))=aCDtγ+m(f(t))(14)

aRLDtγ(dmdtmf(t))=aRLDtγ+m(f(t))−∑k=0m−1f(k)(a)Γ(1+k−γ−m)(t−a)k−γ−m(15)

dmdtm[aCDtγ(f(t))]=aCDtγ+m(f(t))+∑k=nn+m−1f(k)(a)Γ(1+k−γ−m)(t−a)k−γ−m(16)

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 haveL{aRLDtγ(f(t))}=sγF(s)−∑k=0n−1sk⋅aRLDtγ−k−1f(a) 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.

Lemma 2.3 Assumef∈Cn-1([a,t])andn-1<γ≤n, then the Laplace transform of Riemann–Liouville and Caputo’s fractional derivatives are given by

L{ aRLDtγ(f(t)) }=sγF(s)−∑k=0n−1sk⋅aRLDtγ−k−1f(a)(17)

L{aCDtγ(f(t))}=sγF(s)−∑k=0n−1sγ−k−1⋅f(k)(a)(18)

Proof. See [35].

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].

L{tβ−1Eα,β(±atα)}=sα−βsα∓a,R(s)>0,R(α)>0,R(β)>0(19)

Subsequently, we will present the analytic solution for Caputo’s fractional linear time-invariant (LTI) state equation.

2.2 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 domain with x∈[xi-1,xi]and0<γ≤1, as given in Eq. (20), where u(x) is the state variable and v(x) is the input variable.

xi−1CDxγu(x)=Au(x)+Bv(x)(20)

Then, its solution is given by

u(x)=Φ(x−xi−1)u(xi−1)+∫0x−xi−1Φ(x−xi−1−y)Bv^(y)dy(21)

whereΦ(x)=Eγ(Axγ) is the generalized state transition function, Eγ(t) is the one-parameter Mittag–Leffler function, and the fictitious input function v^(y) is obtained by v^(x)=L-1{V(s)s1-γ}.

Proof. Apply Laplace transform on both sides of Eq. (20), it gives

sγU(s)−sγ−1u(xi−1)=AU(s)+BV(s)(22)

Rearrange both sides of Eq. (22) and take the inverse Laplace transform, we haveu(xi)=Φ(Δx)u(xi−1)+B∫0ΔxΦ(y)(Δx−y)γ−1Γ(γ)dy 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 Ωh={xi=iΔx,i=0,1,2,…,N} with step size Δx=L/N. Letting x = xi, we see that two consecutive states have a relation expressed by

u(xi)=Φ(Δx)u(xi−1)+∫0ΔxΦ(Δx−y)Bv^(y)dy(23)

Let us assume that the input function v(t)=1, the fictitious input function is then evaluated by v^(y)=L-1{s-γ}=yγ-1Γ(γ), and considering the commutative property of the convolution integral, then Eq. (23) in this special case can be reformulated as

u(xi)=Φ(Δx)u(xi−1)+B∫0ΔxΦ(y)(Δx−y)γ−1Γ(γ)dy(24)

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]

J0+γf(x)=1Γ(γ)∫0xf(y)(x−y)1−γdy(25)

In the next section, we present discrete approximation formula for the left Riemann–Liouville integral and fractional derivatives.

2.3 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.

Lemma 2.5 Assume0<γ≤1andf∈C4([0,x]) then the following (3+γ)-th order approximation of the left Riemann–Liouville integral can be obtained

J0+γf(x)=(Δx)γ[∑k=0nb2k(γ)f(x2k)+∑k=1nb2k−1(γ)f(x2k−1)]+O(Δx3+γ)(26)

wherexj=jΔx with a positive integer j and step size Δx. Since we let x = x2n, the coefficients bj(γ) satisfies following formula,

b2k−1(γ)=−2(2n−2k+2)2+γ−(2n−2k)2+γΓ(3+γ)+2(2n−2k+2)1+γ+(2n−2k)1+γΓ(2+γ)(27)

b2k(γ)={ (2n)2+γ−(2n−2)2+γΓ(3+γ)−3(2n)1+γ+(2n−2)1+γ2Γ(2+γ)+(2n)γΓ(1+γ), k=0(2n−2k+2)2+γ−(2n−2k−2)2+γΓ(3+γ)−(2n−2k+2)1+γ+6(2n−2k)1+γ+(2n−2k−2)1+γ2Γ(2+γ),360ptk=1,…,n−122+γΓ(3+Γ)−21+γ2Γ(2+γ), k=n (28)

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.

Lemma 2.6 Assume0<γ≤1andf∈C2([0,T]) then the following first-order approximation of the Caputo’s time-fractional derivative can be obtained

0CDtγ(f(t))=1Γ(2−γ)∑m=0kf(tk+1−m)−f(tk−m)Δtγbm,k+1(γ)+O(Δt)(29)

wheretj=jΔt with a positive integer j and step size Δt. Since we let T = tk+1, the coefficients bm,k+1(γ) is determined by

bm,k+1(γ)=(m+1)1−γ−m1−γ(30)

Proof. The discrete approximation in Eq. (29) can be constructed by applying the piecewise quadrature, the error estimates can be derived as follows,φk|∂ΩD=fφ(∂ΩD),k=0,1,…,N where the mean value theorem is applied for f(t) and f′(t) with ξ1, ξ2∈(tm,tm+1), and the continuous function f(2) on a compact domain [0,T] assumes its maximum value, we let M= maxt∈[0,T]f(2)(t).

3  Derivation of the Computational Scheme

Without loss of the generality, we implement the discretization schemes in the two-dimensional 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=Lx/(Nx+1) and in the y direction be Δy=Ly/(Ny+1). The two-dimensional grids are given by xi=iΔx,i=0,1,…,Nx+1 and yj=jΔy,j=0,1,…,Ny+1. In the temporal domain, we let the temporal step size be Δt=T/N, and the temporal girds are then determined by tn=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(tk,xi,yj) and the two components of the electron current density vector (IX(tk,xi+1/2,yj), IY(tk,xi,yj+1/2)) in the discretized electron continuity equation are denoted by ni,jk and (IXi+1/2,jk,IYi,j+1/2k), 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 φi,jk 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.

Figure 1: An illustration of notations in the discrete space for (a) discretized electron continuity equation with electron concentration ni,jk, x-direction component of electron current density IXi,jk, and y-direction component of electron current density IYi,jk; (b) discretized hole continuity equation with hole concentration pi,jk, x-direction component of hole current density JXi,jk, and y-direction component of hole current density JYi,jk

3.1 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

ε~i−1,jφi−1,jk−ε~i,j,1φi,jk+ε~i+1,jφi+1,jkΔx2+ε~i,j−1φi,j−1k−ε~i,j,2φi,jk+ε~i,j+1φi,j+1kΔy2=−q(pi,jk−ni,jk)(31)

where the generalized dielectric coefficients are given by ε̃i-1,j=εi-1/2,j-1/2+εi-1/2,j+1/22, ε̃i+1,j=εi+1/2,j-1/2+εi+1/2,j+1/22, ε̃i,j,1=ε̃i-1,j+ε̃i+1,j, ε̃i,j-1=εi+1/2,j-1/2+εi-1/2,j-1/22, ε̃i,j+1=εi+1/2,j+1/2+εi-1/2,j+1/22, and ε̃i,j,2=ε̃i,j-1+ε̃i,j+1. For i=1,2,…,Nx and j=1,2,…,Ny, at each time step tk, k=1,…,N, rearranging Eq. (31) gives a matrix equation Aφφ(k)=bφ(k), where Aφ∈ℝNxNy×NxNy is a pentadiagonal matrix composed of dielectric permittivity constants, φ(k)=[φ1,1k,φ2,1k,⋯,φNx,Nyk]T is the unknown electrostatic potentials inside the boundary and bφ(k)∈ℝNxNy 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:

φk|∂ΩD=fφ(∂ΩD), k=0,1,…,N(32)

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.

φk|∂ΩN=gφ(φk−1|Ω),k=0,1,…,N(33)

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 Dp, Dn 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 Dp,Dn,μ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):

IXi−1/2,jk+1=−qμnn∂φ∂x+qDn⋅xi−1CDxβn(34)

IYi,j−1/2k+1=−qμnn∂φ∂y+qDn⋅yi−1CDxβn(35)

Referring to Eq. (24) and Theorem 2.4, the solutions to these two fractional-order ODEs can be obtained as follows:

ni,jk+1=Φ2(Δx)ni−1,jk+1+IXi−1/2,jk+1qDnJ0+βΦ2(Δx)(36)

ni,jk+1=Φ1(Δy)ni,j−1k+1+IYi,j−1/2k+1qDnJ0+βΦ1(Δy)(37)

where the generalized state transition functions are defined by

Φ2(t)=Eβ(φi,jk+1−φi−1,jk+1VTtβ−1)(38)

Φ1(t)=Eβ(φi,jk+1−φi,j−1k+1VTtβ−1)(39)

In Eq. (8), the gradient of current density is approximated by ∇⋅I=IXi+1/2,jk+1-IXi-1/2,jk+1Δx+IYi,j+1/2k+1-IYi,j-1/2k+1Δy. By substituting IXi+1/2,jk+1, IXi-1/2,jk+1, IYi,j+1/2k+1 and IYi,j-1/2k+1 into Eq. (8) and approximating the time-fractional derivative using Lemma 2.6, we can get the discrete form of Eq. (8) as follows:

1Γ(2−α)∑m=0kni,jk+1−m−ni,jk−mΔtαbm,k+1(α)=Ci,j−1k+1ni,j−1k+1+Ci−1,jk+1ni−1,jk+1+Ci,jk+1ni,jk+1+Ci+1,jk+1ni+1,jk+1+Ci,j+1k+1ni,j+1k+1+Gn(40)

where bm,k+1(α)=(m+1)1-α-m1-α, and coefficients Ci,j-1k+1, Ci-1,jk+1, Ci,jk+1, Ci+1,jk+1 and Ci,j+1k+1 are given in Eqs. (41)–(45),

Ci,j−1k+1=Φ1(Δy)J0+βΦ1(Δy)DnΔy(41)

Ci−1,jk+1=Φ2(Δx)J0+βΦ2(Δx)DnΔx(42)

Ci,jk+1=[−1J0+βΦ2(Δx)−Φ3(Δx)J0+βΦ3(Δx)]DnΔx+[−1J0+βΦ1(Δy)−Φ4(Δy)J0+βΦ4(Δy)]DnΔy(43)

Ci+1,jk+1=1J0+βΦ3(Δx)DnΔx(44)

Ci,j+1k+1=1J0+βΦ4(Δy)DnΔy(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

Φ3(t)=Eβ(φi+1,jk+1−φi,jk+1VTtβ−1)(46)

Φ4(t)=Eβ(φi,j+1k+1−φi,jk+1VTtβ−1)(47)

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, J0+βΦ1(Δy) is approximated by

J0+βΦ1(Δy)=(Δy)β2Γ(β)[(1−2Φ1(Δy2)+Φ1(Δy))(4β−8β+1+4β+2)−(−3−4Φ1(Δy2)+Φ1(Δy))(2β−2β+1)−4β](48)

The other Riemann–Liouville integrals J0+βΦ2(Δx), J0+βΦ3(Δx), and J0+βΦ4(Δy) can thus be approximated accordingly. As a result, Eq. (40) can then be arranged in a matrix form An(k+1)n(k+1)=bn(k+1), where An(k+1) is a pentadiagonal matrix, n(k+1) is the unknown interior electron concentrations at the (k+1)-th time step, bn(k+1) 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),

1Γ(2−α)∑m=0kpi,jk+1−m−pi,jk−mΔtαbm,k+1(α)=Di,j−1k+1pi,j−1k+1+Di−1,jk+1pi−1,jk+1+Di,jk+1pi,jk+1+Di+1,jk+1pi+1,jk+1+Di,j+1k+1pi,j+1k+1+Gp(49)

where coefficients Di,j-1k+1, Di-1,jk+1, Di,jk+1, Di+1,jk+1 and Di,j+1k+1 are given in Eqs. (50)–(54).

Di,j−1k+1=Φ^1(Δy)J0+βΦ^1(Δy)DpΔy(50)

Di−1,jk+1=Φ^2(Δx)J0+βΦ^2(Δx)DpΔx(51)

Di,jk+1=[−1J0+βΦ^2(Δx)−Φ^3(Δx)J0+βΦ^3(Δx)]DpΔx+[−1J0+βΦ^1(Δy)−Φ^4(Δy)J0+βΦ^4(Δy)]DpΔy(52)

Di+1,jk+1=1J0+βΦ^3(Δx)DpΔx(53)

Di,j+1k+1=1J0+βΦ^4(Δy)DpΔy(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.

Φ^1(t)=Eβ(−φi,jk+1−φi,j−1k+1VTtβ−1)(55)

Φ^2(t)=Eβ(−φi,jk+1−φi−1,jk+1VTtβ−1)(56)

Φ^3(t)=Eβ(−φi+1,jk+1−φi,jk+1VTtβ−1)(57)

Φ^4(t)=Eβ(−φi,j+1k+1−φi,jk+1VTtβ−1)(58)

As remarkably similar in format to Eqs. (40), (49) can also be transformed into a compact matrix equation Ap(k+1)p(k+1)=bp(k+1). The boundary conditions for these linear matrix equations are then specified as follows:

nk|∂ΩD=fn(∂ΩD),nk|∂ΩN=gn(nk−1|Ω),k=0,1,…,N(59)