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

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 fractionalorder 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,JnandJp,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 coefficientsDp,Dnand 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 byJn= −qμnn∇ϕ+qDn∇nandJp= −qμpp∇ϕ−qDp∇p,Einstein’s relations areis the thermal voltage,kis Boltzmann constant,Tis the thermal temperature,qis the charge of an electron,ε is semiconductor’s absolute dielectric permittivity,N+DandN−Aare ionized donor and acceptor concentrations.GnandGpare 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.,However,the long-tail behavior of the photocurrent curve observed in OSCs implies that the MSD in this scenario is given byforα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 carrierspf,nfand trap(localized)charge carrierspt,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),

whereFα(x,t)=τ0cαμpE(x,t)is the anomalous advection coefficient andDα=τ0cαDpis the anomalous diffusion coefficient.The hole mobilityμpand hole diffusion coefficientDpsatisfies the Einstein relation.Eq.(6)coupled with the 1D Poisson equation forms the 1D fractional drift diffusion(Fr-DD)model.1This was initially a simplified Fr-DD model with only time-derivative fractionalized,the order of spatial derivatives remained integer.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),

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

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.3Assume f∈Cn−1([a,t])and n−1<γ≤n,then the Laplace transform of Riemann–Liouville and Caputo’s fractional derivatives are given by

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

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.4Consider 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.

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 functionis obtained by

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

whereId=1 orIdis the identity matrix in caseu,vare 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.Lettingx=xi,we see that two consecutive states have a relation expressed by

Let us assume that the input functionv(t)=1,the fictitious input function is then evaluated byand 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.

2.3 Discrete Approximation of Fractional Integrals and Fractional Derivatives

where the mean value theorem is applied forf(t)andf′(t)withξ1,ξ2∈(tm,tm+1),and the continuous functionf(2)on a compact domain[0,T]assumes its maximum value,we letM=maxt∈[0,T]f(2)(t).

3 Derivation of the Computational Scheme

Figure 1:An illustration of notations in the discrete space for(a)discretized electron continuity equation with electron concentrationx-direction component of electron current density and y-direction component of electron current densitydiscretized hole continuity equation with hole concentration x-direction component of hole current densityand y-direction component of hole current density

3.1 Discretization of Fr-DD Model in Transient State

wherefn,fpare the predefined functions in Dirichlet boundary conditions,andgn,gpare the functions derived from the Neumann or Robin boundary conditions.Furthermore,the initial value conditions are specified in Eq.(61).Given the discrete form of the transient-state Fr-DD model in Eqs.(31),(40)and(49)and the consistency between the initial value and boundary conditions,we propose Algorithm 1 to solve the unknownsϕ,nandpfor each time step.

Algorithm 1:To evaluate the numerical solution of the Fr-DD model in transient state Input:Constant damping parameters ωn,ωp ∈[0,1];Initial guess 0nki,j and 0pki,j for i=1,2,...,Nx;j=1,2,...,Ny; k=1,...,N Output:Unknown interior variables at each time step ϕki,j,nki,j,pki,j for i = 1,2,...,Nx;j=1,2,...,Ny; k=0,1,...,N Step-1.0 Obtain in lectron a itial potentials ϕ0i,j by solving Eq.(31)in presence of the initial conditions of end hole concentrations in Eq.(61)and the boundary conditions of potentials in Eqs.(32)and(33).For each time step k=0,1,...,N −1,do Initialize Gummel iteration counts g=0,old error Err0=1,and divergence counts d=0.While Err>Tol,do Gummel iterations Step-1.1 Generate(g)ϕk+1 i,j by solving Eq.(31)with the initial guess(g)nk+1i,j ,(g)pk+1i,j and the boundary conditions of potentials in Eqs.(32)and(33).Step-1.2 Generate second guess(g+1)nk+1i,j ,(g+1)pk+1i,j by solving Eqs.(40)and(49)with the boundary conditions in Eqs.(59)and(60)and initial conditions(convergent solution from previous time steps) nli,j,pli,j,where l=0,1,...,k.

Step-1.3 Update(g+1)nk+1i,j = ωn ·(g+1)nk+1i,j +(1 −ωn) ·(g)nk+1i,jand(g+1)pk+1i,j =ωp ·(g+1)pk+1 i,j +images/BZ_253_830_419_848_465.png1 −ωp)·(g)pk+1 i,j by damping results to improve the convergence rate.Step-1.4 Compute error Err1 =■■■■■(g+1)nk+1 i,j −gnk+1 i,j gnk+1 i,j■■■■■,Err2 =■■■■■(g+1)pk+1i,j −gpk+1 i,j gpk+1 i,j■■■■■ and Err=max(Err1,Err2).Update iteration counts g=g+1.Step-1.5 If Err>Err0 Update divergence counts d=d+1 and old error Err0=Err.Step-1.6 If d>1000 Update damping parameter ωn=ωn/2 and ωp=ωp/2 to improve the convergence ability,then reset divergence counts d=0.

3.2 Discretization of Fr-DD Model in Steady State

Since Caputo’s fractional derivative of any constant is zero,the time-derivative term with Caputo’s fractional derivatives vanishes in steady state.In contrast to the transient-state Fr-DD model,the discretized steady-state Fr-DD model is formed by Eqs.(31),(62)and(63).

The boundary conditions for Poisson’s equation and the carrier continuity equations are specified in Eqs.(32),(33),(59)and(60).By rearranging Eqs.(31),(62)and(63),three matrix equations,i.e.,Aϕϕ=bϕ,Ann=bnandApp=bp,can be formed for algebraic computations.As a result,we propose Algorithm 2 to solve the numerical solution of the steady-state Fr-DD model.

Algorithm 2:To evaluate the numerical solution of the Fr-DD model in steady state Input:Constant damping parameters ωn,ωp ∈[0,1];Initial guess 0ni,j and 0pi,j for i=1,2,...,Nx;j=1,2,...,Ny Output:Unknown interior variables in steady state ϕi,j,ni,j,pi,j for i = 1,2,...,Nx;j=1,2,...,Ny Initialize Gummel iteration counts g=0,old error Err0=1,and divergence counts d=0.While Err>Tol,do Gummel iterations Step-1.1 Generate(g)ϕi,j by solving Eq.(31)with the initial guess(g)ni,j,(g)pi,j and the boundary conditions of potentials in Eqs.(32)and(33).Step-1.2 Generate second guess(g+1)ni,j,(g+1)pi,j by solving Eqs.(62)and(63)with the boundary conditions in Eqs.(59)and(60).Step-1.3 Update(g+1)ni,j=ωn·(g+1)ni,j+(1 −ωn)·(g)ni,j and(g+1)pi,j=ωp·(g+1)pi,j+images/BZ_253_1872_2583_1891_2629.png1 −ωp)·(g)pi,j by damping results to improve the convergence rate.

Step-1.4 Compute error Err1 =■■■■■(g+1)ni,j −gni,j gni,j■■■■■,Err2 =■■■■■(g+1)pi,j −gpi,j gpi,j■■■■■ and Err =max(Err1,Err2).Update iteration counts g=g+1.Step-1.5 If Err>Err0 Update divergence counts d=d+1 and old error Err0=Err.Step-1.6 If d>1000 Update damping parameter ωn=ωn/2 and ωp=ωp/2 to improve the convergence ability,then reset divergence counts d=0.

3.3 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)degenerate into Eqs.(64)and(65),respectively.

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

5 Numerical Examples

Example 5.1 is a benchmark problem constructed by the method of manufactured solutions[41].The ground truth is known with its solutions att=0.02ssketched 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,Here,the spatial step sizes inxandydimensions are both given byΔx=1/(N+1),whereNis the number of internal grid points in one dimension.

To verify the convergence order in time,we make the spatial step sizeΔxsmall 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.

Figure 2:The contour plots(top)and surface plots(bottom)of the electric potentials(left)and the electron concentrations(right)at t=0.02 s

Table 1:The errors,numerical convergence orders in time and CPU times for different temporal step sizes τ with fixed spatial step size Δx=0.01 and fixed space-derivative order β=1

Table 2:The errors,numerical convergence orders in time and CPU times for different temporal step sizes τ with fixed spatial step size Δx=0.01 and fixed time-derivative order α=1

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

Table 4:The errors,numerical spatial convergence orders and CPU times for different spatial step sizes Δx with fixed temporal step size τ=1e −5 and fixed time-derivative order α=1

Example 5.2Consider the following steady-state single-carrier transport problem in a 2D ptype organic field effect transistor(OFET).

where the effective hole mobilityμpand diffusion coefficientDpare constants for homogeneous materials.The net generation-recombination rateGp≈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 relationDp=VTμ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.

Figure 3:The solution domain of a 2D top-contact bottom-gate(TCBG)OFET device composed of a p-type organic semiconductor layer and a dielectric layer

Table 5:The parameters for OFET simulation

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

Table 6:The parameters for solar cell simulation

In Example 5.3,we let the spatial step size be 1e−8 m and the temporal step size be 1e−6 s.First,we fixβ=1 and setα=0.8,0.6 and 0.4 to solve for solutions att=1e−5 s.As shown in Fig.5a,the hole concentration in this setting displays a decaying trend with the decrease ofα,while the electric potential remains almost the same for differentα.The decay trend in hole concentration can be easily predicted since the hole concentration with smallerαwill reduce 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 att=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β.

Figure 5:The simulated transient-state electric potentials and hole concentrations within the thinner solution domain for the solar cell when(a)space-derivative order β=1 fixed and α=0.8,0.6 and 0.4,respectively;(b)time-derivative order α=0.9 fixed and β=1,0.8 and 0.6,respectively

6 Experimental Validation of the Fractional Drift-Diffusion OFET Model

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.

Figure 6:The solution domain of a 2D top-contact bottom-gate(TCBG)OFET device composed of a p-type organic semiconductor layer,a dielectric layer and encapsulating layers(top and bottom)

Consider that the length of the source electrodeLSand the drain electrodeLDare both 200 μm and the width of the OFET(out-of-plane dimension)Wis 1000 μm,the net current flowing through the drain electrode is evaluated by

whereJyis the y-component of the continuous current density at the Au-DNTT interface,iDis the discrete grid index in the x-direction andjDis 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 fixjD=184 and calculate the summation overiD=1,...,40.The fractional Riemann-Liouville integralis 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,Vgs,Vds).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 onVgsandVds,we can obtain the relationship curves forβ=B(Vgs,Vds)in Figs.7b and 8b by fitting the experimental data.

Whenβis adjusted for a fixedVdsand varyingVgsaccording to Fig.7b,we can utilize Eq.(88)to calculate the drain currentIdsunder differentVgsand 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 varyingVdsand four fixedVgsaccording 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 potentialsVgsandVds.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.

Figure 7:(a)The experimentally measured transconductance at a fixed Vds=−5 V compared with the fitted theoretical transconductance curve obtained from the fractional drift diffusion model;(b)The adjusted β values at different Vgs and a fixed Vds=−5 V

Figure 8:(a)The experimentally measured output curve at a series of fixed Vgs = −5∼−2 V compared with the fitted theoretical output curve obtained from the fractional drift diffusion model;(b)The adjusted β values at different Vds and a series of fixed Vgs=−5∼−2 V

7 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−4to ∼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.

Funding Statement:This work was supported in part by the National Science Foundation through Grant CNS-1726865 and by the USDA under Grant 2019-67021-28990.

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