Drift–diffusion simulation of S-shaped current–voltage relations for organic semiconductor devices

We present an electrothermal drift–diffusion model for organic semiconductor devices with Gauss–Fermi statistics and positive temperature feedback for the charge carrier mobilities. We apply temperature-dependent Ohmic contact boundary conditions for the electrostatic potential and discretize the system by a finite volume based generalized Scharfetter–Gummel scheme. Using path-following techniques, we demonstrate that the model exhibits S-shaped current–voltage curves with regions of negative differential resistance, which were only recently observed experimentally.


Introduction
In inorganic high power transistors and lasers, thermal effects resulting from strong electric and optical fields and from strong recombination processes are of significant importance and have to be included into mathematical models [1].However, electrothermal effects are even more potent in organic semiconductors where the temperature activated hopping transport of charge carriers leads to a strong interplay between electric current and heat flow.They result in interesting phenomena like S-shaped current-voltage relations with regions of negative differential resistance in resistors and organic light-emitting diodes (OLEDs) [2,3] and lead to inhomogeneous luminance in large-area OLEDs.Moreover, electrothermal effects have a strong impact on the performance of organic solar cells and transistors [4][5][6].
As demonstrated in [3], p-Laplace thermistor models are able to capture the positive temperature feedback in OLEDs.Particularly, they can reproduce experimentally observed S-shaped current-voltage relations (see Fig. 1) and inhomogeneous current density and temperature distributions in large-area OLEDs.The p-Laplace thermistor model describes the total current flow and the heat flow in a device.As parameters serve a power law exponent p for the underlying (isothermal) current-voltage relation, an activation energy in an Arrhenius-type temperature law for the hopping transport, and reference electrical and thermal conductivities (see [3]).However, details such as separate electron and hole current flow, space charge regions, generation-recombination and related heat productions as well as energy barriers at material interfaces are neglected.
The aim of the present paper is to investigate the ability of an electrothermal drift-diffusion model for organic semiconductor devices and a suitable numerical approximation to produce S-shaped current-voltage relations.For simplicity, for this proof of concept we use vertically layered device structures.Additionally, we propose a generalization of Ohmic contact boundary conditions for the electrostatic potential to the non-isothermal case.

The PDE system
The electrothermal behavior of organic semiconductor devices is described in a drift-diffusion setting by PDEs for the electrostatic potential , the electrochemical potentials n , p and the temperature T. In the device domain Ω, we consider the following stationary coupled system This system results from the coupling of a generalized van Roosbroeck system and a heat flow equation that includes the Joule heating H J and recombination heat H R as heat sources.The dielectric permittivity is denoted by = 0 r , q is the elementary charge, C represents the density of the charged donors and acceptors, and is the thermal conductivity.
Additionally, we have to take into account the specialities of organic semiconductors: On the one hand, the statistical relation between chemical potentials and charge carrier densities is given by Gauss-Fermi integrals leading to bounded charge carrier densities.On the other hand, the mobility functions n , p depend on temperature, density and electrical field strength.The mobility laws are fitted from a numerical (1) solution of the master equation for the hopping transport in a disordered energy landscape with a Gaussian density of states [7,8].Moreover, the generation/recombination term R (see [9]), the Joule heat H J , the recombination heat H R and the charge carrier densities n and p in (1) are given by w h e r e k B i s B o l t z m a n n's c o n s t a n t a n d G ∶ ℝ × [0, ∞) → (0, 1) is defined by the Gauss-Fermi integral see [10].The LUMO and HOMO levels are denoted by E L , E H , respectively, and 2 n , 2 p are their variances, N n0 , N p0 rep- resent the total densities of transport states.We assume that these parameters are only weakly temperature dependent such that we neglect this weak temperature dependence in our investigations.In Sect.3, we use the abbreviation According to [7], the mobilities of electrons n = n (T, n, |∇ |) and holes p = p (T, p, |∇ |) are tem- perature, density and electrical field strength-dependent functions of the form with with the average hopping distance a.The system (1), ( 2) is closed by mixed boundary conditions on Ω for the station- ary drift-diffusion system combined with Robin boundary conditions for the heat flow equation approximating a heat sink with ambient temperature T a .
Here, Γ D and Γ N denote the Dirichlet and Neumann bound- ary parts, respectively, and denotes the outer unit normal vector on Ω .In [11], the solvability of (1), (2) and (5) (weak solutions of continuity equations and Poisson equation, entropy solution of the heat flow equation) is established, where the Dirichlet data for the potentials are given by H 1 (Ω) ∩ L ∞ (Ω)-functions.For mathematical results concerning the isothermal drift-diffusion model with Gauss-Fermi statistics and field strength-dependent mobility, we refer to [12,13].
We remark that in our model, as well as in [14] for classical semiconductors, additional thermoelectric effects (Peltier, Thomson and Seebeck) are neglected.In [6, Sect.II.D], it is argued that in the case of organic semiconductors, such effects are negligible as the thermal voltages are small compared to the applied voltage.In our simplified model frame, the recombination heat H R is given by the rate of entropy production due to recombination-generation processes, as in [15], multiplied by temperature.For fully thermodynamically consistent energy models for inorganic semiconductors including all these effects, we refer, e.g., to [15][16][17], where [17] also deals with the numerical approximation of the model.

Thermodynamic equilibrium and discussion of boundary conditions
In [11], the consistency of the model with thermodynamic equilibrium is proven.Thermodynamic equilibrium is a physical state with vanishing generation/recombination rate, and vanishing current and heat flow, where we can set 0 = 0 without loss of generality.In this situation, the system (1), ( 2) and ( 5) reduces to the nonlinear Poisson equation In the isothermal case following, e.g.[9], semiconductor-metal contacts, such as Ohmic contacts, are usually modeled by Dirichlet boundary conditions on where V i denotes the constant externally applied voltage at the ith contact.The potential 0i is determined a priori from the condition of local charge neutrality at the contact Γ D i with no applied voltage at ambient temperature T a : Keeping T a in (7) also in the non-isothermal case cor- responds to the unphysical assumption that the influence of temperature increase at the electrical contacts can be neglected for the determination of the equilibrium carrier densities.For the test structure described in the caption of Fig. 2, the simulation using (7) leads to extremely high carrier densities compared to the doping near to the contacts (top left) and to a pinning of the chemical potential v p = p − at both contacts (bottom left).We argue that in the non-isothermal case, the modeling of (ideal) Ohmic contacts requires local charge neutrality at the contact also at the actual temperature-dependent state ( , n , p , T) which leads to the nonlinear relation at the con- tacts Γ D i for the prescribed value 0i = − V i : As a result, the simulated hole densities in Fig. 2 (upper right) are not artificially increased near to the contacts.A straightforward generalization of the standard computational approach for the isothermal case would result in the necessity to update 0i for each modification of the temperature T, leading to an additional iterative loop for the determination of each bias solution.In order to avoid this iteration, we use (8) directly as a nonlinear Dirichlet boundary condition for the electrostatic potential depending on T and treat it with the nonlinear solver along with all other nonlinearities.Note that ( 7) and ( 8) coincide in the isothermal case.

Discretization scheme
We use a finite volume discretization method based on partitioning the computational domain Ω by a Voronoi mesh with m Voronoi cells {V l } l=1,…,m and accompanying collocation points {x l } .The potentials , n , p and the temperature T are evaluated at each node {x l } .The discretized system corresponding to (1) is derived by integrating the equations over each Voronoi cell V l , applying Gauss's theorem to get (9) and then approximating these integrals suitably.Here, N V l stands for the set of Voronoi cells V r which are adjacent to the Voronoi cell V l .We also add the subscript l in all quanti- ties such as potentials, doping density, recombination-generation rate, temperature and recombination heat to denote their corresponding numerical values at the node x l .In the following, we will assume that all material parameters such as the permittivity , the reference mobilities i0 , the densities of state N i0 and the heat conductivity are constant; otherwise, suitable averages have to be used.The surface integrals in ( 9) are split into two parts: integrals over interfaces between two adjacent Voronoi cells and integrals over boundary parts of the device, as follows where mes(K) denotes the measure of a set K. Meanwhile, the corresponding integrals in the continuity equations are approximated with some extra effort where the numerical fluxes J l;r n and J l;r p are determined by a modification of the Scharfetter-Gummel scheme based on averaging of inverse activity coefficients introduced in [18] and discussed with respect to degenerate semiconductors in [9,19].We introduce some notation to define the expressions for J l;r n and J l;r p in ( 12), (13).Let where n and p are defined in (3).Note that according to (4), the mobility over a surface V l ∩ V r depends on the modu- lus of the gradient |∇ | .The finite difference approximation behind the two point flux finite volume ansatz (10) gives only the normal component of the gradient with respect to V l ∩ V r and misses the tangential contribution, allow- ing for weak convergence at best if scaled with the space dimension [20].Therefore, we compute the approximation of |∇ | 2 on V l ∩ V r as the average squared norms of the gradients of the P1 finite element reconstruction over the set T l,r of all simplices (triangles in 2D) in the underlying Delaunay triangulation adjacent to the edge x l x r [21]: n;l ∶= n l , n;l , T l,r , n;r ∶= n r , n;r , T l,r , l,r n ∶= n l,r , n;l,r , T l,r , p;l ∶= p l , p;l , T l,r , p;r ∶= p r , p;r , T l,r , l,r p ∶= p l,r , p;l,r , T l,r , n ∶= n T l,r , n l,r , F l,r , l,r p ∶= p T l,r , p l,r , F l,r , The integrals over interfaces V l ∩ V r must be treated spe- cifically in order to maintain the consistency of the numerical solution (see Sect. 3.1), whereas the surface integrals over V l ∩ Ω are evaluated by quadrature rules after replac- ing the normal flux in the integrand by the corresponding boundary condition (see Sect. 3.2).

Numerical fluxes through interfaces between Voronoi cells @V l ∩ @V r
The integrals of − ∇ ⋅ and − ∇T ⋅ over the interface V r ∩ V l are approximated by the conventional finite differ- ence approximations 1 3 With the above definitions, the numerical fluxes J l;r n and J l;r p have the form Here, B denotes the Bernoulli function, B(x) = x exp (x)−1 .

Numerical treatment of the boundary conditions on @V l ∩ @Ä
The realization of no flux and Robin boundary conditions is based on the evaluation of the corresponding surface integrals by a midpoint quadrature rule.Among several possibilities to implement Dirichlet boundary conditions (e.g., direct elimination, matrix modification), we choose the penalty method for its flexibility and ease of implementation: We replace the Dirichlet boundary conditions for n , p by and treat them like Robin boundary conditions.The penalty parameter Π is a large number which results in marginalizing the normal flux contributions in (14) in the floating point representation, thus giving an accurate implementation of this boundary condition.
In order to approximate the nonlinear Dirichlet boundary condition (8), we use a similar idea.We replace (8) by and treat the resulting boundary condition as a nonlinear Robin boundary condition.Using this approach, the nonlinearity (8) can be treated without any additional iteration along with all the other nonlinearities in the resulting system of equations by the general Newton solver coupled to a parameter embedding scheme.

Volume integrals
The integrals of the charge density C − n + p , the recom- bination-generation rate R and the reaction heat H R are approximated by the midpoint rule, The integral of the Joule heat is approximated by using the numerical fluxes J l;r n and J l;r p , where dim(Ω) denotes the spatial dimension of Ω .Here, we followed the idea proposed in [22] and exploited in [21] allowing to evaluate the Joule heating approximation for electrons and holes by edge contributions.

Path-following technique for the calculation of S-shaped current-voltage curves
As discussed in the introduction, organic semiconductors show a pronounced electrothermal feedback that can lead to S-shaped current-voltage relations.For such situations, a voltage-controlled simulation is unable to cover the full characteristic, since at the lower turning point of the S-curve, one would not find a point on the curve with increased voltage and only slightly increased current and related temperature, see Fig. 1.For such voltage values, only points on the upper branch of the S-curve are J l;r n n;l − n;r + J l;r p p;l − p;r , available, related to very different current and temperature values.In other words, for increasing voltage, if at all the method would converge, one could only jump to the upper part of the S-curve and the (unstable) region of negative differential resistance of the S-curve is impossible to resolve.Formally, a current-controlled simulation would be able to establish this S-shaped relationship.Attempts to perform such simulations failed due to severe convergence problems of the solver (which possibly could be mitigated by a careful embedding strategy).However, there is a deeper reason to choose a voltage-controlled approach: It reasonably can be assumed that in the two-and threedimensional case, the contact voltage is constant at a given metal contact.However, this is not generally true for the current density, and one needs to implement an integral boundary condition for the current which would result in an additional layer of iterations.An alternative to this approach is the implementation of a path-following method to trace the S-curve under voltage-controlled conditions, which we describe in the sequel.We adapt the technique described in [21] which was used in [3] to simulate S-shaped current-voltage relations for organic LEDs resulting from an electrothermal modeling by p-Laplace thermistor models to the drift-diffusion setting.
We consider an organic semiconductor device with two Dirichlet boundary parts Γ D 1 and Γ D 2 .On Γ D 2 , the potential is set to zero and on Γ D 1 to the (spatially constant) externally applied voltage V. We determine the current-voltage relation of the organic device by calculating the current over Γ D 1 .With the equations for the approximation of (9) for all Voronoi cells {V l }, we arrive at a system of 4m coupled nonlinear algebraic equations for u = ( l , n;l , p;l , T l ) l=1,…,m of the form To trace a solution branch, starting from a solution (u 0 , V 0 ) of F(u, V) = 0 we use a predictor-corrector method [23] adapted to PDE calculations implemented in [24] as proposed in [25].For a recent description of this type of methods (based on a different implementation), see [26].The prediction is obtained by moving a step forward along the tangent vector t to the branch.First, we solve F u,V (u 0 , V 0 )t = 0 , t ∈ ℝ 4m+1 , where F u,V (u 0 , V 0 ) stands for the Jacobian of F at (u 0 , V 0 ) .To ensure that t points in the forward direction with respect to the tangent vector t 0 of the last point, we demand t 0 ⋅ t > 0 .In other words, we have to solve Next, we normalize t such that ‖t‖ = 1 .The predictor (u * , V * ) now is obtained by ensures that a step along the branch gives similar proportion to the unknowns and to the parameter, and, by construction, ‖u * − u 0 , V * − V 0 ‖ * = ΔL .The corrector step solves the nonlinear system by Newton's method with the starting value (u * , V * ) .If New- ton's method does not converge since the predictor (u * , V * ) is too far from the desired solution, the step size ΔL (related to the arc length parameter) is locally reduced until the method converges.The convergent Newton procedure yields the next point (u 1 , V 1 ) on the solution branch with a distance of ΔL to (u 0 , V 0 ).

Simulation results
The finite volume method has been implemented in the prototype semiconductor device simulator ddfermi [27] which is based on the PDE solution toolbox pdelib [24].
We provide results of a first investigation intended to be a proof of concept that electrothermal drift-diffusion models from Sect. 2 can exhibit S-shaped current-voltage relations.For simplicity, we restrict our simulations to the unipolar electronic case of a uniformly n-doped layer.Moreover, in the mobility function (4) we only take into account the temperature-dependent factor n0 (T) , describing the posi- tive feedback, but we ignore the density and field dependent factors g 1 and g 2 and set them to 1 to avoid parameter fitting.
We consider a uniformly n-doped, 340-nm-thick organic semiconductor material that is contacted by two metal layers and perform a quasi-1D simulation.Due to the high conductivity of the metal layers, we assume that the potential is The aim of the simulation work in this paper was to find a parameter range leading to a pronounced occurrence of S-shaped current-voltage relations.The resulting essential parameters for the simulation are collected in Table 1.
For the study of the phenomenon of S-shaped current-voltage relations and their appearance in dependence on physical parameters, we performed three types of parameter variations.
In Fig. 3, we investigated the influence of the disorder parameter n on the electrothermal interaction in the device.The resulting current-voltage relations are depicted in the left picture, and the right one shows the maximal device temperature over the applied voltage.We remark that due to the small size of the device, for a given bias voltage, the temperature difference within the device is very small.Whereas for n = 0.05 eV no S-shaped current-voltage relation with region of negative differential resistance occurs, such behavior evolves for higher n .With increasing n , the first turning point of the curve moves to a higher applied voltage, but the related current density decreases, and the 'S' becomes more pronounced.
In Fig. 4, we varied the reference mobility n0 .The resulting current-voltage relations are depicted in the left picture, For all presented calculations, temperatures above ≈ 400 K appear to be out of range for real devices, as the organic semiconductor material would be thermally destroyed.
The exemplary variations of physical parameters show that the complex nonlinear interplay leads to strong variations in the shape of the current-voltage characteristics as well as temperature-voltage relations.In this paper, we were especially looking for parameter constellations leading to complicated characteristics.The simulation of realistic organic semiconductor devices with more adequate physical parameters is planned in future work.

Conclusion and remarks
We presented a discretization scheme for the electrothermal drift-diffusion model (1) for organic semiconductor devices.We formulated temperature-dependent nonlinear Dirichlet boundary conditions for the electrostatic potential (8) at Ohmic contacts, which take into account the shift of the equilibrium potential due to changing device temperature.
We used a finite volume-based generalized Scharfetter-Gummel scheme implemented in the prototype semiconductor device simulator ddfermi [27] on top of the PDE solver toolbox pdelib [24].Implementing a path-following technique, we demonstrated that the model and its discretization for certain parameters exhibit the phenomenon of an S-shaped current-voltage relation with regions of negative differential resistance.To our knowledge, this is the first simulation tool for drift-diffusion-type models mastering this challenge.S-shaped current-voltage relations have been observed experimentally [3] and need to be properly modeled in order to understand and optimize the device behavior.
Besides the characteristics, our model and its discretization are capable to describe spatially resolved the electrothermal behavior of real 3D organic semiconductor devices in terms of charge carrier densities, current densities, potentials, temperature distributions.This in combination with a comparison to measurements has to be done in future work.As a first result, Fig. 6 stems from a 2D simulation for an organic thin-film transistor.It shows the produced Joule heat by using the numerical approximation of the electrothermal drift-diffusion model in Sect. 2 with the temperaturedependent Ohmic contact boundary conditions for the electrostatic potential (8).
In addition, future simulations by means of the model for real organic device structures and realistic physical parameters may help to estimate the region of a stable working regime guaranteeing the absence of material destruction due to overheating.Furthermore, the description of the spatially resolved electrothermal behavior of real devices is very important for understanding the effect of thermal switching, device degradation, device breakdown and local heating (hot spots) in large-area devices.Temperature activated ionic conductivity is observed in solid oxide materials like yttria-stabilized zirconia [28].It would be interesting to apply the methodology developed in this paper to recently derived drift-diffusion models for this type of materials [29].

Fig. 1
Fig. 1 Left: schematic cross section of the OLED stack simulated in [3].Right: simulated and measured S-shaped current-voltage relations with regions of negative differential resistance (NDR) for dif-

3 …
10 5 W m −2 K −1 a 1.5 nm E H 0.0 eV N n0 10 21 cm −3 n 0.05 … 0.08 eV Doping 5 × 10 18 cm −3 Thickness 340 nm constant here and neglect the metal layer entirely.Instead, we prescribe Dirichlet boundary conditions on the parts Γ D 1 and Γ D 2 .On Γ D 2 , the potential is set to zero and on Γ D 1 to the (spatially constant) externally applied voltage V. We determine the current-voltage relation of the organic device (which can be S-shaped) by calculating the current over Γ D 1 .

Fig. 3 Fig. 4
Fig. 3 Simulated current-voltage characteristics using the electrothermal drift-diffusion model for different disorder parameters n (left), the right figure shows the resulting maximal temperature in the device.We choose n0 = 0.8 cm /(Vs), 10 4 W/(m 2 K)