IonMonger: a free and fast planar perovskite solar cell simulator with coupled ion vacancy and charge carrier dynamics

Details of an open-source planar perovskite solar cell simulator, which includes ion vacancy migration within the perovskite layer coupled to charge carrier transport throughout the perovskite and adjoining transport layers in one dimension, are presented. The model equations are discretised in space using a finite element scheme, and temporal integration of the resulting system of differential algebraic equations is carried out in MATLAB. The user is free to modify device parameters, as well as the incident illumination and applied voltage. Time-varying voltage and/or illumination protocols can be specified, e.g. to simulate current–voltage sweeps, or to track the open-circuit conditions as the illumination is varied. Typical simulations, e.g. current–voltage sweeps, only require computation times of seconds to minutes on a modern personal computer. An example set of hysteretic current–voltage curves is presented.


Introduction
Perovskite solar cells (PSCs) are a promising thin-film technology that, due to their rapid rise in power conversion efficiency to 22.7% [15], are attracting a lot of interest and research effort in the photovoltaic community. However, PSCs display unusual transient behaviour (exemplified by current-voltage hysteresis) in the order of seconds to days which affects the power output of the device [24]. The consensus in the literature is that this slow (compared to the timescale of electronic motion) behaviour is due to the motion of mobile ion vacancies within the perovskite layer. The species of ion vacancy deemed most likely to be responsible for the behaviour observed in the order of seconds to minutes is that of the halide (e.g. iodide, I − ) ion vacancy due to its high mean density and high diffusion coefficient (compared to the other ionic species) predicted by DFT calculations [8,26]. Visual evidence of iodide ion migration within a perovskite film has also been obtained experimentally [7]. Recent reviews of the outstanding challenges in the field of perovskite solar cells have been given by Snaith [23], Egger et al. [9] and Phung and Abate [17].
Due to the existence of mobile ion vacancies, the perovskite layer must be treated as a mixed ionic-electronic conductor for the purpose of device simulation. The first works [2,11,25] to apply numerical methods to PSC modelling reported that their simulations suffered from prohibitively long calculation times and inaccuracies in solution for realistic values of the parameters. A combined analytic/numerical method was used by Richardson et al. [5,18] to reveal how iodide ion vacancies accumulate/deplete in very narrow layers (called Debye layers) adjacent to the perovskite boundaries. The associated rapid change in solution across these Debye layers is a significant source of spatial stiffness, while the disparity in timescales between the motion of electronic and ionic charges is a source of temporal stiffness [6], rendering the task of finding solutions to realistic models of PSCs very challenging.
Courtier et al. [6] developed a finite element scheme, implemented in MATLAB [14], that is able to cope with the spatial and temporal stiffness of the problem and obtain Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s1082 5-019-01396 -2) contains supplementary material, which is available to authorised users. accurate solutions to a coupled model for ion migration and charge carrier transport within the perovskite layer of a PSC. Since then, alternative numerical methods have also appeared in [13,27]. Walter et al. [27] have developed a solver that includes the motion of both cations and anions. Their scheme is based on the freely available software Quokka3 [10], originally designed for the solution of models of silicon-based photovoltaic devices. While this provides a thoroughly tested and validated framework for their results, the model being solved only explicitly models the perovskite absorber and there is strong evidence suggesting that the adjacent transport layers play a key role in determining device behaviour [4]. Meanwhile, Jacobs et al. [13] use the proprietary COMSOL package via a MATLAB interface to simulate three layers of a PSC over a range of timescales. However, the details of the modelling and solutions techniques are not given in full, making their results difficult to reproduce and compare with alternatives.
In this work, we present the extension of the finite element code presented in [6] for a single-layer model of a PSC to a model that explicitly describes the three core layers of a PSC: the electron transport layer (ETL), perovskite absorber layer and hole transport layer (HTL). Full details of the charge transport model equations and the implementation of the code are given. The code is freely available on GitHub at https ://githu b.com/Perov skite SCMod ellin g/ IonMo nger (under an AGPL-3.0 copyleft license) and can be used to simulate a variety of different experimental protocols. Current density-voltage ( J-V ) curves are the typical measurement that is performed to assess solar cell performance including, for perovskite solar cells, the extent of J -V hysteresis displayed at a particular scan rate. In addition to J-V curves, the code can be used to simulate photocurrent transients (during which the applied voltage and/or the illumination intensity is varied over time) and photovoltage transients (during which the cell is held at open circuit and the illumination intensity is varied) that occur on timescales of microseconds to minutes. Uncovering the links between model parameters and the results of such simulations will help to improve understanding of the underlying physics of PSCs and hence guide further improvements in their design. One such investigation, into how the extent of observable hysteresis depends on material properties of the two transport layers, has been conducted by Courtier et al. [4]. The investigation is based upon simulations of the same 1 threelayer model as that considered in this work.
In Sect. 2, we show an example set of simulation results obtained using the code. In Sect. 3, we present and discuss the governing equations in each of the three layers as well as the boundary and interface conditions through which the equations couple together. A non-dimensionalisation is also presented which is geared to study the behaviour of the cell on the timescales associated with anion vacancy motion and aids in obtaining uniformly well-resolved solutions by ensuring that the numerical tolerances are applied equally to each of the model variables. In Sect. 4, we detail the finite element discretisation of the system and highlight the differences between how open-circuit and applied voltage protocols are imposed at the discrete level. The focus of Sect. 5 is a discussion of the how the time-stepping is carried out and how the output current density is calculated from the numerical solution. In Sect. 6, we validate the results of the numerical scheme upon which IonMonger is based. Finally, in Sect. 7, we draw our conclusions.

Application example
The main purpose of this paper is to provide the perovskite solar cell community with a high-quality, free and useful tool with which to better understand PSC behaviour, and not to provide a detailed analysis or interpretation of the device physics. The power of the solver in advancing our physical understanding of PSCs is demonstrated in Courtier et al. [4], in which an investigation of the effects of material properties of the transport layers on cell performance is detailed. There it is found that two material properties in particular, namely the permittivity and the effective doping density of the transport layers, have a significant role to play in determining the extent of J-V hysteresis exhibited by a PSC. In addition, characteristics of simulations that can be used to identify the dominant recombination mechanism in a PSC are discussed. Results from two other simulations, also computed using the capabilities offered by IonMonger, have been presented by Idígoras et al. [12] in a study of the role of surface (also often termed interfacial or interface) recombination, which occurs on the interfaces between the perovskite layer and the adjacent transport layers, on PSC performance. However, further work in this area is vital for the future development and optimisation of PSCs. Here, we show how IonMonger can be used to both simulate the most common measurement protocol for assessing the performance of a solar cell, namely a J-V curve, and reveal how each type of recombination included in the model contributes to the observed behaviour. Figure 1a shows the J-V output for an example simulation of a typical J-V measurement protocol performed on a planar, standard-architecture PSC (see Fig. 2 for the cell geometry). The cell is initially preconditioned for 5 s 1 3 before the applied voltage is scanned from 1.2 V to short circuit (0 V) and back to 1.2 V at a scan rate of 100 mV/s. In addition, Fig. 1a displays the corresponding current losses due to the two types of interface recombination included in the model. Note that bulk recombination occurring within the perovskite layer is also included in the model but is not shown. By comparing the shape of the current-loss curves to the J-V curve, it is clear that, for this example, the observed behaviour is controlled primarily by the rate of recombination at the ETL/perovskite interface (the blue line), while the rate of recombination at the perovskite/HTL interface (the red line) has little effect on the performance of the cell. The parameter values used in the simulation are equal to those given in Tables 1 and 2(b) of [4] except that here the effective doping densities d E = d H = 10 24 m −3 and the effective densities of state g E c = g H v = 5 × 10 24 m −3 . The input file for this simulation, along with a GUIDE and documentation to aid in using the code, is provided in the main folder of the IonMonger GitHub repository so that users can utilise these as a starting point for investigations of their own.
In Fig. 1b, the example J-V curve from panel (a) is shown alongside the corresponding results for two other simulations. The only difference in the input parameters between the three simulations is the rate at which the applied voltage is scanned back and forth to measure the J-V curve. The three scan rates are 50, 100 and 200 mV/s. Harvesting the full set of results in Fig. 1b

The charge transport model
In this section, the charge transport model for a planar lead halide perovskite solar cell consisting of a perovskite absorber layer sandwiched between an electron transport layer (ETL) and a hole transport layer (HTL) is presented. Tables of the model variables and parameters along with their definitions are given in the SI. The structure of the cell is shown in Fig. 2. The non-dimensionalisation used by the code is also given.
The purple lines show the current-density output, while the blue and red lines show the current losses due to interface recombination at the ETL/perovskite and perovskite/HTL interfaces, respectively. Losses due to bulk recombination are not shown. The direction of scan is indicated by both the arrows and the style of each line: solid for the reverse scan and dashed for the subsequent forward scan. b A set of three J-V curves. Here, the example in panel (a) is plotted alongside two other J-V curves measured at scan rates of 50 mV/s and 200 mV/s, but for otherwise identical input parameters. These results demonstrate the ability of the model to reproduce the scan rate-dependent J-V hysteresis commonly exhibited by PSCs due to the migration of ion vacancies within the perovskite layer (Color figure online)

Model equations
Perovskite absorber layer (0 < x < b ) A model for the perovskite layer consists of equations for the conservation of conduction (free) electrons, holes and halide ion vacancies coupled with the Poisson equation for the electric potential, (x, t) . We denote the halide ion vacancy density by P(x, t), its flux by F P (x, t) and define N 0 as the mean ion vacancy density. It is assumed that an equal, uniform density N 0 of immobile cation vacancies also exists within the perovskite layer. The electron and hole concentrations are denoted by n and p with current densities j n and j p , respectively. The functions G(x, t) and R(n, p) denote the charge carrier photogeneration and bulk recombination rates, respectively. In the perovskite layer, we thus have with Poisson's equation, Here, D I denotes the diffusion coefficient of the iodide ion vacancies and A is the permittivity of the perovskite absorber layer. These differential equations are supplemented by continuity conditions at the interfaces with the transport layers (given at the end of this section).

Electron transport layer
The majority carriers through the ETL are free electrons. The model for the electrical behaviour of this layer thus consists only of a conservation equation for the free electrons which couples to Poisson's equation. Here, D E denotes the electron diffusion coefficient, E the permittivity and d E the intrinsic free electron density (due to the doping) in the ETL. (1) These equations couple to the perovskite equations via four continuity conditions at the interface (given at the end of this section). On the external boundary with the metal contact, we impose Ohmic boundary conditions. These read where V(t) is the applied voltage and V bi denotes the cell's built-in voltage, which is defined in (17).
The majority carriers in the HTL are holes, and analogously to the case in the ETL, we need specify only two equations, specifically Here, D H is the hole diffusion coefficient, H is the permittivity and d H is the intrinsic hole density (due to the doping) in the HTL. These equations couple to the equations in the perovskite via four continuity conditions at the interface (given at the end of this section) and satisfy the following Ohmic contact conditions on the metal contact.
Continuity conditions on the interfaces (x = 0 and x = b) At the interface between the perovskite and the ETL, (1) the electron flux (and its associated current density) is conserved, (2) the hole flux (and its associated current density) is conserved, (3) there is no flux of halide ion vacancies, (4) and (5) both the electrostatic potential and electric displacement field are continuous, and (6) the majority carrier density (in this case the electrons) at the edge of the ETL is related to the neighbouring carrier density in the perovskite by a factor, k E , which depends upon the relevant band offset and change in effective density of states [11]. Therefore, at the interface between the ETL and the perovskite, the following conditions are applied Analogous conditions are applied at the interface between the perovskite and the HTL (where the holes are the majority carrier). These read Here, the superscripts of ± denote quantities evaluated at either the left-or right-hand side of the perovskite/transport layer interfaces, respectively; R l and R r are the recombination fluxes at the left (ETL/perovskite) and right (perovskite/ HTL) interface, respectively; and k E and k H are constants of proportionality between the charge carrier concentrations on each side of the interfaces according to Boltzmann statistics, given by In these expressions, g c,v denote the effective conduction/ valence band density of states; E c,v are the energies of the conduction/valence band edges; and the superscripts E or H indicate to which transport layer a quantity relates. The validity of each of these expressions relies on the validity of using the Boltzmann approximation to describe the distribution of electrons in a non-degenerate semiconductor. Consequently, users should choose an effective doping density which is less than 20 times smaller than the effective density of states in each transport layer in order to ensure that the equilibrium Fermi level is more than a few thermal voltages away from the band edges, see (15)- (16).

Built-in voltage
The cell's built-in voltage is equal to the difference between the workfunctions of the two metal contacts. Assuming that the contacts form ideal Ohmic contacts with the adjacent transport layer, this difference is equal to the difference between the equilibrium Fermi levels of the two transport layers which are approximated, using Boltzmann statistics, by (12) Hence, the built-in voltage is given by

Carrier generation and recombination rates
For the rate of charge carrier generation within the perovskite layer, we use a simplified Beer-Lambert model of light absorption [16] in which it is assumed that absorption can be characterised by a single absorption coefficient ( ) and photon flux that are independent of the wavelength of light. Taking into account the possibility of choosing either a standard or inverted architecture (i.e. applying the light to either the ETL-or HTL-side of the cell, respectively), this rate can be written as in which F ph denotes the flux of photons incident on the light-facing perovskite surface (after accounting for reflection) under the equivalent of 1 Sun illumination; the function I s (t) is the intensity of the illumination in Sun equivalents; and the parameter l can be set equal to either +1 for light from the left (through the ETL) or −1 for light from the right (through the HTL) by making use of the Inverted option in the parameter input file. We allow bulk recombination to be described by a combination of bimolecular and trap-assisted Shockley-Read-Hall (SRH) recombination mechanisms from, for example, Nelson [16] §4.5.5. Hence, the volumetric bulk recombination rate is in which is the bimolecular rate constant, n i is the intrinsic carrier concentration, n and p are the charge carrier lifetimes and we assume that the trap state energy level lies close to the intrinsic potential energy of the perovskite such that we can apply the approximation that p t = n t = n i .
Similarly, the interfacial recombination fluxes ( R l and R r ) can be chosen as a combination of bimolecular and SRH recombination as follows, noting the use of recombination velocities rather than carrier lifetimes in the SRH recombination rates (see Nelson [16,Sect. 4.5.6]), for consistency with the continuity conditions in (11f) and (12f). Then, in order to keep the number of input parameters to a minimum, we approximate p t = n t = n i in analogy with the approximation made to the bulk recombination rate above. Therefore, the interface recombination rates used by the code are equivalent to

Calculation of the total current density
In order to calculate the total current density from a numerical solution of the drift-diffusion model, we derive an expression that can be evaluated at any point in the domain. The code automatically calculates the current density at the midpoint of the perovskite layer, where the grid spacing is larger and the solution varies more smoothly than in the Debye layers, to minimise numerical error. By subtracting Eq. (2a) from Eq. (1a), we get Then, by substituting the difference in the carrier concentrations ( p − n ) using Poisson's equation for the perovskite layer, given in (4), and multiplying by q, we get (20) Applying the time derivative to the last term in the brackets allows us to eliminate N 0 (the constant and uniform cation vacancy density) and to use the ion vacancy conservation equation in (3) to make a substitution for the time derivative of P, which gives Similarly, for the transport layers, we have After swapping the order of spatial and temporal differentiation, it is possible to integrate these three equations with respect to the spatial variable x. By integrating and applying the continuity conditions and the Ohmic boundary conditions at either metal contact, we get an expression for the total current density which is independent of x and given by The term involving a time derivative in each equation is the displacement current density. It should be noted that the magnitude of the displacement current density is usually even smaller than the magnitude of the expected numerical error in the solution, see Sect. 6.

Boundary conditions for open-circuit conditions
It is possible to simulate open-circuit conditions, rather than a fixed voltage protocol, by changing just two of the model equations. The two conditions that describe the application of a fixed potential difference are (7b) and (10b). The model for simulation of a device at open circuit instead includes boundary conditions to ensure that there is zero flux of electrons across the metal/ETL boundary (and hence no photocurrent) and that the values of the electric potential at each contact are equal and opposite. This amounts to imposing, in place of (7b) and (10b),

Non-dimensionalisation
The code is programmed to solve a non-dimensional form of the model equations. Note, however, that all input parameter values are non-dimensionalised automatically and the output is re-dimensionalised prior to output. Details of the nondimensionalisation are given here to allow the possibility that users can adapt the equations that underlie the model described here. The non-dimensionalisation is given by: where G 0 is a typical value of G (the rate of photogeneration of charge pairs per unit volume) and ion is the characteristic timescale for ion motion into the Debye layers, given by The perovskite (ionic) Debye length is defined as The other input functions and constants are rescaled as follows: The star notation is dropped in the following sections.

Discretisation
The numerical method upon which our code is based was developed by Courtier et al. [6] to solve a simplified model description of a perovskite solar cell in which it was assumed that the transport layers were so highly doped that the potential within them was uniform, thereby reducing the model to equations in the perovskite absorber layer only. In that work, the speed and accuracy of the method are compared against those of two previously used alternatives. It is shown that the method we adopt here is superior to both of these methods for both metrics of performance. Here, we adapt the finite element-based scheme to solve the dimensionless three-layer model set out in the previous section. We do not use the Scharfetter-Gummel scheme [19], commonly used for the solution of drift-diffusion equations, because it is tailored to deal with issues related to charge carrier transport rather than accurately resolving solutions in narrow Debye layers, which is the main difficulty in the present work.

Spatial grid
The discretisation is formulated on a computational grid comprised of N + N E + N H + 1 non-uniformly positioned grid points which partition the (non-dimensional) domain The perovskite layer (including interfaces) contains N + 1 grid points denoted by x = x i for i = 0, … , N with subinterval widths denoted by i+1∕2 = x i+1 − x i . The transport layer domains (excluding interfaces) contain N E and N H grid points, respectively, with grid points denoted by x = x E i for i = 0, … , N E − 1 and x = x H i for i = 1, … , N H , respectively, with corresponding subinterval notation.
It is known that the largest gradients in the solution appear in narrow Debye layers adjacent to the material interfaces [5]. This motivates the use of a grid in which points are concentrated near the domain boundaries and at internal interfaces so that computational resolution is focused there and not wasted where it is not required. One such grid can be created by extending the tanh grid used in [6] to the threelayer cell geometry displayed in Fig. 2. Specifically, we set (34) Our numerical experiments indicate that a good rule of thumb for deciding on a judicious choice for the value for s can be to calculate a value which leads to 20% of the grid points falling within one Debye length of each interface within the perovskite layer (i.e. the intervals x ∈ [0, ] and x ∈ [1 − , 1] ), via numerical solution of The code is set up so that the numbers of grid points N E and N H are chosen based on the input parameter N to give approximately equal spacing immediately either side of the interface.

Finite element scheme
As in [6], we employ a common finite element approach, in which each of the dependent variables is approximated as a linear combination of piecewise linear basis functions with compact support. For a generic dependent variable, u , defined within the perovskite, i.e. for x ∈ (0, 1) , we write where in which i (x) are referred to as the basis functions. Each of the governing equations of interest can be manipulated into the form in which A and B are constants and the function S(x, t, u, v 1 , v 2 , v 3 ) is a source term which depends upon the spatial variable x, the temporal variable t, the generic variable u and a series of other generic dependent variables v i for i = 1, 2, 3 . The electron, hole and halide ion vacancy conservation equations, (1), (2) and (3), are rewritten in this form by eliminating the attendant electron or hole current densities, or the halide ion vacancy flux, respectively. Poisson's equation in the perovskite, (4), is already in this form. The spatially discretised equations in the perovskite are derived by multiplying (40) through each of the test functions j (x) (for j = 0, … , N ), integrating over the domain x ∈ (0, 1) (using integration by parts where appropriate) and substituting form (38) for each of the dependent variables. On doing so, we arrive at Each of the integrals containing expressions that depend solely on the basis functions and/or their derivatives can be computed exactly. Likewise, terms containing quantities evaluated on the boundaries x = 0, 1 can be computed exactly using the continuity conditions (11)- (12), else the relevant equation is replaced by the corresponding Dirichlet condition. The one remaining term that is not so readily computed is the final integral in (41) that depends on the source terms S. For the anion vacancy conservation, S ≡ 0 and so this term is zero. For Poisson's equation, this term is a linear combination of dependent variables and so can be computed exactly. However, for the electron and hole conservation equations, (2) and (1), the source term comprises both the generation and bulk recombination rates, G(x) and R(n, p), which are highly nonlinear, see Sect. 3.2. In order that the integral contained in the final term of (41) can be integrated (at least approximately) regardless of the functional form of the source term, we make a further approximation, that is, to replace the dependent variables in the integrand by functions that are piecewise constant over each subinterval, x ∈ (x i , x i+1 ) , and have a value equal to that of the full series (38) at the midpoint of that interval. In short, we make the additional approximation The additional error incurred as a result of this approximation is comparable to the error associated with the original piecewise linear approximation for the dependent variables. Thus, even though some additional error is introduced, the scheme retains its second-order local accuracy, as demonstrated in Sect. 6. We note that this approach to deal with the nonlinear source terms is a special case of the method used in the work of Skeel and Berzins [22], but we emphasise that in contrast to their method, we only use this additional approximation for treatment of the source terms while the rest of the terms are integrated exactly.
An analogous methodology is used to derive the discrete equations in the transport layers. The only difference is that the basis and test functions are piecewise linear functions with compact support defined within the ETL and HTL, respectively.
For notational convenience, we introduce three discrete operators: a difference operator, i ; an interpolation operator, ℑ i ; and, a linear operator i . For a generic dependent variable u , these three operators are defined as follows, in which the midpoint These discrete operators can be used to obtain discretised versions of the ion flux and carrier current densities in the (42) perovskite layer, given in dimensional form in (3b), (2b) and (1b), as follows: The electron current density in the ETL ( j n,E ) and hole current density in the HTL ( j p,H ) can be expressed in an equivalent way. The carrier generation and bulk recombination terms are approximated to be linear on each interval, and to take the value at the midpoint, hence we define Perovskite absorber layer The discretised equations governing the evolution of the halide ion vacancy density subject to no-flux boundary conditions, corresponding to (3), (11c) and (12c), are The discretisation of Poisson's equation, from (4), becomes The conservation equations for the electrons and holes and the carrier current density boundary conditions, corresponding to (1a), (2a), (11b) and (12b), become

Electron transport layer
The equations and Dirichlet boundary conditions for the electric potential and the electron density in the ETL, from (5a), (6) and (7), become Hole transport layer Similarly, for the electric potential and hole density in the HTL from (8a), (9) and (10), 1∕2 1 3 dp 0 dt + 1 6 dp 1 dt

Continuity conditions
The carrier relations and continuity of the potential across the interfaces from (11d, f) and (12d, f) are applied directly as follows: The continuity of the displacement field and that of the electric potential across the interfaces, from (11e) and (12e), are ensured via Interface recombination is included in the equations for the continuity of carrier current densities, corresponding to (11a) and (12a), as follows: (68) (70)

Boundary conditions for modelling open-circuit conditions
In order to simulate open-circuit conditions using the alternative boundary conditions given in (29), we simply replace (59) and (64) by

Implementation
In this section, we outline the steps performed by the code. We begin by describing the procedure that is used to integrate forward in time. Next, we outline how the parameters, operating protocol and initial conditions are set. Finally, we outline how IonMonger post-processes the results so that quantities of interest, e.g. the dimensional current output, can be extracted and visualised.

Integration in time using MATLAB's ode15s
The system of differential algebraic equations formulated in Sect. 4.2 is evolved forward in time using MATLAB's ode15s [20,21]. A prerequisite for leveraging this algorithm is assembling the state variables into a column vector, (t) . A significant decrease in computational cost (proportional to the length of (t) squared) is available if the size of (t) can be reduced, and so we eliminate superfluous variables, namely F P j n , j p and E between equations (46)-(73) before assembling the 4N + 2N E + 2N H + 4 remaining unknown functions of time, into the column vector (t) as follows: where a superscript T denotes a transpose. In (75), , , and are column vectors of length N + 1 ; E and E are column vectors of length N E ; and H and H are column vectors of length N H . The problem to be solved can now be written in the form in which ( ) is a nonlinear vector function of length 4N + 2N E + 2N H + 4 whose entries are the right-hand sides of (51)-(71) and is a singular diagonal mass matrix whose entries are the coefficients of the time derivative terms in the same equations. Another useful strategy for speeding up computations, and one that we make use of in IonMonger, is to exploit ode15s's jpattern option. This facilitates additional savings in computational cost by specifying entries in the Jacobian of the vector function which are known to always equal zero a priori, thereby preventing the algorithm from having to numerically approximate their value as the integration in time proceeds. The function Jac creates a sparse matrix that indicates which entries of the Jacobian need to be numerically approximated at each time step and which are always equal to zero.

Parameter input and initial conditions
The necessary dimensional parameters, the illumination protocol G(x, t), voltage protocol V(t) and solver options are passed between functions in a MATLAB structure called params. A params structure can most easily be created using the script called parameters.m.
Finding initial conditions which satisfy the requisite boundary conditions is non-trivial. In the code, we opt to supply initial conditions that correspond to a device which has been left to reach a quasi-equilibrium with the applied voltage held equal to either a fixed value or the open-circuit voltage such that there is no output current from the cell. In either case, the task of finding initial conditions amounts to finding a valid steady-state solution to the PSC model and we start by finding initial conditions corresponding to when the applied voltage is set equal to the built-in voltage, as defined in (17). This task is tackled by the function initial_conditions.m. In order to start the simulation protocol from a different value of the applied voltage or from open-circuit conditions, an additional call is then made to either precondition.m or find_Voc.m, respectively. All three of these routines find discrete representations of the (dimensionless) initial conditions which can be written as where the initial profiles P (x) , ̂( x) , n(x) and p(x) satisfy (up to numerical tolerances) the discrete counterpart of the steady-state PDEs. This is achieved by invoking MATLAB's root-finding tools (i.e. fsolve) which act on the nonlinear system ( ) = , defined in (76), with one minor alteration. An additional integral constraint, namely ∫ then save the solution variables ( , , , , E , E , H and H ), spatial vectors ( , xE and xH), time (time), evolution of the applied voltage (V) and evolution of the current densities defined above ( , l and r ). These data are saved into a .mat data file along with the input structure (params).
The saved data can be further analysed and plotted in MATLAB in any way chosen by the user. One example plotting script is included in the IonMonger GitHub repository, and instructions for its use are given in the GUIDE. This script can be used to plot the current density generated by a PSC during the reverse and forward scans of a J-V curve, alongside the current losses due to recombination at each of the perovskite/transport layer interfaces, as shown in Fig. 1a. Such a plot can enable the user to identify the limiting recombination mechanism for a particular set of input parameters as described in Sect. 2.

Validation
For the purpose of verifying the results generated by IonMonger, in Fig. 3a, we plot a measure of the error in eight different solution variables against the number of grid points ( N + N E + N H + 1 ) on a logarithmic scale for five example simulations. The only input parameter that is varied between the five simulations is N, which takes values of 100, 200, 400, 800 and 1600, while all other input parameters are the same as those used for Fig. 1. Due to the lack of an exact solution, the errors are calculated with respect to another simulation performed on an even finer spatial grid consisting of 5613 points (for which the input parameter N is set equal to 3200). The chosen error measure is the sum (an l 1 norm) of the differences between the value of the variable computed by the example simulation and that computed by the 5613-point simulation, averaged over all 300 time points of the simulation protocol after t = 0 . The same error measure was used in [6] to compare the same solution method applied to a single-layer version of the model against two other methods on different spatial grids. Here, the eight solution variables for which the error is calculated are the (dimensionless) ion vacancy density, electric potential, electron concentration and hole concentration at the ETL/perovskite and perovskite/HTL interfaces located at x = 0 and x = 1 , respectively, as listed in the legend. The results demonstrate the expected secondorder pointwise convergence of the finite element scheme on which IonMonger is based [6]. The variation in the magnitude of the error between the eight solution variables is due to differences in the magnitude of the dimensionless variables themselves. Figure 3b shows the computation times associated with each of the five simulations in panel (a), also plotted against the number of grid points on a logarithmic scale. Note that the computation time will also depend upon the length of the simulation protocol. The temporal tolerances for the integration in time performed by MATLAB's ode15s were fixed for all simulations at values of 10 −6 for the relative tolerance and 10 −10 for the absolute tolerance.
A comprehensive verification of the single-layer finite element scheme, from which this code was developed, is provided in Sections 5 and 6 of [6]. This work includes plots of each solution variable across the perovskite layer against corresponding asymptotic results, which show very good agreement between the two approaches for realistic values of the input parameters.

Comparison to asymptotic results
Further validation of the numerical scheme, against results obtained using an alternative (although approximate) solution method, is given in Fig. 4. Here, we compare the simulation results for the typical J-V measurement displayed in Fig. 1 to the equivalent quantities computed using a combined asymptotic/numerical method. This alternative method is described in detail for a single-layer model in [5] and has been used to explain trends in experimental observations in [4]. Excellent agreement is shown between the current densities computed using the two methods.

Conclusions
We have built a fast and robust numerical solver for coupled ionic-electronic charge transport in a realistic three-layer perovskite solar cell architecture. The scheme is able to simulate a variety of relevant device operating regimes, including current-voltage sweeps and open-circuit transients, both Fig. 4 A comparison between the current density calculated using the IonMonger code (purple lines) and an alternative, combined asymptotic/numerical method from [4,5] (green circles) for the same simulation of a J-V measurement as shown in Fig. 1 (Color figure  online) with the possibility of having time-dependent illumination. Simulations of this sort can be carried out in seconds to minutes of computation time on a standard modern personal computer. The only prerequisite for making use of this tool is access to MATLAB and its suite of routines for time integration of ordinary differential equations, specifically the ode15s routine.
This work therefore provides a tool that is capable of playing a major role in guiding the development of perovskite solar cells. Our IonMonger code provides the possibility of independently varying each of the device parameters, so that their roles in determining cell performance can be discerned, something that is difficult, or even impossible, to achieve experimentally. One area of particular practical interest is understanding what can be done to mitigate the amount of parasitic recombination in PSCs, thereby further improving their performance. As demonstrated in [4], it is possible to suppress these losses via careful tuning of the cell's constituent material properties, and we anticipate that further studies in the same spirit will be made possible using the computational tool provided here. A second area where such a simulation tool is surely needed is in understanding the long-term degradation processes that occur within PSCs on timescales of between hours and weeks. While the current version of the code cannot simulate degradation, it can be used to predict the effects of different device parameters on some of the proposed causes of degradation. For example, degradation due to chemical reactions at the perovskite/ transport layer interfaces [3] is likely to be exacerbated by iodide ion accumulation in the Debye layers, while degradation due to the penetration of extrinsic ions, e.g. oxygen, into the perovskite may be enhanced by an accumulation of ion vacancies [1]. The development of IonMonger to include additional physical processes that occur on longer timescales would allow researchers to investigate long-term behaviour and stability much more quickly than is possible via experimentation, and hence it will be the subject of future work.
The authors are committed to maintaining and expanding the functionality of the code, and any updates will be hosted on the GitHub repository which can be accessed at https :// githu b.com/Perov skite SCMod ellin g/IonMo nger. While we cannot promise a high level of technical support to all users, we are happy to receive any suggestions on ways in which the features of the code can be improved and/or expanded, and it is our intention that the code will grow as the research priorities of the perovskite community evolve. Contact details for current code developers can be found in the README file in the main folder of the IonMonger GitHub repository.