A Fractional Diffusion Model for Dye-Sensitized Solar Cells

Dye-sensitized solar cells have continued to receive much attention since their introduction by O’Regan and Grätzel in 1991. Modelling charge transfer during the sensitization process is one of several active research areas for the development of dye-sensitized solar cells in order to control and improve their performance and efficiency. Mathematical models for transport of electron density inside nanoporous semiconductors based on diffusion equations have been shown to give good agreement with results observed experimentally. However, the process of charge transfer in dye-sensitized solar cells is complicated and many issues are in need of further investigation, such as the effect of the porous structure of the semiconductor and the recombination of electrons at the interfaces between the semiconductor and electrolyte couple. This paper proposes a new model for electron transport inside the conduction band of a dye-sensitized solar cell comprising of TiO2 as its nanoporous semiconductor. This model is based on fractional diffusion equations, taking into consideration the random walk network of TiO2. Finally, the paper presents numerical solutions of the fractional diffusion model to demonstrate the effect of the fractal geometry of TiO2 on the fundamental performance parameters of dye-sensitized solar cells, such as the short-circuit current density, open-circuit voltage and efficiency.


Introduction
Dye-sensitized solar cells (DSSCs) were first introduced by O'Regan and Grätzel in their fundamental 1991 paper [1], providing a viable low-cost alternative for renewable solar energy. Functionally, DSSCs operate by a photosensitive dye using absorbed sunlight to inject excited electrons into a nanoporous semiconductor. This approach relives DSSCs of the need for a costly and high purity semiconductor as opposed to the predominant silicon solar cells that have been at the forefront of solar energy since 1954 [2].
Mathematical modelling for DSSCs began shortly afterwards and continues to offer unique insight into the current-voltage relationship and ultimately the efficiency of the DSSC. While some early mathematical models of DSSCs used Maxwell's equations and electric potentials [3], Gregg [4] outlined the prevalence of the photoinduced potential as more influential than electric fields in DSSCs. Södergren et al. [5] were the first to model DSSCs by a simple diffusion equation together with analytical expressions for the short-circuit current density J sc and the open-circuit 2 of 9 voltage V oc . Cao et al. [6] extended the model to include time-dependence, leading to the partial differential equation where n(x, t) is the conduction band electron density at depth x ∈ [0, d] and time t ≥ 0, D 0 is the diffusion coefficient, ϕ is the incident photon flux, α is the absorption coefficient, k R is the recombination constant and n eq is the equilibrium electron density. Anta et al. [7] considered a fully nonlinear diffusion equation for DSSCs and investigated the effect of different power-law diffusivities, which was analysed using Lie symmetry by Maldon et al. 2020 [8]. Andrade et al. 2011 [9] and Gacemi et al. 2013 [10] proposed diffusion equations for modelling the electrolyte concentrations, which were solved analytically by Maldon and Thamwattana in 2019 [11].
Fractional calculus has existed conceptually as long as calculus itself. For close to three centuries, fractional calculus had only one known application in Abel's 1823 tautochrone problem [12] until Nigmatullin [13] suggested a fractional diffusion equation for modelling media exhibiting a fractal geometry. In 2000, Henry and Wearne [14] developed the standard fractional diffusion equation based on continuous-time random-walk (CTRW) models. Henry and Wearne's model [14], together with the CTRW simulation of TiO 2 by Nelson [15], provides motivation for this paper to model DSSCs using fractional partial differential equations.
Modelling DSSCs with fractional calculus is a relatively unexplored research direction. So far, there is only one paper by Sibatov et al. in 2014 [16]. In Sibatov et al. [16], they consider the role of trap states within the TiO 2 network and its effect on the electron hole density. To address the fact that previous mathematical models for DSSCs do not generally consider the effect of the porous network in TiO 2 [6,17], we develop our model based on generalised fractional diffusion-reaction equations, taking this effect into account.
For this study, we use the Caputo fractional derivative, which is given by [18] where λ > 0 is the order of the fractional derivative, n = λ , f (n) (t) denotes the classical derivative of f with respect to its variable of order n ∈ N and Γ is the usual Gamma function The literature employs several definitions for the fractional derivative in modelling diffusion in media exhibiting a fractal geometry. Henry and Wearne [14] suggest the Riemann-Liouville definition, but also remark that the Caputo definition has also seen some use [14,19]. However, Baeumer et al. [20] comment on the use of Caputo derivatives for space-fractional diffusion models that positivity is not preserved under vanishing Neumann boundary conditions. We find that problem is alleviated by the use of a time-fractional derivative on the diffusion term of our equation. Also, the presence of spatially dependent source terms and the combination of Dirichlet and Neumann boundary conditions enjoy a greater level of compatibility with the Caputo fractional derivative [21].

Mathematical Model
In this paper, we adopt a Caputo fractional derivative in time on the diffusion term as shown by Henry and Wearne [14], resulting in the fractional partial differential equation (FPDE) where γ ∈ (0, 1] is the order of the Caputo fractional derivative in time and all other parameters retain their values as in Equation (1). Physically, the parameter γ is the exponent in the mean square displacement of the CTRW simulation [14]. The special case γ = 1 recovers the standard diffusion equation. Lower values for γ correspond to longer path lengths for electron transport through the nanoporous semiconductor [22]. The parameter γ is also strongly influenced by the porosity of the TiO 2 semiconductor. Benkstein et al. [22] found that increasing porosity led to a decrease in γ. We note that Equation (3) does not feature the term where L −1 denotes the inverse Laplace transformation. Though this term is critical for physically meaningful fractional diffusion equations [14], it vanishes in Equation (3) under the Caputo fractional derivative. We prescribe boundary conditions as found in [5], with a Dirichlet boundary condition at x = 0 and a Neumann boundary condition at x = d together with a prescribed initial condition, namely n(x, 0) = n eq e qV mk B T , where q is the standard electron charge, V is the bias voltage, m is the diode ideality factor, k B is Boltzmann's constant and T is the temperature of the DSSC. The diode equation is commonly used to compute the current-density relationship for solar cells, in which the current J as a function of bias voltage V is given by where J 0 is the dark saturation current density, given by Södergren et al. [5] in the form To compute the short-circuit current density J sc we use noting that the standard flux is recovered in the special case of linear diffusion (γ = 1).
Given that the open-circuit voltage V oc satisfies J(V oc ) = 0, we may compute the open-circuit voltage from Maximising the power output P(V) = V J(V) over V, we obtain the maximum power point V max by where W is the Lambert-W function and J max = J(V max ). With P max = V max J max , we compute the efficiency η of the DSSC by where P i is the power of incident light.

Finite Difference Method
Finite difference methods (FDMs) have been used by Hu et al. [23] to solve parabolic FPDEs under the Caputo fractional derivative for fractional time derivatives, and Takeuchi et al. [24] have used finite difference methods to solve FPDEs under fractional spatial derivatives. We refer the reader to Li and Zeng [25] for a finite difference scheme for solving fractional ordinary differential equations over a finite interval under boundary conditions defined at the left boundary.
In this paper, we solve Equation (3) under a FDM scheme using expressions given by Oldham and Spanier [26]. All numerical computations are performed using numerical values of constants provided in Table 1 except for γ, which is given several values in Table 2.  [28] To numerically solve Equation (3) under boundary conditions (4)-(6) with a finite difference scheme, we use the L1 approximation for the fractional derivative given by Oldham and Spanier [26] in which where b

Nodes Determined by Boundary Conditions
To satisfy the initial condition (6), we set for all i ∈ {1, ..., N t }. Finally, for the Neumann boundary condition (5) at x = d we employ a 'ghost node' at x = d + ∆x and a central difference approximation for the first derivative at x = d to set for all i ∈ {1, ..., N t }.

Results and Discussion
Using the finite difference method, we numerically solve Equation (3) for several values of γ using 10 spatial nodes and 5000 temporal nodes over x ∈ [0, d] and t = 100s. To investigate the effect of the parameter γ, we plot numerical solutions for the special cases γ = 0.25, 0.5, 0.75 and 1 in Figure 1. Using the numerical solutions to Equation (3) and the flux estimate given by Equation (12), we are able to compute the short-circuit current density J sc and the open-circuit voltage V oc , leading to the overall efficiency η. Table 2 shows the effect of the parameter γ on these DSSC performance parameters. From Table 2, we see that the short-circuit current density J sc increases when γ increases. Though the open-circuit voltage is not affected to the same extent, the efficiency is notably lower for decreased values of γ. We note the special case γ = 1 is equivalent to the standard diffusion equation without fractional derivatives, and the efficiency η = 7.03% is in agreement with expected efficiencies for DSSCs [8].
As γ decreases, Figure 1 presents two primary trends to the numerical solution to Equation (3). Firstly, the time required to reach steady-state increases as γ decreases by comparison to the numerical solutions for the cases γ = 1 and 0.25. This result is consistent with the observation that lower values for γ imply slower diffusion, based on the CTRW simulations. Benkstein et al. [22] show that γ decreases when porosity increases from p = 0.7 to p = 0.775.
Secondly, the overall electron density is remarkably higher for the cases γ = 0.5 and 0.25 compared to γ = 1 and 0.75. This suggests that the electron density is highly sensitive to the order of the fractional derivative. The standard gradient for flux would consequently produce extreme results, requiring the fractional derivative to rectify this issue. Numerically, the finite difference scheme presents stability problems and is computationally expensive for lower values of γ.
The parameter γ denotes the exponent for the power-law in the mean-square displacement [14]. We conclude from Figure 1 that lower values of γ lead to progressively less realistic behaviour for nanoporous semiconductors used in DSSCs, as the electron density dramatically increases when γ decreases. This observation is consistent with the longer path lengths associated with low values for γ in CTRW simulation. In particular, Ni et al. [27] show that extremely low porosities (such as p = 0.1) show a significant reduction in efficiency (below 1%). From the numerical solution to Equation (3) with γ = 0.612 using 25 spatial nodes, the flux estimate given by Equation (12) and 50,000 temporal nodes to t = 500s, we find J sc = 116.9059Am −2 , V oc = 0.6245V and η = 6.0835%.

Conclusions
We propose a new mathematical model for evaluating the efficiency of dye-sensitized solar cells by using fractional diffusion to incorporate the fractal geometry of the TiO 2 semiconductor. Our results show that lower values of the mean square exponent γ lead to lower efficiencies, a result that is consistent with the literature [11,27]. In particular, figure 11 of Ni et al. [27] shows that efficiency decreases when porosity increases above p = 0.4 or decreases below p = 0.4, which suggests the relationship between porosity and efficiency is nonlinear. We note that the solution profile of the electron density presented here is similar to those obtained from nonlinear diffusion modelling [8], though the orders of magnitude differ significantly. This is due to the longer waiting times associated with lower values for γ, which slows down the diffusion process.
We also develop a finite difference scheme to numerically solve the fractional diffusion equation and provided a tenth-order estimate for obtaining the short-circuit current density. Together, this provides a comprehensive model for incorporating the effect of the random-walk behaviour of the nanoporous semiconductor on the performance of dye-sensitized solar cells.
Future consideration includes incorporating the role of the electrolyte couple by a pair of standard diffusion equations as seen in Maldon and Thamwattana [11].