ARTICLE pubs.acs.org/JPCC
Modeling Nongeminate Recombination in P3HT:PCBM Solar Cells Roderick C. I. MacKenzie,* Thomas Kirchartz, George F. A. Dibb, and Jenny Nelson Department of Physics, Imperial College London, South Kensington Campus, London SW7 2AZ ABSTRACT: By introducing tail states into a continuum drift diffusion model, we are able to self-consistently reproduce the experimental voltage-dependent carrier concentration, the currentvoltage curve, and the nongeminate recombination prefactor as measured by transient photovoltage (TPV) measurements of poly-3-hexylthiophene:phenyl-C61-butyric acid methyl ester (P3HT:PCBM). A mobility edge is introduced into the density of states (DoS) for both electrons and holes thus separating each population into mobile and trapped carriers. The effective mobility for each charge type is calculated from the ratio of free to trapped carriers within the DoS and is thus carrier density dependent. By introducing this effective mobility into Langevin’s formula for charge recombination, we are able to reproduce the experimentally observed quasi-third-order recombination rates. Furthermore, we evaluate three possible DoS shapes, a pure Gaussian, a Gaussian with an exponential superimposed, and a pure exponential. It is found that an exponential DoS is essential to reproduce the experimental data.
’ INTRODUCTION In recent years, organic solar cells have shown great promise as a future low-cost source of low carbon electricity.1 The ability to manufacture organic modules using high volume techniques such as gravure printing2 has attracted considerable attention from both academia and industry. Power conversion efficiencies of 78% have now been reported.3 Despite this success, there is still considerable debate surrounding the fundamental physics of device operation, and there is still no reliable framework for simulation of the opto-electrical behavior. Key device models such as nongeminate recombination,47 mobility,8 and the exact nature of the density of states9 (DoS) are still hotly debated.1012 Part of the problem in finding a correct physical model for organic photovoltaics (OPV) operation is that the current density applied voltage curve (JV curve) itself does not contain enough information to give insight into the physical transport or recombination processes within the device. It is only recently that work using small amplitude transient photovoltage (TPV) and charge extraction (CE) has shown that the relationships between (1) the average photogenerated charge carrier density and photovoltage and (2) current density/charge carrier density can be experimentally determined. The standard drift-diffusion models fail to explain the JV characteristics, often overestimating the fill factor in the light and underestimating the diode ideality factor in the dark.13 The low fill factor has been rationalized in terms of various factors, including a light-dependent shunt resistance,13 electric field dependence of the charge separation rate,14 and a higher-order bimolecular recombination rate.15 In the case of P3HT:PCBM solar cells, several groups have observed a charge recombination rate of the form R = knp where p ≈ 3,10,16,17 lending weight to this r 2011 American Chemical Society
as a possible factor. Shuttle and co-workers have recently demonstrated that the magnitude of the fill factor in P3HT: PCBM solar cells and its dependence on light intensity can be explained by a model in which the only charge loss pathway is via the super second-order recombination route. However, there remains no agreed explanation for the quasi-third-order recombination kinetics observed. Previously, the temperature and light intensity dependence of charge recombination kinetics in polymer:fullerene blends was explained in terms of the trapping of hole polarons in an exponential tail of states lying below the mobility edge, such that the thermal activation from these trap states limits the charge recombination.18,19 This model does not explain the nonideality of the solar cell in the dark, as the recombination takes place between free carriers and does not directly involve carriers in tail states. Recently, an alternative recombination mechanism based on ShockleyReadHall statistics has been proposed whereby recombination takes place between free carriers and carriers trapped in an exponential tail of states.12 In this paper, we show that the dependence of current on voltage and charge density on the open circuit voltage (Voc) can both be explained by incorporating different tail state distributions for electrons and holes. Moreover, by fitting the model to solar cell electrical data, we show that some DoS functions can reproduce the experimental results, while others cannot. This allows us to comment on the probable shape of the DoS in organic solar cells. The net result is that bimolecular losses alone Received: January 9, 2011 Revised: April 12, 2011 Published: April 22, 2011 9806
dx.doi.org/10.1021/jp200234m | J. Phys. Chem. C 2011, 115, 9806–9813
The Journal of Physical Chemistry C
ARTICLE
are sufficient to account for the observed fill factor, so long as a suitable DoS is chosen.
’ THE MODEL The model treats the device as a single layer of bulk heterojunction (BHJ) material. The model uses an effective medium to approximate the BHJ. The electronic levels are defined as the HOMO and LUMO energies, and contacts are considered to consist of planes at x = 0 and x = d where the boundary conditions are applied. The model solves Poisson’s equation to obtain the space charge dependent potential profile across the device, where q is the charge on an electron, ε0 the permittivity of free space, εr the relative permittivity of the medium, n the total electron density (n = nfree þ ntrap), and p the total hole density (p = pfree þ ptrap) d dφ ε0 εr ¼ qðn pÞ dx 3 dx
ð1Þ
The bipolar current continuity equations are solved, where Jn,p are the electron and hole current flux densities, and R and G are the net recombination and generation rates per unit volume, respectively.
full FermiDirac statistics is a better approximation 2 !2 3 N E E e , h e , h 5 Fgauss ðE, Ee, h Þ ¼ pffiffiffiffiffiffiffiffiffiffi exp4 2σ2e, h 2πσ 2
ð8Þ
where Ne,h is the total number of states within the Gaussian; σ is the width of the Gaussian; and Ee,h is the lowest unoccupied molecular orbital (LUMO) in the case of electrons or the highest occupied molecular orbital (HOMO) in the case of holes, respectively. In this work, the difference between these energies is referred to as the band gap. The shape of the density of states is particularly important because it determines the relationship between quasi-Fermi level separation (i.e., applied bias), carrier density, recombination, and transport dynamics. Although Gaussian DoS are often used to approximate the DoS in organic semiconductors,27 it has been shown that a Gaussian DoS cannot explain the spectroscopically observed polaron recombination dynamics as well as an exponential DoS.18 Here, to investigate the impact an exponential tail of states would have on the electrical material properties, three trial DoS functions are used: (1) The pure Gaussian as given in eq 8. (2) A single pure exponential of form
∂Jn ¼ qðR GÞ ∂x
ð2Þ
Fexp ðEÞ ¼ N exp expðE=Eu Þ
∂Jp ¼ qðR GÞ ∂x
ð3Þ
which is used to describe both the traps and the free carriers. (3) A Gaussian as in eq 8 with an exponential function (eq 9) superimposed. Thus, the three DoS functions are
The recombination rate R is given as R ¼ kðn, pÞðnp neq peq Þ
ð4Þ
Fg ðE, Ee, h Þ ¼ Fgauss ðE, Ee, h Þ
ð10Þ
Fe ðE, Ee, h Þ ¼ Fexp ðEÞ
ð11Þ
Fge ðE, Ee, h Þ ¼ Fgauss ðE, Ee, h Þ þ Fexp ðEÞ
ð12Þ
where the subscript eq represents the position-dependent carrier concentration in the dark at equilibrium and the prefactor k is defined in analogy with Langevin’s treatment2022 kðn, pÞ ¼
qðRμe ðnÞ þ βμh ðpÞÞ 2ε0 εr
ð5Þ
where μe and μh are the charge density dependent electron and hole mobilities. Although it is not possible to explicitly include all the effects related to the phase segregation of the polymer and fullerene within this model, the Langevin reduction parameters R and β are included to account for the phase segregation effects in the recombination process.23 It is common to introduce one reduction factor when describing experimental data; however, we found it essential to introduce two reduction parameters to selfconsistently fit all the experimental data. The drift-diffusion equations (momentum conservation equations)24 are given as Jn ¼ qμe n
∂Ec ∂n þ qDn ∂x ∂x
ð6Þ
Jp ¼ q μh p
∂Ev ∂p qDp ∂x ∂x
ð7Þ
where Dn,p are the electron and hole diffusion coefficients. No temperature-related driving terms are considered. The density of states (DoS) is often taken as a parabolic band (E0.5) while assuming MaxwellBoltzmann statistics.25 However, in organic semiconductors, a Gaussian function DoS26 with
ð9Þ
The carrier density is then calculated in the usual way Z ¥ nðTe , Ef n Þ ¼ Fn ðEÞf ðTe , E, Ef n ÞdE ð13Þ ¥
Z pðTh , Ef h Þ ¼
¥ ¥
Fh ðEÞð1 f ðTh , E, Efh ÞÞdE
ð14Þ
where f is the FermiDirac function. Although transient studies18 have shown a slow release of trapped carriers to the free states indicating distinct quasi-Fermi levels for free and trapped carrier populations at the steady state,28 we assume that it is possible to describe both the trapped and free carriers of one species with one quasi-Fermi level. Thus, in total the model contains two quasi-Fermi levels, one for the electrons and one for the holes. To include trap states, we define a carrier mobility edge energy (Eedge e,h ) separating the mobile from the trapped carriers within each population. A diagram showing the DoS model along and Eedge can be seen in with the two mobility edges Eedge e h Figure 1. By assigning a constant mobility μ0e,h to the mobile (or “free”) carriers and assigning zero mobility to the trapped carriers, average carrier mobility is found to vary as a function 9807
dx.doi.org/10.1021/jp200234m |J. Phys. Chem. C 2011, 115, 9806–9813
The Journal of Physical Chemistry C
ARTICLE
Figure 1. Pictorial diagram of the LUMO and HOMO DoS used in the model. The mobility edges for the LUMO and HOMO are marked on the diagram as Eedge and Eedge e h , respectively. The exponential LUMO DoS and the HOMO DoS shapes are marked with red and blue lines. The carrier populations as defined by the multiplication of the DoS function with a FermiDirac function are shown as the solid red and blue areas.
of carrier density.
Z
¥ edge
μn ðTe , Ef n Þ ¼ Z Ee
μ0e Fe ðEÞf ðTe , E, Ef n ÞdE
¥ ¥
¼ μ0e Z
nfree ntrap þ nfree
¥ edge Eh
μh ðTh , Ef p Þ ¼ Z
ð15Þ
μ0h Fh ðEÞf ðTh , E, Ef p ÞdE
¥
¥
¼ μ0h
Fe ðEÞf ðTe , E, Ef n ÞdE
Fh ðEÞf ðTh , E, Ef p ÞdE
pfree ptrap þ pfree
ð16Þ
The introduction of the mobility edge and thus carrierdependent mobility also allows the Langevin prefactors (eq 5) to vary as a function of carrier density. The variation of both carrier mobility and the Langeving prefactor as a function of carrier density have previously been observed experimentally.16,29,30 The contact electron and hole carrier densities are used as Dirichlet boundary conditions and as fitting parameters. From the assumed carrier density on one contact, the equilibrium Fermi level can be calculated. By tracing this Fermi level across the device to the other contact and by using the assumed carrier density on that contact, the built-in potential can be calculated. The above set of equations (eq 1eq 16) is a set of highly nonlinear integro-differential equations such that the only practical solution method is numerical. It is standard practice31,32
to solve such sets of device equations using a fully coupled Newton’s method, which has been shown to offer superlinear convergence33 and a performance far superior to that offered by Gummel-type methods. To force efficient global convergence of the scheme, a back tracking algorithm is implemented.34 The current continuity equations are discretized using a Scharfetter Gummel-type scheme.35 In classical inorganic semiconductor theory, FermiDirac integrals with factors of Ex/2 often occur within the integrand. Therefore, there exist standard computationally efficient methods which avoid the need for full numerical evaluation of eq 13. Unfortunately, there exist no such methods for the evaluation of FermiDirac integrals for the above Gaussian DoS or indeed other trial DoS functions of interest in this work. Therefore, before running the simulation, the FermiDirac integrals for both electrons and holes are evaluated numerically as a function of quasi-Fermi level. Searching the results table with a binary search34 minimizes the impact on computational load. Mobility and the Einstein relations are tabulated in a similar manner. Roichman and Tessler27 have pointed out that it is important to use the full degenerate form Einstein relations Dn = (μen)/ (qdn/dEfn) for organic semiconductors due to the deep trap states making the carriers always degenerate. We found this to be especially important when fitting the model to charge extraction data. If n/dn/dEfn is taken as kT then the quasiFermi levels will not equal the equilibrium Fermi level calculated from the boundary conditions in equilibrium for tail slopes with energies greater than 25 meV. This results in an error in the charge density making it impossible to match the charge extraction data. Both external shunt (Rshunt) and series (Rcontact) resistances are taken into account when calculating the voltage and current. This is done by first calculating the current generated by the diode model when external resistances are neglected. Then, Kirchhoff’s circuit laws are used to relate this internal diode voltage and current to the externally measured voltage and current. The optical problem is treated as the propagation of light within a lossy medium with perfect reflection off the back metallic contact—only absorption in the P3HT:PCBM layer is considered. The model discretizes the incident radiation spectrum in terms of wavelength and propagates each discretized wavelength region independently using eq 17. G(x) is the carrier generation rate per unit volume at a depth x within the device due to the incident photons; I0(λ) is the incident photon flux of wavelength λ on the front of the cell; and η is the wavelength-dependent absorption coefficient of P3HT:PCBM.36 It is difficult to know the exact number of free carriers generated in the cell due to the incident light because the exact absorption spectrum of these binary thin films is not known precisely. Also the number of excitons which disassociate is not known. The fitting factor γ is therefore introduced. This factor was typically varied in the range 0.91.1. Z GðxÞ ¼ γ ðηðλÞI0 ðλÞeηðλÞx þ ηðλÞI0 ðλÞe2ηðλÞd þ ηðλÞx Þdλ ð17Þ The downhill simplex algorithm34,37 was used to perform a multidimensional fit of the model to the experimental data. The algorithm was used to minimize the error function Etotal(ζ) = dark onesun (ζ) þ Fonesun Fdark jv (ζ) þ Fcharge(ζ) þ Fjv charge (ζ) þ Fk(ζ), 9808
dx.doi.org/10.1021/jp200234m |J. Phys. Chem. C 2011, 115, 9806–9813
The Journal of Physical Chemistry C
ARTICLE
where F is the normalized error between the simulated and experimental data and ζ is the vector containing the simulations onesun onesun , Fdark parameters. Fdark jv , Fjv charge, Fcharge , and Fk are the error functions for the dark JV, light JV, dark charge extraction, charge extraction at one sun, and the nongeminate recombination function k(n), respectively. The three optimum DoS functions determined in this way are used to generate JV, charge extraction, and k(n) data for comparison with experimental data below.
’ EXPERIMENTAL DATA TO BE MODELED The experimental method and results have been previously reported.5 A simple charge extraction technique is used, whereby the solar cell is held at a point on the JV curve under constant illumination. The illumination source is turned off, and the device is short circuited simultaneously. The current transient is integrated over time, to obtain the excess charge stored within the device. For P3HT:PCBM solar cells, this charge is far higher than one would expect from a parallel plate capacitor and is attributed to bulk charge carriers within the device. The nongeminate recombination rate is measured using the transient photovoltage (TPV) technique,16 where the device is connected to a high impedance oscilloscope and a light bias applied to the device to set the charge density. Then a short optical pulse from a 620 nm nitrogen laser is applied to the device which results in a voltage transient being measured by the oscilloscope. The resulting voltage transient has been shown to be identical in form to the transient absorption obtained using reflection mode transient absorption spectroscopy (TAS) on the same devices. From this signal, the recombination rate is recovered.5 The cells consisted of blends of 1:1 P3HT:PCBM and were produced in accordance to ref 38. The experimental results are plotted in the next section along with the theoretical predictions. ’ COMPARISON OF THEORY AND EXPERIMENT The LUMO and HOMO densities of states obtained from minimizing Etotal when using the DoS models Fe, Fge, and Fg are plotted in Figure 2. The DoS Fge and Fe were both able to reproduce the experimental JV, charge extraction, and k(n) data. To achieve a good fit to the experimental data when using the Fge DoS, it was found essential that the exponential component dominated the Gaussian component. Although the magnitude of the Gaussian component was small (see Figure 2(c)), its presence was found to improve the fit to the k(n) data. It was however not possible to self-consistently fit all the experimental data with the pure Gaussian density of states (Fg) alone. It was found that either a good fit could be obtained to the JV/charge extraction data, or a relatively good fit could be obtained to the k(n) data, but not all. The slope of the JV curve as described by the ideality factor has been linked to the slope of the trap states within the band39 and thus recombination. The inability of the Gaussian DoS to reproduce both the experimental ideality factor and the k(n) data simultaneously suggests that the pure Gaussian is thus an inappropriate model for the P3HT:PCBM DoS. A consequence of choosing an exponential DoS is that at high energies far above the mobility edge the number of states will become very large. These high-energy states will however remain mostly unoccupied as they are far away from the quasi-Fermi levels. The fitting parameters for the best performing (Fge) DoS are given in Table 1, and the corresponding comparisons with
Figure 2. DoS obtained when minimizing Etotal for DoS functions Fg, Fe, and Fge, respectively. (a) The best fit Gaussian DoS function Fg from eq 10. (b) The best fit exponential DoS function Fe from eq 11. (c) The best fit DoS of the form of an exponential DoS with small Gaussians superimposed on Fge from eq 12. The magnitude of the Gaussian was small but however improved the fit to the k(n) data.
experimental data are given below. Parameters held constant such as electron affinities and temperature are given in Table 2. In this work, we fixed the mobility edges at the LUMO/HOMO 9809
dx.doi.org/10.1021/jp200234m |J. Phys. Chem. C 2011, 115, 9806–9813
The Journal of Physical Chemistry C
ARTICLE
Table 1. Parameters Forming the Multidimensional Fitting Space description
parameter
value
left contact electron density
nl
2.108 1024 m3
right contact hole density Langevin reduction parameter
pr R
1.980 1027 m3 1.088 1004
Langevin reduction parameter
β
8.491 1004
effective electron trap density
Nexp e Nexp h Eue
8.150 1024 m3 eV1
characteristic energy for hole
Euh
39 meV
exponential tail free electron mobility
μ0e
3.564 1005 m2 V1 s1
free hole mobility
μ0h
1.572 1006 m2 V1 s1
effective hole trap density characteristic energy for
1.534 1027 m3 eV1 45 meV
electron exponential tail
width of LUMO Gaussian
σe
0.032 eV
width of HOMO Gaussian
σh
0.027 eV
effective DoS in Gaussian LUMO
Ne
3.564 1023 m3
effective DoS in Gaussian HOMO
Nh
4.179 1025 m3
shunt resistance
Rshunt
9.972 1005 Ω
Table 2. Constant Model Parameters description
parameter
value
h Elumo
230 nm 3.8 eV
HOMO electron affinity
Ehomo
4.8 eV
relative permittivity of active layer
εr
3.8
temperature
T
300.0 K
contact resistance
Rcontact
10.0 Ω cm2
electron mobility edge energy
Eedge e
3.800 eV
hole mobility edge energy
Eedge h
4.800 eV
active layer thickness LUMO electron affinity
energies. However, the exact position of the mobility edges are to a large extent arbitrary; if for example an electron mobility edge below the LUMO level is chosen, the free electron mobility would have to be decreased to compensate for the higher number of free electrons. Thus, the mobility averaged over all carrier energies would remain the same. It should be noted that the charge extraction technique can only probe the DoS around the quasi-Fermi levels, and because of the low charge carrier density in operating solar cells, this means that only the tail of the density of states is probed. However, the low-energy part of the external quantum efficiency EQE data recently published by Presselt and coworkers40 can be approximately reproduced with the DoS presented in Figure 2c. Before going further, it is worth discussing the choice of the fitted parameters in Table 1. Although there are 15 free parameters in this table, our choice of parameters is far more constrained than it would first appear. First, as can be seen from Figure 2, the impact of the Gaussian DoS superimposed on the exponential is minimal, and thus we can effectively discard the four terms associated with the Gaussian. Although the characteristic tail energies are a fit to experiment, they are closely related to the ideality factor as was proved theoretically by Berkel et al.39 If these tail slopes are moved by more than 5 meV, the dark ideality cannot be reproduced. The Langevin reduction parameters are a
Figure 3. (a) Comparison of numerical and experimental5 light JV curves at 0.5, 1.0, 2.4, and 3.5 suns. The crosses represent the experimental data, and the solid lines represent the simulation using a DoS formed from the Fge DoS eq 12, Figure 2(c). (b) A comparison of numerical and experimental5 dark JV curves. The crosses represent the experimental data, and the solid lines represent the simulation using a DoS formed from the Fge DoS eq 12, Figure 2(c). Also plotted on this graph is the dark curve resimulated with no trap states included in the model (dashed line).
pure fit to the experimental data; however, they are very close to the values reported from experimental measurements10 for P3HT:PCBM. The shunt resistance is associated with imperfections in the device, and it is relatively independent of the other device parameters and can be accurately estimated by examining the low voltage (
