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

### Keywords

- organic photovoltaics
- bimolecular recombination
- Langevin dynamics
- bulk heterojunction
- reaction-diffusion
- kinetic Monte Carlo

## 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

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

where *M* and *M*, randomly splits the space into D domains wherever

The diffusion and recombination of charge carriers are simulated using the KMC technique where the charge carrier transfer rate *i* to one of the six nearest neighboring sites, *j*, is described by the Miller-Abrahams model [46]:

where *i*, in our simulations ^{−1} denotes the intrinsic attempt frequency, ^{−1} indicates the inverse of localization radius, *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, *k*th step in the KMC simulation from

where *N* represents the total number of sites, *i* is occupied by a charge carrier and zero otherwise, and

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

and (ii) an exponential distribution of states:

where the DoS concentration and distribution width are given by

## 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^{17} cm^{−3} 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 1F compares the KMC simulation results (solid lines) with analytical results (color dashed lines) for the decay of the charge density, *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

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 (

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,

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.

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

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 ^{−3} and exponential distribution with

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, ^{−3}. Comparing Figure 5A and B suggests that BMR kinetics is highly influenced by exponential DoS compared to Gaussian DoS. Specifically,

Looking at high DoS concentrations (100% of energies of lattice sites taken from a DoS distribution) shown in Figure 6 reveals that Gaussian DoS (^{−3}) does not change BMR exponent, leading to Langevin behavior of charge carriers at long time (^{−3}) 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

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 ^{−3} 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.

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 ^{−3} to ^{−3}, from red to blue, respectively, for isoenergetic DoS with ^{−3} with

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.

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 ^{−3}, for different distribution width of ^{−3} and distribution width of

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 ^{−3}, with the total DoS concentration being constant, ^{−3}, 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 (

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 ^{−2}, respectively, we perform simulations with various initial charge densities. We utilize our theoretical model with initial charge densities of ^{−3}, respectively. To match the experimental data, we apply the model to an energetically disordered D/A network with Gaussian ^{−3} and ^{−3}, respectively.

Figure 11 demonstrates that simulation results (solid lines) for decay of charge density, ^{–1}, ^{−3}, and ^{−3}, 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, ^{−3}, 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 (

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]

where *df* is the fractal dimension of the medium and *dw* 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 ^{−3}. Transition to Langevin-type behavior occurs at about 10^{−6} s at which point the density starts to decay as

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

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), *x* for time interval *t* drawn from the PDF,

where *k* and *u* represent Fourier and Laplace variables, respectively,

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

where

**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

where

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

with a finite mean jump time

where *u*, one can write

With initial condition

where

with initial condition

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

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

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

into Eq. (A7) renders

where

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

to yield *u* is [71]

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

where

where the condition of *u* being small has been applied. Eq. (A19) can be generalized to three dimensions, illustrating that

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

where

The choice of

where

Although charge carrier lifetime can be seen as the mean survival time, which is directly related to the survival probability

for which substituting Eq. (B3) into Eq. (B4) leads to *t* tends to zero,

To resolve this instantaneous recombination at short time, we modify the above time-dependent rate coefficient to

for which BMR rate can be defined as

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

which can very well be approximated by *t* tends to zero, giving rise to a non-zero initial 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