Analytical expression for the current-voltage characteristics of organic bulk heterojunction solar cells

Organic bulk heterojunction (BHJ) solar cells have a huge potential as a source of clean energy in the future.^1,2 In order to improve the performance, understanding of the current–voltage (J-V) characteristics is important. Organic solar cells differ from conventional inorganic solar cells. In organic solar cells, the absorbed light generates excitons first instead of free electrons and holes as in conventional inorganic solar cells. The excitons must be dissociated into free electrons and holes in order to produce useful electric current. Two materials (donor and acceptor) with different lowest unoccupied molecular orbitals (LUMOs) and different highest occupied molecular orbitals (HOMOs) are needed to assist the dissociation of excitons into free charge carriers. In order to generate and extract free charge carriers efficiently, the donor and the acceptor materials are blended together to create donor-acceptor interfaces throughout the device’s active layer (called bulk heterojunction).

Drift-diffusion model with metal–insulator–metal picture³ is a widely used model for organic BHJ solar cells.⁴ This model solves the continuity and the Poisson’s equations, accounting for the drift and diffusion of the charge carriers, the effect of space charge, and the generation and loss of the charge carriers.^4,5 A numerical procedure has been widely used in order to obtain the J-V characteristics from this model.^4,6 There are a number of analytical studies attempting to describe the J-V characteristics of organic BHJ solar cells that can be found in the literature.^7–10 In previous analytical studies,^7–9 the electric field is assumed constant with negligible recombination. Hence, the J-V characteristics are obtained by adding the photogenerated and the dark current densities instead of solving the above mentioned drift-diffusion model. Recently, Chowdhury and Alam¹⁰ proposed an analytical model based on the drift-diffusion model but did not incorporate the intrinsic bimolecular recombination.

The aim of this work is to derive an expression to describe the J-V characteristics of organic BHJ solar cells by solving the drift-diffusion model analytically, where the effect of bimolecular recombination is considered. Compared to previous analytical works, this expression should be more appropriate for analyzing the J-V characteristics of organic BHJ solar cells when the bimolecular recombination is not negligible.

First, let us present the operational mechanism for organic BHJ solar cells that is used in this study. After the light is absorbed, the generated excitons must reach the donor-acceptor interface before they can dissociate into free electrons and holes. At the interface, the exciton splits into an electron in the acceptor material and a hole in the donor material, where the electron and hole are presumed to form a bound geminate electron-hole pair (called polaron pair). The polaron pair needs to further dissociate into a free electron and a free hole in order to produce photocurrent. At the interface, the polaron pair can also decay, which is a geminate monomolecular (first order) recombination.¹¹ Hence, the probability for a polaron pair to dissociate into free charge carriers can be defined as the ratio of the dissociation rate of a polaron pair to the sum of the dissociation rate and the decay rate of a polaron pair. The free electron in the acceptor and the free hole in the donor can also meet at the interface and recombine into a polaron pair again. This is a non-geminate bimolecular (second order) recombination where the recombination rate is the rate for free electrons and free holes finding each other, which can be predicted by the Langevin theory.¹² 

There are also other possible recombination mechanisms (see. Ref. 11 for example) and their importance is still being studied and debated. It is generally believed that trap-assisted recombination can be significant if trap densities are high. For example, Cowan et al.¹³ concluded that the recombination in organic BHJ solar cells evolves from being monomolecular at short circuit to being bimolecular at open circuit. This finding suggests that non-geminate monomolecular recombination (i.e. trap-assisted recombination) should be considered. However, Kniepert et al.¹⁴ recently demonstrated that the recombination in their as-prepared and annealed poly(3-hexylthiophene) and [6,6]-phenyl-C₇₁-butyric acid methyl ester (P3HT:PCBM) solar cells is strictly bimolecular and trap-assisted recombination cannot account for their results. In our study here, we assume the contribution of trap-assisted recombination in the studied organic BHJ solar cells is negligible.

Next, we present the basics of the drift-diffusion model for organic BHJ solar cells based on the operational mechanism described above. Details of the model can be referred to Refs. 4 and 5 for example. As usual, the transport, generation, and loss of the charge carriers are assumed to be 1-dimensional. The Poisson’s equation is given by

where F is the electric field, q is the elementary charge, ε is the effective permittivity of the organic active layer, p is the hole density, and n is the electron density.

Based on the above operational mechanism, the continuity equations for electrons and holes in steady-state condition can be written as^4,15

where J_n is the electron current density, J_p is the hole current density, P is the dissociation probability of polaron pairs into free charge carriers, G is the generation rate of polaron pairs per unit volume, and R_L is the non-geminate bimolecular recombination rate.

The electron and the hole current densities consist of the drift and the diffusion components, and are given by¹⁶ 

where μ_n is the electron mobility, μ_p is the hole mobility, D_n is the electron diffusion coefficient, and D_p is the hole diffusion coefficient. For simplicity, μ_n and μ_p are assumed to be constant, which is needed in order to solve the drift-diffusion model analytically. In a more realistic model, the carrier mobility is electric field, temperature, and carrier density dependent.¹⁷ The diffusion coefficients are given by

where k is the Boltzmann constant, and T is the absolute temperature. Regarding the non-geminate bimolecular recombination in organic BHJ solar cells, it was found that the recombination rate determined experimentally is reduced compared to the Langevin’s theory.¹⁸ Taking this into account, a reduction factor ς is often included, and the bimolecular recombination rate reads

where $γ=q μ n + μ p /ε$ is the Langevin recombination constant.

Regarding the polaron pair dissociation in organic semiconductor blends, the exact mechanism is still uncertain. The Onsager-Braun theory^19,20 is a commonly used quantitative model for considering the polaron pair dissociation at the interface (see Refs. 4, 5, and 21 for instance). In this study, we use the Onsager-Braun model where the dissociation probability of polaron pairs P is defined as

where k_d is the polaron pair dissociation rate, a is the polaron pair separation, and k_f is the polaron pair decay rate. Besides the Onsager-Braun model, the role of charge delocalization in the dissociation of polaron pairs into free charge carriers has also been studied and proposed.²² 

In order to numerically solve the continuity and the Poisson’s equations, boundary conditions at the contacts between the electrodes and the active layer are required. For the electrical potential ψ (where ψ/q = − ∫Fdx), the boundary condition is

where V_bi is the built-in voltage and V_a is the applied voltage. Here, x = 0 and x = L are the contacts at the anode and the cathode, respectively. The built-in voltage V_bi is defined as

where χ_a and χ_c are the work functions of the electrodes at anode and at cathode, respectively. The carrier densities at the contacts (x = 0, L) can be simply given by the Boltzmann statistics with

where N_c and N_v are the effective density of states for the LUMO of the acceptor and the HOMO of the donor, respectively, E_g is the effective band gap (i.e. the difference between the LUMO of the acceptor and the HOMO of the donor), ϕ_a is the hole injection barrier at the anode (i.e. the difference between the HOMO of the donor and the work function of the electrode at anode), and ϕ_c is the electron injection barrier at the cathode (i.e. the difference between the LUMO of the acceptor and the work function of the electrode at cathode). For simplicity, the energetic disorder in the organic active layer is not considered in this study.

In order to obtain the J-V characteristic, the continuity and the Poison’s equations are numerically solved using the Scharfetter–Gummel discretization²³ and an iteration approach based on the work of Gummel.²⁴ 

In order to solve the drift-diffusion model analytically, several approximations and assumptions are required. The electric field F inside the active layer is assumed to be uniform. This should be valid since drift-diffusion simulations show that F is almost constant inside the active layer.⁴ The uniform electric field F is given by (see Ref. 7 for example)

where L is the thickness of the active layer. The polaron pair dissociation probability P is also assumed to be uniform inside the device. For the Onsager-Braun model, a uniform F together with constant charge mobilities will automatically lead to a uniform P. Furthermore, we assume that the generation rate of polaron pairs per unit volume G is constant and uniform inside the device. For organic BHJ solar cells, it has been shown that a constant generation profile gives identical electrical characteristics compared to the generation profile obtained from optical model as long as the thickness of the active layer is less than 250 nm.²⁵ Therefore, this assumption should be valid when the thickness of the active layer is less than 250 nm.

The bimolecular recombination rate R_L [Eq. (5)] is dependent on the position inside the active layer x because the product of n and p is different at different x. In order to obtain an analytical solution, let us also assume that the bimolecular recombination rate inside the device is uniform (denoted by R_C) and equal to the average R_L. Therefore, R_C is independent of x and its value is given by

The recombination current density J_R is given by

where R is any recombination rate in general. From Eq. (12), for R = R_C (where R_C is independent of x) and integrating over x from 0 to L, we get

Equation (13) gives the total recombination current density in the device for R = R_C. Applying Eq. (11) to Eq. (13), we get

Therefore, the total recombination current density using R = R_C is the same as the one using R = R_L. We can argue that for the current-voltage (J-V) characteristic of a solar cell, it is not important where inside the active layer the charge carriers are generated and lost, as long as the net charge carriers that are collected by the electrodes are the same. For example, by integrating Eq. (2a) over x from 0 to L (with P and G are independent of x), it can be seen that the net flow of the electron current density inside the device [i.e. J_n(0) − J_n(L)] depends on the total recombination current density. This means that we can replace R_L with R_C in Eq. (2a) and still be able to get the same net electron current density. Therefore, the basic idea of assuming uniform bimolecular recombination rate R_C is to obtain an analytical solution that should be able to describe the J-V characteristics of organic BHJ solar cells. However, this assumption may not lead to very accurate descriptions of carrier densities, electron current density, and hole current density as functions the position inside the active layer x.

Figure 1 shows the ratio of the recombination current density for R = R_L (denoted as J_RL) to the polaron pair dissociation probability P as a function of applied voltage V_a obtained using drift-diffusion numerical simulation. In this study, all numerical simulations are performed using a commercial program called SETFOS,²⁶ which is specifically made for simulating organic BHJ solar cells based on the drift-diffusion model as described in Sec. II. The values of the parameters used in the simulation are shown in Table I. The simulation produces a short circuit current density J_sc of -102.6 A/m², an open circuit voltage V_oc of 0.59 V, and a fill factor FF of 0.578, which are typically measured values for a P3HT:PCBM solar cell.

The behavior seen in Fig. 1 can be easily understood. As the applied voltage V_a is increased, the electric field decreases. This makes it harder for the electrodes to extract the generated free charge carriers (given by the product of P, G, and L). As a result, this increases the carrier densities per generated free charge carriers, thus increases the recombination per generated charge carriers. Therefore, it is expected that the ratio of the bimolecular recombination to the generated charge carriers (or to the dissociation probability P, since G and L are constants) increases with V_a. As seen in Fig. 1, we find that the natural logarithmic of J_RL/P is somewhat nonlinear as a function of applied voltage V_a. Varying the values of various parameters, we find that the general behavior seen in Fig. 1 is the same at least between zero applied voltage and the open circuit voltage. From this, let us approximate and write

where J₀, g, and s are (fitting) constants. Equation (15) can be rewritten as

If we assume that the bimolecular recombination rate is uniform inside the device as given by Eq. (11), the total recombination current densities using R_C and R_L must be equal (i.e. J_RC = J_RL), as shown before. Since Eq. (13) equal to Eq. (16), we get

where R₀ = J₀/qL.

Inserting Eqs. (3a) and (3b) into Eqs. (2a) and (2b), and replacing R_L [Eq. (5)] with R_C [Eq. (17)], the continuity equations become

Equations (18a) and (18b) can be solved using

where A_n, B_n, A_p, and B_p are independent of x. Applying the boundary conditions at the contacts [i.e. Eqs. (9a) to (9d)] to Eqs. (19a) and (19b), we can solve for A_n, B_n, A_p, and B_p, which are given by

It should be noticed that the values of A_n, B_n, A_p, and B_p are different at different applied voltages V_a. Inserting Eqs. (19a) and (19b), and their differentiations with respect to x into Eqs. (3a) and (3b), we can get the expressions for J_n and J_p, where

The total current density J is the sum of J_n and J_p, which is given by

It can be seen that the total current density J is independent of x, which must be the case for a steady current density. Using Eq. (10), we can rewrite the term R_C/F in Eq. (22) as $R C L/ V a − V bi$⁠. Applying Eq. (11), this equal to $∫ 0 L R L dx/ V a − V bi$⁠. Hence, it can be seen that although we assume uniform bimolecular recombination rate R_C to obtain Eq. (22), it should still give a J-V characteristic as if we consider the actual bimolecular recombination rate R_L. It is worth noting that Eq. (22) can be derived using any other form of bimolecular recombination rate [i.e. not limited to the form of Eq. (17)] as long as the recombination rate used is independent of x and equal to the average of the actual bimolecular recombination rates R_L. Similarly, any other form of polaron pair dissociation (i.e. not limited to the Onsager-Braun model) can also be used as long as it is independent of x.

In order to test this analytical expression, let us compare Eq. (22) with the J-V curve obtained using drift-diffusion simulation, which is widely used to analyze the experimental J-V characteristics of organic BHJ solar cells. Figure 2 compares the J-V characteristic obtained using drift-diffusion simulation with the one obtained using the analytical expression. The parameters used for the simulation and the analytical expression are given in Table I. It can be seen that the J-V characteristic obtained using the analytical expression is in good agreement with the J-V characteristic obtained using drift-diffusion simulation. In order to obtain the constants R₀, g, and s [see Eq. (17)], let us rewrite Eq. (22) as

where the values of m and c at any given V_a can be determined from the values of the parameters in Table I. By inserting the simulated current density at short circuit into Eq. (23), we can easily solve for R₀. To obtain g and s, we can insert the simulated current densities at two other applied voltages to obtain two equations that enable us to solve for g and s. However, this requires numerical techniques where numerical solvers such as MATLAB can be helpful. Alternatively, we can insert the open circuit voltage and current density (i.e. J = 0 A/m²) into Eq. (23), and guess the value of s in order to get the value of g. Then, the obtained g and s, together with R₀ determined earlier are used to obtain the J-V curve. This process of guessing the values of s and g is repeated until a good fit to the simulated J-V curve is obtained. The dots in Fig. 2 are obtained with R₀ = 3.69 × 10²⁹m⁻³s⁻¹, g = 17.65V⁻¹, and s = 5.

Figure 3 compares the carrier densities as functions of x obtained using drift-diffusion simulation with the ones obtained using Eqs. (19a) and (19b) with the values of R₀, g, and s as mentioned above. As expected, the carrier densities obtained using the analytical approach are only in moderate agreement with the ones obtained using drift-diffusion simulation. As mentioned in Sec. III, the assumption of uniform bimolecular recombination rate may lead to somewhat inaccurate descriptions of carrier densities, electron current density, hole current density, and others as functions of x.

There are a few advantages of the analytical expression over numerical simulation of the drift-diffusion model. For example, it is practically impossible for the simulation to exactly reproduce the experimental J-V curve of an organic BHJ solar cell. Moreover, the quality of fittings to experimental J-V curves obtained using simulations may not be very good for some solar cells. This means that the recombination current densities obtained from simulations may differ quite significantly from the actual recombination current densities of the studied solar cells. Using the analytical expression, it is possible to directly extract the recombination current densities from experimental J-V characteristics. This can be done by inserting the experimental current densities and the associated applied voltages into Eq. (22) [or Eq. (23)], and then directly determine the value of R_C at each applied voltage. Then, the recombination current density can be obtained by multiplying the extracted R_C with the thickness of the solar cell and the elementary charge. Therefore, this analytical expression enables us to directly extract the overall recombination loss of an organic BHJ solar cell from the experimental J-V curve. It is worth noting that series and shunt resistances in solar cells should be minimized and maximized, respectively. If the effects of these resistances are significant, then the experimental J-V curve should be adjusted accordingly before the analytical expression can be used to analyze it.

Furthermore, we also find that our numerical program fails to produce the desired output when simulating an organic BHJ solar cell with a relatively large band gap and a high injection barrier. A large band gap organic BHJ solar cell has a useful application as a sub-cell of a tandem solar cell. Hadipour et al.²⁷ showed that the J-V characteristic of a tandem solar cell can be constructed from the J-V characteristics of individual sub-cells. Hence, it is important to analyze the J-V characteristic of each sub-cell. The values of the parameters used in the simulations are shown in Table II. The simulation with ϕ_a = 0 eV produces J_sc = − 28.92 A/m², V_oc = 1.06 V, and FF = 0.72. However, for a tandem solar cell with series configuration, it is unlikely that both ϕ_a and ϕ_c are zero since the anode (cathode) of a large band gap sub-cell is connected to the cathode (anode) of a smaller band gap sub-cell, and the effective band gaps of the sub-cells almost certainly overlap with each other. When ϕ_a = 0.9 eV is used in the drift-diffusion simulation, we get J_sc = − 26.08 A/m², V_oc = 0.367 V, and FF = 0.513. However, when ϕ_a = 1.0 eV is used, the drift-diffusion simulation starts to become unstable, producing V_oc = 0.267V but J_sc = 0 A/m². When ϕ_a = 1.1 eV, we get J_sc = 0 A/m² and V_oc = 0 V. When ϕ_a = 1.2 eV is used, the numerical simulation simply fails to produce any result. Since a commercial program is used for the simulations, we are only able to vary the numerical mesh size and the residual error. Increasing those iteration parameters decrease the accuracy, but increase the possibility of getting the desired output in general. The exact reason for the failure is uncertain since we are unable to get the desired output even if those two iteration properties are increased and decreased. It is possible that when certain values of certain parameters are used, the drift-diffusion simulation would fail to fully converge (which is unsurprising for a numerical approach). However, by using the analytical expression, we should be able to, for example, extract and analyze the recombination current density as a function of applied voltage if the experimental J-V curve of the large band gap solar cell with properties as shown Table II is available.

In order to clearly demonstrate the advantage of this analytical expression, let us try to estimate J-V curves of the solar cell with properties as shown in Table II when ϕ_a = 1.0eV and ϕ_a = 1.1eV. Since we do not have the experimental J-V data, let us use the simulated J-V curve of the solar cell when ϕ_a = 0.9 eV in order to estimate the R_C [refer Eq. (17)] when ϕ_a = 1.0eV and ϕ_a = 1.1 eV. In order to satisfactorily reproduce the simulated J-V curve when ϕ_a = 0.9 eV by using the analytical expression, we set V_oc = 0.35V instead of V_oc = 0.367V. Therefore, it is worth noting that in order to satisfactorily reproduce a given J-V curve by using the analytical expression with R_C as given by Eq. (17), we may need to set the V_oc and/or J_sc to be slightly different compared to the V_oc and J_sc of the given J-V curve. We find that R₀ = 2 × 10²⁸m⁻³s⁻¹, g = 3.423V⁻¹, and s = 1 can satisfactorily reproduce the simulated J-V curve for the solar cell when ϕ_a = 0.9eV. From this, let us write R_C [refer Eq. (17)] for the solar cell when ϕ_a = 1.0eV and ϕ_a = 1.1eV as

Since R₀, g, and s are obtained from the simulated J-V curve when ϕ_a = 0.9eV, V_diff is introduced to account for the difference between the built-in voltage when ϕ_a = 0.9eV and the built-in voltages when ϕ_a = 1.0eV and ϕ_a = 1.1eV. Here, we simply take V_diff = 0.1V when ϕ_a = 1.0eV and V_diff = 0.2V when ϕ_a = 1.1eV. Figure 4 shows the J-V curves estimated using the analytical expression with R_C given by Eq. (24) when ϕ_a = 1.0eV and ϕ_a = 1.1eV. When ϕ_a = 1.0eV, we get J_sc = − 24.42 A/m², V_oc = 0.245V, and FF = 0.53. When ϕ_a = 1.1eV, we get J_sc = − 22.12 A/m², V_oc = 0.145 V, and FF = 0.482. Therefore, we demonstrate that when the numerical simulation fails to produce the J-V curve of a specific organic BHJ solar cell, the analytical expression is useful in order to estimate the J-V curve and analyze the organic BHJ solar cell.

In conclusion, we have proposed an analytical approach to describe the current-voltage (J-V) characteristics of organic BHJ solar cells. The derivation is obtained by analytically solving the drift-diffusion model for organic BHJ solar cells with the assumption of uniform bimolecular recombination rate that equal to the average of the actual bimolecular recombination rates. We show that when the uniform bimolecular recombination rate is assumed, the derivation should still produce an expression for the current as a function of applied voltage as if the actual bimolecular recombination rate is used. Compared to previous analytical expressions, this analytical expression should be more suitable for analyzing the J-V characteristics of organic BHJ solar cells when the bimolecular recombination is not negligible. There are a few advantages of this analytical expression compared to numerical simulation of the drift-diffusion model. As discussed in Sec. IV, applying the analytical expression to experimental J-V data enable us to directly extract and analyze the overall recombination loss of an organic BHJ solar cell. Furthermore, as shown in Sec. IV, we find that when the injection barrier of a large band gap organic BHJ solar cell is high, our numerical program fails to give the desired output. We then demonstrated that this analytical approach can be used to analyze the large band gap organic BHJ solar cell. Therefore, when the numerical simulation fails to produce the desired results for a specific organic BHJ solar cell, the analytical expression is useful in order to analyze the organic BHJ solar cell.

The authors thank University of Malaya for financial support through BKP grant (BK025-2013) and the Malaysian Ministry of Higher Education for funding under the Long Term Research Grant Scheme (LR003-2011A).