Impact of Active Layer Morphology, Density of States, Charge Carrier Concentration, and Local Charge Density Fluctuations on Bimolecular Recombination of Bulk Heterojunction Solar Cells: A Theoretical Perspective

Abstract

We study the merits of a reaction-diffusion model to unravel the effects of active layer morphology and donor-acceptor interfacial roughness, density of states, charge carrier concentration, and local charge density fluctuations on the bimolecular recombination kinetics in bulk heterojunction organic semiconductors. We consider the cases of a single and composite electronic density of states (DoS) that consists of a superposition of a Gaussian and an exponential DoS. Using kinetic Monte Carlo (KMC) simulations, we apply the reaction-diffusion model in order to investigate the factors impacting bimolecular recombination (BMR) kinetics and rates at short and long time scales. We find that morphology, donor-acceptor interfacial roughness, and charge carrier concentration only affect BMR time, whereas DoS characteristics as well as local charge density fluctuations can significantly impact BMR kinetics and rates.

1. Introduction

Early organic photovoltaic (OPV) devices based on bilayer architectures suffered from limited interfacial area between the donor (D) and acceptor (A) components, which resulted in low charge separation efficiency. To overcome this issue, the bulk heterojunction (BHJ) architecture was introduced in the mid-1990s and consists of a D/A blend. Compared to the bilayer architecture, the blend structure provides a much larger D/A interfacial area that is distributed throughout the active layer, facilitating exciton dissociation and thereby charge separation. Recent data [1, 2, 3, 4, 5] using high-resolution temporal techniques underlie that exciton in high-performance devices can readily dissociate within less than 100 fs. In such cases, exciton formation, diffusion, and dissociation are efficient. However, it is the competition between charge migration toward the electrodes and bimolecular charge recombination [6], which appears to play the leading role in decreasing the performance of OPV devices. Thus, a detailed understanding of the loss mechanism emerging from BMR is critical.

While the dominant loss mechanism is still a matter of debate [7], there is increasing evidence [1, 8, 9] that in BHJ-based solar cells, non-geminate (bimolecular) recombination represents the primary factor limiting device performance [10, 11, 12, 13]. However, establishing the most relevant physical model of BMR remains a challenge. Previously, a number of experimental data show that, depending on the mobility of charge carriers [14], recombination evolves from trap-assisted first-order (monomolecular) dynamics under short-circuit current conditions to second-order (bimolecular) Langevin dynamics under open-circuit voltage conditions [1, 15, 16, 17, 18, 19]. Though, many recent studies [20, 21, 22, 23, 24, 25, 26, 27, 28, 29] indicate that BMR in organic solar cells significantly deviates from a traditional Langevin description.

In the case of polymer-fullerene-based devices, the charge decay dynamics at open-circuit voltage exhibit approximately a third-order dependence on charge density [24, 25]. Thus, in order to be consistent with experimental data, the proportionality constant of the Langevin model must depend on charge density. This third-order dependence of the BMR rate on charge density has been suggested to arise as a result of either a carrier lifetime dependence on charge density [30], recombination via an exponential tail of states [31, 32], or carrier trapping in an inhomogeneous distribution of localized states [20, 21, 26]. Also, it was suggested that although traps can sometimes enhance the dissociation of geminate pairs into free carriers [33], they also act as recombination centers leading to a Shockley-Read-Hall (SRH) recombination dynamics [34, 35, 36, 37, 38].

While experimental efforts have been extensive, there are few theoretical studies to date [26, 39, 40, 41] that have been conducted to understand the charge carrier loss mechanisms in organic solar cells. As a result, the mechanisms underlying BMR in OPV devices remain poorly understood. Thus, it is desirable to develop a thorough understanding of the BMR loss mechanisms that currently limit the enhancement of the efficiency of organic solar cells.

The purpose of the present work is to investigate how the kinetics of non-geminate recombination is affected by (i) the detail of the morphology of the active layer; (ii) energetic disorder and various distributions of densities of electronic states, (iii) charge carrier density, and (iv) local charge density fluctuations. Utilizing KMC simulations, we develop a reaction-diffusion model to explain the role that abovementioned factors play in BMR of BHJ OPV solar cells.

2. Theoretical model

To study the impact of interfacial roughness, DoS, and charge carrier density on BMR in a BHJ device at open-circuit voltage, we develop a theoretical three-dimensional reaction-diffusion lattice model. In this model, we only look at electron-hole recombination after exciton dissociation.

We represent electrons and holes as particles and antiparticles that can randomly diffuse through the space, restricted to a cubic lattice, and react (recombine) once they collide (occupy same lattice site). In our model, we restrict holes and electrons to move exclusively within the D/A domains of a randomly generated morphology, which are discussed next in this section. Based on the reaction-diffusion model, we perform KMC simulations on a cubic lattice of $100\times 100\times 100$ sites with lattice spacing equivalent to 1 nm, unless otherwise mentioned. We ignore additional charge generation or charge extraction processes at the electrodes. Since we are mainly interested in BMR kinetics at low charge carrier concentrations, we also neglect the effect of electrostatic interactions.

To evaluate the kinetics of BMR at the device scale, unless otherwise mentioned, we apply a reflective boundary condition on all sides such that charge carriers cannot exit from the simulation box until they entirely recombine. We perform KMC simulations on disordered morphologies. In creating bulk heterojunction OPV using solution processing, a disordered morphology, also called random morphology, is formed. This random morphology is comprised of two organic semiconductors, intermixed randomly within the active (light-absorbing) layer. We computationally generate the random morphologies utilizing the potential model [42, 43, 44, 45] in which random morphologies of D/A domains in space are represented by means of a random potential $\psi \left(x,y,z\right)$. $\psi$ is a random scalar function represented as a sum of several propagating plane waves in random directions given by

$$\psi \left(x,y,z\right)=\sum _{i=1}^{M}sin\left(\frac{2\pi }{\lambda }\left[k_{i1}x+k_{i2}y+k_{i3}z\right]+{\theta }_i\right),$$E1

where $k_{i1}=cos\left({\alpha }_i\right)sin\left({\phi }_i\right)$, $k_{i2}=sin\left({\alpha }_i\right)sin\left({\phi }_i\right)$, $k_{i3}=cos\left({\phi }_i\right)$, $M=25$, ${\alpha }_i$ and ${\theta }_i$ are uniformly distributed random numbers between $\left[0,2\pi \right]$, and ${\phi }_i$ is another uniformly distributed random number between $\left[0,\pi \right]$. $\lambda$ in Eq. (1) is directly proportional to the average size of D and A domains. $\psi$, whose value fluctuates between –M and M, randomly splits the space into D domains wherever $\psi \ge 0$ and A domains wherever $\psi <0$. In this way, $\psi$ allows the representation of a random morphology within the simulation box.

The diffusion and recombination of charge carriers are simulated using the KMC technique where the charge carrier transfer rate ${\nu }_{\mathrm{ij}}$ from any site i to one of the six nearest neighboring sites, j, is described by the Miller-Abrahams model [46]:

$${\nu }_{\mathit{ij}}={\nu }_{0}exp\left(-2\gamma r_{\mathit{ij}}\right)exp\left[-\beta \left(\mathrm{\Delta }{E}_{\mathit{ij}}+\left|\mathrm{\Delta }{E}_{\mathit{ij}}\right|\right)/2\right],$$E2

where ${\overrightarrow{r}}_{\mathit{ij}}={\overrightarrow{r}}_j-{\overrightarrow{r}}_i$, $r_{\mathit{ij}}=\left|{\overrightarrow{r}}_{\mathit{ij}}\right|$, and ${\overrightarrow{r}}_i$ are the position vectors of site i, in our simulations $r_{\mathit{ij}}=a=1$ nm represents the lattice site spacing, ${\nu }_0=7\times {10}^{12}$ s⁻¹ denotes the intrinsic attempt frequency, $\gamma =3\times {10}^7$ cm⁻¹ indicates the inverse of localization radius, $\beta ={\left(k_{B}T\right)}^{-1}$ is the Boltzmann constant, and $\mathrm{\Delta }{E}_{\mathit{ij}}=E_j-E_i$ is the energy level difference between adjacent sites i and j.

We make use of the first reaction method [47, 48] in the KMC simulations. At each step of the KMC simulations, we compute the hopping rates of all possible hops of all charge carriers and make a hopping list. The next hop is randomly selected from the hopping list. We find the hopping time, ${\tau }_k$, for the kth step in the KMC simulation from

$${\tau }_k=\frac{-ln\left(X\right)}{{\sum }_{i,j,i\ne j}^N{\nu }_{\mathit{ij}}{L}_i},$$E3

where N represents the total number of sites, $L_i$ is equal to unity if site i is occupied by a charge carrier and zero otherwise, and $X$ is a random number uniformly distributed between 0 and 1 (note that the ${\nu }_{\mathit{ij}}$ values are the rates of possible hops updated at each time step). We run each simulation until all charge carriers recombine.

To study the effect of DoS on the kinetics of BMR, we consider the diffusion of electrons and holes in energetically disordered systems, where we separately describe the energetic spectrum of the system by (i) a Gaussian distribution of states, centered at $E=0$, with a total DoS concentration $N_G$ and a distribution width ${\delta }_G$:

$${\rho }_G\left(E\right)=\frac{N_G}{\sqrt{2{\mathit{\pi \delta }}_G^2}}exp\left(\frac{-E^2}{2{\delta }_G^2}\right),$$E4

and (ii) an exponential distribution of states:

$${\rho }_E\left(E\right)=\frac{N_E}{{\delta }_E}exp\left(\frac{E}{{\delta }_E}\right),E\le 0,$$E5

where the DoS concentration and distribution width are given by $N_E$ and ${\delta }_E$, respectively. We also study the BMR kinetics in the presence of a composite DoS for which we superpose a Gaussian and an exponential distribution. To generate an energetic disorder in the KMC simulations, a random Gaussian or exponentially distributed value taken from Eq. (4) or (5), respectively, is assigned to each lattice site in the beginning of each simulation. We start each simulation with a random distribution of equal numbers of electrons and holes, whose initial densities are given by ${\rho }_0$, within the A and D domains. We allow single occupancy on each lattice site only.

3. Results and discussions

It has been shown experimentally that morphology of the active layer is indeed one of the primary factors affecting bimolecular recombination and so the device efficiency [49, 50, 51]. To demonstrate the role of morphology in BMR kinetics, we apply our reaction-diffusion model to four different morphologies. Random morphologies are obtained using the potential model according to Eq. (1). These morphologies are depicted in Figure 1A–D where hollow and red regions represent D and A domains, respectively. Figure 1A and B represents a bilayer solar cell with flat and rough interfaces, correspondingly. The structures in Figure 1C and D are representative of BHJs with large and small domain sizes of D and A, generated using Eq. (1) with λ =100 nm and 20 nm (with domains of size 50 and 10 nm), respectively. Figure 1E illustrates an initial distribution of 600 electrons and holes illustrated by red and blue dots, respectively, in the simulation box of 100 × 100 × 100 lattice sites representing a charge density of 6 × 10¹⁷ cm⁻³ in a BHJ device. In Figure 1E the domain size of D and A is smaller than the lattice size (1 nm); thus, for this structure, charge carriers are allowed to occupy any lattice sites during diffusion.

Figure 1.

Regardless of morphology, in a bilayer with (A) flat interface and (B) rough interface, a BHJ with (C) coarse and (D) fine donor (hollow regions)-acceptor (red regions) network, and a BHJ with D/A cluster size of less than 1 nm (E) long-time behavior of charge carriers during BMR follows Langevin dynamics ($\rho \sim t^{-1}$), i.e., variation in D/A interfacial area does not change BMR exponent. Solid lines are from KMC simulations, and color dashed lines are analytical form (F).

Figure 1F compares the KMC simulation results (solid lines) with analytical results (color dashed lines) for the decay of the charge density, $\rho$, with time, $t$, using morphologies in Figure 1A–E. Details of the analytical form of charge carrier density decay are given in Appendix B. In these simulations, we assume an isoenergetic energy spectrum, where DoS is theoretically a delta function peaked at E = 0, to consider solely the effect of morphology and the interfacial roughness on BMR kinetics. According to Figure 1F, density decay is only delayed as the D/A interfacial area increases from Figure 1A–E. The quantities of interest are the impact of morphology and interfacial roughness on BMR exponent, α, which is defined by

$$\rho \left(t\right)\sim t^{-\alpha },\mathrm{as}t\to \infty ,$$E6

and charge carrier lifetime. Figure 1F suggests that any variation in D/A interfacial area via morphological change does not influence BMR exponent. This is clearly seen from the charge density decay curves in Figure 1F at long time, exhibiting Langevin dynamics ($\rho \sim t^{-1}$) regardless of the morphology.

Details of how to obtain charge carrier lifetime from density decay are given in Appendix B (Eq. (B7)). As shown in Figure 1F, the time it takes to reach a certain density decreases from morphology 1A–E, suggesting that charge carrier lifetime decreases with the D/A interfacial area from Figure 1A–E. Time shifts in the density decay illustrated in Figure 1F indicate that each sequential increase in the D/A interfacial area from the morphology in Figure 1A–E, which corresponds to a decrease in D and A average domain sizes, enhances the charge carrier lifetime.

Figure 2 illustrates the effect of morphology shown in Figure 1A–E on charge carrier lifetime, $\tau$ (see Appendix B for details). Figure 2A illustrates that morphology only affects charge carrier lifetime at short time, whereas at long time, charge carrier lifetime linearly increases with time, independent of the morphology. This long-time behavior demonstrates the charge carrier diffusion-limited recombination, which does not depend on the morphology. Figure 2B demonstrates that for a fixed charge carrier density, charge carrier lifetime always decreases with D/A interfacial area. This decrease can be as large as two orders of magnitude between morphologies in Figure 1A and E. Figure 2B also shows that in isoenergetic DoS, charge carrier lifetime decreases with carrier density, proportional to ${\rho }^{-1/\alpha }$, independent of the morphology.

Figure 2.

Variation of charge carrier lifetime with (A) time and (B) charge carrier density for different morphologies shown in Figure 1A–E.

We now look at how BMR rate is affected by morphology and D/A interfacial area. Utilizing the analytical form (see Appendix B) which approximates KMC simulation results for charge density decay shown in Figure 1F, we find that BMR rate at short time is intimately correlated with the D/A domain size. As shown in Figure 3A, at short time, BMR rate is inversely proportional to the D/A interfacial area and domain size, i.e., BMR rate is the highest for morphology in Figure 1E and the lowest for the morphology in Figure 1A. However, at long time, BMR rate becomes proportional to the D/A interfacial area and domain size, meaning that BMR rate is the highest for morphology in Figure 1A and the lowest for the morphology in Figure 1E. Results shown in Figure 3A suggest that BMR rate at long time in Figure 1A can be two orders of magnitude larger than that in Figure 1E. Such a difference in BMR rate will be explained in the following paragraphs.

Figure 3.

Variation of BMR rate with (A) time and (B) charge carrier density for different morphologies shown in Figuere 1A–E.

Initially, due to very close proximity of counter charge carriers at short time when charge carrier density is high, BMR rate is large for morphologies with small domain size. However, at longer time, when charge carrier density decreases, morphologies with small domain size restrict the motion of charge carriers (holes within D and electrons within A domains) more than those with large domain size. This restriction lowers BMR rate in morphologies with smaller domain size, and so diffusion-limited recombination will be the dominant effect at low charge carrier density.

The abovementioned BMR kinetics is only observed when we look at the effect of BMR on BMR rate in time. However, if we express BMR rate in terms of charge carrier density, as shown in Figure 3B, for a given carrier density, morphologies with small domain size (large D/A interfacial area shown in Figure 1E) always exhibit larger BMR rate than those with large domain size (small D/A interfacial area depicted in Figure 1A). This trend in BMR rate is consistently observed regardless of charge carrier density shown in Figure 3B.

To justify that variation in morphology and D/A interfacial area of an organic solar cell does not affect BMR rate at long time, we look at charge density decay using the morphologies in Figure 1A–E but now in the presence of two different DoS, Gaussian and exponential distributions. Figure 4 compares KMC simulation results of charge density decay due to BMR in the presence of DoS with energies of all lattice sites taken from a Gaussian distribution with ${\delta }_G=3k_{B}T$ (Figure 4A) and an exponential distribution with ${\delta }_E=3k_{B}T$ (Figure 4B). From KMC simulation results in Figure 4, we observe that for a given DoS distribution and different morphologies, charge density decay curves stay parallel at long time. Our observation suggests that morphology impacts the delay in BMR kinetics by shifting the density decay curves in time as the D/A interfacial area increases, whereas DoS affects the BMR rate (${{\rho }_0}^{-1}\mathit{d\rho }\left(t\right)/\mathit{dt}$) and BMR exponent in particular.

Figure 4.

Impact of DoS, with (A) Gaussian distribution using $N_G={10}^{21}$ cm⁻³ and ${\delta }_G=3k_{B}T$ and (B) exponential distribution using $N_E={10}^{21}$ cm⁻³ and ${\delta }_E=3k_{B}T$, on BMR kinetics. Different morphologies in Figure 1A–E, corresponding to red to blue, only shift the density decay in time, whereas DoS impacts the BMR exponent and thus rate at long time.

Figure 5 compares the effect of DoS distribution at low DoS concentration on the charge carrier density decay during BMR. Figure 5A and B illustrates charge density decay in the presence of Gaussian DoS distribution with $N_G={10}^{19}$ cm⁻³ and exponential distribution with $N_E={10}^{19}$ cm, respectively. For these simulations, 1% of lattice site energies are taken from Gaussian or exponential distributions with different distribution widths ranging from $1k_{B}T$ to $4k_{B}T$. The energies for the rest of lattice sites are set to $E=0$ (isoenergetic spectrum). Our KMC simulation results of BMR illustrated in Figure 5A in the presence of Gaussian DoS demonstrate that increasing the distribution width, from ${\delta }_G=1k_{B}T$ to $4k_{B}T$, only delays BMR, while it does not affect BMR exponent. These results suggest that the behavior of charge carriers at long time within a Gaussian DoS follows Langevin dynamics, $\rho \sim t^{-1}$, as shown by long-time asymptotes in dashed lines.

Figure 5.

KMC simulations of density decay due to BMR in 1% DoS ($N_G=N_E={10}^{19}$ cm⁻³) in an otherwise isoenergetic spectrum with a random morphology generated by potential model ($\lambda =10$ nm, see inset in panel A) using (A) Gaussian and (B) exponential DoS. For Gaussian DoS, the results manifest a constant BMR exponent (−1) and thus decay rate at long time regardless of distribution width. However, exponential DoS exhibits BMR exponent (varying from −1 to −0.25) and decay rates depending on the width of the distribution.

In contrast to Gaussian DoS, 1% exponential DoS (low concentration) in an otherwise isoenergetic spectrum not only delays BMR but also decreases BMR exponent. This is clearly demonstrated in Figure 5B, where KMC simulation results (solid lines) agree with the theoretical long-time asymptotic predictions, $\rho \sim t^{-\left(T/T_0\right)}$ [52, 53], $T_0={\delta }_E/k_B$, illustrated by dashed lines for the same range of distribution width, ${\delta }_E$, and DoS concentration, $N_E={10}^{19}$ cm⁻³. Comparing Figure 5A and B suggests that BMR kinetics is highly influenced by exponential DoS compared to Gaussian DoS. Specifically, $3 k_{B}T$ variation in distribution width slows BMR by 10 orders of magnitude to reduce carrier density to 0.001 of initial carrier concentration if exponential DoS is present (compare abscissa in Figure 5B for ${\delta }_E=1k_{B}T$ and $4k_{B}T$). This is comparable to the Gaussian DoS in which BMR is only slowed to about two orders of magnitude by the same variation in the distribution width (compare abscissa in Figure 5A for ${\delta }_G=1k_{B}T$ and $4k_{B}T$). It is noteworthy that results in Figure 5 are rendered when DoS concentration is low. Next, we will consider cases where DoS concentration is high.

Looking at high DoS concentrations (100% of energies of lattice sites taken from a DoS distribution) shown in Figure 6 reveals that Gaussian DoS ($N_G={10}^{21}$ cm⁻³) does not change BMR exponent, leading to Langevin behavior of charge carriers at long time ($\rho \sim t^{-1}$). However, when exponential DoS with high concentration ($N_E={10}^{21}$ cm⁻³) is employed, BMR exponent, BMR rate, and the long-time kinetics of charge carriers during BMR are highly influenced by the width of the distribution. Particularly, there is a significant change in the BMR exponent and thus the rate for ${\delta }_E>1k_{B}T$ as shown in Figure 6B. Comparing Figures 5 and 6 suggest that DoS concentration only delays the BMR, whereas DoS distribution (Gaussian or exponential) affects the BMR exponent and rate at long time.

Figure 6.

KMC simulation results (solid lines) of density decay of charge carriers in a BHJ with a random morphology generated by potential model (see Figure 5A inset) using (A) Gaussian and (B) exponential DoS with high concentration ($N_G=N_E={10}^{21}$ cm⁻³) and different distribution widths varying from $1k_{B}T$ to $4k_{B}T$. Dashed lines indicate asymptotes at long time.

To further study the effect of DoS concentration on the BMR kinetics, we now look at different DoS concentrations for two cases: a Gaussian and an exponential DoS distribution with fixed distribution width ${\delta }_G={\delta }_E=3k_{B}T$. Figure 7 compares the impact of DoS concentration on BMR kinetics. For a Gaussian DoS, increasing the DoS concentration, $N_G$, from $0.1\times {10}^{20}$ to $9\times {10}^{20}$ cm⁻³ only delays BMR (shift the density decay curve in time) without any impact on the BMR exponent and rate as shown in Figure 7A. A similar trend is observed when an exponential DoS is utilized as shown in Figure 7B. Note that in our KMC simulations, BMR exponent, which is constant as shown in Figure 7A and B, at long time only depends on the width of the DoS distribution for exponential DoS but is independent of the width for Gaussian DoS. The BMR exponent and thus the rate can be readily obtained from asymptotic dashed lines shown in Figure 7A and B. These results indicate that DoS concentration only delays BMR and does not affect BMR exponent and rate at long time.

Figure 7.

KMC simulations of charge density decay (solid lines) due to BMR for different DoS concentrations with (A) Gaussian DoS, $N_G=0.1\times {10}^{20}$, $1\times {10}^{20}$, $3\times {10}^{20}$, and $9\times {10}^{20}$ cm⁻³ and (B) exponential DoS, $N_E=0.1\times {10}^{18}$, $1\times {10}^{18}$, $3\times {10}^{18}$, and $9\times {10}^{18}$ cm⁻³, in an otherwise isoenergetic spectrum. Dashed lines illustrate long-time asymptotes to the density curves.

We now turn to the effect of initial charge carrier density on the kinetics of BMR. Figure 8 illustrates the KMC simulation results of charge density decay due to BMR for different initial charge carrier concentrations, ranging from ${\rho }_0={10}^{17}$ cm⁻³ to $27\times {10}^{17}$ cm⁻³, from red to blue, respectively, for isoenergetic DoS with $E=0$ (Figure 8A), Gaussian DoS (Figure 8B), and exponential DoS (Figure 8C) with DoS concentration of ${10}^{21}$ cm⁻³ with $3k_{B}T$ distribution width. Results in Figure 8 reveal that charge carrier density at long-time approaches a unique asymptote regardless of the initial charge concentration. We observe that the only influence of initial charge carrier concentration is to delay the transition time at which BMR exponent and rate transition from nearly zero at very short time to the asymptotic value at long time. We find that this shift strongly correlates with the DoS distribution, which is larger in exponential (Figure 8C) than Gaussian (Figure 8B) DoS. Overall, the results illustrated in Figure 8 suggest that initial charge carrier density does not affect BMR exponent and rate. It also does not render a delay in BMR, as all charge density decay curves with different initial charge densities merge at long time.

Figure 8.

KMC simulation results of charge carrier density decay (solid lines) during BMR for different initial charge carrier densities ranging from ${\rho }_0={10}^{17}$ to $27\times {10}^{17}$ cm⁻³ with (A) isoenergetic DoS, (B) a Gaussian DoS with concentration $N_G={10}^{21}$ cm⁻³, and (C) an exponential DoS with concentration $N_E={10}^{21}$ cm⁻³ and fixed $3k_{B}T$ distribution width. Dashed lines are the long-time asymptotes. These results indicate that initial charge carrier concentration neither changes BMR exponent or rate nor delays BMR.

Thus far, we have looked at the kinetics of BMR in a BHJ device in the presence of single DoS distribution. In reality, however, DoS distributions may better be explained by a composite DoS, representing two distributions for band (conducting) and band-tail (trap) states. To account for the effect of a composite DoS, we now consider a DoS consisting of Gaussian band states and exponential band-tail states, illustrated in Figure 9.

Figure 9.

Energetic disorder described by superimposing two distributions, a Gaussian DoS (black), ${\rho }_G\left(E\right)$, centered at $E=0$, and an exponential DoS (blue), ${\rho }_E\left(E\right)$, positioned below $E=0$.

KMC simulation results for BMR and charge density decay in the presence of a composite DoS shown in Figure 9 are presented in Figure 10. Figure 10A illustrates the density decay (solid lines) during BMR using Gaussian band states with DoS concentration of $N_G=0.99\times {10}^{21}$ cm⁻³, for different distribution width of ${\delta }_G=$ 1, 2, 3, and 4 $k_{B}T$, corresponding to red to blue, respectively, and exponential band-tail states with DoS concentration of $N_E=0.01\times {10}^{21}$ cm⁻³ and distribution width of ${\delta }_E=3k_{B}T$. Results in Figure 10 demonstrate that as the Gaussian distribution of band states broadens, BMR kinetics is delayed, exhibiting the dominance of Langevin dynamics ($\rho \sim t^{-1}$) at intermediate time with BMR exponent of 1. However, the long-time BMR exponent is always governed by the width of band-tail state distribution ($\rho \sim t^{-033}$), which is 0.33, depicted by dashed line in Figure 10A.

Figure 10.

KMC simulation results of density decay due to BMR for a composite DoS shown in Figure 9, with (A) different ${\delta }_G$ varying from $1k_{B}T$ to $4k_{B}T$, ${\delta }_E=3k_{B}T$, $N_G=0.99\times {10}^{21}$ cm⁻³, $N_E=0.01\times {10}^{21}$ cm⁻³; (B) different ${\delta }_E$ varying from $1k_{B}T$ to $4k_{B}T$, ${\delta }_G=2k_{B}T$, with the same DoS concentration described in Figure 10A; and (C) band-tail state DoS concentration varying from $N_E=0.05\times {10}^{19}$ to $4\times {10}^{19}$ cm⁻³, subject to $N_G+N_E={10}^{21}$ cm⁻³, ${\delta }_G=2k_{B}T$, and ${\delta }_E=3k_{B}T$.

The dependence of long-time BMR density decay on the band-tail state distribution width can be further justified by Figure 10B, depicting the impact of exponential band-tail state distribution width on the BMR density decay. The width of distribution ranges from ${\delta }_E=$ 1 to 4 $k_{B}T$, and its DoS concentration is the same as that utilized in Figure 10A. Results in Figure 10B suggest that BMR decay rate is directly correlated with the width of band-tail state distribution. As the band-tail state distribution broadens (the band-tail state distribution width increases), BMR exponent and rate decrease. In this case, the charge density decay approaches that of analytical form for exponential DoS shown by dashed lines ($\rho \sim t^{-\left(T_0/T\right)}$, $T_0={\delta }_E/k_B$). Elevating the DoS concentration of band-tail state distribution (increasing DoS concentration of exponential distribution in the composite DoS shown in Figure 9), from $N_E=0.05\times {10}^{19}$ to $4\times {10}^{19}$ cm⁻³, with the total DoS concentration being constant, $N_G+N_E={10}^{21}$ cm⁻³, only delays the BMR (shifts charge density decay in time) without changing long-time BMR rate as shown in Figure 10C.

Overall, Figure 10 indicates that (i) Gaussian band states only widen the intermediate time window at which charge carrier behavior can be described by Langevin dynamics ($\rho \sim t^{-1}$), (ii) long-time behavior of charge carriers is governed by the distribution of band-tail states, and (iii) the ratio of band to band-tail state DoS concentrations only induces a time delay to the instant at which initial BMR regime with small BMR rate transitions to the long-time BMR regime with large BMR rate. Charge carriers’ behavior during this transition state is governed by band state DoS distribution.

Thus far, our theoretical treatment provides a detailed prediction of the physical behavior underlying BMR in bilayer and BHJ solar cells. Specifically, our theoretical model predicts the effect of morphology, DOS, trap (i.e., band-tail states) concentration on the BMR, and its rate in such devices. We now demonstrate that our theoretical model reproduces experimental data [25] taken from 170 nm thick 1:1 blend film of P3HT and PCBM at various light-pulse intensities at open-circuit voltage. Using the potential model in Eq. (1) with $\lambda =50$ nm, we generate random morphologies reminiscent of BHJs, as shown in the inset of Figure 11. To simulate various light-pulse intensities as shown in Figure 11, from top to bottom corresponding to 60, 24, 6, and 3.6 μ J cm⁻², respectively, we perform simulations with various initial charge densities. We utilize our theoretical model with initial charge densities of $n=6.3\times {10}^{17}$, $4\times {10}^{17}$, $1.3\times {10}^{17}$, and $0.58\times {10}^{17}$ cm⁻³, respectively. To match the experimental data, we apply the model to an energetically disordered D/A network with Gaussian $\left({\delta }_G=2k_{B}T\right)$ and exponential $\left({\delta }_E=2.2k_{B}T\right)$ DOS for band and band-tail states, respectively. We set the density of “conducting” (band) and “trapping” (band-tail) states, $N_G$ and $N_E$ to $0.97\times {10}^{21}$ cm⁻³ and $2.6\times {10}^{19}$ cm⁻³, respectively.

Figure 11.

Comparison of charge density decay with time from experimental data (black dots) at open-circuit voltage and various light-pulse intensities and from simulation results (solid lines) at various initial charge densities utilizing Gaussian band states and exponential band-tail states.

Figure 11 demonstrates that simulation results (solid lines) for decay of charge density, $\rho$, at open-circuit voltage are in reasonably good agreement with experimental data [25] (black dots) taken at various light-pulse intensities. The question that arises is whether the superposition of two DOS distributions for band and band-tail states is necessary to reproduce experimental data shown in Figure 11. We attempted to fit experimental data using only exponential distribution for band-tail states with parameters ${\nu }_0=5.9\times {10}^{11}$ s^–1, $N_E=2.3\times {10}^{19}$ cm⁻³, and ${\delta }_E=2.2k_{B}T$. Our simulation results (not shown) suggest that using only exponentially distributed band-tail states requires initial charge concentration $n={10}^{18}$, $6\times {10}^{17}$, $1.4\times {10}^{17}$, and $0.65\times {10}^{17}$ cm⁻³, from high to low light intensities, respectively, which seems to be systematically higher than those with two DOS distributions. For instance, given the simulation results, we suspect that charge concentration at the highest light-pulse intensity, ${10}^{18}$ cm⁻³, is beyond the physical values; thus, making it difficult to justify that only exponential distribution for band-tail states without the need for band states is enough.

Altogether, our results, illustrating the impact of DoS on BMR exponent and rate, demonstrate that for narrow exponential distribution of DoS (${\delta }_E\le k_{B}T$), charges move diffusively; therefore, recombination exhibits a Langevin-type dynamics. As distribution becomes broader, electrons and holes transition from a diffusive to sub-diffusive motion, leading to a sub-diffusion-limited reaction process.

Sub-diffusive motion of polarons comes from anomalous diffusion in disordered media. This is due to long waiting times occurring during diffusion through traps. There are physical mechanisms of sub-diffusion including random walks in dynamically disordered medium. Essentially, sub-diffusive motions and thus sub-diffusive limited reactions (electron-hole collision) arise when the waiting time between hops follows a power-law distribution with finite variance at long time [54, 55, 56, 57]. This is true provided that reactants are distributed homogenously throughout the medium. Given this condition, utilizing thermally activated or Miller-Abrahams hopping mechanisms, one will find (see Appendix A) that reacting random walkers, such as electrons and holes, in quench disordered exponentially distributed DoS sampled from a power-law hopping time distribution, exhibit a sub-diffusive motion. Consequently, if electrons and holes move sub-diffusively in a medium and are allowed to recombine, they manifest a sub-diffusive limited reaction. This is reminiscent to what is observed in transient absorption spectroscopy and transient photovoltage experimental data [25] taken from P3HT/PCBM BHJ device.

Another mechanism by which sub-diffusive limited reactions can occur is where reactants are inhomogenously dispersed; thus there will be local relative density fluctuations (density difference) in the medium. The role of local relative density fluctuations on the charge density decay can be better understood by invoking of binary reactions with two species. In Langevin dynamics, the relative density is zero due to the uniform density assumption. However, because of initial fluctuations in the relative density of the reactants, experimental observations [58, 59, 60, 61], theoretical predictions, numerical simulations [62, 63, 64, 65, 66, 67, 68] of reaction-diffusion of two species with periodic boundary conditions in one and two dimensions [63, 69], and reflective boundary conditions in three dimensions [65] demonstrate deviation from Langevin dynamics. These relative density fluctuations lead to the density decay of species with a power law given by Ref. [69]

$$\rho \sim t^{-\frac{d_f}{2d_w}},$$E7

where d_f is the fractal dimension of the medium and d_w is the fractal dimension of the random walk. Here, it is assumed that charges diffuse in the medium similar to a random walk.

Similar to binary reactions, in random distribution of electrons and holes in BHJs, the uniform charge density assumption may not be valid, especially in a D/A network, where initial charge separation is indicative of non-zero local relative density fluctuations. To find the extent of influence of relative density fluctuations in BHJs on BMR, we perform KMC simulations. We employ our theoretical model in a cubic cell of size $200\times 200\times 200$ with morphology given in Figure 1D. Our simulation results (solid lines) shown in Figure 12 indicate deviation from Langevin dynamics at charge concentration greater than $4\times {10}^{18}$ cm⁻³. Transition to Langevin-type behavior occurs at about 10⁻⁶ s at which point the density starts to decay as $t^{-1}$. Such deviation is not observed for lower initial charge densities. We also find that our observation is consistent with previous calculations for binary reactions [64].

Figure 12.

KMC simulation results (solid lines) illustrating the effect of local relative density fluctuation on the dynamics of bimolecular recombination applied to BHJ morphology shown in Figure 1D. At high density ($\ge 4\times {10}^{18}$ cm⁻³), deviation from Langevin dynamics is apparent with charge density decay ∼t^−3/4. Effect of local relative density fluctuations is not present in low initial densities ($<4\times {10}^{18}$ cm⁻³).

4. Conclusions

Developing a reaction-diffusion model, we have studied the effect of morphology, DoS, charge carrier density, and local charge density fluctuations on BMR kinetics in OPV cells at short and long time scales. We have looked at single as well as composite DoS that consists of a superposition of Gaussian and exponential distributions. In all of our KMC simulations, we always find two regimes during BMR: a slow charge density decay (small BMR rate) followed by a fast electron-hole annihilation (large BMR rate). Our KMC simulation results obtained from bilayer as well as BHJs with random morphology indicate that D/A interfacial roughness and morphology in organic solar cells only affect the total BMR time and do not influence BMR exponent and rate at short or long time scales. Our simulation results also illustrate that charge carrier concentrations in BHJ devices neither change BMR rate nor time. It only shifts the time at which transition from slow to fast BMR occurs.

Furthermore, we find that the effect of morphology, interfacial roughness, and charge carrier density is in contrast to the effect of DoS, whose characteristics can completely modify BMR kinetics and rate. We find that kinetics of charge carriers due to BMR in a Gaussian DoS can be explained by Langevin dynamics, demonstrating a density decay inversely proportional to time. However, dynamics of charge carriers in the presence of exponential DoS deviates from Langevin dynamics, illustrating a slower dynamics for distribution width higher than the thermal energy. Particularly, our KMC simulations show that the BMR rate at long time is controlled by the width of DoS distribution and that the DoS concentration changes BMR time without changing the BMR rates at short or long time. Extending our findings from single to composite DoS indicates that the BMR rate is controlled by the characteristics of trap state distribution and that the ratio of band to tail state concentration dictates the BMR time, specifically the transition time from slow to fast BMR regime.

We also demonstrate that decrease in recombination rate can be due to local fluctuation in the electron and hole density difference and sub-diffusive motion of polarons. Sub-diffusivity arises from a power-law distribution of jump time. Utilizing Miller-Abrahams charge transfer relation, we find that such power-law distribution is present in the hopping mechanism provided by our theoretical model. The power-law distribution necessitates an exponentially distributed dynamically disordered system, to which Miller-Abrahams or thermally activated hopping is applied. We also find that existence of a quenched material is not necessary for a sub-diffusive motion of charges. Altogether, we find a very rich predicted behavior on which bimolecular loss mechanisms can depend.

Future work will address the effect of positional disorder and its coupling effect with energetic disorder on charge transfer rate and transport properties of OPV cells using hopping mechanisms and lattice model.

Acknowledgments

This work has been supported by the seed funding from the College of Engineering at University of Illinois, Chicago.

The long-time behavior of charge carriers in bulk heterojunction OPVs is approximated by the behavior of reacting continuous time random walk (CTRW) in Euclidean medium, representing a sequence of jump-trap-release events. For reacting CTRW, $\rho \left(t\right)\sim R\left(t\right)$, where $\rho \left(t\right)$ and $R\left(t\right)$ are the carrier density and survival probability at time $t$, respectively. It turns out that for a reacting CTRW, $R\left(t\right)\sim 1/S\left(t\right)$ for three dimensions [70], where $S\left(t\right)$ is the number of distinct sites visited by a CTRW. $S\left(t\right)$ can be visualized as the volume that a carrier sweeps during diffusion. One may also note that in regular lattices, $S\left(t\right)\sim ⟨\overrightarrow{r}{\mathrm{}}^2\left(t\right)⟩$ [67], where $⟨\overrightarrow{r}{\mathrm{}}^2\left(t\right)⟩$ is the mean square displacement (MSD) of CTRW at time $t$. From the above argument, one concludes that

$$\rho \left(t\right)\sim \frac{1}{⟨\overrightarrow{r}{\mathrm{}}^2\left(t\right)⟩}.$$EA1

Without loss of generality, in the one-dimensional analogue of CTRW model, we describe the transport in terms of succession of jumps of length x, which is drawn from a probability distribution function (PDF), $\lambda \left(x\right)$, followed by a trapping event at the same position x for time interval t drawn from the PDF, $w\left(x\right)$. In this model, the probability that a carrier is found at position $x$ in space at time $t$ is given by propagator $W\left(x,t\right)$, whose Fourier-Laplace transform obeys the relation [54]

$$W\left(k,u\right)=\frac{1-w\left(u\right)}{u}\frac{W_0\left(k\right)}{1-w\left(u\right)\lambda \left(k\right)},$$EA2

where k and u represent Fourier and Laplace variables, respectively, $W_0\left(k\right)$ indicates the Fourier transform of the initial condition $W_0\left(x\right)$, and $w\left(u\right)$ and $\lambda \left(k\right)$ are the Laplace and Fourier transform of $w\left(t\right)$ and $\lambda \left(x\right)$, respectively.

In our KMC simulations, jump length is equal to the lattice spacing; thus for small k values.

$$\lambda \left(x\right)=\delta \left(x-a\right)\text{and}\lambda \left(k\right)=exp\left(-\mathit{ika}\right)=1-a^2{k}^2+O\left(k^4\right),$$EA3

where $\delta \left(x\right)$ is the delta function. Note that in our reaction-diffusion model, waiting time PDF is governed by the DoS distribution. Here, we consider Gaussian, uniform, and exponential DoS as follows.

Gaussian DoS. Using Eq. (4), a thermally activated process [53] with hopping frequency from a site with energy E to any other site is expressed by

$$\nu ={\nu }_{0}exp\left(\frac{E}{k_{B}T}\right),E\le 0,$$EA4

where ${\nu }_0$ is the intrinsic attempt frequency. Since the jump time, $\tau$, is inversely proportional to the hopping frequency; thus,

$$\tau =\frac{1}{{\nu }_0}exp\left(-\frac{E}{k_{B}T}\right)$$EA5

The randomness in jump time comes from the randomness in energies sampled from DoS distribution. Thus, finding E from Eq. (A5) and its derivative with respect to $\tau$, one arrives at the PDF of the jump time

$$w\left(\tau \right)\mathit{d\tau }=-{\rho }_G\left(E\right)\mathit{dE}$$EA6

with a finite mean jump time $\overline{\tau }$ equal to

$$\begin{array}{c}\overline{\tau }=\underset{0}{\overset{\infty }{\int }}\tau w\left(\tau \right)\mathit{d\tau }=\underset{-\infty }{\overset{0}{\int }}{\rho }_G\left(E\right)\tau \mathit{dE}\\ =\frac{1}{\sqrt{2\pi }{\delta }_G{\nu }_0}\underset{-\infty }{\overset{0}{\int }}exp\left(\frac{-E^2}{2{\delta }_G^2}-\frac{E}{k_{B}T}\right)\mathit{dE}\\ =\frac{1}{\sqrt{2\pi }{\delta }_G{\nu }_0}exp\left[\frac{2{\delta }_G^2}{{\left(k_{B}T\right)}^2}\right]\underset{-\infty }{\overset{0}{\int }}exp\left[-{\left(\frac{E}{\sqrt{2}{\delta }_G}+\frac{\sqrt{2}{\delta }_G}{k_{B}T}\right)}^2\right]\mathit{dE}\\ =\frac{1}{2{\nu }_0}\left[1+erf\left(\frac{\sqrt{2}{\delta }_G}{k_{B}T}\right)\right].\end{array}$$EA7

where $w\left(t\right)=\gamma exp\left\{-{\gamma }^2{\left[\mathit{log}\left({\nu }_0\tau \right)\right]}^2/2\right\}/\left(\sqrt{2\pi }\tau \right)$ and $\gamma =k_{B}T/{\delta }_G$. For any DoS with finite mean jump time, using Taylor series for small u, one can write

$$w\left(u\right)=1-u\overline{\tau }+O\left({\overline{\tau }}^2\right)$$EA8

With initial condition $W_0\left(x\right)=\delta \left(x\right)$, Eqs. (A2), (A3), (A7), and (A8) imply

$$W\left(k,u\right)=\frac{1}{u+{\mathit{Dk}}^2}$$EA9

where $D=a^2/\overline{\tau }$ is the diffusion constant. Taking the inverse Fourier-Laplace transform from Eq. (A9), we arrive at the diffusion equation

$$\frac{\partial W\left(x,t\right)}{\partial t}=D\frac{{\partial }^{2}W\left(x,t\right)}{\partial x^2}$$EA10

with initial condition $W\left(x,0\right)=W_0\left(x\right)$. Solution to Eq. (A10) is the well-known Gaussian propagator

$$W\left(x,t\right)=\frac{1}{\sqrt{4\mathit{\pi Dt}}}exp\left(-\frac{x^2}{4\mathit{Dt}}\right)$$EA11

Note that the long-time form of the propagator was considered through the assumption of k being very small. Recall that for a Gaussian propagator, MSD is given by

$$⟨\overrightarrow{r}{\mathrm{}}^2\left(t\right)⟩\sim t.$$EA12

Using Eqs. (A1) and (A12), one finds that for a Gaussian DoS in a thermally activated process, charge carrier density decays as

$$\rho \left(t\right)\sim t^{-1},$$EA13

which is indicative of Langevin dynamics. This relation holds true for any DoS distribution, whose jump time distribution has a finite mean value.

Uniform DoS. Plugging the uniform distribution

$${\rho }_U\left(E\right)=\frac{N_U}{{\delta }_U},-E_U\le E\le 0,$$EA14

into Eq. (A7) renders

$$\overline{\tau }=\frac{1}{E_U{\nu }_0}\underset{-E_U}{\overset{0}{\int }}exp\left(-\frac{E}{k_{B}T}\right)\mathit{dE}=\frac{k_{B}T}{E_U{\nu }_0}\left[exp\left(\frac{N_U}{k_{B}T}\right)-1\right],$$EA15

where $w\left(t\right)=k_{B}T{\left(E_{U}t\right)}^{-1}$. Note that $\overline{\tau }$ is a finite value. Thus, following the same procedure performed for Gaussian DoS, we conclude that $\rho \left(t\right)\sim t^{-1}$.

Exponential DoS. For an exponential DoS given in Eq. (5), one can rewrite Eq. (A7)

$$\overline{\tau }={\mathit{\alpha \nu }}_0\underset{{\nu }_0^{-1}}{\overset{\infty }{\int }}{\left({\nu }_0\tau \right)}^{-\alpha }\mathit{d\tau }$$EA16

to yield $w\left(t\right)={\mathit{\alpha \nu }}_0{\left({\nu }_{0}t\right)}^{-\left(\alpha +1\right)}$, where $\alpha =k_{B}T/{\delta }_E$ and $0<\alpha <1$. Note that $\overline{\tau }$ diverges which is indicative of a dispersive process. The Laplace transform of $w\left(t\right)$ for small u is [71]

$$w\left(u\right)\approx 1-{\left(\frac{u}{{\nu }_0}\right)}^{\alpha }+O\left(u\right).$$EA17

Substitution of Eqs. (A3) and (A17) into Eq. (A2) leads to the propagator

$$W\left(k,u\right)=\frac{1}{u\left(1+D_{\alpha }{u}^{-\alpha }{k}^2\right)}$$EA18

where $D_{\alpha }=a^2/{\nu }_0^{\alpha }$. From Eq. (A18), MSD is found by inverse Laplace transform of $⟨x^2\left(u\right)⟩={lim}_{k\to 0}\left\{-d^{2}W\left(k,u\right)/{\mathit{dk}}^2\right\}$, which leads to

$$⟨x^2\left(t\right)⟩=\frac{2D_{\alpha }}{\mathrm{\Gamma }\left(\alpha +1\right)}{t}^{\alpha },$$EA19

where the condition of u being small has been applied. Eq. (A19) can be generalized to three dimensions, illustrating that $⟨\overrightarrow{r}{\mathrm{}}^2\left(t\right)⟩\sim t^{\alpha }$. Using this relation and Eq. (A1), we conclude that $\rho \left(t\right)\sim t^{-\alpha }$.

Assuming the general form for non-Langevin dynamics in BMR, we express the time variation of charge carrier density by

$$\frac{\mathit{d\rho }\left(t\right)}{\mathit{dt}}=-\beta \left(t\right)\rho {\left(t\right)}^2$$EB1

where $\rho \left(t\right)$ is the charge carrier density and $\beta \left(t\right)$ is the time-dependent rate coefficient [72]. Considering boundary condition $\rho ={\rho }_0$ at $t=0$, one finds that

$$\frac{\rho \left(t\right)}{{\rho }_0}={\left[1+{\rho }_0\underset{0}{\overset{t}{\int }}\beta \left(t\right)\mathit{dt}\right]}^{-1}.$$EB2

The choice of $\beta \left(t\right)$ is semiempirical. Note that the charge density decay follows a power law at long time ($\rho \left(t\right)\sim t^{-\alpha }$), so it has been suggested [72] that $\beta \left(t\right)={\beta }_0{t}^{\alpha -1},0<\alpha <1$, where ${\beta }_0$ is some constant. For this particular choice of $\beta \left(t\right)$, solving BMR Eq. (B1) leads to

$$\frac{\rho \left(t\right)}{{\rho }_0}={\left[1+{\left(\frac{t}{{\tau }_0}\right)}^{\alpha }\right]}^{-1},$$EB3

where ${\tau }_0={\left[\alpha /\left({\rho }_0{\beta }_0\right)\right]}^{1/\alpha }$ is called the effective bimolecular lifetime [72].

Although charge carrier lifetime can be seen as the mean survival time, which is directly related to the survival probability $S\left(t\right)$ [73, 74] via $\tau =\underset{0}{\overset{\infty }{\int }}S\left(t\right)\mathit{dt}$, calculating carrier lifetime using the concept of survival probability seems to be impractical as the abovementioned integral is diverging due to long-lived charge carriers trapped in the band-tail states. However, neglecting those long-lived carriers, we can intuitively find carrier lifetime using [75]

$$\tau =\frac{\rho \left(t\right)}{-\frac{\mathit{d\rho }\left(t\right)}{\mathit{dt}}},$$EB4

for which substituting Eq. (B3) into Eq. (B4) leads to $\tau ={\left({\rho }_0{\beta }_0\right)}^{-1}\left[1+{\left(t/{\tau }_0\right)}^{\alpha }\right]t^{1-\alpha }$. However, in this equation, as t tends to zero, $\mathit{d\rho }/\mathit{dt}$ diverges; thus, carrier lifetime diminishes, meaning that charge carriers annihilate immediately after charge separation. Such a fast recombination mechanism can also be seen for the other choice of time-dependent rate coefficient, $\beta \left(t\right)$, mentioned before.

To resolve this instantaneous recombination at short time, we modify the above time-dependent rate coefficient to $\beta \left(t\right)={\beta }_0{\left(1+t/\eta \right)}^{\alpha -1}$ where $\eta =\alpha /\left({\rho }_0{\beta }_0\right)$ is a time scaling constant. Experimental data [72] shows that $\eta$ is on the order of $\mathit{\mu s}$ at room temperature for MDMO-PPV and PCBM blend. For this choice of $\beta \left(t\right)$, from Eq. B1, one finds

$$\frac{\rho \left(t\right)}{{\rho }_0}={\left(1+\frac{t}{\eta }\right)}^{-\alpha },$$EB5

for which BMR rate can be defined as

$$\frac{1}{{\rho }_0}\frac{\mathit{d\rho }\left(t\right)}{\mathit{dt}}=-\frac{\alpha }{\eta }{\left(1+\frac{t}{\eta }\right)}^{-\alpha -1}=-\frac{\alpha }{\eta }{\left(\frac{\rho }{{\rho }_0}\right)}^{1+\frac{1}{\alpha }}.$$EB6

Using Eqs. (B5) and (B6), one can express carrier lifetime in Eq. (B4) in terms of the carrier density

$$\tau =\frac{\eta }{\alpha }\left(1+\frac{t}{\eta }\right)=\frac{\eta }{\alpha }{\left(\frac{\rho \left(t\right)}{{\rho }_0}\right)}^{-\frac{1}{\alpha }},$$EB7

which can very well be approximated by $\tau \approx t/\alpha$ at long time. In contrast to the density decay in Eq. (B3), Eq. (B5) leads to a finite value for $\mathit{d\rho }/\mathit{dt}$ in Eq. (B6) as t tends to zero, giving rise to a non-zero initial carrier lifetime $\tau =\eta /\alpha$. Eq. (B7) also demonstrates that carrier lifetime decreases with carrier density (see Figure 2) which suggests that ultimately, long-lived carriers in the deep band-tail states dominate the carrier lifetime.

Figure B1 illustrates the results taken from analytical form in Eq. (B6), indicating the impact of BMR exponent on BMR rate for a fixed $\eta ={10}^{-6}{\mathrm{s}}^{-1}$. Figure B1A indicates that there is a transition at $t=\eta$ where BMR rate transitions from a constant to a rate decreasing with time. Furthermore, the effect of BMR exponent on the BMR rate is reversed at this point, suggesting that BMR rate increases with BMR exponent at short time but decreases at long time. Figure B1B demonstrates that at a given charge carrier density, as BMR exponent decreases (e.g., band-tail state distribution becomes deeper), BMR rate also decreases.

Figure B1.

Impact of BMR exponent, $\alpha$, on BMR rate with (A) time and (B) charge carrier density utilizing Eq. (B6). In panel A, vertical dotted line corresponds to $t=\eta ={10}^{-6}{\mathrm{s}}^{-1}$ using which these results were obtained.