Systematic derivation of a surface polarization model for planar perovskite solar cells

Increasing evidence suggests that the presence of mobile ions in perovskite solar cells can cause a current-voltage curve hysteresis. Steady state and transient current-voltage characteristics of a planar metal halide CH$_3$NH$_3$PbI$_3$ perovskite solar cell are analysed with a drift-diffusion model that accounts for both charge transport and ion vacancy motion. The high ion vacancy density within the perovskite layer gives rise to narrow Debye layers (typical width $\sim$2nm), adjacent to the interfaces with the transport layers, over which large drops in the electric potential occur and in which significant charge is stored. Large disparities between (I) the width of the Debye layers and that of the perovskite layer ($\sim$600nm) and (II) the ion vacancy density and the charge carrier densities motivate an asymptotic approach to solving the model, while the stiffness of the equations renders standard solution methods unreliable. We derive a simplified surface polarisation model in which the slow ion dynamic are replaced by interfacial (nonlinear) capacitances at the perovskite interfaces. Favourable comparison is made between the results of the asymptotic approach and numerical solutions for a realistic cell over a wide range of operating conditions of practical interest.


Introduction
Since the first use of methylammonium lead tri-halide perovskite as a sensitizer in a dyesensitized solar cell [16], and its subsequent incorporation into a novel thin film solar technology as a bulk solar absorber [18,15], the efficiency of perovskite solar cells (PSCs) has increased extremely rapidly from around 3% to above 20% [7], a level that is comparable to the standard crystalline silicon devices. This increase, along with advances in the material properties and stability of PSCs, makes this area of photovoltaic research a very hot topic [22,38]. transient behaviour has also been observed in dark current transients (whereby the cell is first held in the dark, then the applied voltage is suddenly changed and the resulting current measured) [23]. More recently, very long timescale transients lasting many hours have been observed in cell efficiency [10]. These decays in PCE can be reversed by allowing the cell to recover in the dark. Various explanations have been proposed for these transient behaviours, including (a) large trap state densities close to the interfaces with the transport layers, (b) slow ferroelectric polarisation of the perovskite material and (c) the motion of iodide (I − ) vacancies within the perovskite material [35]. As discussed in Richardson et al. [28], it is now widely accepted that the only one of these mechanisms capable of explaining the data is iodide vacancy motion.
Various approaches may be used to model PSCs ranging from atomistic density functional theory (DFT) simulations, to drift-diffusion models of charge carrier and ion motion, to lumped parameter device models (equivalent circuits). DFT calculations, while perhaps the most fundamental approach, are so computationally intensive that they are incapable of describing the behaviour of a full cell. In practice they are used to obtain estimates of macroscopic quantities, such as ion vacancy densities and mobilities, from the atomistic structures of the materials forming the device [11]. In contrast, drift-diffusion models, which are applicable on the nanometre length scale and upward, describe the motion of electrons, holes and ion vacancies. Such models have been presented and solved in a number of works [28,30,23,34,6,10,14,40,21]. However it is notable that, with the exception of two [23,28], all of these works use parameter values that are very far from realistic. This may be ascribed to the extreme numerical stiffness of the problem owing to very narrow (≈ 2 nm) Debye (boundary) layers that form as a result of ion accumulation/depletion at the edges of the perovskite layer. In order to overcome this difficulty, Richardson et al. [28] adopted a combined numerical and asymptotic approach, in which the electrical properties of the Debye layers are modelled by a nonlinear surface (Debye layer) capacitance, based on estimates for the equilibrium ion vacancy density and mobility obtained from DFT calculations performed by Eames et al. [11]. The purpose of that work was to demonstrate that experimental J-V hysteresis data could be explained by the motion of ion vacancies in the perovskite layer and so the derivation of the asymptotic solution was not given there.
The aim of this paper is to systematically derive the asymptotic approach used in the earlier work by Richardson et al. [28] and validate it against numerical solutions to the full model. A similar approach has been used for (i) asymptotic derivations of equivalent circuit models from drift-diffusion models [12,13,32] (in the context of organic solar cells, PSCs and bipolar silicon devices, respectively); (ii) a matched asymptotic analysis of np-diodes [24]; (iii) asymptotic derivations of the standard 'regional' models of semiconductors from a drift-diffusion model [31]; (iv) multidimensional models of bulk heterojunction solar cells [4,29]; and (v) the asymptotic analysis of quantum drift-diffusion models [3]. Subsequent to Richardson et al. [28], Ravishankar et al. [26] published a heuristic model similar to the surface capacitance model used in this earlier work, which they term a surface polarization model. We argue that the systematic derivation of such models from the underlying driftdiffusion equations, as here, has the significant advantage of directly relating the surface capacitances to the device physics.
This work is set out as follows. In §2, we formulate the drift-diffusion model for a PSC, non-dimensionalise and estimate the model parameters. In §3, we use formal asymptotic methods, based on the parameter estimates made in §2, to derive a hierarchy of simplified models to the full PSC model including the surface polarization model of Ravishankar et al. [26]. In §4, the results of the simplified models are compared to numerical solutions of the full PSC model and finally, in §5, we draw our conclusions.

Problem formulation
Here we consider a perovskite absorber layer, sandwiched between an ETL and an HTL (typically TiO 2 and spiro, respectively). We make the assumption that the transport layers are sufficiently highly doped that they are effectively equipotential across their width and take the same potential as their respective contacts. In the perovskite, in line with DFT calculations on its chemical structure [11], we assume there exists a high density of mobile anion vacancies, in addition to the charge carriers. The resulting dimensional model for the perovskite layer (0 < x < b), following earlier work [28], is outlined below.
Dimensional model. Conservation of holes (density p) and conduction electrons (density n) is described by where G is the photo-generation rate; R(n, p) is the bulk recombination and thermal generation rate (henceforth abbreviated to recombination rate); φ is the electric potential; j n and j p are electron-and hole-currents, respectively; and V T = kT /q is the thermal voltage. Similar equations for the conservation of positively-charged anion vacancies (density P ) and negatively-charged cation vacancies (density N) take the form where F p (and F n ) are the fluxes of the positive (and negative) ion vacancies (as opposed to the current fractions carried by these species). Both sets of equations couple to Poisson's equation for the electric potential Boundary conditions at the edges of the perovskite, x = 0 (the interface with the ETL) and x = b (the interface with the HTL) take the form where V ap is the applied voltage; V bi is the built-in potential; R l and R r are the interfacial charge recombination rates on x = 0 and x = b, respectively; and the carrier densities on the interfaces are given by the expressions (see e.g. Nelson [20]) Here, g c and g v are the effective density of states in the conduction and valence bands of the perovskite, respectively;μ n andμ p are the perovskite conduction and valence band energies, respectively. In addition, we model the highly doped ETL and HTL as metals in which E F d , the HOMO energy level of the HTL (Spiro), and E Fa , the conduction band energy of the ETL (TiO 2 ), play the roles of the Fermi levels in these materials. These equations are supplemented by initial conditions, which we choose as follows to ensure charge neutrality, The built-in voltage. This quantity can be found from (1) with boundary conditions (4) by noting that, at equilibrium, the photo-generation rate, applied voltage and electron-and hole-currents are all zero (G = 0, V ap = 0 and j p = j n = 0). The equilibrium solutions for n and p have the form in which the constants A and B are determined by the boundary conditions such that Furthermore, since the rate of thermal generation and recombination must be equal (R = 0) at equilibrium (see e.g. (7)), we require np = n 2 i . It follows that which, with parameter estimates in Table 1, turns out to be 1V ≈ 39V T .
Recombination and photo-generation. At the radiation intensities associated with sunlight, the bulk recombination rate within the perovskite, R, is believed to be predominantly trap assisted (although at higher radiation intensities bimolecular recombination becomes significant) [37]. It is therefore appropriate to model bulk recombination by the Shockley-Read-Hall rate equation (see e.g. Nelson [20] §4.5.5) where τ n and τ p are the pseudo-lifetimes of conduction electrons and holes, respectively, and k 3 is a constant related to the pseudo-lifetimes and trap state energy level (typically negligible to the other terms in the denominator of (7) when the cell is under illumination). Furthermore, Stranks et al. [37] suggest that bulk recombination is hole dominated (τ p ≫ τ n ), an assumption which is in line with that made in Richardson et al. [28]. There is still no consensus on the relative importance of interfacial recombination (at the interfaces between the perovskite and the transport layers) in comparison to bulk recombination although we note that this may be sample dependent. For example, de Quilettes et al. [9] note that recombination within the perovskite occurs predominantly at crystal boundaries, which implies the magnitude of bulk recombination is strongly dependent upon the perovskite structure. The photo-generation rate, G, is assumed to follow the Beer-Lambert law of light absorption; with light entering the device through the ETL (TiO 2 ). This has the form where F ph is the incident photon flux and α is the light absorption coefficient of the perovskite.

Non-dimensionalisation
Dimensionless variables (denoted by a star) are introduced by rescaling (i) space with the width of the perovskite layer; (ii) voltages with the thermal voltage; (iii) charge carrier densities with the typical photo-generated charge density, Π 0 (see (10)); (iv) current densities with the typical photo-generated current density, qF ph ; and (v) ion densities with the typical ion density, N 0 . The rescaling reads Here, L d is the Debye length calculated on the basis of the ion vacancy density and τ ion is the characteristic timescale for ion motion defined, respectively, by Furthermore, we take Π 0 to be the characteristic charge carrier density required to remove the photo-generated charge in the absence of an electric field whereD is a typical carrier diffusivity. The non-dimensionalisation gives rise to the following dimensionless quantities that characterise the system: The dimensionless problem. The system of equations obtained by applying the rescaling (9) to the variables in (1)- (8) is The dimensionless recombination and generation rates (for a device under constant illumination) are given by Henceforth, we drop the star superscript from the dimensionless variables.

Parameter estimates for real devices
A list of parameter estimates obtained from the literature is supplied in Table 1. Note that F ph , τ n , τ p and D + are in line with the range of values found in the literature but have been specifically chosen to give good agreement to the experimental J-V curves presented by Richardson et al. [28]. Based on this data, the dimensionless parameters, corresponding to a cell with perovskite width b = 600nm, are While, for a cell with perovskite width b = 150nm, they remain unchanged except that For the range of possible perovskite layer thicknesses considered it always holds that δ ≪ λ ≪ 1 and this observation motivates the asymptotic solution to the model considered in the next section.
3 Asymptotic simplification of the model (δ ≪ λ ≪ 1) Here we assume dimensionless parameter sizes consistent with (15) and in particular require that δ ≪ λ ≪ 1. In this scenario, the problem for the anion vacancy density and potential (P and φ) decouples from that for the charge carrier densities (n and p) so that a very good estimate of φ can be obtained by ignoring the contributions of n and p in (11). We shall further assume that the cation vacancies are effectively immobile on the timescales of interest, reflected in the choice of ∆ = D − /D + = 0. This assumption, coupled to equations (11) and initial conditions (13), imply that the cation vacancy density remains constant with N ≡ 1.

Sym.
Description Hole pseudo-lifetime 9 × 10 −10 s, [28] Table 1: Parameters for the device described in 2.2, where ε 0 is the permittivity of free space. Here α is calculated from Loper et al. [19] based on light wavelength of 585nm (close to the peak absorption of the perovskite layer). Unless stated otherwise, the parameters are for the perovskite layer.

The ion problem
A good approximation to the potential can be obtained from the ion vacancy dependent equations in (11)-(13) at leading order, i.e. with Since λ ≪ 1, these equations can be further approximated by using asymptotic boundary layer theory, in a similar vein to Richardson et al. [27]. In the limit λ → 0, the solution can Debye Layer Figure 2: Schematic representation of the Debye layers and the solution for the electric potential, φ.
be subdivided into three regions consisting of a bulk (or outer) region which is separated from the two boundaries by boundary layers of width O(λ), see figure 2. As is usual in this type of problem, these boundary layers are termed either Debye layers or double layers (we opt for the former usage).
Bulk Region. Away from the boundaries (i.e. for x ≫ λ and 1 − x ≫ λ) the variables P , F p and φ can be expanded, in powers of λ and δ, as follows: Substituting these expansions into (16), and assuming δ/λ ≪ 1, gives, at leading order, Note that correction terms in the expansions of P and φ are O(δ) and O(δ/λ), respectively. These arise from the presence of the O(δ) charge carrier terms in Poisson's equation (last of (11)) and this is why the expansion breaks down if the value of either n or p becomes comparable to O(λ/δ). It follows that φ for arbitrary functions of time W − (t) and W + (t). It follows, on substituting into (19), that the leading order ion flux is given by The Debye layers. The asymptotic solution in the Debye layer about x = 0 is obtained by rescaling space in the governing equations (16)-(17) via and substituting the asymptotic expansions into the rescaled equations to obtain the leading order problem. The solution to which is given in Appendix A and can be summarised as follows: (A) the leading order potential, φ where W − is the potential at the left-hand side of the bulk, to which φ Similarly, the asymptotic solution in the Debye layer about x = 1 is obtained by rescaling space in the governing equations (16)-(17) using the transformation and substituting the asymptotic expansions into the resulting equations and solving at leading order. Once again this process is described in detail in Appendix A. As in the other Debye layer, the leading order potential φ (D) 0 (ξ, t) and vacancy distribution P where W + is the potential at the right-hand side of the bulk (to which φ (D) 0 matches as ξ → +∞), V + (t) is the potential drop across this Debye layer (see figure 2) and θ(ξ, V + (t)) is once again a solution to the problem (24).
In order to fully determine the leading order solutions in both Debye layers and the bulk region it is necessary to solve for the time-dependent functions V − , W − , V + and W + . The requirement that the leading order solutions in the Debye layers, (23) and (27), satisfy the potential boundary conditions in (17) gives Charge conservation within the Debye layers. A further two conditions on these four functions can be obtained by matching the flux of vacancies into the Debye layers with the leading order expansion of the vacancy conservation equations in the Debye layers. This leads to solvability conditions (described in Appendix A) which can be interpreted in terms of global conservation of charge within the Debye layers. Since the leading order solutions for the vacancy densities within the Debye layers are quasi-steady, the total (dimensionless) charge per unit area within each Debye layer, Q, can be related to the potential drop across the layer, V, in the form of a nonlinear capacitance relation. Here the charges per unit area contained within each Debye layer (Q − in that about x = 0 and Q + in that about x = 1) are defined, in terms of the local Debye layer variables ζ and ξ, by and, as shown in Appendix A, are related to the potential drops across the Debye layers (V − and V + , respectively) via the capacitance relations where the function Q(V) is defined by This relation is plotted in figure 3. Furthermore, since vacancies (and hence charge) are conserved, the rate of change of the total charge per unit area within the Debye layers must equal the flux of (positively charged) vacancies flowing into each layer from the bulk region.
Since the vacancy flux in the bulk region, F p,0 , is spatially independent, and given by (20), this observation corresponds to the conditions Alternatively, on eliminating W + and W − in favour of V + and V − , we have the equivalent conditions which can be solved in conjunction with (30)- (31). Adding these two equations together and integrating with respect to t implies that the total charge in the Debye layers is conserved, i.e. Q − (t) + Q + (t) is constant. This is to be expected given that the predominant mobile charge carriers are the positive vacancies which cannot leave the perovskite region. Furthermore, since the net charge arising from both positive and negative vacancies will initially always be zero, it remains so for all time, i.e.
At this stage we can choose either to solve an ODE for V + (t) or one for Q + (t). Since neither of these problems admit exact solution we opt to solve for Q + because this is preferable from a numerical point of view. We do this by noting that the inverse of (31) is where LambertW 0 (·) and LambertW −1 (·) are the 0'th and −1'st branch of the Lambert W function. On substituting the above functional relation in (32), together with (33), we obtain a single ODE for Q + (t): The solution to (34) may be used to obtain the leading order bulk potential via (19), that is Remark. The dimensional surface charge density (in the Debye layers), Q (dim) , is related to its non-dimensional counterpart, Q, by The uniformly valid approximation to φ. We can now write down a uniformly valid approximation to φ that is valid throughout the bulk and both Debye layers: where the function θ(z, V) is defined implicitly in (68) in Appendix A. The corresponding uniformly valid asymptotic approximation for the anion vacancy density, P , is

Asymptotic approximation to the charge carrier equations
As we demonstrate in §4, the potential is well-approximated by the solution to the ion problem (16)- (17) and is almost entirely unaffected by the carrier distributions. Furthermore, since the Debye layers are extremely thin, the effects of both photo-generation and recombination within these layers are negligible so that, from (11), the electron and hole currents are to a good approximation spatially independent across these layers, Furthermore, in these narrow regions, electron and hole densities are in approximate quasithermal equilibrium. In particular, in the Debye layers close to x = 0 and x = 1 respectively Referring to the boundary conditions (12), we find that For the purposes of predicting the output current of the device, we need only determine the carrier concentrations within the bulk region. Matching conditions on the bulk carrier problems (for n and p) are obtained from the far-field behaviour of the Debye layer solutions, namely The appropriate boundary conditions on the bulk carrier densities are thus The corresponding equations for the carrier densities in the bulk, as obtained from (11), and are, on taking the physically appropriate limit ν → 0, where E Hence the asymptotic approximation to the charge carrier problem can be found from the solution of (38)- (39) in which the electric field term, E  (34)). Usually the solution will have to be obtained numerically because of the nonlinearity of the recombination term. Nonetheless, numerically solving this reduced problem is considerably less challenging than directly tackling (11)-(13) because it excludes the Debye layers, over which the solution varies very rapidly. Finally, we note that the net current density j (o) (t) = j p (x, t) is independent of the spatial variable x and so can be found simply by evaluating the sum of the electron and hole current densities at any point in the domain.

An analytic solution in the limit ǫ → 0 with zero interfacial recombination
It is notable that the parameter ǫ = τ n /τ p is typically small (we estimate, on the basis of earlier work [28], ǫ ≈ 3.3 × 10 −3 ) while the other parameters in the SRH recombination term (14), N i and K 3 , are both very small. These observations lead us to set N i ≡ 0, K 3 ≡ 0 and to investigate the small ǫ limit. In which case, provided that p/n is not large, R(n, p) can be approximated by If we restrict our interest to the case where interfacial recombination is negligible (i.e. if we take R l ≡ 0 and R r ≡ 0), it follows that the equation for the hole density decouples from that for the electron density (see (39)) and can be reformulated as the following linear equation for p (o) : These may be solved by eliminating j This can be rewritten in the form where On noting that E (o) 0 (t) = β 1 (t) + β 2 (t), the boundary conditions (38) can be stated as The solution to (42) and (44) is where , , , An expression for the total current in the device. An expression for the hole current density j The total current J(t) = j n (x, t) is determined from the condition that j   (1, t). It follows that Asymptotic solution for the bulk electron density. In order to monitor whether this asymptotic solution breaks down, it is useful to derive an asymptotic expression for the bulk electron density, n (o) , while recalling that we require p (o) /n (o) ≫ ǫ in order for the validity of the expansion. The equations and boundary conditions for n (o) are, at leading order, in which we once again write E where time-dependent functions D, G and H are given by , , . (51)

Comparison between numerical and asymptotic solutions to the model
In this section, we compare the results obtained from (i) a numerical solution to the full model, (11)- (13), to those obtained from (ii) a combined asymptotic/numerical approach, in which the ion problem is solved asymptotically as in §3.1, and from (iii) the special case described in §3.3 which is entirely based on asymptotic approximations. In particular, we show that the results from (ii) the combined asymptotic/numerical approach, adopted in an earlier work [28], compare extremely favourably to (i) numerical solution of the model.

Numerical methods
In approach (i), we use the method of lines. A detailed description of the numerical scheme is given by Courtier et al. [8], here we restrict ourselves to a brief outline. The spatial derivatives in equations (11) are treated using a finite difference approach that is secondorder accurate in space, both on the internal and boundary points, and chosen in such a way that conservation of species is also exact up to second order. After application of the finite difference approximations, the problem is reduced to a system of differential algebraic equations (DAEs) in which the ODEs arise from the evolution equations for P , n, and p, in (11), and the algebraic equations are a result of Poisson's equation for the potential. Solving systems of DAEs presents a challenging numerical problem which we tackle using the ode15s routine in Matlab [1]. Owing to rapid changes of the solution curves within the narrow Debye layers, we find that the problem is sufficiently stiff to require non-uniform grid spacing and the additional precision offered by Advanpix's Multiprecision Computing Toolbox [2]. In approach (ii), the system of equations requiring numerical treatment is that for the charge carriers in the bulk, (38)- (39). Having taken the asymptotic limits δ, λ and ν → 0, the remaining problem is a second-order boundary value problem (BVP). Crucially, since asymptotic expressions have been derived for the narrow Debye layers, only the solution in the bulk needs to be resolved numerically. This problem exhibits significantly reduced stiffness and, as a result, a straightforward application of the bvp4c routine in Matlab [1] suffices.

Results
In figures 4-7, we show results for a device characterised by the parameters given in Table  1 with the perovskite layer width equal to 600 nm, corresponding to the set of dimensionless parameters given in (15). All numerical calculations are performed on a spatial grid consisting of 800 points. Figures 4-6 show the internal state of the cell at five equally-spaced values of time during a variation of the applied voltage, in a scenario in which the cell is abruptly illuminated at t = 0s having been preconditioned in the dark with V ap = V bi . For figure 4, the applied bias is varied smoothly from V ap = V bi at t=0 s and V ap = 0.8 V at t=10 s (precisely, V ap = V bi −0.2 tanh (t)/ tanh (10)). For figure 5, the applied bias is instantaneously decreased from V ap = V bi to V ap = 0V at t = 0s and held there for 4 seconds. Plots show solutions at t = 0.8, 1.6, 2.4, 3.2, 4.0 s. Finally, for figure 6, the applied bias is varied linearly from V ap = V bi at t=0 s to V ap =1.25V at t=10s. Plots are for t = 2, 4, 6, 8, 10 s. In figure 7, comparison is made between current-voltage (J-V) curves calculated using all three approaches and which model the experimental data presented by Richardson et al. [28]. Here, the cell is preconditioned for 5 seconds at 1.2V in the light before the J-V curve is measured. The current is calculated at equally-spaced intervals in time as the applied voltage is varied at a constant rate from 1.2V (forward bias) to 0V (short-circuit) and back; the four different scan rates are 50mVs −1 , 100mVs −1 , 250mVs −1 and 500mVs −1 . The colour scheme has been chosen for consistency with Fig. 7 (b) from Richardson et al. [28]. In panel (a), solutions calculated using (i) the fully numerical (solid lines) and (ii) the combined asymptotic/numerical approach (dashed lines) are shown. Note that both of these methods calculate currents based on the full SRH recombination rate, (7). While in panel (b), solutions from (iii) the fully asymptotic approach are shown.

Discussion
We have looked at three approaches to solving the drift-diffusion model (11)- (13). Approach (i) is fully numerical and involves solution of the full problem. In contrast, in approach (ii) (used previously [28]), we formally take the limits λ → 0 (small Debye length) and δ → 0 (charge carrier concentration negligible in comparison to ion vacancy concentration). The comparison between the results of these two approaches is extremely favourable, as illustrated by the very small discrepancies in the J-V curves calculated using both approaches, for a range of scan rates, in figure 7(a).
The other main approximations to the drift-diffusion model that we make use of are the quasi-steady carrier limit, ν → 0 and the approximation of SRH recombination by hole dominated monomolecular recombination, ǫ → 0. The former limit (ν → 0) and its use, or otherwise, makes negligible difference to the results obtained. The latter, however, is frequently problematic, despite the very small value of ǫ (= 1/300) we use in the simulations. This slightly surprising result is best illustrated by the significant differences between J-V curves calculated using the fully numerical method (solid curves in figure 7(a)) and those calculated using the fully asymptotic method in the limit ǫ → 0 ( figure 7(b)). Where there are significant differences between the two approaches this can be ascribed to strong spatial variations in charge carrier concentrations across the cell, resulting in regions where n ≤ O(ǫ p) so that the approximation of R(n, p) in (14) by R(n, p) ≈ γp no longer holds.

Conclusion
In this work we outlined a model for charge carrier transport and ion vacancy motion in a tri-layer planar perovskite solar cell (previously discussed in Richardson et al. [28]). Using parameters extracted from the literature, we were able to identify two key small dimensionless parameters that characterise the model: λ, which gives the ratio of the Debye length in the perovskite to the width of perovskite layer, and δ, the ratio of the typical charge carrier (electron and hole) densities to the typical ion vacancy density. Based on the small size of these parameters, we performed an asymptotic analysis of the model which showed that: (a) the problem for the ion vacancy density and the electric potential is almost completely independent of the charge carrier densities and (b) the decoupled problem for ion vacancies and electric potential is well-approximated by the solution to a single first order ODE that describes the evolution of charge in the Debye layers (at the edge of the perovskite) in terms of the current through a resistor and a nonlinear capacitor in series. In dimensional form, this simplified model states that the charge (per unit area) in the right-hand Debye layer, Q + , evolves according to the equation where the term in the brackets is the (uniform) electric field in the perovskite bulk (away from the Debye layers) and V(Q) is the inverse to the nonlinear capacitance relation A good approximation to the full model can then be obtained by solving this much simplified problem for ion vacancy density and electric potential and using the resulting electric potential as an input into the charge carrier equations. The resulting model can sensibly be termed a surface polarization model of charge transport because it describes the effect on the current in a cell of the polarization of the perovskite layer, as ionic charge is transported from one of its surfaces to the other. In general, the simplified problem that we are left to solve for the charge carrier densities is nonlinear and so requires numerical solution. However, in contrast to the problem for the ion vacancies and potential, it is non-stiff and so this is not usually problematic. Moreover, parameter estimates suggest that the Shockley-Read-Hall recombination term in the charge carrier equations can be well approximated by monomolecular hole dominated recombination (R(n, p) ≈ γp). This allows the charge carrier equations to be linearised and, in turn, solved analytically. Where this is the case, an asymptotic solution to the entire model can be obtained from the solution to the single first order ODE discussed above. In order to test the validity of the asymptotic method used to solve this model, we compared our asymptotic results to the results of a numerical solution of the full model. The latter was conducted using a recently developed numerical procedure [8] that is able to accurately solve the full model in realistic parameter regimes. Where we used a combined asymptotic/numerical approach (solving for the ion vacancy distribution and electrical potential using the asymptotic model and solving for the charge carrier densities and currents numerically), we found extremely good agreement to the full numerical solution. In the case where we additionally linearised the charge carrier equations and solved them analytically, the comparison to the full numerical solution, while still good, was less impressive.
The physics of perovskite solar cells is still far from fully understood and in order to improve this situation it is vital that drift-diffusion models and their solution techniques continue to be developed. One obvious, and important, extension to the model discussed here is the explicit inclusion of charge transport in the electron-and hole-transport layers on either side of the perovskite. Such an extension will be able to elucidate how the choice of these layers affects the cell's transient behaviours. In particular, this extended model could be used to investigate cell architectures giving rise to so-called low hysteresis behaviour and would also be better able to account for interfacial recombination, see for example [6,39]. Here we assume cation vacancies are immobile, which is justified by the relatively short timescales. However, it is believed that mobile methylammonium vacancies can lead to slow (over the timescale of many hours) but reversible changes in efficiency [10].

Author Contributions
NEC and JMF coded the numerical solver and produced the plots. GR performed the asymptotic analysis. GR and NEC wrote the manuscript. ABW, SEJO'K and GR conceived the project and formulated the model.

A Solution for P and φ in the Debye layers
In the bulk region, we obtain a solution for the leading order vacancy density P = W − (t)(1 − x) + W + (t)x. These expressions satisfy the potential boundary conditions but in general cannot satisfy the flux boundary conditions, see (17). In order to resolve this seeming paradox, we need to account for narrow boundary layers (Debye layers) of width O(λ) about x = 0 and x = 1.
Debye layer about x = 0. Considering first the Debye layer about x = 0, we use the rescaling (21) to rewrite the governing equations (16)- (17) in terms of the rescaled spatial variable ζ, yielding the boundary layer equations: The expansions for P , φ and F p proceed as in (22) so that to leading order in (54) we obtain the following equation for P This has the solution for some as yet undetermined function of time, W (t).
Matching to the outer. In order for the leading order Debye layer solution to match to the leading order outer solution, through (18)- (19), we require Applying the matching condition (58) to the solution (57) determines a relation between the arbitrary functions W (t) = W − (t) motivating us to eliminate one of them by writing On substituting this expression into (55) balanced at leading order, we find which satisfies boundary conditions obtained from the leading order terms in (56) and from (58), namely The corresponding expansion for the total charge per unit area in the left-hand Debye layer, Q − (defined in (29)) is and, by substituting this into (29), we obtain We can reformulate the problem for φ (d) 0 , given by (60)-(61), in a generic form by writing φ where V − (t) represents the potential gained across the Debye layer, i.e.
). It is then straightforward to show that the function θ(z, V) must satisfy the generic modified Poisson-Boltzmann problem Furthermore, in order that φ (d) Thus if we are able to determine V − (t), we can determine the unknown function W − (t) in the leading order outer solution for the potential in (19). It is straightforward to obtain a first integral to the autonomous equation (64) in the standard fashion by multiplying by θ z and integrating with respect to z. This yields, on applying the far-field condition (65), the expression where the −sign(V) is to account for the fact that if V < 0 (V > 0) the gradient of θ must be negative (positive). We can integrate (67) once more to obtain a relation for z as a function of θ which reads z = Z 1 (θ) − Z 1 (−V) for V > 0 with θ < 0, z = Z 2 (θ) − Z 2 (−V) for V < 0 with θ > 0, where Z 1 (θ) = 1 Characterising the capacitance of the Debye layer. We now seek to relate Q −,0 , as given in (62), to V − . We note that Q −,0 = ∞ 0 e −θ − 1 dζ, which we can rewrite as On substituting for θ ζ from (67) and writing θ = −V this integral transforms to This relation is plotted in figure 3, from which it can be seen that Q −,0 is a single valued function of V − . Hence, given the Debye layer charge density, Q −,0 , we can invert to find the potential jump across the Debye layer, V − . This motivates us to consider the evolution of Q −,0 (t) as charge (in the form of positively charged vacancies) flow out into (or in from) the bulk region.
A solvability condition on Q −,0 (t). It remains to determine the evolution of V − (t). This can be done by tracking the charge build up in the Debye layer through the leading order expansion of the positively charged vacancy conservation equation (54), and the boundary conditions p,0 | x=0 as ζ → +∞.
These conditions are obtained from the leading order expansion of (56b) and from matching to the leading order outer solution as ζ → +∞, respectively. By writing ∂P The Debye layer about x = 1. The analysis of this right-hand layer proceeds in a similar fashion to the left-hand Debye layer presented above. We introduce the rescaled spatial variable ξ, defined in (25), and then expand as follows.
Once again, a solvability condition may be derived from the problem for the leading order anion vacancy density, namely p,0 | x=1 as ξ → +∞.
The solvability condition we obtain on integrating this system is the following evolution equation for Q +,0 (t):