Coupled 3D master equation and 1D drift‐diffusion approach for advanced OLED modeling

A novel simulation approach for excitonic organic light‐emitting diodes (OLEDs) is established by combining a continuous one‐dimensional (1D) drift‐diffusion (DD) model for the charge carrier dynamics with a three‐dimensional (3D) master equation (ME) model describing the exciton dynamics in a multilayer OLED stack with an additional coupling to a thin‐film optics solver. This approach effectively combines the computational efficiency of the 1D DD solver with the physical accuracy of a discrete 3D ME model, where excitonic long‐range interactions for energy transfer can be taken into account. The coupling is established through different possible charge recombination types as well as the carrier densities themselves. We show that such a hybrid approach can efficiently and accurately describe steady‐state and transient behavior of optoelectronic devices reported in literature. Such a tool will facilitate the optimization and characterization of multilayer OLEDs and other organic semiconductor devices.


| INTRODUCTION
A lot of attention has been drawn by organic semiconducting materials and devices in the recent years as they enable disruptive lighting and flat panel displays and even flexible electronics applications with superior performance. [1][2][3][4] The exact mechanisms and interactions between charge carriers, excitons, and photons in these materials are complex, and no analytical solutions describing the dynamics are available in general; therefore, numerical methods have to be applied to investigate the processes in organic light-emitting diodes (OLEDs) and to support the development process.
Two widely used modeling frameworks are continuum approaches such as deterministic drift-diffusion (DD) models and atomistic approaches like the probabilistic kinetic Monte Carlo (kMC) method. A master equation (ME) model is situated in between these two approaches as it is atomistic like kMC yet deterministic like DD models, describing the evolution of the states' probability mass function over time. The computational cost of solving such a system is consequently also to be found in between those of the highly efficient DD and the computationally highly demanding kMC approaches. 5 Although there exist advanced DD models for the charge carrier dynamics, excitons are generally described by simple local reaction-diffusion models. An ME approach to describe the exciton dynamics provides certain advantages in the physical accuracy over the former as transport can be accurately described using established expressions for Förster and Dexter energy transfer rates where long-range interactions are taken into account. Furthermore, the atomistic nature of ME models renders the modeling of complex host-guest systems as well as multilayer structures with various interfaces straightforward.

| MODEL
The model presented consists of two parts: a 1D DD model describing the electronic part, that is, the charge carrier dynamics of electrons and holes as well as the coupled electric potential distribution inside the device. The second part consists of a 3D ME model describing the exciton dynamics. These two models are solved selfconsistently, where the two are coupled through the carrier recombination rate, carrier density for polaron quenching, and the potential in the device.
The DD model 6,7 for the carrier dynamics consists of the two charge continuity equations for electrons and holes as well as the Poisson equation for the potential (see Equations 1a-1c). Furthermore, the electron and hole currents are modeled as the sum of a drift and a diffusion part as shown in Equations 1d and 1e.
∂n ∂t −r Á J n = −R n, p ð Þ−R nt n, n t ð Þ, ð1bÞ J n ≔nμ n rψ −D n rn, ð1dÞ The variables to be found are the potential ψ as well as n and p, the electron and hole densities, respectively. Furthermore, C (m −3 ) denotes the density of immobile charges due to doping and trap states, R (m −3 s −1 ) the carrier recombination rate, R nt,pt (m −3 s −1 ) the net trapping rates for electrons and holes, J n,p (m −2 s −1 ) the electron and hole current fluxes, and μ n,p (m 2 V −1 s −1 ) the electron and hole mobilities. D n,p are the electron and hole diffusion constants, related to the mobilities through a generalized Einstein relation. 8 The ME model was adapted from Zhou et al. 9 It is assumed that the excitons in the material are strongly localized Frenkel excitons due to the strong Coulomb interaction of the electron-hole pair (low permittivity) in organic semiconductors. This allows the definition of occupation numbers for each molecule in the material matrix. Furthermore, as relaxation from higher excited states to the lowest excited state happens on timescales much faster than other processes, 1 0 only the first singlet and triplet states S 1 and T 1 , respectively, are considered in the model. Using these assumptions, the ordinary differential equations (ODEs) describing the local change in occupation numbers for the two exciton species can be formulated as follows 9 : where p s,t i (−) denote local singlet and triplet occupation numbers, respectively, ω s,t,F,D ij (s −1 ) the Förster and Dexter transfer rates from site i to site j, and τ i (s) the local characteristic time constants for the different mechanisms. E b, i and E g,i (J) denote the local exciton binding energies and LUMO-HOMO gap energies, respectively, where LUMO denotes the lowest unoccupied molecular orbital whereas HOMO denotes the highest occupied molecular orbital. Furthermore, the factors g rec,i , g opt,i (−) describe the efficiencies for the singlet generation by charge carrier recombination and optical absorption, respectively.
In both Equations 2a and 2b the first line represents Förster and Dexter transfer, the second line radiative and nonradiative decay, the third line (reverse) intersystem crossing, the fourth line polaron quenching and triplettriplet annihilation, and the last line generation due to carrier recombination and photon absorption. The exponentials in the (non-)radiative decay rates are included because of the detailed balance considerations of these processes. It is assumed the lifetimes for the intersystem crossing processes are proportional with a Boltzmann factor Furthermore, it is approximately assumed that triplet-triplet annihilation results in generation of a single singlet, that is, T 1 +T 1 ! S 1 . In Equations 2a and 2b the quantities n,p,R, and G are directly taken as computed by the DD model. As they represent spatial densities coming from a continuum formalism, they need to be converted to the discrete ME formalism by multiplying with the unit cell volume a 3 0 . Depending on the type of molecule (host or guest), the recombination term can be taken differently, for example, from Langevin recombination for host molecules and trap-trap recombination for guest molecules. This leads to spatially nonuniform generation throughout the yz cross sections. A similar procedure holds for the charge carrier densities, where for the guest molecules the corresponding trap densities from the DD model are taken. Note that as an approximation, polaron quenching and annihilation processes are modeled purely locally. Alternatively, these processes can also be modeled using long-range Förster type interactions. [11][12][13] Coupling to a thin film optics simulation furthermore provides the quantity F i , which is defined as known as Purcell factor, which describes the radiative lifetime modulation of dipoles due to microcavity effects, which results in a spatial dependency of said lifetime. This factor depends on the intrinsic photoluminescence (PL) spectrum of the dipole and can therefore vary between materials and exciton type. The calculated radiative decay events from Equations 2a and 2b are then again used as input to a thin film optics emission model in order to calculate (among others) the outcoupled photons and finally the luminance and spectrum of the device.

| Auxiliary models
The transfer rates ω F,D are modeled using the canonical expressions for Förster and Dexter energy transfer, namely, 9 where R F 0 denotes the Förster radius, that is, the donoracceptor separation where hopping due to Förster transfer is equally likely as spontaneous decay on the donor in absence of the acceptor molecule. ΔE et ij denotes the difference in exciton energies between donor and acceptor. The function f penalizes the jump from lower to higher energies by introducing The Förster radius R F 0 is taken constant as an approximation; however, the energy penalty term f ΔE et ij strongly prohibits transfers from guest to matrix molecules while it favors the reverse process. Alternatively, R F 0 can be taken differently according to the donor-acceptor material combination involved in the transfer process. Similarly, the Dexter transfer rates are modeled as 9 k D 0 the Dexter hopping rate in the zero-distance limit and R D 0 the spatial extension of the electron wave function. As the absolute LUMO and HOMO levels appear in these expressions, a weak electric field dependency emerges for the Dexter hopping rates.
The HOMO and LUMO levels of the different molecules are taken from the DD model and consist of intrinsic material HOMO and LUMO levels superposed with the contribution of the electric potential. Furthermore, there is a random contribution added in the ME formalism. The random energy contribution is used to model energetic disorder in the organic material and can be either directly sampled from a Gaussian density of states (DOS): in which case the disorder is spatially uncorrelated or by using random dipole interactions 14 resulting in a spatially correlated Gaussian distribution with a standard deviation given by 14 where d is the magnitude of the dipole moment which is fixed for a specific material.

| RESULTS AND DISCUSSION
In this section, three application examples of the hybrid DD and master equation approach are presented, which illustrate the feasibility of this approach for comprehensive electronic, excitonic, and optical (OLED) device simulation.

| Steady-state electroluminescence
We simulate with the present hybrid model the steadystate behavior of a green OLED structure illustrated in Figure 1, as analyzed previously by Zhou et al. 9 using a full ME approach (i.e., charge carriers and excitons modeled with an ME model). For the DD model, the host-guest complex in the emitting layer (EML) is modeled as a layer with electron traps using trap state densities and energy depths corresponding to the dye dopants. In the electron and hole transport layers (ETL and HTL, respectively) the recombination rate is given by the expression for Langevin recombination with n i being the intrinsic carrier concentration in thermal equilibrium. In the EML, it is assumed that hole transport only happens over the Ir(ppy) 3 dopants (hole transport level in the layer is given by the HOMO of these dopants), and therefore, only recombination over these trap states is possible. For this reason, the recombination in such layer was modeled using the Shockley-Read-Hall (SRH) expression with E t the trap energy level, N t , N 0,n , and N 0,p the trap, mobile LUMO and mobile HOMO states densities, respectively, and c n and c p electron and hole capture rates. The carrier densities and potential and recombination rate calculated by the 1D DD model are then interpreted as layer averaged quantities and distributed over the 3D grid accordingly for the ME model, where in the EML, R SRH was used as recombination term on the dye molecules (i.e., no recombination on matrix molecules). The resulting electric field distributions and charge carrier profiles for both cases as defined in Figure 1 are given in Figure S1 and Figure S2 of the supplementary information. Figure S2 also shows a comparison of the F I G U R E 1 Green organic lightemitting diode (OLED) structure used for the simulation with Alq 3 and NPB as electron transport layer (ETL) and hole transport layer (HTL), respectively, and a PH1/Ir(ppy) 3 host-guest complex as emitting layer (EML); energies in electronvolts electric field distribution for case (b) in Figure 1 using different mesh sizes for the 1D-DD part where it is shown that the resulting electric field is not sensitive on the mesh size if it is chosen to be small enough (h < 0.1 nm). Excitonic parameters for the different materials were used as also used by Zhou et al.; see table 2 of Zhou et al. 9 . The OLED analyzed here uses Ir(ppy) 3 as phosphorescent emitter, and therefore, the triplet density profile is of main interest as the luminescence will mainly result from radiative decay events on the dye dopants. The resulting yz-averaged triplet density profiles are shown in Figure 2 and compared with the profiles obtained using the full ME approach. The results are shown for different dye dopant density profiles (molar fraction of total number of sites ρ Ir ppy ð Þ 3 ) at a bias of 4 V. Each profile calculated using the present hybrid approach has been averaged over three runs where the random energies were sampled anew. Thanks to the computationally efficient DD model implementation in SETFOS, the electronic part can be simulated in 1s-3s, whereas the more demanding ME model for the exciton dynamics takes around 3-5 min to compute on a standard personal computer (Intel ® Core™ i5-6600 CPU @3.3GHz) for a simulation grid of 120 × 50 × 50 sites. Even though only a 1D DD model was used for the charge carrier dynamics, the resulting exciton densities match well with the profiles obtained by Zhou et al. using a computationally demanding full ME approach.
The peak in triplet density at the EML-HTL interface results from the fact that electrons and holes accumulate at said interface resulting in high local recombination. When the Ir(ppy) 3 doping is increased, hole mobility is effectively increased in the EML because of denser transport states, resulting in a broader recombination profile inside the EML. Specifically, the increase in hole mobility results from the fact that the doping molecules, over which hole transport occurs, are closer to each other on average. If it is assumed transport is dominated by hopping between such localized states, the microscopic hopping rates are given by Miller-Abrahams hopping 15 where ν 0 is the hopping rate in the zero-distance limit, γ the inverse localization radius of the electron wavefunction, and r ij the spatial separation between states. f(ΔE) is defined as given in Section 2. Therefore, because of the exponential dependency on the spatial separation between the localized states on the guest molecules, the hopping rates are strongly increased, resulting in a higher macroscopic hole mobility. Such a spatial mobility profile is achieved in the DD formalism by subdividing the EML into sublayers, for each of which a specific mobility was defined. The linearly decreasing doping profile prevents at the same time too high recombination close to the EML-ETL interface leading to exciton quenching into the Alq 3 layer. The resulting hole mobility distribution is exemplarily shown for case (b) of Figure 1 in Figure S3 of the supplementary information.

| Transient electroluminescence
The model is further assessed by simulating the transient turn-on and turn-off behavior of a simple single layer OLED. A short voltage pulse is applied on the electrodes (turn-on at t = 0μs at a bias of 10 V, turn-off at t = 150μs, see also Figure 3A) with a total simulation time of 300 μs. As a simplification on the excitonics only generation by charge recombination, migration through the layer and radiative decay were enabled (see Table 1 for the relevant excitonic parameters used). The resulting singlet densities were again averaged over the yz dimensions and are plotted against time in Figure 3B. The results represent a single realization of the random energy distribution, which leads to a pattern of spatially varying singlet density where excitons are drawn to sites with lower binding energy leading to higher exciton densities in planes with a high share of such low-energy sites. The high hole mobility in the layer leads to an accumulation of holes close to the cathode at x = 120 nm, which leads to strong recombination and consequently a high singlet density (Purcell effect was neglected here). The very sharp decrease in singlet numbers at turn-off (t = 150 μs) happens because of the sudden decrease in electric field leading to lower carrier mobilities and consequently in a sharp decrease of Langevin recombination (see Figure 3C).

| Combined transient electrophotoluminescence experiment
Wehrmeister et al. 16 investigated the influence of triplettriplet annihilation as well as polaron quenching processes in phosphorescent OLEDs on the efficiency roll-off at high current densities. The device under investigation was a multilayer OLED featuring a Ir(MDQ) 2 (acac):α-NPD blend EML with 8 wt% of the red phosphorescent dye Ir(MDQ) 2 (acac). The device was subjected to a voltage pulse with an intermittent optical laser pulse to introduce an additional nonequilibrium exciton distribution, where the subsequent decay of luminescence after the optical pulse was recorded for various voltage levels.
A similar case was investigated using the present hybrid model, which lends itself for this study of exciton-  16 Energies given in electronvolts polaron interactions. As the electric characteristics of the device are governed by the hole and electron blocking layers (HBL and EBL, respectively) as well as the EML, these three layers were considered for the electrical DD simulation with energy levels as shown in Figure 4. The host-guest complex in the EML was again modeled as a layer with trap states; however, here the recombination on guest sites was modeled as a direct interaction of the electron and hole trap states (trap-trap recombination) instead of using a Shockley-Read-Hall expression. Following Wehrmeister et al., only triplets were considered in the ME model as a simplification due to the phosphorescent nature of the emitter material.
As a first step, the steady-state of the device under different applied voltages was simulated ranging from 0 up to 3.5 V using the DD model from which the background recombination (and hence the exciton density due to electrical excitation in the ME model) and the background charge carrier densities were calculated. The resulting JV-curve along with some of the calculated recombination profiles is depicted in Figure 5A and Figure 5B. On top of these results, a transient simulation was run with a 50ns short laser pulse after t = 1 μs leading to the optical generation of nonequilibrium excitations on the emitter molecules. The absorption profile was calculated using a thin-film optics simulation for the entire stack (with different layers L1-L8 consisting of glass substrate, ITO, HTL (p-NHT5), EBL (α-NPD), EML (Ir(MDQ) 2 (acac):α-NPD), HBL (NET5), ETL (n-NET5), and the Ag electrode, respectively), showing strong absorption in the HTL as can be seen in Figure 5C. However, because of the fairly large distance to the EML, it was assumed that these excitons will decay nonradiatively before being able to reach the EML and that they therefore do not contribute to the triplet densities in the EML. As also taken into account by Wehrmeister et al., the microcavity effects lead to a spatial dependency of the radiative exciton lifetime in the EML as shown in the inset of Figure 5C.
Using these inputs from the DD model as well as the thin-film optics, the ME model was then used to calculate the transient luminance of the device under study for the different voltage levels, as the optically excited triplets decayed over time and with them the luminance as the device approached again the steady-state electroluminescence caused by the electrical injection of charge carriers.
The relevant excitonic parameters for Ir(MDQ) 2 (acac) were taken as reported by Wehrmeister et al. and are given again for completeness in Table 2.
In order to compare the different transient luminance for different current densities, the resulting data were normalized where the electrical background luminance was set to zero and the optically induced nonequilibrium part was normalized to one. The resulting curves are plotted in Figure 6 with the corresponding lifetimes according to a monoexponential fit. From Figure 6A it is clearly visible (the staircase pattern stems from the discrete time stepping of the simulation) how the effective lifetimes of the emitter triplet states are reduced for higher voltages (and consequently higher current densities). This is due to the stronger quenching of triplets by polarons and other triplets in the vicinity. In order to check the relative impact, triplet-triplet annihilation was artificially turned off and the simulation was repeated with the results depicted in Figure 6B. Comparison with Figure 6A makes it clear that TTA is only a substantial quenching process for high voltages as the resulting effective lifetimes stay approximately the same except for the case of V = 3.5 V, despite the lower numerical value used for the TTA lifetime. If the simulation is run once more but with additionally tripletpolaron quenching deactivated, the transient luminescence decay changes dramatically; see Figure 6C. There is no dependency of the effective lifetimes on voltage bias anymore, as all the remaining decay processes are independent of charge carrier and triplet densities. Differences in the lifetimes arise probably purely out of the probabilistic nature of parts of the ME model. The dominance of TPQ as quenching process at higher current densities is in accordance with the findings of Wehrmeister et al. We see, however, that toward very high current densities, the role of TTA as quenching process starts to become more important.
The effective lifetimes reported here are in slight disagreement with the ones found by Wehrmeister et al. This is to be expected, however, as the rates used in the present ME model were taken as reported in the reference, which were obtained by fitting zero-dimensional models and might therefore not be equivalent to exact rates in a discrete 3D transport problem where a combination of various processes is taken into account together. Still, the present hybrid method shows clearly its ability to model such strongly coupled problems involving optics, carrier dynamics, and excitonics.

| CONCLUSION
A hybrid model for OLED devices consisting of a 1D DD model for the electronic part (charge carrier densities and potential) and a 3D ME model for the excitonics has been presented. The two models are coupled through the charge carrier densities (through polaron quenching), recombination rates, and the electric field (weak electric field dependency of Dexter transfer rates). Furthermore, the 3D ME model can be directly coupled to an optical thin-film simulation providing the necessary information of cavity effects and light outcoupling. The 1D DD model can be solved very efficiently and fast. At the same time, the 3D ME model allows for physically sound transport models using Förster and Dexter transfer rates taking long-range interaction into account. Such long-range interactions could be additionally introduced in order to describe nonlocal quenching processes for even higher physical accuracy. Furthermore, the discrete nature of the model allows for very easy treatment of layer interfaces as well as interspersed host-guest complexes. Despite the lower dimensionality of the DD model, the subsequent ME simulation for the exciton dynamics yields physically accurate results with the benefit of lower computation times compared with full 3D ME or kMC approaches for both the charge carrier and the exciton dynamics.
The scope of applicability of the present model includes OLEDs, quantum dot solids, and other optoelectronic devices featuring materials with strongly discrete structures, and it can be used for steady-state as well as transient characterization of such devices. The modeling of host-guest complexes is also possible through trapping dynamics in the 1D DD model, which is then translated into different densities and rates for host and guest molecules in the 3D ME formalism. The ME model is then capable of taking into account a multitude of physical processes at still moderate computational cost.

ACKNOWLEDGMENTS
We wish to thank Ralf Hiptmair and Stéphane Altazin for fruitful discussions as well as Yannick Masson for Dr. Christoph Kirsch specializes in applied mathematics. He got his master's degree from ETH Zurich in 2001 and his PhD from the University of Basel in 2005. He has used mathematical methods and developed numerical simulation tools in various fields such as plasma physics, thermal process engineering, and organic photovoltaics during postdoctoral stays in France, Germany, and in the United States. Since 2012, he has been working as a research associate in the Organic Electronics and Photovoltaics group at the Institute of Computational Physics of the Zurich University of Applied Sciences in Winterthur, Switzerland.
Dr. Urs Aeberhard received his PhD in theoretical solid state physics in 2008 from ETH Zurich, Switzerland, for a thesis on a microscopic theory of quantum well solar cells and performed in the Condensed Matter Theory group at Paul Scherrer Institute (PSI). From 2009 to 2018, he was a postdoc, senior scientist, and group leader at the Institute of Energy and Climate Research-Photovoltaics (IEK-5), Forschungszentrum Jülich (Germany), heading the activities in theory and multiscale simulation. In 2013, he stayed as visiting research scholar at the National Renewable Energy Lab (NREL) in Golden, Colorado. Since 2018, he has been employed as an R&D scientist at Fluxim and as a guest lecturer for the simulation of photovoltaic devices at ETH Zurich. Dr. Aeberhard is the developer of pioneering quantum-mechanical simulation approaches for nanostructure-based solar cells and the author of numerous publications on this topic.
Dr. Balthasar Blülle received his PhD in Physics in 2016 from ETH Zurich, where he investigated charge dynamics in organic molecular crystals by means of numerical simulations and charge transport measurements. During a postdoctoral stay at the University of Tokyo, he further focused on monolayer-thin solution-processed organic FETs, before joining the Institute of Computational Physics, Zurich University of Applied Sciences ZHAW, as a research associate in 2017. In parallel, he worked at Fluxim Inc. as an R&D engineer in the hardware and software development for OLED and emerging PV research. Since 2018, he has been product manager of Fluxim's optoelectronic simulation software SETFOS.
Dr. Sandra Jenatsch obtained her PhD from the EPFL (CH) in 2017. Her work on organic solar cells and light-emitting device which was carried out at Empa, was awarded with the Prof. René Wasserman price in 2018. Currently, she contributes to a variety of research projects of Fluxim and is also coaching PhD and MSc students in collaboration with several European universities. Some of her responsibilities include the training of business partners and customers using the scientific tools of Fluxim.
Prof. Dr. Beat Ruhstaller is a lecturer at the Zurich University of Applied Sciences ZHAW and founder of Fluxim. After a Diploma in Physics from ETH Zürich, he obtained his PhD in Physics at the University of California, Santa Cruz (USA), in 2000. He was a postdoc at the IBM Zurich Research Laboratory in the display technology group before joining ZHAW, where he headed the Institute of Computational Physics from 2007 to 2010. In 2006 he founded Fluxim, which he has managed as CEO since 2011. Fluxim has successfully brought R&D tool innovations from the lab to the OLED display and lighting as well as solar cell industry worldwide.

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section at the end of this article.