Disorder compensation controls doping efficiency in organic semiconductors

Conductivity doping of inorganic and organic semiconductors enables a fantastic variety of highly-efficient electronic devices. While well understood for inorganic materials, the mechanism of doping-induced conductivity and Fermi level shift in organic semiconductors remains elusive. In microscopic simulations with full treatment of many-body Coulomb effects, we reproduce the Fermi level shift in agreement with experimental observations. We find that the additional disorder introduced by doping can actually compensate the intrinsic disorder of the material, such that the total disorder remains constant or is even reduced at doping molar ratios relevant to experiment. In addition to the established dependence of the doping-induced states on the Coulomb interaction in the ionized host-dopant pair, we find that the position of the Fermi level and electrical conductivity is controlled by disorder compensation. By providing a quantitative model for doping in organic semiconductors we enable the predictive design of more efficient redox pairs.

Supplementary Figure 6. The evolution of the HOMO peak depending on the layer thickness. a-e , Peak of the HOMO distribution (HOMOpeak) as a function of the distance to the center of the "single-layer device" at five different layer thicknesses from 10 nm (a) to 50 nm (e). Dopant molar ratio is 0.1%. EF denotes Fermi level, bulk value of the HOMO peak and the HOMO peak in the center of the smallest, 10 nm sample are denoted as HOMO peak center (bulk) and HOMO peak center (10 nm), respectively. Shift of the Fermi level with respect to the HOMO peak is denoted by Δ F peak . For small thickness, the HOMO peak in the center of the simulated layer may look flat, although an organic film is fully depleted. Only for the thickness of about 40...50 nm (see (d) and (e)), the HOMO in the center of the system stops changing its position, indicating that the bulk value of the HOMO is reproduced in the center of the layer.

Supplementary Note 1: Polaron dynamics in doped organic semiconductors
First, we consider factors that determine ionization of the dopant. We start with the disorder-free doped organic semiconductor that consists of the host molecules (H) and dopant molecules (D). In this scenario, all host (dopant) molecules have the same electron affinity EAh (EAd) and ionization energy IPh (IPd) such that IPh -EAd is small in comparison with the gap. An adjacent pair of D and H molecules can be found either in a neutral or ionized microstate (see Supplementary Figure 1) according to the local equilibrium of the reaction H + D ⇆ H + + D − . Neglecting electrostatic interaction with other molecules, the energy difference between these two microstates would be: where Δ off = EA d − IP h is the molecule specific frontier orbital energy difference, and C = −e 2 /4π 0 is the electrostatic interaction energy between the two molecules with a being the host-dopant separation, electric constant, 0 relative material permittivity and e elementary charge. Because this energy is always negative, charge separation can even occur if EA D < IP H .
In organic thin films and devices, the electrostatic interaction between the ionized pair (state "I") and positive ( ) and negative ( ) charges in the environment cannot be neglected, which introduces a third term to Δ ext : , , h , d are the positions of positive and negative charges in the environment and the positions of the host and dopant molecule, respectively. Due to the random dopant distribution in an amorphous thin film, ext behaves as a stochastic quantity that can be either positive or negative. Moreover, this quantity changes dynamically depending on the position of all charges during the simulation. Note that ext is computed using Ewald summation 1 as described in reference 2 .
Intrinsic energy disorder results from slightly different molecular conformations and the unique electrostatic environment of each molecule, which is inherent to all disordered organic semiconductors. To include this effect, we employ the standard model that IP of host and EA of dopant molecules are distributed according to Gaussian distributions with the expectation values IPh and EAd and standard deviation IP and EA . Therefore the energy difference between the neutral and ionized state has an additional time- Including all these contributions, the central equation of the dopant ionization process takes on a form: Supplementary Equation 1 shows that the dopant ionization energy is defined by the following contributions: energy difference between IP h and EA d , host-dopant Coulomb interaction (always negative, for a dielectric permittivity of = 4 and a host-dopant distance of 1 nm equal to 0.36 eV), intrinsic disorder (positive or negative), (doping-induced) environmental disorder (positive or negative), the local potential due to the applied electric field ( and are vector quantities; signs "+" and "-" are for electrons and holes, respectively). In spite of their simplicity, the last two or even three terms of Supplementary Equation 1 to date not been considered in explicit mesoscopic models, mainly because of the associated numerical complexity. The hopping process between two host molecules and injection from electrodes if any is described with similar rates.
We note that IP and EA levels computed with Supplementary Equation 1 can only be compared directly to measurements of bulk materials and not gas or liquid phase measurements due to environment specific shifts of the energy levels 3,4 .

Supplementary Note 2: Microscopic processes
The following kinds of microscopic processes are considered: hopping of holes/electrons; dopant ionization/neutralization and electrons/holes injection/ejection.

Hopping of holes and electrons.
It is assumed that molecules can be in the neutral or first oxidation state. If a molecule is positively charged, this is interpreted as a missing electron (hole) in the HOMO; if it is negatively charged, it is interpreted as an additional electron in the LUMO. Changing of the oxidation state of the molecule due to charge transfer (for example in the reaction between two host molecules): H + H → HH + is interpreted as a hole hop. The hopping rate of an electron/hole from site to , is described by the Miller-Abraham (MA) expression 5 : where and energies of the final and initial state, 0 is the attempt-to-escape frequency, is the localization radius of a charge carrier, B is Boltzmann constant, T is the temperature. The usage of MA expression in this work has been preferred to the Marcus hopping rate 6 to avoid one more free parameter (reorganization energy), which is not important to study the general doping-induced phenomena. Apart of this, it has been shown that both methods lead to similar results in bulk material simulations 7 .
Corresponding energy difference is described by the last three terms of Supplementary Equation 1. The energy difference for the hopping process of holes or electrons corresponds to the last three terms of Supplementary Equation 1: where = − is a hopping vector, where E ext and E int are differences in the energy of the site and due to extrinsic and intrinsic disorder.

Dopants activation.
Hole dopant is by definition the molecule, whose EA is close to the IP of the host material (in terms of the molecular orbitals, HOMOh ≈ LUMOd). To enable dopant activation (ionization of the dopant and generation of the hole on the host material) we allow the redox process between the host H and the dopant D to occur: The generated charge in the host is now treated as a hole. The reciprocal process, a dopant deactivation, is also allowed. The rate of the dopant ionization (or activation) is also described by the MA expression (Supplementary Equation 2), where and are now states before and after dopant activation, is the dopant-host distance. It must be noted that the energy of the ionized dopant-host pair is reduced due to attractive Coulomb energy of the generated charge pair. For each kMC simulation performed in this work, we compute the time-averaged fraction of the ionized dopants as: where is a dwelling time of the i-th microstate, Σ = ∑ ; and d − are number of all and ionized dopants. The averaging is performed over the last 50% of the kMC trajectory. In this manuscript, we have simulated a very strong dopant with EA higher than the host IP by 1.5 eV. Therefore, the computed value of is 1 within a machine accuracy. This makes the molar ratio of ionized dopant defined as:

Injection rate is heuristically described by an Arrhenius-like equation:
where 0 is an effective metal-organic coupling analogously to the attempt to jump frequency in the MA expression, is the distance between the metal surface and the given molecule, is an empirical parameter that accounts for the distance dependence of the metal-organic coupling.
For electron/hole injection, el / hole are defined as where is the interaction energy between a charge and its mirror charge.

Supplementary Note 3: Extraction of the HOMO and LUMO + distribution
The distributions of the HOMO and LUMO + levels of the host molecules have been computed from kMC trajectory "measuring" site energies of neutral hosts as − IP and site energies of ionized hosts as − EA + . The resulting DOS is the sum of these two quantities. Note that hereinafter IP, EA, HOMO and LUMO + though used without subscript "h" relate to host molecules.
In more details, at each i-th frame of the kMC trajectory we compute a binned distribution of the IP of the hosts, which are neutral at the i-th frame, IP ( ), and IP of the hosts, which are positively charged at the i-th frame, EA + ( ), as well as the dwelling time of the system at the i-th frame. Then, we find the time-average of these distributions: We note that this definition of the density of states takes into account that the energy of a given particle depends on the position of other − 1 particles with being the total number of particles. To our understanding, this is the maximum information we can retrieve from the kMC trajectory to DOS. To our knowledge, although the density of states of doped organic semiconductors has been recently reported 8 , this type of the DOS averaging has never been done before. This limited the insights gained from kMC simulations, in particular about reliable information about the position of the doping-induced trap states and the DOS modification, which is essential to understand what governs the doping efficiency.
There is another possible way that one can use to compute DOS from the kMC trajectory. In this approach, one first computes time-averaged charge density at each site during some number of iterations (and thus, time-average external potential on each host site). After that, one computes the distribution of IP and EA + of the sites subjected to this average potential.
In this case, valuable information about the correlations contained in the kMC trajectory is lost. Such a method is essentially a mean-field method. It is likely that in already mentioned work 8 authors used this "mean-field" method, as they observe no signs of EA + (LUMO + ).
We have implemented this latter method (referred to as "mean-field") and compared it to the one used throughout this work (referred to as "many-body"). The results are shown in Supplementary Figure 2 for low (a) and high (b) dopant concentration. One can conclude that the LUMO + states are indeed missing. At low dopant concentration, HOMO distribution computed by "many-body" method and the total DOS computed by "mean-field" method are almost the same. At higher dopant concentration, however, even this fails: the "meanfield" DOS is shifted towards higher energies. Such "mean-field" method obviously prevent reasonable interpretation of the doping in organic semiconductors, which was fixed by our many-body approach.

Supplementary Note 4: Extraction of the bulk Fermi level shift
In case of the "bulk semiconductor" The Fermi level EF is extracted from kMC simulations of the "bulk" system a posteriori via the charge neutrality equation: Here DOS( ) is the density of states (namely the number of host states per molecule of a material per energy unit) and ( , ) is the Fermi-Dirac distribution function. Polarons in the organic semiconductors are assumed to obey Fermi statistics.
The left part of this equation, DMR − is the molar ratio of ionized dopants and the right part is the molar ratio of positively charged host molecules. As far as the material is neutral, the Fermi level can be found that turns this equation into equality. For the energetics of hosts and dopants used in this manuscript, all dopants are ionized and in Supplementary Equation 5: DMR − = DMR.
As mentioned in the manuscript, the Fermi level is always the intersection energy of LUMO and HOMO associated DOS. This is because holes associated with LUMO distribution and electrons associated with HOMO distribution obey Fermi statistics, in which the Fermi level corresponds to the energy where the state (if exists) is occupied with 50% probability. Therefore, the density of LUMO + and HOMO states is always the same at this energy. Supplementary Figure 3  Note that in this work we have computed the position of the LUMO + /HOMO distribution with respect to the Fermi level using two independent methods. The first method is based on the charge neutrality equation (Supplementary Equation 5). The second method is based on the equilibration with the fixed chemical potential of electrodes (see below).
Although the second method does not require explicit postulation that the particles are Fermi distributed, both provide the same results and put the Fermi level exactly at the HOMO/LUMO intersection, which means that holes obey Fermi statistics.

In case of "single layer devices"
We simulate the "single layer device" system with initial energy levels as shown in Supplementary Figure 4a. To obtain a Fermi level of the bulk semiconductor, we use the following system specifications: the length of the system for each dopant concentration is chosen in a way that the space charge region is significantly shorter than the half of the layer thickness. This ensures that the properties in the center of the layer are the same as in a bulk semiconductor; to efficiently reach convergence of the charge distribution, hopping from the electrodes to each site of the left/right half of the organic system is allowed. The parameter in Supplementary Equation 4 is assumed to be infinite, which means that each molecule of the semiconductor is connected to electrodes with the same electronic coupling strength. This influences the time scale of the kMC simulation but fulfills detailed balance and thus does not change the equilibrium distribution. In other words, each material site is connected directly to the electron reservoir (electrodes), and thus equilibration is achieved in the fastest possible way. In practice, all simulated systems reach equilibrium after 10 5 kMC steps that has been shown by longer simulations of 10 6 and 10 7 steps. The DOS results presented in the paper were time averaged over the10 3 last steps of a kMC trajectory of 10 6 steps.
When extracting Fermi level in the centre of organic thin films shorter than double the width of the depletion layer (this quantity has a physical sense of a "hole injection barrier"), the thickness of the layer is set to the value we are interested in (Supplementary Note 5 provides more details). The rest of the described procedure holds.

Supplementary Note 7: Effect of the trap states filling
In the ultralow doping regime, realistic semiconductors show another important effect, the passivation of traps 9 . Different trap states concentration have been reported ranging from 10 16 cm -3 in pentacene 10 to 10 18 ... 10 19 in C60. These effects are observed as the prompt increase of the conductivity upon doping in the ultralow doping regime, normally below DMR < 10 -3 as documented in the literature 9,11,12 . In order to show the applicability and the relevance of the disorder compensation effects to realistic organic semiconductors, we made a comparison of a pure materials simulated in the manuscript and the same materials, but with traps. Figure 8 shows the effect of the traps on the conductivity and the Fermi level position. We have added exponentially distributed trap level with parameters taken from the reference 9 (molar ratio of traps is 0.007) to the doped materials we have simulated in Figure 4. We have also extended and shifted the interval of the simulated dopant molar ratios down to DMR = 10 -4 to better see the effect of traps passivation (this additional interval of DMR is filled in grey in Supplementary Figure 8). Conductivity of the materials with and without traps are plotted in the same color, but for materials with traps lines are semitransparent. As reported in the literature and as expected from traps concentration, the main effect of the trap states vanishes shortly after dopant molar ratio exceeds 10 -3 . After that, at moderate doping concentrations, conductivity increases either linearly or superlinearly depending on the intrinsic material disorder, which is explained in this work by means of the disorder compensation effect. The slope of the Fermi level shift upon doping slightly increases mainly below DMR = 10 -3 , but for bulk samples it is still 10 times less than 1 eV/dec. To explain UPS measurements that reports such a prompt (1 eV/dec) Fermi level shift, one needs to involve effects due to finite layer thickness as described above and in the manuscript.

Supplementary
It is noticeable that highly disordered materials are less sensitive to the traps that we have introduced (for example, the conductivity of the material with intrinsic disorder 0.15 eV doesn't visually change upon introducing traps, Supplementary Figure 8a). This is because exponentially distributed shallow traps, which we used in this case, are largely overlapped with the tail of the HOMO states, if HOMO intrinsic disorder is high. Thus, in highly disordered materials, relevant traps must be deeper to cause pronounced trap-filling effects.