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

Abstract

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.

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

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 different thermal outcoupling regimes realized by different values of the outcoupling coefficients. Both pictures are taken from [3]

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.

2 Electrothermal drift–diffusion description of organic semiconductor devices

2.1 The PDE system

The electrothermal behavior of organic semiconductor devices is described in a drift–diffusion setting by PDEs for the electrostatic potential $\psi$, the electrochemical potentials $\varphi _n$, $\varphi _p$ and the temperature T. In the device domain $\Omega,$ we consider the following stationary coupled system

$$\begin{aligned} \begin{aligned}&-\nabla \cdot (\varepsilon \nabla \psi ) = q(C-n+p),\\&-\nabla \cdot j_n = -qR, \quad j_n = - q n \mu _n \nabla \varphi _n, \\&\nabla \cdot j_p = -qR, \quad j_p = - q p \mu _p \nabla \varphi _p, \\&-\nabla \cdot ( \lambda \nabla T) = H_\mathrm{{J}} + H_\mathrm{{R}}. \end{aligned} \end{aligned}$$

(1)

This system results from the coupling of a generalized van Roosbroeck system and a heat flow equation that includes the Joule heating $H_\mathrm{{J}}$ and recombination heat $H_\mathrm{{R}}$ as heat sources. The dielectric permittivity is denoted by $\varepsilon =\varepsilon _0\varepsilon _r$, q is the elementary charge, C represents the density of the charged donors and acceptors, and $\lambda$ 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 $\mu _n,\,\mu _p$ depend on temperature, density and electrical field strength. The mobility laws are fitted from a numerical 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_\mathrm{{J}}$, the recombination heat $H_\mathrm{{R}}$ and the charge carrier densities n and p in (1) are given by

$$\begin{aligned} \begin{aligned} R&=r_0(\cdot ,n,p,T)\,n \, p \, \Big (1-\exp {\frac{q(\varphi _n- \varphi _p)}{k_{\mathrm {B}}T}}\Big ),\\ H_\mathrm{{J}}&=qn \mu _n |\nabla \varphi _n|^2 + qp \mu _p |\nabla \varphi _p|^2 =-j_n\cdot \nabla \varphi _n -j_p\cdot \nabla \varphi _p,\\ H_\mathrm{{R}}&= q R(\varphi _p-\varphi _n),\\ n&=N_{n0}(T)G\Big (\frac{q(\psi - \varphi _n) - E_\mathrm{{L}}(T)}{k_{\mathrm {B}}T};\frac{\sigma _n(T)}{k_{\mathrm {B}}T}\Big ), \\ p&=N_{p0}(T)G\Big (\frac{E_\mathrm{{H}}(T) - q( \psi - \varphi _p)}{k_{\mathrm {B}}T};\frac{\sigma _p(T)}{k_{\mathrm {B}}T}\Big ), \end{aligned} \end{aligned}$$

(2)

where $k_{\mathrm {B}}$ is Boltzmann’s constant and $G:{{\mathbb {R}}}\times [0,\infty )\rightarrow (0,1)$ is defined by the Gauss–Fermi integral

$$\begin{aligned} G(\eta ,s):= \frac{1}{\sqrt{2 \pi }} \int \limits _{-\infty }^{\infty } \exp \left( -\dfrac{\xi ^2}{2}\right) \dfrac{1}{\exp \left( s \xi - \eta \right) + 1} \,{\text {d}}\xi , \end{aligned}$$

see [10]. The LUMO and HOMO levels are denoted by $E_\mathrm{{L}}$, $E_\mathrm{{H}}$, respectively, and $\sigma _n^2$, $\sigma _p^2$ are their variances, $N_{n0}$, $N_{p0}$ represent 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

$$\begin{aligned} \eta _n= & {} \eta _n(\psi ,\varphi _n,T):= \frac{q(\psi -\varphi _n)-E_\mathrm{{L}}}{k_{\mathrm {B}}T},\quad \nonumber \\ \eta _p= & {} \eta _p(\psi ,\varphi _p,T):= \frac{E_\mathrm{{H}} - q( \psi - \varphi _p)}{k_{\mathrm {B}}T}. \end{aligned}$$

(3)

According to [7], the mobilities of electrons $\mu _n=\mu _n(T,n,|\nabla \psi |)$ and holes $\mu _p=\mu _p(T,p,|\nabla \psi |)$ are temperature, density and electrical field strength-dependent functions of the form

$$\begin{aligned} \mu _n(T,n,F)=\mu _{n0}(T)\times g_1(n,T)\times g_2(F,T) \end{aligned}$$

(4)

with

$$\begin{aligned} \begin{aligned} \mu _{n0}(T)&=\mu _{n0}c_1\exp \big \{{-}c_2s_n^2\big \},\quad s_n=\frac{\sigma _n}{k_{\mathrm {B}}T},\\ g_1(n,T)&=\exp \Big \{\frac{1}{2}( s_n^2- s_n)(2na^3)^\delta \Big \}, \quad \delta =2\frac{\ln ( s_n^2- s_n)-\ln (\ln 4)}{ s_n^2},\\ g_2(F,T)&=\exp \Big \{0.44( s_n^{3/2}-2.2)\Big ( \sqrt{1+0.8\Big (\frac{Fqa}{\sigma _n}\Big )^2}-1\Big )\Big \} \end{aligned} \end{aligned}$$

with the average hopping distance a. The system (1), (2) is closed by mixed boundary conditions on $\partial \Omega$ for the stationary drift–diffusion system combined with Robin boundary conditions for the heat flow equation approximating a heat sink with ambient temperature $T_\mathrm{{a}}$.

$$\begin{aligned} \begin{aligned}&\psi = \psi ^D,\quad \varphi _n= \varphi _n^D,\quad \varphi _p= \varphi _p^D \quad {\text {on}}\quad \Gamma _{\mathrm{{D}}}, \\&\varepsilon \nabla \psi \cdot \nu = j_n\cdot \nu = j_p\cdot \nu = 0\quad {\text {on}}\quad \Gamma _\mathrm{{N}}, \\&\lambda \nabla T\cdot \nu +\kappa (T- T_\mathrm{{a}})=0\quad {\text {on}} \quad \partial \Omega . \end{aligned} \end{aligned}$$

(5)

Here, $\Gamma _{\mathrm{{D}}}$ and $\Gamma _{\mathrm{{N}}}$ denote the Dirichlet and Neumann boundary parts, respectively, and $\nu$ denotes the outer unit normal vector on $\partial \Omega$. 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(\Omega )\cap L^\infty (\Omega )$-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_\mathrm{{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.

2.2 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,

$$\begin{aligned} j_n=j_p=0,\quad T=\text {const}=T_\mathrm{{a}}, \quad \varphi _0:=\varphi _n=\varphi _p=\text {const}, \end{aligned}$$

where we can set $\varphi _0=0$ without loss of generality. In this situation, the system (1), (2) and (5) reduces to the nonlinear Poisson equation

$$\begin{aligned} -\nabla \cdot (\varepsilon \nabla \psi )= & {} q\Big (C- N_{n0}G\Big (\frac{q\psi - E_\mathrm{{L}}}{k_{\mathrm {B}}T_\mathrm{{a}}};\frac{\sigma _n}{k_{\mathrm {B}}T_\mathrm{{a}}}\Big )\\&+N_{p0}G\Big (\frac{E_\mathrm{{H}} - q\psi }{k_{\mathrm {B}}T_\mathrm{{a}}};\frac{\sigma _p}{k_{\mathrm {B}}T_\mathrm{{a}}}\Big )\Big ). \end{aligned}$$

In the isothermal case following, e.g. [9], semiconductor–metal contacts, such as Ohmic contacts, are usually modeled by Dirichlet boundary conditions on $\Gamma _\mathrm{{D}}=\cup _{i=1}^I \Gamma _{\mathrm{{D}}_i}$,

$$\begin{aligned} \psi =\psi _{0i}+V_i,\quad \varphi _n=V_i,\quad \varphi _p=V_i\quad {\text {on}}\quad \Gamma _{\mathrm{{D}}_i}, \end{aligned}$$

(6)

where $V_i$ denotes the constant externally applied voltage at the ith contact. The potential $\psi _{0i}$ is determined a priori from the condition of local charge neutrality at the contact $\Gamma _{\mathrm{{D}}_i}$ with no applied voltage at ambient temperature $T_\mathrm{{a}}$:

$$\begin{aligned}&C- N_{n0}G\Big (\frac{q\psi _{0i} - E_\mathrm{{L}}}{k_{\mathrm {B}}T_\mathrm{{a}}};\frac{\sigma _n}{k_{\mathrm {B}}T_\mathrm{{a}}}\Big )\nonumber \\&\quad +N_{p0}G\Big (\frac{E_\mathrm{{H}} - q\psi _{0i}}{k_{\mathrm {B}}T_\mathrm{{a}}};\frac{\sigma _p}{k_{\mathrm {B}}T_\mathrm{{a}}}\Big )=0. \end{aligned}$$

(7)

Keeping $T_\mathrm{{a}}$ in (7) also in the non-isothermal case corresponds 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=\varphi _p-\psi$ at both contacts (bottom left).

Fig. 2

Simulation of a 200-nm-thick structure of three (50 nm, 100 nm, 50 nm) pip-doped layers for different applied voltages. The p-doped layers are doped by $10^{18}\,{\text {cm}}^{-3}$, the intrinsic region by $10^{10}\,{\text {cm}}^{-3}$. Upper row: hole density. Lower row: chemical potential $v_p=\varphi _p-\psi$. The plots in the left column have been obtained with the boundary conditions resulting from (7) and show unphysical increased hole densities in the p-doped regions. This is avoided by the contact description using (8) in the right column

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 $(\psi ,\varphi _n,\varphi _p,T)$ which leads to the nonlinear relation at the contacts $\Gamma _{\mathrm{{D}}_i}$ for the prescribed value $\psi _{0i}=\psi -V_i$:

$$\begin{aligned} C_{\mathrm{{D}}_i}(\psi ; V_i, T):= \,& {} C- N_{n0}G\Big (\frac{q(\psi - V_i)-E_\mathrm{{L}}}{k_{\mathrm {B}}T};\frac{\sigma _n}{k_{\mathrm {B}}T}\Big )\nonumber \\&+\,N_{p0}G\Big (\frac{E_\mathrm{{H}} -q(\psi -V_i)}{k_{\mathrm {B}}T};\frac{\sigma _p}{k_{\mathrm {B}}T}\Big )=0. \end{aligned}$$

(8)

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 $\psi _{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 $\psi$ 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.

3 Discretization scheme

We use a finite volume discretization method based on partitioning the computational domain $\Omega$ by a Voronoi mesh with m Voronoi cells $\{V_l\}_{l=1,\dots ,m}$ and accompanying collocation points $\{x_l\}$. The potentials $\psi$, $\varphi _n$, $\varphi _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

$$\begin{aligned} \begin{aligned} \int \limits _{\partial V_l} -\varepsilon \nabla \psi \cdot \nu \,{\text {d}}\Gamma&= \int \limits _{V_l} q \left( C - n + p\right) \,{\text {d}}x, \\ \int \limits _{\partial V_l} -j_n \cdot \nu \,{\text {d}}\Gamma&= \int \limits _{V_l} -qR \,{\text {d}}x, \\ \int \limits _{\partial V_l} j_p \cdot \nu \,{\text {d}}\Gamma&= \int \limits _{V_l} -qR\,{\text {d}}x, \\ \int \limits _{\partial V_l} -\lambda \nabla T \cdot \nu \,\text {d}\Gamma&= \int \limits _{V_l} (H_\mathrm{{J}}+H_\mathrm{{R}})\,{\text {d}}x \end{aligned} \end{aligned}$$

(9)

and then approximating these integrals suitably. Here, $\mathcal {N}\left( V_l\right)$ 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 quantities 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 $\varepsilon$, the reference mobilities $\mu _{i0}$, the densities of state $N_{i0}$ and the heat conductivity $\lambda$ 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

$$\begin{aligned} \begin{aligned} \int \limits _{\partial V_l} -\varepsilon \nabla \psi \cdot \nu \,{\text {d}}\Gamma&= \sum _{V_r \in \mathcal {N}\left( V_l\right) } \int \limits _{\partial V_l \cap \partial V_r}-\varepsilon \nabla \psi \cdot \nu \,{\text {d}}\Gamma + \int \limits _{\partial V_l \cap \partial \Omega } -\varepsilon \nabla \psi \cdot \nu \,{\text {d}}\Gamma , \\ \int \limits _{\partial V_l} -j_n \cdot \nu \,{\text {d}}\Gamma&= \sum _{V_r \in \mathcal {N}\left( V_l\right) } \int \limits _{\partial V_l \cap \partial V_r} -j_n \cdot \nu \,{\text {d}}\Gamma + \int \limits _{\partial V_l \cap \partial \Omega } -j_n \cdot \nu \,{\text {d}}\Gamma , \\ \int \limits _{\partial V_l} j_p \cdot \nu \,{\text {d}}\Gamma&= \sum _{V_r \in \mathcal {N}\left( V_l\right) } \int \limits _{\partial V_l \cap \partial V_r} j_p \cdot \nu \,{\text {d}}\Gamma + \int \limits _{\partial V_l \cap \partial \Omega } j_p \cdot \nu \,{\text {d}}\Gamma , \\ \int \limits _{\partial V_l} -\lambda \nabla T \cdot \nu \,{\text {d}}\Gamma&= \sum _{V_r \in \mathcal {N}\left( V_l\right) } \int \limits _{\partial V_l \cap \partial V_r}-\lambda \nabla T \cdot \nu \,{\text {d}}\Gamma + \int \limits _{\partial V_l \cap \partial \Omega } -\lambda \nabla T \cdot \nu \,{\text {d}}\Gamma . \end{aligned} \end{aligned}$$

The integrals over interfaces $\partial V_l \cap \partial V_r$ must be treated specifically in order to maintain the consistency of the numerical solution (see Sect. 3.1), whereas the surface integrals over $\partial V_l \cap \partial \Omega$ are evaluated by quadrature rules after replacing the normal flux in the integrand by the corresponding boundary condition (see Sect. 3.2).

3.1 Numerical fluxes through interfaces between Voronoi cells ${\varvec{\partial V_l \cap \partial V_r}}$

The integrals of $-\varepsilon \nabla \psi \cdot \nu$ and $-\lambda \nabla T \cdot \nu$ over the interface $\partial V_r \cap \partial V_l$ are approximated by the conventional finite difference approximations

$$\begin{aligned} \int \limits _{\partial V_r \cap \partial V_l} -\varepsilon \nabla \psi \cdot \nu \,{\text {d}}\Gamma&\approx \frac{{\mathrm {mes}}\left( \partial V_r \cap \partial V_l\right) }{\left| x_l - x_r\right| } \varepsilon \left( \psi _l - \psi _r\right) , \\ \int \limits _{\partial V_r \cap \partial V_l} -\lambda \nabla T \cdot \nu \,{\text {d}}\Gamma&\approx \frac{{\mathrm {mes}}\left( \partial V_r \cap \partial V_l\right) }{\left| x_l - x_r\right| } \lambda \left( T_l - T_r\right) , \nonumber \end{aligned}$$

(10)

where ${\mathrm {mes}}(K)$ denotes the measure of a set K. Meanwhile, the corresponding integrals in the continuity equations are approximated with some extra effort

$$\begin{aligned} \begin{aligned} \int \limits _{\partial V_r \cap \partial V_l} j_n \cdot \nu \,{\text {d}}\Gamma&\approx \frac{{\mathrm {mes}}\left( \partial V_r \cap \partial V_l\right) }{\left| x_l - x_r\right| } J_n^{l; r},\\ \int \limits _{\partial V_r \cap \partial V_l} j_p \cdot \nu \,{\text {d}}\Gamma&\approx \frac{{\mathrm {mes}}\left( \partial V_r \cap \partial V_l\right) }{\left| x_l - x_r\right| } J_p^{l; r}, \end{aligned} \end{aligned}$$

where the numerical fluxes $J_n^{l; r}$ and $J_p^{l; r}$ 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_n^{l; r}$ and $J_p^{l; r}$ in (12), (13). Let

$$\begin{aligned} \psi _{l,r}:= & {} \frac{\psi _l + \psi _r}{2}, \ \varphi _{n; l,r} := \frac{\varphi _{n;l} + \varphi _{n;r}}{2}, \ \varphi _{p; l,r} := \frac{\varphi _{p;l} + \varphi _{p;r}}{2}, \ T_{l,r} := \frac{T_l + T_r}{2}, \\ \overline{\eta }_{n;l}:= & {} \eta _n\left( \psi _l, \varphi _{n;l}, T_{l,r}\right) , \ \overline{\eta }_{n;r} := \eta _n\left( \psi _r, \varphi _{n;r}, T_{l,r}\right) , \ \overline{\eta }_n^{l,r} := \eta _n\left( \psi _{l,r}, \varphi _{n; l,r}, T_{l,r}\right) , \\ \overline{\eta }_{p;l}:= & {} \eta _p\left( \psi _l, \varphi _{p;l}, T_{l,r}\right) , \ \overline{\eta }_{p;r} := \eta _p\left( \psi _r, \varphi _{p;r}, T_{l,r}\right) , \ \overline{\eta }_p^{l,r} := \eta _p\left( \psi _{l,r}, \varphi _{p; l,r}, T_{l,r}\right) , \\ U_T^{l,r}:= & {} \frac{k_{\mathrm {B}} T_{l,r}}{q}, \ s_n^{l,r} := \frac{\sigma _n}{k_{\mathrm {B}} T_{l,r}}, \ s_p^{l,r} := \frac{\sigma _p}{k_{\mathrm {B}} T_{l,r}}, \\ n^{l,r}:= & {} N_{n0} G\left( \overline{\eta }_n^{l, r}; s_n^{l,r}\right) ,\ p^{l,r} := N_{p0} G\left( \overline{\eta }_p^{l, r}; s_p^{l,r}\right) ,\\ \mu _n^{l,r}:= & {} \mu _n\left( T_{l,r}, n^{l,r}, F^{l,r}\right) , \ \mu _p^{l,r} := \mu _p\left( T_{l,r}, p^{l,r}, F^{l,r}\right) , \end{aligned}$$

where $\eta _n$ and $\eta _p$ are defined in (3). Note that according to (4), the mobility over a surface $\partial V_l\cap \partial V_r$ depends on the modulus of the gradient $|\nabla \psi |$. The finite difference approximation behind the two point flux finite volume ansatz (10) gives only the normal component of the gradient with respect to $\partial V_l\cap \partial V_r$ and misses the tangential contribution, allowing for weak convergence at best if scaled with the space dimension [20]. Therefore, we compute the approximation of $|\nabla \psi |^{2}$ on $\partial V_l\cap \partial V_r$ as the average squared norms of the gradients of the P1 finite element reconstruction $\psi _\tau$ over the set $\mathcal T_{l,r}$ of all simplices $\tau$ (triangles in 2D) in the underlying Delaunay triangulation adjacent to the edge $\overline{x_lx_r}$ [21]:

$$\begin{aligned} \big |\nabla \psi \big |^2 |_{\partial V_l\cap \partial V_r} \approx \frac{\sum _{\tau \in {\mathcal {T}}_{l,r}} {\mathrm {mes}}(\tau ) |\nabla \psi _\tau |^2}{\sum _{\tau \in {\mathcal {T}}_{l,r}} {\mathrm {mes}}(\tau )}=:(F^{l,r})^2. \end{aligned}$$

(11)

With the above definitions, the numerical fluxes $J_n^{l;r}$ and $J_p^{l;r}$ have the form

$$\begin{aligned} J_n^{l;r}= & {} - q N_{n0} \mu _n^{l,r} U_T^{l,r} \frac{G\left( \overline{\eta }_n^{l,r}; s_n^{l,r}\right) }{\exp \left( \overline{\eta }_n^{l,r}\right) } \left[ \exp \left( \overline{\eta }_{n; l} \right) B\left( \frac{\psi _l - \psi _r}{U_T^{l,r}}\right) - \exp \left( \overline{\eta }_{n; r}\right) B\left( -\frac{\psi _l - \psi _r}{U_T^{l,r}}\right) \right] , \end{aligned}$$

(12)

$$\begin{aligned} J_p^{l;r}= & {} q N_{p0} \mu _p^{l,r} U_T^{l,r} \frac{G\left( \overline{\eta }_p^{l,r}; s_p^{l,r}\right) }{\exp \left( \overline{\eta }_p^{l,r}\right) } \left[ \exp \left( \overline{\eta }_{p; l} \right) B\left( -\frac{\psi _l - \psi _r}{U_T^{l,r}}\right) - \exp \left( \overline{\eta }_{p; r}\right) B\left( \frac{\psi _l - \psi _r}{U_T^{l,r}}\right) \right] . \end{aligned}$$

(13)

Here, B denotes the Bernoulli function, $B\left( x\right) = \frac{x}{\exp \left( x\right) - 1}$.

3.2 Numerical treatment of the boundary conditions on ${\varvec{\partial V_l \cap \partial \Omega }}$

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 $\varphi _n, \varphi _p$ by

$$\begin{aligned} j_n\cdot \nu + \Pi (\varphi _n -V) =0, \quad j_p\cdot \nu + \Pi (\varphi _p -V)=0 \end{aligned}$$

(14)

and treat them like Robin boundary conditions. The penalty parameter $\Pi$ 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

$$\begin{aligned} -\varepsilon \nabla \psi \cdot \nu +\Pi C_{\mathrm{{D}}_i}(\psi ; V, T)=0 \end{aligned}$$

(15)

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.

3.3 Volume integrals

The integrals of the charge density $C - n + p$, the recombination–generation rate R and the reaction heat $H_\mathrm{{R}}$ are approximated by the midpoint rule,

$$\begin{aligned} \begin{aligned} \int \limits _{V_l} (C - n + p) \,{\text {d}}x&\approx {\mathrm {mes}}\left( V_l\right) (C_l - n_l + p_l), \\ \int \limits _{V_l} R \,{\text {d}}x&\approx {\mathrm {mes}}\left( V_l\right) R_l, \\ \int \limits _{V_l} H_\mathrm{{R}} \,{\text {d}}x&\approx {\mathrm {mes}}\left( V_l\right) R_l(\varphi _{p;l}-\varphi _{n;l}). \end{aligned} \end{aligned}$$

The integral of the Joule heat is approximated by using the numerical fluxes $J_n^{l;r}$ and $J_p^{l;r}$,

$$\begin{aligned} \int \limits _{V_l} H_\mathrm{{J}}\,{\text {d}}x \approx \sum \limits _{V_r \in \mathcal {N}\left( V_l\right) } \dfrac{{\mathrm {mes}}\left( \partial V_l \cap \partial V_r\right) }{2{\mathrm {dim}}\left( \Omega \right) } \Big ( J_n^{l;r} \left( \varphi _{n;l} - \varphi _{n;r}\right) +J_p^{l;r} \left( \varphi _{p;l} - \varphi _{p;r}\right) \Big ), \end{aligned}$$

(16)

where ${\mathrm {dim}}(\Omega )$ denotes the spatial dimension of $\Omega$. 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.

3.4 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 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 three-dimensional 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 $\Gamma _{\mathrm{{D}}_1}$ and $\Gamma _{\mathrm{{D}}_2}$. On $\Gamma _{\mathrm{{D}}_2},$ the potential is set to zero and on $\Gamma _{\mathrm{{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 $\Gamma _{\mathrm{{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=(\psi _l,\varphi _{n;l},\varphi _{p;l},T_l)_{l=1,\dots ,m}$ of the form

$$\begin{aligned} F(u,V)=0, \quad F:{{\mathbb {R}}}^{4m}\times {{\mathbb {R}}}\rightarrow {{\mathbb {R}}}^{4m}. \end{aligned}$$

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\in {{\mathbb {R}}}^{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\cdot t>0$. In other words, we have to solve

$$\begin{aligned} \begin{pmatrix} F_{u,V}(u_0,V_0)\\ t_0^\top \end{pmatrix}t =\begin{pmatrix} 0\\ 1 \end{pmatrix}. \end{aligned}$$

Next, we normalize t such that $\Vert t\Vert =1$. The predictor $(u^*,V^*)$ now is obtained by

$$\begin{aligned} \begin{pmatrix} u^*\\ V^* \end{pmatrix} =\begin{pmatrix} u_0\\ V_0 \end{pmatrix} +\frac{\Delta L}{\Vert t\Vert _*}t,\quad {\text {where}} \Vert t\Vert _*^2=\frac{1}{4m}\sum _{i=1}^{4m}t_i^2 + t_{4m+1}^2 \end{aligned}$$

ensures that a step along the branch gives similar proportion to the unknowns and to the parameter, and, by construction, $\Vert u^*-u_0,V^*-V_0\Vert _*=\Delta L$. The corrector step solves the nonlinear system

$$\begin{aligned} \begin{pmatrix} F(u,V)\\ \Vert u-u_0,V-V_0\Vert ^2_*-(\Delta L)^2 \end{pmatrix}=0 \end{aligned}$$

by Newton’s method with the starting value $(u^*,V^*)$. If Newton’s method does not converge since the predictor $(u^*,V^*)$ is too far from the desired solution, the step size $\Delta 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 $\Delta L$ to $(u_0,V_0)$.

4 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 $\mu _{n0}(T)$, describing the positive 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 constant here and neglect the metal layer entirely. Instead, we prescribe Dirichlet boundary conditions on the parts $\Gamma _{\mathrm{{D}}_1}$ and $\Gamma _{\mathrm{{D}}_2}$. On $\Gamma _{\mathrm{{D}}_2},$ the potential is set to zero and on $\Gamma _{\mathrm{{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 $\Gamma _{\mathrm{{D}}_1}$. 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.

Table 1 Simulation parameters

Fig. 3

Simulated current–voltage characteristics using the electrothermal drift–diffusion model for different disorder parameters $\sigma _n$ (left), the right figure shows the resulting maximal temperature in the device. We choose $\mu _{n0}=0.8$ cm$^2$/(Vs), $\kappa =10^4$ W/(m$^2$ K)

Fig. 4

Simulated current–voltage characteristics using the electrothermal drift–diffusion model for different reference mobilities $\mu _{n0}$ and $\kappa =10^4$ W/(m$^2$ K), $\sigma _n=0.08$ eV. The right figure shows the maximal temperature in the device

Fig. 5

Simulated current–voltage characteristics (left) and maximum temperature (right) using the electrothermal drift–diffusion model for different thermal outcoupling factors $\kappa$, $\mu _{n0}=0.8$ cm$^2$/(Vs) and $\sigma _n=0.08$ eV

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 $\sigma _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 $\sigma _n=0.05$ eV no S-shaped current–voltage relation with region of negative differential resistance occurs, such behavior evolves for higher $\sigma _n$. With increasing $\sigma _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 $\mu _{n0}$. The resulting current–voltage relations are depicted in the left picture, and the right one shows the maximal device temperature over the applied voltage. Lower mobilities enforce a more pronounced S-shape which is additionally shifted to the right. Higher mobilities lead to a stronger increase in the current density.

Finally, Fig. 5 contains a variation of the thermal outcoupling factor $\kappa$ corresponding to the Robin boundary conditions $\lambda \nabla T\cdot \nu +\kappa (T-T_\mathrm{{a}})=0$ on $\partial \Omega$. (left: current–voltage relations, right: the maximal device temperature over the applied voltage). Better cooling broadens the ‘S’; for the two turning points, the applied voltage as well as the current density increase.

For all presented calculations, temperatures above $\approx$ 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.

5 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 temperature-dependent Ohmic contact boundary conditions for the electrostatic potential (8).

Fig. 6

Simulated Joule heat density [W/cm$^{3}$] in an organic thin-film transistor using the numerical approximation of the electrothermal drift–diffusion model (1)

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