Exploring Nonlinear Diffusion Equations for Modelling Dye-Sensitized Solar Cells

Dye-sensitized solar cells offer an alternative source for renewable energy by means of converting sunlight into electricity. While there are many studies concerning the development of DSSCs, comprehensive mathematical modelling of the devices is still lacking. Recent mathematical models are based on diffusion equations of electron density in the conduction band of the nano-porous semiconductor in dye-sensitized solar cells. Under linear diffusion and recombination, this paper provides analytical solutions to the diffusion equation. Further, Lie symmetry analysis is adopted in order to explore analytical solutions to physically relevant special cases of the nonlinear diffusion equations. While analytical solutions may not be possible, we provide numerical solutions, which are in good agreement with the results given in the literature.


Introduction
Given the increasing demand for affordable renewable energy, solar cell research seeks ways to alleviate the relatively high cost of materials while maintaining high efficiency. In 1991, O'Regan and Grätzel [1] introduced dye-sensitized solar cells (DSSCs) and demonstrated their potential to provide an alternative solution for solar cell technology. In particular, expansive research areas on determining new low cost materials, optimal structures, and configurations of the device have emerged to maximise the performance and efficiency of DSSCs.
A DSSC is composed of four primary components; a photosensitive dye, a nano-porous semiconductor, an electrolyte couple, and a counter electrode. Traditionally, the Ruthenium (II) photosensitive dye, the nano-porous semiconductor TiO 2 , the iodide-triiodide electrolyte couple I − -I − 3 , and a platinised counter electrode are typically used in DSSCs [1,2]. As depicted in Figure 1, sunlight excites dye molecules to a high energy state. This provides energy for the dye to donate an electron to the nano-porous semiconductor (a process known as electron injection), which leaves the DSSC to power the load. Electrons are reintroduced through the counter electrode, which returns electrons back to the photosensitive dye through the redox electrolyte couple.
To determine the key performance parameters of a DSSC (including its efficiency), it is necessary to derive the current-voltage relationship. To obtain this, one modelling approach is to use a drift-diffusion model based on studies of junction solar cells [3,4]. However, the electric field that is usually present in typical junction solar cells is considered to have an insignificant influence on electron transport in DSSCs. As such, it can be argued that drift-diffusion models may not be appropriate for modelling DSSCs [5]. Alternatively, Södergren et al. proposed a diffusion equation for the electron density in the Since the study of Södergren et al. in 1994 [2], the diffusion model for DSSCs has undergone significant development. In particular, in 1996, Cao et al. added time-dependence to the diffusion equation and introduced a nonlinear electron diffusion term into their partial differential equation (PDE) [6]. While the model proposed by Cao et al. [6] was only concerned with the electron density in the nano-porous semiconductor, Papageorgiou et al. [7] considered the effect of electrolyte concentrations, to allow understanding of the distribution of electron density in DSSCs beyond the nano-porous semiconductor.
In 2006, Anta et al. [8] revisited the electron diffusion equation given in [6] and introduced a nonlinear term for loss of electrons at material interfaces, known as electron recombination. Later, Barnes and O'Regan [9] introduced modified recombination kinetics for the linear diffusion equation. In 2010, Bisquert and Mora-Seró studied [10] the electron diffusion equation with a constant diffusion coefficient and a modified recombination term to gain insight into the electron diffusion length. In 2011, Andrade et al. [11] combined the linear PDE for the electron density in the nano-porous semiconductor given by Cao et al. [6] and the electron density model in the electrolyte solution given by Papageorgiou et al. [7] to produce a new system of linear PDEs. These equations were also numerically investigated by Gacemi et al. [12].
In the following subsection, we briefly state the nonlinear PDE used to model the diffusion of electrons in the nano-porous semiconductor of DSSCs. Furthermore, the relationships between the electron density, the current-voltage profile, and the efficiency are given in Section 1.2.

Diffusion Equation
Given a DSSC of thickness d, let n(x, t) be the electron density at depth x ∈ [0, d] and time t ≥ 0.
Here, x = 0 denotes the location of the transparent conductive oxide (TCO), and x = d is the location of the counter electrode of the device. From [8], the electron density satisfies: where D 0 is the diffusion coefficient, n eq is the dark equilibrium electron density, ϕ is the incident photon flux, α ab is the absorption coefficient of the Ruthenium (II) dye, k R is the recombination coefficient, and m is the diffusion order. In this paper, we refer to Equation (1) as the electron diffusion equation.
To prescribe boundary conditions to Equation (1), we assume a Boltzmann approximation for the density at the TCO substrate (which is valid for nanoporous semiconductors with fast electron transfer rates [8,10,13]) and that the current density vanishes at the counter electrode [14,15]. This yields the boundary conditions: n(x, 0) = n eq e qV m I k B T , n(0, t) = n eq e qV m I k B T , where q is the standard electron charge, V is the applied bias voltage of the DSSC, m I is the diode ideality factor, k B is Boltzmann's constant, and T is the temperature of the DSSC. Short-circuit conditions may therefore be found under the special case V = 0. For open-circuit conditions, we replace the Dirichlet boundary condition at x = 0 with the Neumann boundary condition: To nondimensionalise Equation (1), we use the same scaling parameters as of Cao et al. [6], namely: with the additional nondimensional parameter ω = qV m I k B T , which depends on the value of applied voltage V. The nondimensional equation is therefore given by: where our nondimensional parameters are µ = d 2 ϕα ab D 0 n eq , ν = α ab d and ξ = k R d 2 D 0 . Dropping the bar notation, the boundary and initial conditions for V = V oc become: n(x, 0) = e ω , n(0, t) = e ω , ∂n ∂x x=1 = 0, and the conditions for open-circuit are given by: where ω oc = qV oc m I k B T . In Table 1, the values of physical parameters of DSSCs found in Anta et al. [8] and Gacemi et al. [12] are presented. These values are used to compute our nondimensional parameters µ, ν, and ξ, the values of which are also given in Table 1. We comment that this paper only considers m ≥ 0. For the special case of m = 0, Equation (4) reduces to a linear PDE for which analytical solutions can be found, as shown in Section 2. Further, Anta et al. [8] studied the special cases of m = 1, 2, 3, and Cao et al. [6] considered the case of m = 1.

Current-Voltage Characteristics
Here, we give the relationship between the electron density, the current, and the voltage and the formula for calculating the efficiency of DSSCs.
First, the current density J as a function of the electron density n is given by: Similarly, the short-circuit current density J sc is given by: where n sc is the electron density under short-circuit conditions (V = 0). Given the short-circuit current density J sc , we apply the diode equation for the voltage-dependent current density J(V) [2]: where J 0 is the dark saturation current density [2], given by: By definition of the open-circuit voltage V oc , we have J(V oc ) = 0. Therefore, the open-circuit voltage V oc is obtained from [18]:

Efficiency
To calculate efficiency, we first compute the power output of the DSSC by the standard formula P = JV. Writing J = J(V), maximising P is traditionally done over V to obtain the maximum power point P max = J(V max )V max . The efficiency of the DSSC, η, is then calculated from: where P i is the total available power from sunlight, a value typically taken as 1000 Wm −2 [9,16,17]. In the following section, we consider a special case of Equation (4) where we assume m = 0. This assumption leads to a linear diffusion equation for which exact analytical solutions are obtained as shown in Sections 2.1 and 2.2 for bias voltage and open-circuit models, respectively. In Section 3, we investigate the case of m ≥ 0 and use Lie symmetry techniques to explore the possibility of finding new analytical solutions to the nonlinear diffusion equations. Finally, Section 4 presents numerical results of the nonlinear model together with some discussion of the results and potential research directions in this area.

Analytical Solution for the Linear Electron Diffusion Equation
Assuming m = 0, Equation (4) reduces to a linear PDE of the form: This special case readily admits an analytical solution, using a standard technique as shown in [19]. Since Equation (9) is linear, the solution can be expressed as: where f is a particular solution to: Using Equation (9), we find that u(x, t) must satisfy: The boundary conditions applied for f are imposed to allow separation of variables for finding u, according to whether the model is under open-circuit conditions (resulting in two Neumann boundary conditions) or a given voltage bias V (resulting in a Dirichlet boundary condition and a Neumann boundary condition).

Short-Circuit Conditions
Under the standard bias voltage conditions (for which V = V oc ), we solve Equation (9) under the conditions: Here, we find the particular solution f v to be given by: where the constants of integration A and B are given by: Finally, using a separation of variables approach, u is given by: where the constants C k are given by: In summary, the analytical solution for Equation (9) with a bias voltage V = V oc is given by: In Figure 2, we plot the exact solution (11) under short-circuit conditions (by setting ω = 0) using the data given in Table 1 for the first 100 terms of the infinite series in MAPLE. From the figure, we can see that the electron density rises quickly from its dark equilibrium due to the influence of the exponential source term of electron generation. The electron density continues to increase until it reaches the steady-state, as expected for a linear diffusion problem.

Open-Circuit Model
For the open-circuit model, we solve Equation (9) under the following conditions: Using the same approach as before, we find the particular solution f oc for open-circuit conditions as: where the constants of integration are given by: and the nondimensionalised open-circuit voltage ω oc is given by: Using the separation of variables, we find that u oc is given by: where the constants V k are of the form: In summary, the analytical solution for the open-circuit model is given by: Using the values of the parameters given in Table 1, we plot the exact solution (14) for the first 100 terms of the infinite series in MAPLE (see Figure 3). This solution has a different profile to that of the short-circuit case shown in Figure 2. Here, the electron density is initially at an extremely high magnitude (of the order of 10 10 ). Then, it shifts slightly to match the exponential source term in accordance with the Neumann boundary conditions. We note that the solution profile of our model is in agreement with that of Gómez and Salvador [17]. We note that the discrepancy between the solutions of the short-circuit and open-circuit is due primarily to the Dirichlet boundary condition for the short-circuit model and the parameter ξ. The short-circuit electron density prescribed at the TCO electrode is small relative to the electron density found for the open-circuit model. Next, we use k R = 400 s −1 , so that ξ = 10 5 (µ and ν have the same values as in Table 1) to plot the short-circuit and open-circuit electron densities, as shown together in Figure 4. As we increase the value of ξ, which results in a higher level of recombination, we find that the discrepancy between the open-circuit and short-circuit electron densities is much smaller. Additionally, the solutions do not vary much from their initial states in comparison to Figures 2 and 3. In the following section, we consider the nonlinear electron diffusion equation for which m > 0.

Classical Lie Symmetry for the Nonlinear Electron Diffusion Equation
In this section, we seek analytical solutions for particular forms of the nonlinear diffusion equation by using classical Lie symmetry techniques. First, we rewrite Equation (4) in a general form as: The method of Lie symmetry was first developed by Sophus Lie in the late 19th Century [20]. This framework offered a systematic approach to finding transformations of variables in a differential equation that lead to a reduction of order (or integrating factors for a first order ODE). Finding these transformations amounts to solving a system of linear differential equations, which is referred to in the literature as classical symmetries. Classical symmetry analysis has been used to analyse diffusion equations for decades. The studies have led to a vast knowledge of particular structures of the diffusion equation that admit classical Lie symmetry (see, for example, [21][22][23]). Nevertheless, many special cases involving spatially dependent source terms have not yet been considered. In this paper, we look for special cases of Equation (15) that may admit analytical solutions by the classical Lie symmetry method.

Classical Lie Symmetry
To find other classical symmetries of diffusion equations that are not covered in [21], we use an algebraic computer package DIMSYM [24]. For Equation (15) in its current form, the only symmetry operator is of the form: which corresponds to an invariant under translation in time. Furthermore, DIMSYM suggests six special cases that may potentially yield further symmetries:
Given that this model is based on the assumption of electron diffusing within the conduction band of DSSCs, the diffusion coefficient must not be zero. Assuming {x dG dx , dG dx , G, 1} is linearly dependent or R(n) = 0 does not admit any new symmetries without also making assumptions about the diffusion coefficient. Note that we consider R(n) = 0 as part of other special cases.

Constant Diffusion Coefficient
Here, we suppose that D(n) is a constant. Without loss of generality, we may assume D(n) = 1. Equation (15) is therefore of the form: We note that if R(n) = an + b for some constants a and b, then Equation (16) reduces to a linear partial differential equation, which can be solved by using standard techniques. Generally, Equation (16) does not immediately admit symmetries other than Γ 1 above. Two special cases are: 1.
G is a constant.
If G is a constant, Equation (15) reduces to the case that has already been classified by Dorodnitsyn [22]. We therefore assume that G is not a constant and consider recombination terms for R that ensure {n ∂R ∂n , ∂R ∂n , R, n, 1} is linearly dependent. That is, R(n) is of the form: where c 1 , c 2 , c 3 , c 4 , and c 5 are arbitrary constants. Though there are several special cases giving rise to symmetries (found in Table 2), none of these yield analytical solutions unless R(n) is linear.

Particular Diffusion Coefficient
Recall the third order ordinary differential equation for the diffusion coefficient given in Section 3.1, We notice that only the power-law and exponential functions for D(n) are able to satisfy Equation (17) (as shown by Edwards [25], for example). Under a power-law diffusivity, the special cases from DIMSYM are that m = −1 and − 4 3 and {x dG dx , dG dx , G, 1} is linearly dependent. Since m ≥ 0 for the study of DSSCs, we assume that {x dG dx , dG dx , G, 1} is linearly dependent. We note that if G(x) is constant, the equation is of the form n t = (D(n)n x ) x + R(n), which has already been extensively studied (in [21,23,26], for example). Clarkson and Mansfield [27] give an account of a wide range of analytical solutions found for this special case. Otherwise, G(x) is of the form: for some constants A, B, and C. Table 2 contains special cases for (4) and (17) admitting symmetries. Though these special cases for G(x) admit nontrivial symmetries for both power-law and exponential functions of D(n), none of these special cases reduce to analytical solutions unless G is a constant.

Particular Source Term
From the fourth special case in Section 3.1, the electron generation G(x) and the recombination R(n) can be shown to have the form: for some constant λ. In this case, Equation (15) reduces to: for some constant A. This special case admits the symmetry Γ 2 = x ∂ ∂x + 2t ∂ ∂t , an invariant of which is the Boltzmann similarity variable s = xt − 1 2 . In this case, the above partial differential equation reduces to the ordinary differential equation: However, we find that there are no symmetries that can reduce this equation further under power-law or exponential diffusivities. Table 2 summarises all classical Lie symmetries found with DIMSYM [24]. Under these special cases, we obtain invariants of the form F(x, t, n) by solving X ∂F ∂x + T ∂F ∂t + N ∂F ∂n = 0. For example, we consider the PDE:

Summary
From Table 2, we find that the PDE admits the symmetry: Upon solving Γ(F) = 0, we find the invariants s = x + (m+1) Bm ln(t) and ϕ = nt 1 m . Using these invariants, we may reduce (18) to the ODE: This ODE may also be analysed by classical Lie symmetry, and we find that the ODE does not admit more symmetries.

Nonclassical Lie Symmetry for the Nonlinear Electron Diffusion Equation
Developed in 1969, Bluman and Cole [28] added an additional invariant surface condition, which transformed the linear equations of classical symmetry to nonlinear equations. This approach is now known as nonclassical symmetry or Q-conditional symmetry in the literature. Results on nonclassical symmetry are far more recent, but nevertheless extensive [26,27,[29][30][31].
A particular study of nonclassical symmetry developed in 1993 by Nucci [32] is the method of Heir equations, which assumes particular forms for the invariant surface conditions. Recently, Bradshaw-Hajek [33] investigated the partial differential equation: for the special cases R(x, u) = Q(u) and R(x, u) = r(x)Q(u). This study was motivated by problems in mathematical biology, and we refer to Joshi and Morrison [34] for a generalised form for R(x, u) with the additional assumption of constant diffusion coefficient (D(u) = D).

Nonclassical Symmetry Analysis for D(n) = n m
To find nonclassical symmetries for a partial differential equation of the form: we use MAPLE to compute and expand the second prolongation Γ (2) , which is the differential operator given by: where subscripts for n denote the usual partial derivatives and N (k) x i x j are the extended infinitesimals. We refer the reader to [25] for details of the differential operator Γ (2) .
This step leads to the Lie symbol: Following the nonclassical technique given by Bluman and Cole [28], we also include the invariant surface condition: Adopting a common convention for nonclassical symmetry analysis for evolution equations, we set T = 1 so the invariant surface condition (20) can be rewritten as: leading to the following series of determining equations: Solving Equation (21) as an ODE for X(n) yields: for some arbitrary functions f (x, t) and g(x, t). Upon substituting X(x, t, n), Equation (22) simplifies to: Given that the parameter m is non-negative, the general solution to this equation is given by: where h and k are arbitrary functions of x and t. Splitting (23) and (24) as a polynomial in n leads only to the trivial symmetry Γ = ∂ ∂t for general source terms G(x) and R(n), as we assume m ≥ 0. For the special case when m = 1, G(x) = x + A, and R(n) = 0, the PDE can be written as: which admits the symmetry: for some constant µ. We note that the symmetry is equivalent to the classical symmetries found for this special case (as seen in Table 2). Nevertheless, this symmetry produces the invariant variables ψ = n(t + 3µ) −3 and s = (x + A)(t + 3µ) −2 , leading to the reduction: There are no classical symmetries that may reduce this equation further.

Nonclassical Symmetry Analysis for D(n) = e mn
To find nonclassical symmetries for the partial differential equation, we use MAPLE to compute and expand the second prolongation Γ (2) with the same setup as with the power law diffusion coefficient. The determining equations are: The general solution for Equation (27) is given by: for arbitrary functions f (x, t) and g(x, t). Equation (28) implies: for some arbitrary functions h, k of x and t. Splitting each term for n in (29) and (30) leads to the trivial symmetry Γ = ∂ ∂t for general source terms G(x) and R(n). We find that the special cases yielding a nontrivial symmetry are precisely those in which G is constant, which has already been studied (in [35], for example).

Numerical Results and Discussion
In this paper, we present a nonlinear diffusion model for electron density in the conduction band of the nano-porous semiconductor TiO 2 in a DSSC. This equation admits an exact analytical solution under linear diffusion (m = 0). For the general case of m > 0, both classical and nonclassical Lie symmetry analyses indicated no new analytical solutions for the physically relevant special cases of the partial differential equations studied here. As a result, we proceeded to obtain solutions numerically.
Using a standard forward time central space finite difference scheme in MATLAB (detailed in Maldon et al. [36]), we provided numerical solutions to Equation (4) for both power-law and exponential diffusion coefficients under short-circuit conditions. We simulated for dimensional time t ∈ [0, 1].
In Figure 5, we plot the numerical solution to the following partial differential equation, where we assume m = 1 and subject to short-circuit boundary conditions. We note that D(n) = n was found to best resemble the nonlinear diffusion of electron density present in the semiconductor conduction band [6,8,37]. The solution profile in Figure 5 is similar to that of the linear case (m = 0) shown in Figure 2. However, the magnitude of the electron density differs by a factor of O(10 2 ). Since m in the diffusion coefficient D(n) = n m determines the strength of electron trap in the TiO 2 structure, we expected a lower value of the electron density for the case of m = 1 as this corresponds to a higher effect of the trap, preventing electron transport throughout the DSSC. In Figure 6, we elucidate the influence of m on the electron density as given in Equation (1). In particular, we consider m = 0, 1 4 , 1 2 , 3 4 , and m = 1. From the figure, we see that for lower power-law diffusivities, the electron density reached the steady-state slower. This case could cause potential stability issues in numerical calculations, owing to the higher overall electron densities. Conversely, the higher power-law diffusivities (stronger traps) in the TiO 2 nano-porous structure led to significantly lower electron density. For higher values of m, the solutions also reached the steady-state faster. To investigate the effect of an exponential diffusivity, in Figure 7, we plot the numerical solution n(x, t) for the following equation, for which the diffusion coefficient is D(n) = e mn . We note that the recombination term also took the same exponential coefficient (as shown in [8]). To better resemble the power-law diffusivity and avoid stability issues, m was assumed to take the form m = ln(n eq ) n eq . Under the data in Table 1, this meant that m = O(10 −21 ). The numerical solution shown in Figure 7 fundamentally resembles the power-law diffusivities presented in Figure 6. In Figure 8, we use parameter values consistent with Cao et al. [6] to solve Equation (4) and find good agreement with their results. In comparison to our studies, their diffusion coefficient was significantly higher, and the recombination term was omitted. These densities nevertheless compared favourably well as the recombination parameter ξ was very small (ξ = 10 −5 in our case).

Efficiency Calculations
Using the formulae in Sections 1.2 and 1.3, we computed the efficiency using different diffusivity functions D(n) as shown in Table 3. For each D(n), we simulated the numerical solution to time t = 10 3 and computed the short-circuit current density J sc , the open-circuit voltage V oc , and efficiency η. From Table 3, we see that linear diffusion outperformed all other diffusion coefficients. For the diffusivity, D(n) = n, √ n, n 2 , we see that they hindered the electron density so strongly that the resulting short-circuit current density J sc and the efficiency were adversely affected. From the values of J sc , V oc , and η presented in Table 3, we found that the linear diffusion model with a constant D(n) was in agreement with data of DSSCs available in the literature [1].
In summary, this paper explored the use of nonlinear diffusion equations to model electron density within the nano-porous semiconductor of DSSCs. We showed that classical and nonclassical Lie symmetry analysis indicated no physically relevant analytical solutions for the general form of the nonlinear diffusion Equation (1). However, using a finite difference method to solve (1) for particular cases of diffusion coefficients and recombination terms, we saw good agreement between our solutions and those given in the literature.
The role of nonlinear diffusion is shown by Figure 6, in which the electron density is greatly affected by the power-law. Given the depth of active traps within the TiO 2 network [8], as well as a nonlinear recombination term [10], we see the importance of using nonlinear PDEs to model diffusion in DSSCs.
Finally, we comment that in order to develop more realistic mathematical models for DSSCs, the role of the electrolytes and the efficiency of the counter electrode should be included in the model to better describe the photoelectrochemical activity in DSSCs.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: DSSC Dye-sensitized solar cell IPCE Incident photon to current efficiency PDE Partial differential equation