Open access peer-reviewed chapter

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

By Daniel Christiansen and Shafigh Mehraeen

Submitted: October 8th 2018Reviewed: February 8th 2019Published: April 26th 2019

DOI: 10.5772/intechopen.85074

Downloaded: 249


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.


  • 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 100×100×100sites 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 ψxyz. ψis a random scalar function represented as a sum of several propagating plane waves in random directions given by


where ki1=cosαisinφi, ki2=sinαisinφi, ki3=cosφi, M=25, αiand θiare uniformly distributed random numbers between 02π, and φiis another uniformly distributed random number between 0π. λin Eq. (1) is directly proportional to the average size of D and A domains. ψ, whose value fluctuates between –M and M, randomly splits the space into D domains wherever ψ0and A domains wherever ψ<0. In this way, ψ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 νijfrom any site i to one of the six nearest neighboring sites, j, is described by the Miller-Abrahams model [46]:


where rij=rjri, rij=rij, and riare the position vectors of site i, in our simulations rij=a=1nm represents the lattice site spacing, ν0=7×1012s−1 denotes the intrinsic attempt frequency, γ=3×107cm−1 indicates the inverse of localization radius, β=kBT1is the Boltzmann constant, and ΔEij=EjEiis 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, τk, for the kth step in the KMC simulation from


where N represents the total number of sites, Liis equal to unity if site i is occupied by a charge carrier and zero otherwise, and Xis a random number uniformly distributed between 0 and 1 (note that the νijvalues 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 NGand a distribution width δG:


and (ii) an exponential distribution of states:


where the DoS concentration and distribution width are given by NEand δ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 ρ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 × 1017 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 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 (ρ∼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, ρ, 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


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 (ρt1) 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, τ(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 ρ1/α, 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 δG=3kBT(Figure 4A) and an exponential distribution with δE=3kBT(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 (ρ01t/dt) and BMR exponent in particular.

Figure 4.

Impact of DoS, with (A) Gaussian distribution using NG=1021 cm−3 and δG=3kBT and (B) exponential distribution using NE=1021 cm−3 and δE=3kBT, 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 NG=1019cm−3 and exponential distribution with NE=1019cm, respectively. For these simulations, 1% of lattice site energies are taken from Gaussian or exponential distributions with different distribution widths ranging from 1kBTto 4kBT. 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 δG=1kBTto 4kBT, 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, ρt1, as shown by long-time asymptotes in dashed lines.

Figure 5.

KMC simulations of density decay due to BMR in 1% DoS (NG=NE=1019 cm−3) in an otherwise isoenergetic spectrum with a random morphology generated by potential model (λ=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, ρtT/T0[52, 53], T0=δE/kB, illustrated by dashed lines for the same range of distribution width, δE, and DoS concentration, NE=1019cm−3. Comparing Figure 5A and B suggests that BMR kinetics is highly influenced by exponential DoS compared to Gaussian DoS. Specifically, 3 kBTvariation 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 δE=1kBTand 4kBT). 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 δG=1kBTand 4kBT). 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 (NG=1021cm−3) does not change BMR exponent, leading to Langevin behavior of charge carriers at long time (ρt1). However, when exponential DoS with high concentration (NE=1021cm−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 δE>1kBTas 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 (NG=NE=1021 cm−3) and different distribution widths varying from 1kBT to 4kBT. 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 δG=δE=3kBT. Figure 7 compares the impact of DoS concentration on BMR kinetics. For a Gaussian DoS, increasing the DoS concentration, NG, from 0.1×1020to 9×1020cm−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.

Figure 7.

KMC simulations of charge density decay (solid lines) due to BMR for different DoS concentrations with (A) Gaussian DoS, NG=0.1×1020, 1×1020, 3×1020, and 9×1020 cm−3 and (B) exponential DoS, NE=0.1×1018, 1×1018, 3×1018, and 9×1018 cm−3, 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 ρ0=1017cm−3 to 27×1017cm−3, 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 1021cm−3 with 3kBTdistribution 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 ρ0=1017 to 27×1017 cm−3 with (A) isoenergetic DoS, (B) a Gaussian DoS with concentration NG=1021 cm−3, and (C) an exponential DoS with concentration NE=1021 cm−3 and fixed 3kBT 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), ρGE, centered at E=0, and an exponential DoS (blue), ρEE, 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 NG=0.99×1021cm−3, for different distribution width of δG=1, 2, 3, and 4 kBT, corresponding to red to blue, respectively, and exponential band-tail states with DoS concentration of NE=0.01×1021cm−3 and distribution width of δE=3kBT. Results in Figure 10 demonstrate that as the Gaussian distribution of band states broadens, BMR kinetics is delayed, exhibiting the dominance of Langevin dynamics (ρt1) 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 (ρt033), 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 δG varying from 1kBT to 4kBT, δE=3kBT, NG=0.99×1021 cm−3, NE=0.01×1021 cm−3; (B) different δE varying from 1kBT to 4kBT, δG=2kBT, with the same DoS concentration described in Figure 10A; and (C) band-tail state DoS concentration varying from NE=0.05×1019 to 4×1019 cm−3, subject to NG+NE=1021 cm−3, δG=2kBT, and δE=3kBT.

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 δE=1 to 4 kBT, 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 (ρtT0/T, T0=δE/kB). Elevating the DoS concentration of band-tail state distribution (increasing DoS concentration of exponential distribution in the composite DoS shown in Figure 9), from NE=0.05×1019to 4×1019cm−3, with the total DoS concentration being constant, NG+NE=1021cm−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 (ρt1), (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 λ=50nm, 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−2, respectively, we perform simulations with various initial charge densities. We utilize our theoretical model with initial charge densities of n=6.3×1017, 4×1017, 1.3×1017, and 0.58×1017cm−3, respectively. To match the experimental data, we apply the model to an energetically disordered D/A network with Gaussian δG=2kBTand exponential δE=2.2kBTDOS for band and band-tail states, respectively. We set the density of “conducting” (band) and “trapping” (band-tail) states, NGand NEto 0.97×1021cm−3 and 2.6×1019cm−3, 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, ρ, 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 ν0=5.9×1011s–1, NE=2.3×1019cm−3, and δE=2.2kBT. Our simulation results (not shown) suggest that using only exponentially distributed band-tail states requires initial charge concentration n=1018, 6×1017, 1.4×1017, and 0.65×1017cm−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, 1018cm−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 (δEkBT), 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]


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 200×200×200with 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×1018cm−3. Transition to Langevin-type behavior occurs at about 10−6 s at which point the density starts to decay as t1. 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 (≥4×1018 cm−3), 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×1018 cm−3).

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.


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, ρtRt, where ρtand Rtare the carrier density and survival probability at time t, respectively. It turns out that for a reacting CTRW, Rt1/Stfor three dimensions [70], where Stis the number of distinct sites visited by a CTRW. Stcan be visualized as the volume that a carrier sweeps during diffusion. One may also note that in regular lattices, Str2t[67], where r2tis the mean square displacement (MSD) of CTRW at time t. From the above argument, one concludes that


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, followed by a trapping event at the same position x for time interval t drawn from the PDF, wx. In this model, the probability that a carrier is found at position xin space at time tis given by propagator Wxt, whose Fourier-Laplace transform obeys the relation [54]


where k and u represent Fourier and Laplace variables, respectively, W0kindicates the Fourier transform of the initial condition W0x, and wuand λkare the Laplace and Fourier transform of wtand λx, respectively.

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


where δxis 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


where ν0is the intrinsic attempt frequency. Since the jump time, τ, is inversely proportional to the hopping frequency; thus,


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 τ, one arrives at the PDF of the jump time


with a finite mean jump time τ¯equal to


where wt=γexpγ2logν0τ2/2/2πτand γ=kBT/δG. For any DoS with finite mean jump time, using Taylor series for small u, one can write


With initial condition W0x=δx, Eqs. (A2), (A3), (A7), and (A8) imply


where D=a2/τ¯is the diffusion constant. Taking the inverse Fourier-Laplace transform from Eq. (A9), we arrive at the diffusion equation


with initial condition Wx0=W0x. Solution to Eq. (A10) is the well-known Gaussian propagator


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 wt=kBTEUt1. Note that τ¯is a finite value. Thus, following the same procedure performed for Gaussian DoS, we conclude that ρtt1.

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


to yield wt=αν0ν0tα+1, where α=kBT/δEand 0<α<1. Note that τ¯diverges which is indicative of a dispersive process. The Laplace transform of wtfor small u is [71]


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


where Dα=a2/ν0α. From Eq. (A18), MSD is found by inverse Laplace transform of x2u=limk0d2Wku/dk2, which leads to


where the condition of u being small has been applied. Eq. (A19) can be generalized to three dimensions, illustrating that r2ttα. Using this relation and Eq. (A1), we conclude that ρttα.

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


where ρtis the charge carrier density and βtis the time-dependent rate coefficient [72]. Considering boundary condition ρ=ρ0at t=0, one finds that


The choice of βtis semiempirical. Note that the charge density decay follows a power law at long time (ρttα), so it has been suggested [72] that βt=β0tα1,0<α<1, where β0is some constant. For this particular choice of βt, solving BMR Eq. (B1) leads to


where τ0=α/ρ0β01/α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 St[73, 74] via τ=0Stdt, 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]


for which substituting Eq. (B3) into Eq. (B4) leads to τ=ρ0β011+t/τ0αt1α. However, in this equation, as t tends to zero, /dtdiverges; 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, βt, mentioned before.

To resolve this instantaneous recombination at short time, we modify the above time-dependent rate coefficient to βt=β01+t/ηα1where η=α/ρ0β0is a time scaling constant. Experimental data [72] shows that ηis on the order of μsat room temperature for MDMO-PPV and PCBM blend. For this choice of βt, from Eq. B1, one finds


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/αat long time. In contrast to the density decay in Eq. (B3), Eq. (B5) leads to a finite value for /dtin Eq. (B6) as t tends to zero, giving rise to a non-zero initial carrier lifetime τ=η/α. 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 η=106s1. Figure B1A indicates that there is a transition at t=η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, α, on BMR rate with (A) time and (B) charge carrier density utilizing Eq. (B6). In panel A, vertical dotted line corresponds to t=η=10−6s−1 using which these results were obtained.

© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Daniel Christiansen and Shafigh Mehraeen (April 26th 2019). 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, Solar Cells, Majid Nayeripour, Mahdi Mansouri and Eberhard Waffenschmidt, IntechOpen, DOI: 10.5772/intechopen.85074. Available from:

chapter statistics

249total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Nanoplasmonic for Solar Energy Conversion Devices

By Samy K.K. Shaat, Hussam Musleh, Jihad Asad, Nabil Shurrab, Ahmed Issa, Amal AlKahlout and Naji Al Dahoudi

Related Book

First chapter

EU Energy Policies and Sustainable Growth

By Carlo Andrea Bollino and Silvia Micheli

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us