Review of diffusion models for charge-carrier densities in dye-sensitized solar cells

In essence, there are seven main chemical reactions used to describe DSSCs in operation [1, 6]. Four of these chemical reactions enhance the performance of DSSCs, whereas the other three mechanisms are detrimental to the generation of electricity.

In DSSCs, photosensitive dye molecules absorb photons and enter an excited energy state:

$$\begin{eqnarray}&&{\rm{S}}+h\nu \to {{\rm{S}}}^{* },\end{eqnarray} \tag{ 1 }$$

where S is the dye molecule, h is Planck's constant and ν is the frequency of the photon. The photosensitive dye is designed to use this excited energy state to donate (or inject) electrons into the semiconductor's conduction band:

$$\begin{eqnarray}&&{{\rm{S}}}^{* }\to {{\rm{S}}}^{+}+{e}_{{CB}}^{-}.\end{eqnarray} \tag{ 2 }$$

The photosensitive dye is regenerated by the electrolyte couple:

$$\begin{eqnarray}&&{{\rm{S}}}^{+}+{\rm{M}}\to {{\rm{M}}}_{{ox}}+{\rm{S}},\end{eqnarray} \tag{ 3 }$$

where M is the electrolyte and M_ox is the oxidized electrolyte. This process enables the chemical reaction (1) to repeat. Finally, regeneration of the electrolyte is carried out by the counter electrode:

$$\begin{eqnarray}&&{{\rm{M}}}_{{ox}}+{e}_{{CE}}^{-}\to {\rm{M}}.\end{eqnarray} \tag{ 4 }$$

As a result of chemical reactions (1) to (4), a DSSC is able to generate electricity provided that there is sufficient sunlight as an energy input. However, DSSCs are also subject to loss mechanisms, which generally involve recombination of electron at interfaces. Excited dye molecules do not necessarily inject electrons into the conduction band as described by chemical reaction (1), instead they may emit the photon as radiation:

$$\begin{eqnarray}&&{{\rm{S}}}^{* }\to {\rm{S}}+h\nu .\end{eqnarray} \tag{ 5 }$$

Alternatively, injected electrons may recombine with the DSSC before leaving the device to power a load. Active electrons in the semiconductor's conduction band may be used for regeneration either in the dye or in the electrolyte, which are respectively described by

$$\begin{eqnarray}&&{{\rm{S}}}^{+}+{e}_{{CB}}^{-}\to {\rm{S}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{\rm{M}}}_{{ox}}+{e}_{{CB}}^{-}\to {\rm{M}}.\end{eqnarray} \tag{ 7 }$$

Construction of DSSCs must therefore enhance reactions (1) to (4) and reduce the potential loss of electrons through the processes (5) to (7).

1.2.1. Photosensitive dyes

The unique feature of DSSCs is its photosensitive dyes which are used to absorb sunlight [6]. Even though several types of dye molecules have already been proposed, the development of new dyes continues an active area of research.

One of the dyes that saw continual use since its introduction is the Ruthenium(II) N3 dye [7, 8]. The relationship between the molar absorption coefficient of N3 dye and the wavelength is given in figure 2 [8]. In their paper, Wu et al [8] show that the molar absorption coefficients of CYC-B1 and CYC-B5 dyes are generally higher than that of N3 dye, despite their inferior performances at lower wavelengths. However, the Ruthenium based N3 dye has continued its popularity, owing to its broad light absorption spectrum, increased duration in the excited state and high thermal stability [9, 10].

Zoom In Zoom Out Reset image size

Another well-known Ruthenium based dye is N719 [3]. The absorption coefficient of N719 relative to the wavelength is shown in figure 3 [11]. From the figure, we see that the absorbance of the N719 dye is greatly improved when anchored to the ZnO semiconductor as opposed to ${\mathrm{TiO}}_{2}$.

Zoom In Zoom Out Reset image size

1.2.2. Nanoporous semiconductor

Given that the nanoporous semiconductor is responsible for the bulk of the electron transport in DSSCs, the choice of semiconductor is extremely important. A nanoporous semiconductor requires a wide band gap, which allows electrons to easily receive energy from excited dye molecules. In table 1, we summarise commonly used semiconductors together with their band gaps and efficiencies [12].

Table 1. Nanoporous semiconductors used in DSSCs [12].

Semiconductor	Band Gap	Efficiency
${\mathrm{TiO}}_{2}$	3.23	12.30%
ZnO	3.30	5.60%
Nb2O5	3.49	5.00%
Zn2SnO4	3.70	3.70%
SrTiO3	3.20	1.80%
WO3	2.60–3.00	0.75%
SnO2	3.60	0.30%

From table 1 we see that ${\mathrm{TiO}}_{2}$ vastly surpasses other semiconductors by way of efficiency. In their paper, Lau et al [12] show that ${\mathrm{TiO}}_{2}$ is generally regarded as superior for its wide band gap, low cost and resistance to photo-corrosion. Figure 4 shows the porous structure of ${\mathrm{TiO}}_{2}$ under electron microscope [13].

Zoom In Zoom Out Reset image size

1.2.3. Electrolyte couple

1.2.4. Counter electrode

One approach to modelling DSSCs is based on the study of traditional junction solar cells. In particular, the drift-diffusion equations are used to describe the electron density in the semiconductor's conduction band while the hole density in the semiconductor's valence band and the strength of the electric field are modelled by Poisson's equation [2, 15, 30]. These equations are inherited from the models for the silicon solar cells pioneered by Chapin et al [31] and are based on Shockley's mathematical models for semiconductors [32]. Although this approach may benefit from the studies in the past decades, Gregg [33] mentions that it is not necessarily the best modelling approach for DSSCs. This is because the electric field is considered to have a negligible effect on the overall operation of DSSCs and is therefore ignored in subsequent models of DSSCs [22]. In this paper, we review the development of the diffusion models for DSSCs.

Due to the structure of DSSCs, only one spatial dimension is considered. Given a DSSC of thickness d, we orient the device so that the counter electrode is located at x = d. From this, the conduction band electron density n(x, t) at position x and time t in the DSSC obeys the following partial differential equation (PDE) [34]:

$$\begin{eqnarray}&&\displaystyle \frac{\partial n}{\partial t}=\displaystyle \frac{\partial }{\partial x}\left[D(n)\displaystyle \frac{\partial n}{\partial x}\right]+G(x)+R(n),\end{eqnarray} \tag{ 8 }$$

where D(n) is the density-dependent diffusion coefficient, G is the spatially dependent electron generation term and R(n) is the density dependent recombination term. For DSSCs, we assume no electron flux at the counter electrode, giving rise to the boundary condition [35]

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial n}{\partial x}\right|}_{x=d}=0.\end{eqnarray} \tag{ 9 }$$

The electron density at the TCO interface (x = 0) is dependent on the bias voltage V on the DSSC [36], and is given by

$$\begin{eqnarray}&&n(0,t)={n}_{{eq}}{e}^{\tfrac{{qV}}{{{mk}}_{B}T}},\end{eqnarray} \tag{ 10 }$$

where n_eq is the dark equilibrium electron density, m is the ideality factor, q is the electron charge, k_B is Boltzmann's constant and T is the temperature of the DSSC. For short-circuit conditions, we set V = 0 and use the resulting boundary condition $n(0,t)={n}_{{eq}}$. Under open-circuit conditions (J = 0), an alternative boundary condition for x = 0 is given by [37]

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial n}{\partial x}\right|}_{x=0}=0,\end{eqnarray} \tag{ 11 }$$

which allows us to compute the electron density if V_oc is unknown. The DSSC is assumed to be initially at equilibrium, giving the simple initial condition

$$\begin{eqnarray}&&n(x,0)={n}_{{eq}}{e}^{\tfrac{{qV}}{{{mk}}_{B}T}}.\end{eqnarray} \tag{ 12 }$$

This model was developed in 1994 by Södergren et al [36], in which the steady-state special case was considered under the assumptions of linear diffusion and recombination with the boundary conditions based on short-circuit conditions (V = 0). This diffusion model is a popular choice for modelling DSSCs (as seen in [20, 36–39]). In particular, the linear diffusion equation is coupled with electrolyte equations in some models (which are discussed in more detail in section 2.5)

In 1996, Cao et al [22] extended this model to include time dependence, resulting in the first PDE for electron diffusion. This paper also suggested the first nonlinear version of the equation, with the diffusion coefficient D(n) = n.

The diffusion coefficient D(n) in equation (8) is an important component, as it reconciles the boundary condition $n(0,t)={n}_{{eq}}$ with the source term G(x) that is typically at its maximum at x = 0. In the diffusion model of Södergren et al [36], the diffusion coefficient is assumed to be a constant,

$$\begin{eqnarray*}&&D(n)={D}_{0}.\end{eqnarray*}$$

Though many papers adopt a constant value for D₀ [35, 36, 40], a study in 2006 introduced a model incorporating the effect of the porosity of the nanoporous semiconductor into the diffusion coefficient [38].

In particular, Ni et al [38] gave the following expression for the diffusion coefficient that is dependent on the porosity p of the nanoporous semiconductor (${\mathrm{TiO}}_{2}$ in this case):

$$\begin{eqnarray}&&{D}_{0}(p)=1.69\times {10}^{4}\left(-17.48{p}^{3}+7.39{p}^{2}-2.89p+2.15\right).\end{eqnarray} \tag{ 13 }$$

The electron generation term G(x) in (8) is a source term designed to incorporate the increase of conduction band electrons induced by sunlight exposure. The Beer–Lambert model is used exclusively, and is given by

$$\begin{eqnarray}&&G(x)={\int }_{{\lambda }_{\min }}^{{\lambda }_{\max }}{\eta }_{{inj}}\phi (\lambda )\alpha (\lambda ){e}^{-\alpha (\lambda )x}d\lambda ,\end{eqnarray} \tag{ 14 }$$

where photons absorbed are in the wavelength range $[{\lambda }_{\min },{\lambda }_{\max }]$, η_inj is the electron injection efficiency (the proportion of electrons successfully injected into the semiconductor conduction band, also referred to as Φ_inj and η_inj ), ϕ is the wavelength dependent incident photon flux and α is the wavelength dependent absorption coefficient of the photosensitive dye.

For ease of computation, several papers omit wavelength dependence in (14) (such as [22, 37, 41, 42]), leading to a simplified Beer–Lambert model of the form

$$\begin{eqnarray}&&G(x)={\eta }_{{inj}}\varphi \alpha {e}^{-\alpha x},\end{eqnarray} \tag{ 15 }$$

where φ is an integrated version of ϕ and wavelength dependence of α is omitted. This generalized form of the Beer–Lambert law is also used in [43–45], for example.

We comment that equation (15) was first introduced by Södergren et al [36], who considered two simplified forms of (14). The first form, which is 15, assumes that substrate-electrode (SE) illumination is located at x = 0. The second form assumes that electrolyte-electrode (EE) illumination is located at x = d. In this case, the Beer–Lambert law is of the form [20, 46]

$$\begin{eqnarray}&&{G}_{{EE}}(x)=\varphi \alpha {e}^{\alpha (x-d)}.\end{eqnarray} \tag{ 16 }$$

Finally, an analytical approach to model the electron generation term was developed by der Maur et al [16]. Their analytical model is motivated by relaxing the assumption of a constant absorption coefficient in DSSC modelling. By applying a rate equation to the density of photoactive dye molecules and a continuity equation on the electron generation term, their model for electron generation term, referred to here as g(x), is given implicitly by

$$\begin{eqnarray}&&g(x){e}^{{\tau }_{r}\sigma \left(g(x)-g(0)\right)}=g(0){e}^{-{D}_{0}\sigma x},\end{eqnarray} \tag{ 17 }$$

where τ_r is the dye regeneration time, σ is the attenuation cross section of the absorbing dye molecule and D₀ is the total density of photoactive dye molecules. If τ_r  = 0, we recover the classical Beer–Lambert model $g(x)=g(0){e}^{-{D}_{0}\sigma x}$. This extended model allows the DSSC modelling to incorporate the effect of light intensity on electron generation within an operating DSSC.

2.3.1. Incident photon flux

Central to the Beer–Lambert equation (14) is the incident photon flux ϕ, which is the input from sunlight. Figure 5 is a plot of ϕ, the measured AM1.5 solar spectrum against wavelength λ [47].

Zoom In Zoom Out Reset image size

The literature models the input of sunlight in solar cells by two expressions, namely the incident photon flux and the sunlight intensity. Though equivalent, the incident photon flux considers the photons per unit area provided by sunlight, which has unit ${{\rm{m}}}^{-2}{{\rm{s}}}^{-1}$. On the other hand, sunlight intensity is a measurement of the AM1.5 solar spectrum, which considers sunlight as energy per unit area, so that the unit is given by Wm⁻². Dividing the sunlight intensity (Wm⁻²) by energy (Joules), we can convert from the sunlight intensity to photon flux (m⁻²s⁻¹). Here, the energy we use is defined as the energy of a photon,

$$\begin{eqnarray*}&&E=\displaystyle \frac{{hc}}{\lambda },\end{eqnarray*}$$

where h is Planck's constant, c is the speed of light and λ is the wavelength of the photon.

In the Beer–Lambert expression G(x), ϕ is usually considered as a photon flux and assumed to be a constant φ. The majority of the literature agrees upon φ = 10¹⁷ cm⁻²s⁻¹ [37–39, 48, 49]. Nevertheless, other values used for φ include 1.5 × 10¹⁷ [1], 4.5 × 10¹⁵ [20], 9.87 × 10¹⁶ [41], 5 × 10¹⁶ [50] and 8.19 × 10¹⁶ cm⁻²s⁻¹ [51].

Instead of assuming ϕ as a constant φ, we can use Planck's law, or an approximation to Planck's law to prescribe ϕ as a function of wavelength [52], namely

$$\begin{eqnarray}&&\phi (\lambda )=\displaystyle \frac{2\pi {{hc}}^{2}}{{\lambda }^{5}}{\left[{e}^{\tfrac{{hc}}{\lambda {k}_{B}T}}-1\right]}^{-1}.\end{eqnarray} \tag{ 18 }$$

Integrating ϕ over all wavelengths, we have

$$\begin{eqnarray*}\begin{array}{rcl}\varphi & = & {\displaystyle \int }_{0}^{\infty }\phi (\lambda )d\lambda ={\displaystyle \int }_{0}^{\infty }\displaystyle \frac{2\pi {{hc}}^{2}}{{\lambda }^{5}}{\left[{e}^{\tfrac{{hc}}{\lambda {k}_{B}T}}-1\right]}^{-1}d\lambda \\ & = & (2\pi {c}^{2}h){\left(\displaystyle \frac{{k}_{B}T}{{hc}}\right)}^{4}{\displaystyle \int }_{0}^{\infty }{y}^{3}{\left({e}^{y}-1\right)}^{-1}{dy}\quad \ \mathrm{with}\ y=\displaystyle \frac{{hc}}{\lambda {k}_{B}T}\\ & = & \displaystyle \frac{2\pi {\left({k}_{B}T\right)}^{4}}{{h}^{3}{c}^{2}}{\rm{\Gamma }}(4)\zeta (4),\end{array}\end{eqnarray*}$$

where Γ is the gamma function defined by ${\rm{\Gamma }}(z)={\int }_{0}^{\infty }{x}^{z-1}{e}^{-x}{dx}$ and ζ is the Riemann zeta function, given in its integral definition by [53]

$$\begin{eqnarray*}&&\zeta (s)=\displaystyle \frac{1}{{\rm{\Gamma }}(s)}{\int }_{0}^{\infty }{x}^{s-1}{\left({e}^{x}-1\right)}^{-1}{dx}.\end{eqnarray*}$$

As Γ(4) = 3! = 6 and ζ(4) = π⁴/90 [54], we obtain

$$\begin{eqnarray}&&\varphi =\displaystyle \frac{2{\pi }^{5}{k}_{B}^{4}{T}^{4}}{15{c}^{2}{h}^{3}}.\end{eqnarray} \tag{ 19 }$$

Under the standard temperature T = 5800 K, we have φ ≈ 6.4168 × 10⁷ Wm⁻². This value justifies the assumption of a constant φ = 10¹⁷ cm⁻²s⁻¹ as used in [37–39, 48, 49].

Alternatively, Anta et al [40] define the incident photon flux as

$$\begin{eqnarray}&&\phi (\lambda )=\displaystyle \frac{2{Fc}}{{\lambda }^{4}}{e}^{-\tfrac{{hc}}{\lambda {k}_{B}T}},\end{eqnarray} \tag{ 20 }$$

where $F\in [0,1]$ is an irradiance coefficient used to measure the intensity of the sunlight. Integrating this form of ϕ over all wavelengths gives

$$\begin{eqnarray*}\begin{array}{rcl}\varphi & = & {\displaystyle \int }_{0}^{\infty }\phi (\lambda )d\lambda ={\displaystyle \int }_{0}^{\infty }\displaystyle \frac{2{Fc}}{{\lambda }^{4}}{e}^{-\tfrac{{hc}}{\lambda {k}_{B}T}}d\lambda \\ & = & {\displaystyle \int }_{0}^{\infty }2{Fc}{\left(\displaystyle \frac{{k}_{B}T}{{hc}}\right)}^{4}{x}^{4}{e}^{-x}\left(\displaystyle \frac{{hc}}{{k}_{B}T}\right){x}^{-2}{dx}\quad \ \mathrm{with}\ x=\displaystyle \frac{{hc}}{\lambda {k}_{B}T}\\ & = & 2{Fc}{\left(\displaystyle \frac{{k}_{B}T}{{hc}}\right)}^{3}{\displaystyle \int }_{0}^{\infty }{x}^{2}{e}^{-x}{dx}=2F\displaystyle \frac{{k}_{B}^{3}{T}^{3}}{{h}^{3}{c}^{2}}{\rm{\Gamma }}(3)=4F\displaystyle \frac{{k}_{B}^{3}{T}^{3}}{{h}^{3}{c}^{2}}.\end{array}\end{eqnarray*}$$

Using F = 1 and T = 5960 K [40], we find φ ≈ 8.5237475580 × 10²⁵m⁻²s⁻¹.

Finally, Malyukov et al [55] propose the integral equation for ϕ as

$$\begin{eqnarray}&&\varphi ={\int }_{{\lambda }_{\min }}^{{\lambda }_{\max }}C{\left(\displaystyle \frac{{R}_{s}}{{r}_{0}}\right)}^{2}\left(\displaystyle \frac{2\pi {cF}}{{\lambda }^{4}}\right){\left({e}^{\tfrac{{hc}}{\lambda {k}_{B}T}}-1\right)}^{-1}d\lambda ,\end{eqnarray} \tag{ 21 }$$

where C is a coefficient used to model the atmosphere's influence on available photons for the DSSC, R_s is the radius of the Sun and r₀ is the average distance from Earth to the Sun. The parameter F has the same meaning as in the model by Anta et al [40]. Integrating this expression over all possible wavelengths, we find

$$\begin{eqnarray*}\begin{array}{rcl}\varphi & = & {\displaystyle \int }_{0}^{\infty }C{\left(\displaystyle \frac{{R}_{s}}{{r}_{0}}\right)}^{2}\left(\displaystyle \frac{2\pi {cF}}{{\lambda }^{4}}\right){\left({e}^{\tfrac{{hc}}{\lambda {k}_{B}T}}-1\right)}^{-1}d\lambda \\ & = & 2\pi {cCF}{\left(\displaystyle \frac{{R}_{s}}{{r}_{0}}\right)}^{2}{\displaystyle \int }_{0}^{\infty }{\lambda }^{-4}{\left({e}^{\tfrac{{hc}}{\lambda {k}_{B}T}}-1\right)}^{-1}d\lambda \\ & = & 2\pi {cCF}{\left(\displaystyle \frac{{R}_{s}}{{r}_{0}}\right)}^{2}{\left(\displaystyle \frac{{k}_{B}T}{{hc}}\right)}^{3}{\displaystyle \int }_{0}^{\infty }{x}^{2}{\left({e}^{x}-1\right)}^{-1}{dx}\quad \ \mathrm{with}\ x=\displaystyle \frac{{hc}}{\lambda {k}_{B}T}\\ & = & 2\pi {cCF}{\left(\displaystyle \frac{{R}_{s}}{{r}_{0}}\right)}^{2}{\left(\displaystyle \frac{{k}_{B}T}{{hc}}\right)}^{3}{\rm{\Gamma }}(3)\zeta (3)\\ & = & 4\pi {cCF}{\left(\displaystyle \frac{{R}_{s}}{{r}_{0}}\right)}^{2}{\left(\displaystyle \frac{{k}_{B}T}{{hc}}\right)}^{3}\zeta (3).\end{array}\end{eqnarray*}$$

Though there is no analytical expression for ζ(3), the result is a finite number known as Apéry's constant [56]. With R_s  = 6.9951 × 10³ m [57], r₀ = 1.495978707 × 10¹¹ m [57] and F = 1, we have φ ≈ 4.8144702590 × 10¹¹m⁻²s⁻¹.

The values of φ from the literature and from the above derivations are summarised in table 2 (given in units m⁻²s⁻¹).

Table 2. Numerical values for the constant φ used in DSSC modelling (in m⁻²s⁻¹).

Value	Author	References
9.87 × 1020	Lee et al 2004	[41]
1021	Gómez and Salvador 2005	[37]
8.52 × 1025	Anta et al 2006	[40]
4.5 × 1019	Barnes et al 2009	[20]
1.5 × 1021	Andrade et al 2011	[1]
5 × 1020	Bertoluzzi and Ma 2013	[50]
8.19 × 1020	Shi et al 2013	[51]
4.81 × 1011	Malyukov et al 2014	[55]

2.3.2. Absorption coefficient

The wavelength-dependent absorption coefficient α(λ) (sometimes referred to as α_ab to avoid ambiguity with other parameters with the same name) is usually assumed to be a constant, despite a clear relationship with wavelength as shown in figures 2 and 3.

Other model for α is given by Ni et al [38] where the absorption coefficient is assumed to be dependent on porosity p of the nanoporous semiconductor, namely

$$\begin{eqnarray}&&\alpha (p)=2.97\times {10}^{4}{p}^{2}.\end{eqnarray} \tag{ 22 }$$

In contrast to the electron generation, electron recombination has seen a great variety of modelling approaches. The first such model is shown in Södergren et al [36] which is given by

$$\begin{eqnarray*}&&R(n)=-\displaystyle \frac{n-{n}_{{eq}}}{{\tau }_{0}},\end{eqnarray*}$$

where n_eq is the dark equilibrium density and τ₀ is the electron lifetime. Referred to as the first order recombination [35, 36], this model is used frequently for its simplicity (used in [20, 37–39, 41, 46, 49, 55]).

Simple extension to this recombination term includes introducing nonlinear recombination rate (as considered in Fisher et al [58] and Barnes et al [48]):

$$\begin{eqnarray*}&&R(n)=-{k}_{R}{\left(n-{n}_{{eq}}\right)}^{\beta },\end{eqnarray*}$$

where k_R  = 1/τ₀ is a recombination constant and β is the order of recombination. In particular, Fisher et al [58] consider second order recombination (β = 2) and Barnes et al [48] consider 0 < β < 1.

A slight modification of this model is given in Anta et al [40] as

$$\begin{eqnarray*}&&R(n)=-{k}_{R}{\left(\displaystyle \frac{n}{{n}_{{eq}}}\right)}^{\beta }\left(n-{n}_{{eq}}\right).\end{eqnarray*}$$

This recombination term is used in [40] in conjunction with the diffusion coefficient $D(n)={D}_{0}{\left(n/{n}_{{eq}}\right)}^{\beta }$. The extra factor (n − n_eq ) is introduced to facilitate the effectiveness of recombination based on the electron density and its relative distance from equilibrium.

Finally, in Nithyanandam et al [2] and Tanaka [59], a nonlinear recombination model is proposed to incorporate the role of the electrolytes in the DSSCs:

$$\begin{eqnarray*}&&{R}_{{e}^{-}}(x)=-{k}_{e}\left({n}_{{e}^{-}}(x)\sqrt{\displaystyle \frac{{n}_{{I}_{3}^{-}}(x)}{{n}_{{I}^{-}}(x)}}-\overline{{n}_{{e}^{-}}}\sqrt{\displaystyle \frac{\overline{{n}_{{I}_{3}^{-}}}}{{\overline{{n}_{{I}^{-}}}}^{3}}}{n}_{{I}^{-}}\right),\end{eqnarray*}$$

where the bar notation denotes the equilibrium values. This equation together with the model proposed by Nithyanandam et al [2] are not popularly used for DSSCs. This is because the model is based on a drift-diffusion equation, and since the effect of drift is negligible in DSSCs [22, 33], pure diffusion models are more widely adopted for DSSCs.

The diffusion equation given by (8) is primarily concerned with the electron density in the conduction band. Here, we summarise the extended model of DSSCs that also include the effect of electrolyte concentrations. This model allows insights into the distribution of charged particles in DSSCs beyond the nanoporous semiconductor.

The earliest model in this area was published in 1996 by Papageorgiou et al [60]. Although this model does not consider the influence of electrons in the conduction band on the electrolytes, it provides the framework for such extension, which is investigated by Andrade et al [1]. In [1], the iodide-triiodide electrolyte couple is modelled by the following diffusion equations:

$$\begin{eqnarray}&&{D}_{{I}^{-}}\displaystyle \frac{{\partial }^{2}{n}_{{I}^{-}}}{\partial {x}^{2}}-\displaystyle \frac{3}{2p}\left[{\eta }_{{inj}}\alpha {\phi }_{0}{e}^{-\alpha x}-\displaystyle \frac{{n}_{e}(x,t)-{n}_{{eq}}}{{\tau }_{e}}-\displaystyle \frac{\partial {n}_{e}}{\partial t}\right]=\displaystyle \frac{\partial {n}_{{I}^{-}}}{\partial t},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{D}_{{I}_{3}^{-}}\displaystyle \frac{{\partial }^{2}{n}_{{I}_{3}^{-}}}{\partial {x}^{2}}+\displaystyle \frac{1}{2p}\left[{\eta }_{{inj}}\alpha {\phi }_{0}{e}^{-\alpha x}-\displaystyle \frac{{n}_{e}(x,t)-{n}_{{eq}}}{{\tau }_{e}}-\displaystyle \frac{\partial {n}_{e}}{\partial t}\right]=\displaystyle \frac{\partial {n}_{{I}_{3}^{-}}}{\partial t},\end{eqnarray} \tag{ 24 }$$

where ${n}_{{I}^{-}}$ and ${n}_{{I}_{3}^{-}}$ are the densities of the iodide and triiodide, respectively, ${D}_{{I}^{-}}$ and ${D}_{{I}_{3}^{-}}$ are the diffusion coefficients of the iodide and triiodide and p is the porosity of ${\mathrm{TiO}}_{2}$. The accompanying initial and boundary conditions result from assuming a given initial concentration, a net flux of zero at the TCO electrode and a constant total concentration throughout the operation [1]. Mathematically, we can express these conditions as

$$\begin{eqnarray}&&{n}_{{I}^{-}}(x,0)={n}_{{I}^{-}}^{{init}},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{n}_{{I}_{3}^{-}}(x,0)={n}_{{I}_{3}^{-}}^{{init}},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial {n}_{{I}^{-}}}{\partial x}\right|}_{x=0}=0,\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial {n}_{{I}_{3}^{-}}}{\partial x}\right|}_{x=0}=0,\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{\int }_{0}^{d}{n}_{{I}^{-}}(x,t){dx}={{dn}}_{{I}^{-}}^{{init}},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{\int }_{0}^{d}{n}_{{I}_{3}^{-}}(x,t){dx}={{dn}}_{{I}_{3}^{-}}^{{init}},\end{eqnarray} \tag{ 30 }$$

where ${n}_{{I}^{-}}^{{init}}$ and ${n}_{{I}_{3}^{-}}^{{init}}$ are respectively the initial concentrations of the iodide and triiodide.

So far, this is the only model for the electrolyte couple under pure diffusion, and is used in a number of papers, such as [1, 25, 35, 55].

The purpose of modelling electron density by diffusion is to use the computed conduction band electron density to determine the current-voltage relationship in DSSCs. From there, we may calculate efficiency. Included with the diffusion model given by Södergren et al [36], the paper also provides the derivation of the J − V characteristics of a DSSC based on the diode equation.

The diode equation for current J as a function of voltage V is given by

$$\begin{eqnarray}&&J(V)={J}_{{sc}}-{J}_{0}\left({e}^{\tfrac{{qV}}{{{mk}}_{B}T}}-1\right),\end{eqnarray} \tag{ 31 }$$

where J_sc is the short-circuit current density, J₀ is the reverse saturation coefficient, q is the electron charge, m is the ideality factor, k_B is Boltzmann's constant and T is the temperature of the DSSC. By construction J(0) = J_sc . Likewise, we define the open-circuit voltage V_oc to be the solution to J(V_oc ) = 0.

We obtain the short-circuit current density J_sc and the open-circuit voltage V_oc from the electron density n(x), which is the solution to the diffusion equation:

$$\begin{eqnarray}&&{D}_{0}\displaystyle \frac{{d}^{2}n}{{{dx}}^{2}}+\varphi \alpha {e}^{-\alpha x}-\displaystyle \frac{n(x)-{n}_{0}}{\tau }=0,\end{eqnarray} \tag{ 32 }$$

under the boundary conditions given in section 2.1. The photocurrent density J is given by [61]

$$\begin{eqnarray*}&&J={\left.{{qD}}_{0}\displaystyle \frac{\partial n}{\partial x}\right|}_{x=0},\end{eqnarray*}$$

where q is the electron charge. As the short-circuit current density J_sc is found by this equation under short-circuit conditions (in which the bias voltage V = 0 for the electron density n_sc ), this yields the expression [35]

$$\begin{eqnarray}&&{J}_{{sc}}={\left.{{qD}}_{0}\displaystyle \frac{\partial {n}_{{sc}}}{\partial x}\right|}_{x=0},=\displaystyle \frac{q\phi L\alpha \left(-L\alpha \cosh (\tfrac{d}{L})+\sinh (\tfrac{d}{L})+\alpha {{Le}}^{-\alpha d}\right)}{\left(1-{L}^{2}{\alpha }^{2})\cosh (\tfrac{d}{L}\right)},\end{eqnarray} \tag{ 33 }$$

where $L=\sqrt{{D}_{0}\tau }$ is the electron diffusion length. The open-circuit voltage is found by solving equation (8) under the alternative boundary condition ${\left.\displaystyle \frac{\partial n}{\partial x}\right|}_{x=0}=0$ to find the open-circuit electron density n_oc . Therefore, the open-circuit voltage is found by

$$\begin{eqnarray}&&{V}_{{oc}}=\displaystyle \frac{{{mk}}_{B}T}{q}\mathrm{ln}\left[\displaystyle \frac{{n}_{{oc}}(0)}{{n}_{{eq}}}\right],=\ \displaystyle \frac{{{mk}}_{B}T}{q}\mathrm{ln}\left[1+\displaystyle \frac{\phi \alpha \tau \left(\alpha L\cosh (\tfrac{d}{L})-\alpha {{Le}}^{-\alpha d}-\sinh (\tfrac{d}{L})\right)}{{n}_{0}\left({L}^{2}{\alpha }^{2}-1\right)\sinh (\tfrac{d}{L})}\right].\end{eqnarray} \tag{ 34 }$$

Using the expressions for J_sc and V_oc , we obtain the dark saturation current density as given by Södergren et al [36] as

$$\begin{eqnarray*}&&{J}_{0}=\displaystyle \frac{{J}_{{sc}}}{{e}^{\tfrac{{{qV}}_{{oc}}}{{{mk}}_{B}T}}-1}=\displaystyle \frac{{{qD}}_{0}{n}_{{eq}}}{L}\tanh \left(\displaystyle \frac{d}{L}\right).\end{eqnarray*}$$

The performance of DSSCs can be determined from the efficiency of the devices. The efficiency η is defined as [62]

$$\begin{eqnarray}&&\eta =\displaystyle \frac{{J}_{{sc}}{V}_{{oc}}\mathrm{FF}}{{P}_{{in}}}=\displaystyle \frac{{P}_{\max }}{{P}_{{in}}}=\displaystyle \frac{{V}_{\max }{J}_{\max }}{{P}_{{in}}},\end{eqnarray} \tag{ 35 }$$

where P_in is the total power available from sunlight. For the AM1.5 solar spectrum, P_in is taken to be 1000 Wm⁻² [63]. We calculate power by the standard formula P = JV, which takes the form

$$\begin{eqnarray}&&P(V)={VJ}(V)=V\left[{J}_{{sc}}-{J}_{0}\left({e}^{\tfrac{{qV}}{{{mk}}_{B}T}}-1\right)\right].\end{eqnarray} \tag{ 36 }$$

To find the maximum power point $({V}_{\max },{J}_{\max })$, we optimise P over V to obtain V_max and then compute ${J}_{\max }=J({V}_{\max })$. Differentiating both sides of equation (36) with respect to V, we obtain

$$\begin{eqnarray*}\begin{array}{rcl}\displaystyle \frac{{dP}}{{dV}} & = & {J}_{{sc}}+{J}_{0}-{J}_{0}{e}^{\tfrac{{qV}}{{{mk}}_{B}T}}\left[\displaystyle \frac{{qV}}{{{mk}}_{B}T}+1\right],\\ \displaystyle \frac{{d}^{2}P}{{{dV}}^{2}} & = & -\displaystyle \frac{{{qJ}}_{0}}{{{mk}}_{B}T}\left[\displaystyle \frac{{qV}}{{{mk}}_{B}T}+2\right].\end{array}\end{eqnarray*}$$

As all parameters in the expression for $\displaystyle \frac{{d}^{2}P}{{{dV}}^{2}}$ are positive, we get $\displaystyle \frac{{d}^{2}P}{{{dV}}^{2}}\lt 0$ for all voltages $V\in [0,{V}_{{oc}}]$. Thus, all critical points will be maximum values in this domain. Using the Lambert-W function, we obtain an explicit expression for the maximum power point, given by

$$\begin{eqnarray*}\begin{array}{rcl}{V}_{\max } & = & \displaystyle \frac{{{mk}}_{B}T}{q}\left[W\left(e\displaystyle \frac{{J}_{0}+{J}_{{sc}}}{{J}_{0}}\right)-1\right],\\ {J}_{\max } & = & \displaystyle \frac{({J}_{{sc}}+{J}_{0})\left[W\left(e\tfrac{{J}_{{sc}}+{J}_{0}}{{J}_{0}}\right)-1\right]}{W\left(e\tfrac{{J}_{{sc}}+{J}_{0}}{{J}_{0}}\right)}.\end{array}\end{eqnarray*}$$

Therefore, the maximum power output of the DSSC is given by

$$\begin{eqnarray*}&&{P}_{\max }=P({V}_{\max })=\displaystyle \frac{{{mk}}_{B}T({J}_{{sc}}+{J}_{0}){\left[W\left(e\tfrac{{J}_{{sc}}+{J}_{0}}{{J}_{0}}\right)-1\right]}^{2}}{{qW}\left(e\tfrac{{J}_{{sc}}+{J}_{0}}{{J}_{0}}\right)}.\end{eqnarray*}$$

We compute the short-circuit current density J_sc and the open-circuit voltage V_oc under each model by using the expressions given in equations (33) and (34). In section 4, we calculate efficiencies for DSSCs using these equations.

In addition to efficiency η, the incident photon to current efficiency (IPCE) curve can be used to gauge the performance of DSSCs with regard to photons of a particular wavelength. While it is used as a performance indicator, the short-circuit current density J_sc may be calculated using the IPCE(λ) by [6]

$$\begin{eqnarray}&&{J}_{{sc}}=\displaystyle \frac{{hc}}{q}\int \lambda \mathrm{IPCE}(\lambda )\phi (\lambda )d\lambda ,\end{eqnarray} \tag{ 37 }$$

where h is Planck's constant, c is the speed of light, q is the electron charge and ϕ is the incident radiative flux. There are several definitions for IPCE(λ). The first definition is directly from O'Regan and Grätzel [4],

$$\begin{eqnarray*}&&\mathrm{IPCE}(\lambda )=\mathrm{LHE}(\lambda ){\phi }_{{inj}}{\eta }_{e},\end{eqnarray*}$$

where LHE(λ) is the wavelength dependent light harvesting efficiency, ϕ_inj is the electron injection efficiency (also denoted by η_inj ) and η_e is the electron collection efficiency.

An alternative definition for IPCE $(\lambda )$ is given in Tanaka [59],

$$\begin{eqnarray}&&\mathrm{IPCE}(\lambda )=\displaystyle \frac{{{hcJ}}_{{sc}}}{q\lambda \phi (\lambda )},\end{eqnarray} \tag{ 38 }$$

where J_sc is the short-circuit current density and ϕ is the incident radiative flux.

We note that the quantity hc/q is often replaced by the numerical value 1240 [3, 14] or 1239 [64]. This value is obtained from converting h, c and q to SI units, i.e.

$$\begin{eqnarray*}&&\displaystyle \frac{{hc}}{q}\approx 1240\ {\mathrm{kgm}}^{2}{\mathrm{nms}}^{-3}{{\rm{A}}}^{-1}.\end{eqnarray*}$$

Finally, Barnes and O'Regan [48] define IPCE(λ) as [48]

$$\begin{eqnarray*}&&\mathrm{IPCE}=-{\left.\displaystyle \frac{{D}_{0}}{\phi }\displaystyle \frac{\partial {n}_{c}}{\partial x}\right|}_{x=0},\end{eqnarray*}$$

where D₀ is the diffusion constant for electron density in the DSSC, ϕ is the incident radiative flux and n_c is the conduction band electron density of the DSSC.

Consider a special case of equation (8) as studied by Anta et al [40]. That is, suppose $D(n)={D}_{0}{\left[\tfrac{n}{{n}_{{eq}}}\right]}^{\beta }$, $G(x)=\varphi \alpha {e}^{-\alpha x}$ and $R(n)={k}_{R}{\left[\tfrac{n}{{n}_{{eq}}}\right]}^{\beta }\left(n-{n}_{{eq}}\right)$. Using the change of variables $\bar{n}=\tfrac{n}{{n}_{{eq}}}$, $\bar{x}=\tfrac{x}{d}$ and $\bar{t}=\tfrac{{D}_{0}t}{{d}^{2}}$, we can rewrite the diffusion equation as

$$\begin{eqnarray}&&\displaystyle \frac{\partial \bar{n}}{\partial \bar{t}}=\displaystyle \frac{\partial }{\partial \bar{x}}\left[{\bar{n}}^{\beta }\displaystyle \frac{\partial \bar{n}}{\partial \bar{x}}\right]+\mu {e}^{-\nu \bar{x}}-\xi {\bar{n}}^{\beta }\left[\bar{n}-1\right],\end{eqnarray} \tag{ 39 }$$

where $\mu =\tfrac{{d}^{2}\varphi \alpha }{{D}_{0}{n}_{{eq}}}$, ν = α d and $\xi =\tfrac{{k}_{R}{d}^{2}}{{D}_{0}}$ are our nondimensional constants. Under short-circuit conditions, the boundary conditions become

$$\begin{eqnarray}&&\bar{n}(0,\bar{t})={\left.{e}^{\omega },\displaystyle \frac{\partial \bar{n}}{\partial \bar{x}}\right|}_{\bar{x}=1}=0,\quad \bar{n}(\bar{x},0)={e}^{\omega },\end{eqnarray} \tag{ 40 }$$

where $\omega =\tfrac{{qV}}{{{mk}}_{B}T}$.

Using the values of parameters given in Anta et al [40] and Gacemi et al [35], we obtain the values of nondimensional constants μ, ν and ξ for the nonlinear diffusion equation as shown in table 3.

Assuming β = 0, (39) reduces to a linear PDE which has an analytical solution (provided by Maldon et al in their 2020 paper [5]). The analytical solution is given by

$$\begin{eqnarray}&&n(x,t)=1+{{Ae}}^{\sqrt{\xi }x}+{{Be}}^{-\sqrt{\xi }x}-\displaystyle \frac{\mu }{{\nu }^{2}-\xi }{e}^{-\nu x}+\displaystyle \sum _{k=0}^{\infty }{C}_{k}\sin \left(\displaystyle \frac{(2k+1)\pi }{2}x\right){e}^{-\left[{\left(\tfrac{(2k+1)\pi }{2}\right)}^{2}+\xi \right]t},\end{eqnarray} \tag{ 41 }$$

where A, B, and C_k are constants given respectively as

$$\begin{eqnarray*}\begin{array}{rcl}A & = & -\displaystyle \frac{\mu \nu {e}^{\sqrt{\xi }-\nu }+{\xi }^{\tfrac{3}{2}}\left({e}^{\omega }-1\right)-\sqrt{\xi }\left({e}^{\omega }{\nu }^{2}-{\nu }^{2}+\mu \right)}{\sqrt{\xi }\left({\nu }^{2}-\xi \right)\left({e}^{2\sqrt{\xi }}+1\right)},\\ B & = & \displaystyle \frac{{e}^{\sqrt{\xi }-\nu }\left[{e}^{\sqrt{\xi }+\nu +\omega }\left(\sqrt{\xi }{\nu }^{2}-{\xi }^{\tfrac{3}{2}}\right)+\mu \nu +{e}^{\sqrt{\xi }+\nu }\left({\xi }^{\tfrac{3}{2}}+\sqrt{\xi }(\mu -{\nu }^{2})\right)\right]}{\sqrt{\xi }\left({\nu }^{2}-\xi \right)\left({e}^{2\sqrt{\xi }}+1\right)},\\ {C}_{k} & = & -2{\displaystyle \int }_{0}^{1}\sin \left(\displaystyle \frac{(2k+1)\pi }{2}x\right)\left[{{Ae}}^{\sqrt{\xi }}+{{Be}}^{-\sqrt{\xi }}-\displaystyle \frac{\mu }{{\nu }^{2}-\xi }{e}^{-\nu x}\right]{dx}.\end{array}\end{eqnarray*}$$

Unlike other existing solutions in the literature (such as that of [22]), equation (41) not only presents the nondimensional solution but also incorporates bias voltage V despite the majority of papers considering only short-circuit conditions (such as [1] and [65]).

In figure 6, we plot this solution under short-circuit conditions using 100 terms of the infinite series for x ∈ [0, 1], t ∈ [0, 1] and the parameters given in table 3.

Zoom In Zoom Out Reset image size

Based on the Lie symmetry analysis carried out by Maldon et al in their 2020 paper, we find that physically relevant versions of the nonlinear PDE (39) do not exhibit symmetry solutions [5]. Thus, we solve this equation numerically, as shown in [40, 66]. In [67], we solve equation (39) with the boundary conditions (40) using a forward time continuous space finite difference scheme. All simulations use 100 spatial nodes and 100, 000 temporal nodes with constant values taken from table 3.

In figure 7, we see that the numerical solution greatly resembles the exact solution for the linear case, suggesting that nonlinear diffusion only affects the overall magnitude of the solution. Nevertheless the nonlinear diffusion coefficient significantly lowers the overall electron density, corresponding to stronger traps in the ${\mathrm{TiO}}_{2}$ semiconductor [68].

Zoom In Zoom Out Reset image size

We investigate the effect of the parameters μ, ν and ξ on the solution profile for the linear case β = 1 as given by (39).

In figure 8, we plot the numerical solution for t = 1 for the electron density to elucidate the effect of β. We see that higher values of β lead to an overall decrease in electron density. Furthermore, numerical solutions reach an equilibrium faster with increased values for β. Given that β governs the density of trap states in a DSSC (higher values of β leading to deeper traps [40]), this result shows that the nonlinear diffusion mechanism is functioning in equation (39) as expected.

Zoom In Zoom Out Reset image size

Figure 9 shows the relationship between μ and the electron density. By varying μ by several powers of 10, we see that the increase in values of μ leads to a significant overall increase in electron densities. This is as expected since μ is the main coefficient in the source term $G(x)=\mu {e}^{-\nu x}$.

Zoom In Zoom Out Reset image size

Next, we consider the effect of ν on the electron density. For this study, we run simulations varying values of ν by different multiples of the original value in table 3. Since $G(x)=\mu {e}^{-\nu x}$, an increase in ν leads to a decrease in the electron generation, which in turn decreases the electron density as we see in figure 10.

Zoom In Zoom Out Reset image size

Finally, figure 11 demonstrates the effect of the recombination parameter ξ on the electron density. The increased values of ξ correspond to higher levels of recombination, which is reflected in figure 11 as the electron density decreases for higher values of ξ.

Zoom In Zoom Out Reset image size