Numerical simulation of organic semiconductor devices with high carrier densities

Organic semiconducting materials have been the subject of much attention in recent years due to the prospect of tailoring their photoelectric properties to specific applications.^{1; 2; 3} Doping organic semiconductors, by embedding electron- or hole-donor species in the material, provides a particularly effective way of tuning device properties.⁴ Dopant concentrations can be easily and precisely controlled, although the effect on device performance is not always known a priori.

Simulation is a powerful tool for understanding and optimizing device performance, allowing systematic studies of device design that are experimentally infeasible. Computational approaches taken by other authors include Monte Carlo simulation,⁵ master equation approaches,⁶ and solution of the one-dimensional drift-diffusion-Poisson system of equations,⁷ which is the approach we follow here. This last approach is highly suited to the extraction of model parameters from experimental data and device optimization, since it is less computationally-expensive than methods dealing with more dimensions. One-dimensional drift-diffusion models originate from the theory of conventional crystalline semiconductors. However, when combining models for disordered materials with high doping intensities, numerical approaches originally developed for inorganic crystalline materials may fail. We derive a generalization of existing methods, namely the Scharfetter-Gummel discretization⁸ and the Gummel-iteration map⁹, to overcome these problems and then compare our results with a state-of-the-art simulation approach presented in Ref. ¹⁰. The numerical scheme that we develop here improves the numerical stability and computational efficiency of the approach.

The rest of this paper is structured as follows: In Sec. II, we introduce a model for steady-state charge transport in disorded semiconductors. In Sec. III, we describe a scaling scheme to improve numerical behavior. We present in Sec. IV the proper Scharfetter-Gummel discretization of the model, accounting for the generalized Einstein relation. In Sec. V, we derive the appropriately generalized Gummel iteration for the model. Numerical data for three example calculations are given in Sec. VI, comparing the proposed method with that of Ref. ¹⁰. Finally, we draw our conclusions in Sec. VII. The procedure is written for hole-transporting and bipolar devices in Appendix A, and the scaling factors we use are tabulated in Appendix B. Throughout this article, we use the symbol $e$ for the elementary charge.

The starting point for the derivations of most macroscopic charge transport models in semiconductors is the drift-diffusion-Poisson system introduced for inorganic crystalline semiconductors by van Roosbroeck,¹¹

Here, $J$ is the particle current, $\rho$ is the density, $\psi$ is the electrostatic potential, $C$ is the dopant density (assumed to be positive in case of $n$-type doping), $\mu$ is the electron mobility, and $D$ is the diffusion constant. The subscripts $n$ and $p$ refer to electrons and holes, respectively, and the total charge current is given by $-e(J_{n}-J_{p})$. To simplify the notation, we write the equations henceforth only for electrons using the electron density $n=\rho_{n}$. The following derivation is equally applicable for holes, since the equations in this case only differ by the sign of the drift term in the continuity equation and the carrier density in the Poisson equation. Appendix A gives the extension of our results to bipolar and hole-transport devices.

The model of Eqs. 1–4 was derived for crystalline silicon-based semiconductors¹² and has been applied in that area for many years with great success. Although there is no ab initio derivation, application of Eqs. 1–4 to organic semiconductors seems promising since the charge continuity and Poisson equations intuitively apply in some form. Similar models are used for organic electronics in, e.g., Refs. ¹⁰; ¹³.

In fact, Eqs. 1–4 must be extended in order to accurately describe disorded organic semiconductors.¹⁴ Available models include the extended Gaussian disorder model (EGDM) and the extended correlated disorder model (ECDM).^{15; 16} Both models assume a hopping mechanism for charge transport with a Gaussian density of states. The EGDM and ECDM introduce density- and field-dependence into the mobility and diffusion coefficients $\mu(T,\rho,\frac{\operatorname{d}\!\psi}{\operatorname{d}\!x})$ and $D(T,\rho,\frac{\operatorname{d}\!\psi}{\operatorname{d}\!x})$, respectively. The methods that we present in Secs. III, IV and V are applicable to arbitrary $\mu$ and $D$; here we discuss the EGDM and ECDM since they are recent and popular models.

It should be noted that in the context of so-called deep traps, introduced e.g. in Refs. ¹⁰; ¹⁷, the trapped carriers act as a stationary background charge and their treatment is completely analogous to the doping. Our method should therefore also lead to better performance in situations where deep traps are considered.

We restrict the discussion to one-dimensional steady-state current simulations. It is useful to introduce the functions $F_{1}$ and $F_{2}$, writing the steady-state system as

where $T$ is the temperature. When the system is accurately described by Maxwell-Boltzmann statistics, the classical Einstein relation holds,

with the thermal voltage $V_{\text{th}}=k_{B}T/e$. Note that in this equation the Boltzmann constant $k_{B}$ is taken to have units $\mathrm{JK}^{-1}$. In order to align notation with other authors in the field, we will frequently use $k_{B}$ in units of $\mathrm{eVK}^{-1}$. To still be able distinguish these two conventions we will use $k$ in this case.

When facing high charge-carrier densities, the effects of the Pauli exclusion principle become important and one is forced to resort to Fermi-Dirac statistics. In the context of organic semiconductors, it was first pointed out by Roichmann and Tessler¹⁸ that one should then use a generalized version of Eq. 7,

where $\varphi_{n}$ is the electron quasi-Fermi level (sometimes referred to as the chemical potential; we use the term quasi-Fermi level to avoid confusion and to stress that it is a non-equilibrium quantity). Within the EGDM/ECDM framework, Eq. 8 is usually written as

and $g_{3}$ is termed the diffusion enhancement. We adopt this convention for clarity and compatibility with other authors.

For a given density, the quasi-Fermi level is defined implicitly by the identity

where DOS($E$) represents the density of states.

Equation 8 is called the generalized Einstein relation¹⁹ and accounts for the fact that charge carriers near the quasi-Fermi level are more likely to contribute to diffusion. The derivative occurring in Eq. 8 can be computed analytically by differentiating Eq. 10. The differentiated expression can be found in e.g., the Appendix of Ref. ²⁰.

In SI units, the absolute values of the carrier density $[\mathrm{m}^{-3}]$ and electrostatic potential $[\mathrm{V}]$ can differ by over $20$ orders of magnitude. To avoid numerical problems we thus rescale the equations. Several scaling schemes have been proposed by other authors for both inorganic crystalline and organic devices.^{21; 22} We find that the scaling leading to the most stable behavior is state-dependent.

We suggest scaling the density by its maximum value (taking doping into account). We refer to the density scaling factor as $n_{\mathrm{scal}}$. Since we use a thermionic injection model ²³, which can be coupled with barrier lowering ²⁴, the boundary values and therefore maxima of the density are not known a priori. The scaling is therefore applied at every iteration. The potential is scaled by $\psi_{\mathrm{scal}}$, the maximum value of (the absolute value of) the applied voltage $V_{\mathrm{appl}}$ and the thermal voltage. This ensures that the scaled potential is reasonably small even for high applied voltages. On the other hand, we force the scaling parameter to be at least the thermal voltage to avoid scaling by zero. Since the mobility may differ strongly between materials, we scale all mobilities by their zero-field, zero-density limit $\mu_{0}$. Finally, space is scaled by $L$, the total thickness of the device. The full list of the scaling parameters including the resulting scaling of the current density is given in Table 2.

The scaled system of equations reads (using the same symbols for the scaled quantities)

where we have used the constants

Here, $\lambda$ is the scaled Debye length.

This scaling ensures that the potential and density vary between $0$ and $1$ for most physical situations. Failing to implement a suitable scaling scheme can result in numerical instability and loss of accuracy. Furthermore, it reduces the number of constants used in the actual computations. Except where explicitly stated, we henceforth use the scaled quantities.

To solve Eqs. 5 and 6 (note that these are the unscaled equations), Scharfetter and Gummel⁸ proposed a scheme for the case of constant mobility and diffusion. Instead of using central or upwind differences to approximate the current, a weighted difference depending on the field is used. It can be derived as a solution of the boundary value problem (BVP)

which describes the physics on a sub-interval for given carrier densities at the boundaries. The BVP of Eq. 14 is first order, but we can prescribe two boundary values because $J_{i+\frac{1}{2}}$ is a free parameter of the system, and there is only one value of $J_{i+\frac{1}{2}}$ that admits a solution. This value can be used as an approximation of the current, assuming that the mobility, diffusion coefficient and electric field are locally constant. Expressed in terms of $n$ and $\psi$, it reads

The scheme can be viewed as an upwind scheme, where the amount of upwind difference depends on the local strength of the field. This is easily seen by considering the asymptotic behavior of $J_{i+\frac{1}{2}}$ when the field is large or small,

so that for small fields we use a classic central difference approximation for the gradient of the density, i.e., the diffusive part of the current. When the magnitude of the field is large, we obtain

and

which are exactly upwind approximations for the density.

Upwind schemes introduce artificial diffusion into the solution,²⁵ which is a source of error. In the case of constant mobility and the simple Einstein relation of Eq. 7, it can be proven that the discretization introduced in Eq. 15 yields the optimal artificial diffusion (i.e., the minimum required for numerical stability). ²⁶

Within high-density capable models such as EGDM/ECDM, the generalized Einstein relation of Eq. 8 increases physical diffusion. Use of Eq. 15 to approximate the current is then suboptimal in terms of the accuracy of the solution because more artificial diffusion is added than is required for stability. Using the generalized Einstein relation Eq. 8 and applying the scaling introduced in Sec. III, we obtain the following generalization of the BVP of Eq. 14,

This corresponds to the current in Eq. 11.

Assuming $\mu$, $g_{3}$ and $\nabla\psi$ to be constant on each interval, or put equivalently, approximating them by their value at $i+\frac{1}{2}$, we solve Eq. 19 as before to get a properly-adjusted discretization scheme,

This scheme accounts for the additional physical diffusion in our system by decreasing the artificial diffusion, improving the accuracy of the current.

In the EGDM/ECDM, $\mu=\mu_{0}(T)g_{1}(n,T)g_{2}(\nabla\psi,T)$. In this case $g_{1}$ and $g_{3}$ at half-integer gridpoints can be obtained by various averaging methods. Throughout this work we use the geometric average

which is an exact interpolant for exponentially varying functions. Note that since $g_{2}$ depends on the gradient of the potential, it is already known on half-integer gridpoints only.

To solve the system of discretized equations, Eqs. 5 and 6, using the discretization of Sec. IV, one can either use the Newton algorithm or the so-called Gummel iteration, a physically-motivated variant of the Gauss-Seidel relaxation scheme.⁹

Numerical comparisons as in the work of Knapp et al.²⁷ and additional theoretical considerations show that Newton’s method converges in fewer iterations than the Gummel method. However, Newton iterations are generally more expensive to compute. Furthermore, Newton’s method often requires initial values in the vicinity of the solution in order to obtain convergence.

This is a severe problem for optimization and parameter extraction, where the solver is routinely called several hundred times. In these cases a single failure to converge can be catastrophic for the optimizer. We therefore derive an iterative solution scheme in the spirit of Gummel.⁹ As a starting point we will use the continuous (non-discretized) problem to avoid working in spaces of grid functions.

The idea is to alternately solve $F_{1}$ for $n$ at a given $\psi$ and $F_{2}$ for $\psi$ and given $n$ until self-consistency is achieved.

Since the drift-diffusion-Poisson system of equations can be nonlinear (even when $\mu$ and $D$ are constant; there is implicit dependency of e.g., $\psi$ on $n$), we must linearize first. Our iteration starts with the Poisson equation.

To highlight the advantages of the schemes described above, we have simulated a simple single-layer structure using the ECDM model. We consider three different situations: First, we simulate a low density situation without doping — here we show that our approach reproduces known results with a small gain in computational efficiency. Second, we simulate a device with high carrier densities due to a low injection barrier — here we observe a significant performance gain with the proposed method. Finally, we simulate a strongly-doped device, with approximately every $100$th molecule replaced by a dopant. In this third case, conventional methods fail but our procedure converges. The parameters used for the simulations are collected in Table 1. The resulting charge carrier densities and potentials at $3\mathrm{V}$ are shown in Fig. 2. All of the device simulations, including those using the conventional iteration, used the discretization described in Sec. IV and had identical initial values. Our calculations did not include the barrier-lowering effect of image-charges in the metal contacts. Its inclusion would have the effect of increasing the carrier injection, further enhancing the computational advantage of our method.

We have introduced a numerical method consisting of a scaling scheme, a discretization scheme, and a fixed point iteration for the van Roosbroeck drift-diffusion-Poisson system with arbitrary density- and field-dependent mobility functions. The method properly accounts for the implications of Fermi-Dirac statistics by incorporating the generalized Einstein relation. The iteration is derived by linearizing the Poisson equation with explicit treatment of the nonlinear dependence of the density on the potential.

Failure to take consequences of the generalized Einstein relation into account in the iteration results in numerical instability when the density is large, e.g., in the case of strong doping or large trap density, whereas the generalized method converges for arbitrary doping. In other cases, we find that the generalized method is significantly more efficient than the conventional method while reproducing converged current values.