Parameter-free continuous driftdiffusion models of amorphous organic semiconductors

Continuous drift–diffusion models are routinely used to optimize organic semiconducting devices. Material properties are incorporated into these models via dependencies of diffusion constants, mobilities, and injection barriers on temperature, charge density, and external field. The respective expressions are often provided by the generic Gaussian disorder models, parametrized on experimental data. We show that this approach is limited by the fixed range of applicability of analytic expressions as well as approximations inherent to lattice models. To overcome these limitations we propose a scheme which first tabulates simulation results performed on small-scale off-lattice models, corrects for finite size effects, and then uses the tabulated mobility values to solve the drift–diffusion equations. The scheme is tested on DPBIC, a state of the art hole conductor for organic light emitting diodes. We find a good agreement between simulated and experimentally measured current–voltage characteristics for different film thicknesses and temperatures.


Introduction
Macroscopic and mesoscopic properties of organic semiconductors, such as charge carrier mobility or the width of the density of states, are often extracted by fitting the solution of the drift-diffusion equation [1][2][3][4][5] to the experimentally measured current-voltage (I-V) characteristics.The mobility and diffusion constant of charge carriers, which enter these equations, depend on charge carrier density, r, electric field, F, and temperature, T. 6,7 For onedimensional transport and specific rate expressions these dependencies can be obtained analytically. 8,91][12] The extended Gaussian disorder model (EGDM), 6 for example, provides a parametrization of the mobility, m(r,F,T), for uncorrelated, Gaussian-distributed site energies, while the extended correlated disorder model (ECDM) 11 additionally accounts for spatial site energy correlations due to long-range charge-dipole interactions.
The aforementioned approach has become a standard tool for analyzing experimental data. 13It has, however, several issues: (i) Gaussian disorder models are parametrized only for materials with moderate energetic disorder, s o 0.15 eV at room temperature, while many amorphous materials have a higher s.(ii) The spatial correlation of site energies in the ECDM is material-independent and has an (approximate) 1/r decay, where r is the intermolecular distance, but recent studies show that this decay may be different. 14(iii) Due to the non-Gaussian shape of the density of states, 15 the energetic disorder and the lattice constant are different from those provided by microscopic calculations, 14 thus making them merely fitting parameters without a comprehensive link between macroscopic properties and the chemical composition of the material.
In this paper we propose an approach which does not have these limitations.In a nutshell, the mobility dependence on charge density, field, and temperature is first tabulated by combining quantum mechanical, classical atomistic and coarsegrained stochastic models for charge transfer and transport.These tables, corrected for finite-size effects, are then used to solve the drift-diffusion equations.
To illustrate the advantages of the method, we compare it to the ECDM and the Mott-Gurney model 16 as well as to experimental measurements performed on amorphous layers of Tris[(3-phenyl-1H-benzimidazol-1-yl-2(3H)-ylidene)-1,2-phenylene]Ir (DPBIC), a hole-conducting material used in organic light emitting diodes (OLEDs) 17 and organic photovoltaic cells (OPVs). 18he paper is organized as follows.In the Methods section we describe the coarse-grained, off-lattice transport model, the procedure used to tabulate the charge carrier mobility, the algorithm used to solve drift-diffusion equations, and the parametrization of the extended correlated Gaussian disorder model.The entire workflow is summarized in Fig. 1.We also recapitulate the main results of the Mott-Gurney model and provide details of experimental measurements.The I-V curves, electrostatic potential, and charge density profiles are then compared in Section 3, where we also validate the transferability of the method by studying different layer thicknesses and temperatures.A short summary concludes the paper.

Tabulated mobility
To tabulate charge carrier mobility as a function of temperature, field, and charge density, we first simulate amorphous morphologies of N = 4000 molecules using molecular dynamics simulations in the NPT ensemble with a Berendsen barostat and thermostat. 19The simulation box is equilibrated at 700 K for 1 ns, which is well above the glass transition temperature, and then quenched to 300 K during 1.3 ns.The force-field is tailored for the DPBIC molecule as described elsewhere 20 by performing potential energy scans using density functional theory (DFT) calculations with the B3LYP functional and the 6-311g(d,p) basis set.The Gaussian package 21 was used for all energy calculations.
The charge transport network is then generated as follows.A list of links is constructed from all molecules with adjacent conjugated segments closer than 0.7 nm.For each link a charge transfer rate is calculated using Marcus theory, i.e., in the hightemperature limit of the non-adiabatic charge transfer theory, 22 where h is the Planck constant, k B the Boltzmann constant, and T the temperature.Electronic couplings J ij are evaluated for each dimer by using the dimer projection method, 23 the PBE functional and the def2-TZVP basis set.These calculations were performed using the TURBOMOLE package. 24Note that the values of electronic couplings can deviate by up to 50%, depending on the functional and the basis set size. 25,26This deviation is, however, systematic and will result in a constant prefactor for the mobility, i.e., we do not expect any changes in functional dependencies on the external field, charge density, or temperature.
The hole reorganization energy, 27 l ij = 0.068 eV, was evaluated in the gas phase using the B3LYP functional and 6-311g(d,p) basis set.Site energy differences, DE ij = E i À E j were evaluated using a perturbative scheme 28 with the molecular environment modeled by a polarizable force-field, parametrized specifically for these calculations.In this approach, the site energy + qFÁr i is the sum of the gas phase ionization potential, E int i = 5.87 eV, an electrostatic part, E el i , an induction contribution, E pol i , and the contribution due to an external electric field, qFÁr i .The mean value of these energies gives an ionization potential of E IP = 5.28 eV.
The electrostatic contribution was evaluated using the Ewald summation technique 29,30 adapted for charged, semi-periodic systems 31,32 and distributed multipole expansions. 33,34Note that using an interaction cutoff would yield a shifted energetic landscape with an underestimated spatial correlation of energies. 35he induction contribution, E pol i , was calculated self-consistently using the Thole model 36,37 with a 3 nm interaction range.Note that the set of Thole polarizabilities were scaled in order to match the volume of the polarizability ellipsoid calculated using the B3LYP functional and 6-311g(d,p) basis set.This step is required to account for larger polarizabilities of conjugated, as compared to biological, molecules.
The resulting charge transport network is used to parametrize the coarse-grained model, by matching characteristic morphological and transport properties of the system, such as the radial distribution function of molecular positions, the list of neighboring molecules, the site energy distribution and spatial correlation, and the distance-dependent distribution of transfer integrals. 14,38The coarse-grained model allows to study larger systems, here of 4 Â 10 4 and 4 Â 10 5 sites, which are required to perform simulations at low charge carrier densities, in our case from 0.025 down to 10 À5 carriers per site.
Charge transport is modeled using the kinetic Monte Carlo (KMC) algorithm.Note that charge carriers interact only via the exclusion principle, i.e., a double occupation of a molecule is forbidden.Charge mobility is evaluated by averaging the carrier velocity along the field, m = hviÁF/F 2 .KMC simulations are repeated for eight different temperature values, from 220 K to 992 K, and twelve field values, in the range of 2.5-30 Â 10 7 V m À1 .
To avoid finite size effects, an extrapolation procedure 39,40 is used for small charge carrier densities.The mobility is simulated at a range of higher temperatures, where mobilities are non-dispersive and hence system-size independent.The extrapolation to lower temperatures is performed by parametrizing the analytic mobility versus temperature dependence available for one-dimensional systems 8 or, alternatively, using the box-size scaling relation. 39he tabulated mobility is finally interpolated and smoothed by the scattered data interpolation method using radial base functions, 41 which can treat many-dimensional, unstructured data.

Drift-diffusion modeling
Macroscopic dynamics of electrons/holes (n/p) is modeled using one-dimensional drift diffusion equations @r n=p @t ¼ Àr Á J n=p ; (3) coupled to the Poisson equation Here c denotes the electrostatic potential, D is the diffusion constant, e 0 the vacuum permittivity and e r the relative permittivity.J = I/A is the current density, where I is the current and A the electrode area.In case of DPBIC we are interested in hole transport only, hence the electron current density, J n , and density, r n , are set to zero and only the hole equations need to be solved.To simplify the notation we omit the index n/p.Here we are interested only in the steady state, i.e., qr/qt = 0 in eqn (3).
Since charge carriers occupy energetic levels according to Fermi-Dirac statistics, the carrier density is related to the quasi-Fermi level, Z, as where g(E) is the density of states and V the box volume.The diffusion coefficient and mobility in eqn (2) are related via the generalized Einstein relation 42 Eqn ( 2)-( 6) are solved using an iterative scheme, until a selfconsistent solution for electrostatic potential, c, density, r, and current, I, is found. 43First the equations are rescaled to ensure numerical stability, which is necessary since carrier density and electrostatic potential vary by several orders of magnitude.Then they are discretized according to a scheme proposed by Scharfetter and Gummel, 44 linearized, 45 and solved by using the Gummel iteration method, 46 adapted to organic semiconductors at finite carrier density.This method is less sensitive to the initial value than a Newton algorithm and thus is the method of choice despite its slower convergence 47 in terms of iteration steps.The tabulated mobility values, m(F,r,T), computed in Section 2.1, are used while solving eqn ( 2)-( 6).
We use Dirichlet boundary conditions for the electrostatic potential, c, by setting the potential difference at the boundaries to c eff = V app À V int , where V app is the applied potential and V int the built-in potential, defined as the difference of the materials' work functions.For ITO and Aluminum we use experimental values: the work function of ITO is reported to lie in the range from 4.15 eV to 5.3 eV, [48][49][50][51] and for Aluminum from 4.06 eV to 4.26 eV. 52Here we assume average values of 4.73 eV for ITO and 4.16 eV for Aluminum.In combination with the calculated DPBIC solid-state ionization potential (IP) of 5.28 eV, which is the mean value of the site energies, E i , that are calculated as described before, this yields injection barriers of DE ITO = 0.55 eV and DE Al = 1.12 eV.
The charge density at the electrodes is fixed to the density resulting from inserting DE ITO/Al into eqn (5).To model the doped interlayers (see Section 2.5) within a five nanometer range from both electrodes, an additional charge concentration of 3 Â 10 À4 carriers per site, estimated from previous calculation 53 is added in these regions when solving the Poisson eqn (4), leading to high hole densities in the doped regions even without space-charge limited effects.

Lattice model
To test the validity of lattice models, we have also parametrized the extended correlated disorder model (ECDM) 11 by fitting the simulation results to the ECDM expression for mobility.The fit was performed for charge densities, r, in the range of 8.7-140 Â 10 23 m À3 , including an extrapolated, non-dispersive zerodensity mobility, 40 and electric fields in the range of 3-9 Â 10 7 V m À1 .For the extrapolation temperatures from 1200 K to 50 000 K, giving non-dispersive transport in the small system, were used.All simulations were performed at 300 K.
The fit to the ECDM model yields a lattice constant of a = 0.44 nm, an energetic disorder of s = 0.211 eV, and a zero-field zero-density mobility of m 0 (300 K) = 1.8 Â 10 À13 m 2 V À1 s À1 . 20hese values serve mainly for providing a fitting and extrapolation function as they differ from the values observed in microscopic simulations (a = 1.06 nm, s = 0.176 eV, m 0 (300 K) = 3.4 Â 10 À12 m 2 V À1 s À1 ).

Mott-gurney model
The Mott-Gurney, or trap-free insulator model, 16 predicts a current density of where d is the thickness of the sample and e r the material's relative permittivity (here we have chosen e r = 3).This expression is only valid under the assumptions of (i) hole-only (or electron-only) transport, (ii) no doping, (iii) constant mobility and relative permittivity and (iv) no injection barriers.The electrostatic potential and hole density throughout the sample are then given by where 0 o x o d.A mobility of m = 3 Â 10 À22 m 2 V À1 s À1 has been chosen to provide the best match of the experimental data.

Experimental measurements
I-V curves were measured for three film thicknesses: 203 nm, 257 nm and 314 nm, including two interlayers of DPBIC doped with molybdenum trioxide (MoO 3 ) of 5 nm thickness on both sides of the DPBIC film.These serve to enhance the injection efficiency, which has been taken into account in our model by the previously mentioned additional charge in these regions.The hole-conducting DPBIC layer was sandwiched between a 140 nm indium tin oxide (ITO) anode and a 100 nm aluminum cathode.To control the temperature, the samples are placed into the oil reservoir of a cryostat, which allows for a variation between 220 K and 330 K.The voltage was varied between 0 V and 20 V.All films were fabricated by vacuum thermal evaporation of DPBIC on a glass substrate, patterned with the ITO layer.Thicknesses were determined by optical ellipsometry after a simultaneous deposition of the same amount of DPBIC on a silicon wafer.

Validation
We first compare the current-voltage characteristics, the electrostatic potential and the hole density profiles calculated using tabulated and ECDM mobilities, which are shown in Fig. 2  together with the experimentally measured current-voltage characteristics.One can see that the experimentally measured currentvoltage characteristics are well reproduced using the tabulated mobilities.The ECDM underestimates the current by an order of magnitude and there is a clear mismatch of the slope, as it can be seen in Fig. 2(a).It also predicts a negative electrostatic force at the beginning of the slab, Fig. 2(b), and a very steep charge accumulation at the injecting anode, Fig. 2(c).The disagreement is due to the high energetic disorder obtained from the fit, s = 0.211 eV, which is outside the range used to parametrize the ECDM expression.In addition, the ECDM does not reproduce the spatial correlation of site energies well.Finally, the Mott-Gurney model does not reproduce experimental results even qualitatively: it neither takes into account doped layers nor field-or density-dependence of the mobility.
To illustrate the transferability of the proposed method we also compare current-voltage characteristics for different temperatures and different film thicknesses.Fig. 3 shows that for high temperatures the agreement between theory and experiment is excellent.At 233 K deviations are significant and can be attributed to the breakdown of the drift-diffusion description, since at low temperature and large energetic disorder charge transport becomes dispersive, showing anomalous diffusion. 53Its description using equilibrium distributions, mobility and diffusion constant cannot be justified in this situation.Moreover, Marcus theory only applies to sufficiently high temperatures.The crossover temperature below which Miller-Abrahams rates 54 become a more appropriate description has been estimated to be about 250 K. 55

Conclusions
To conclude, we have proposed a parametrization scheme for drift-diffusion equations which is based on evaluation of charge transfer rates, simulation of charge transport in a coarse-grained charge transport network, and tabulation of charge carrier mobility as a function of field, charge density and temperature.The method is rather general, in part because it is not limited to functional dependencies build into the ECDM and EGDM models and, therefore, allows to treat systems with large energetic disorder and material-specific spatial site energy correlation functions.
Using this scheme, we have simulated I-V characteristics of a single-layer device, and found them to be in a good agreement with the experimentally measured I-V curves, whereas significant deviations have been observed for the ECDM and Mott-Gurney models.

Fig. 1
Fig. 1 Overview of the method.(a) Chemical structure used to parametrize atomistic force field.(b) Amorphous morphology obtained using molecular dynamics simulations.(c) A coarse-grained model for the charge transport network.(d) Kinetic Monte Carlo simulations are used to tabulate the mobilities.(e) Solution of drift-diffusion equations. ,

Fig. 2
Fig. 2 (a) Current-voltage characteristics, (b) electrostatic potential profiles, and (c) hole density profiles.Slab thickness 314 nm, temperature 300 K. (b and c) are plotted for an external voltage of 4 V.