IonMonger 2.0: software for free, fast and versatile simulation of current, voltage and impedance response of planar perovskite solar cells

The second generation of the open-source MATLAB-based software tool IonMonger, for solving drift–diffusion models of charge transport in planar perovskite solar cells, is presented here. This version is based upon a generalisation of the original drift–diffusion model of charge carrier and ion motion in the perosvkite cell, as described in Courtier (J Comput Electron 18:1435–1449, 2019). The generalised model has the flexibility to capture (1) non-Boltzmann statistics of charge carriers in the transport layers, (2) steric effects for the ions in the perovskite layer, (3) generation of charge carriers from light made up of a spectrum of different wavelengths and, (4) Auger recombination. The updated software is significantly more stable than the original version and also adds the ability to simulate impedance spectroscopy measurements as well as transient voltage and/or illumination protocols. In addition, it is fully backwards compatible with the original version and displays improved performance through refinement of the underlying numerical methods. Furthermore, the software has been made accessible to a wider user base by the addition of IonMonger Lite, a version that leverages MATLAB’s live scripts and eliminates the need for a detailed knowledge of MATLAB’s syntax.


Introduction
Over the past decade, charge transport modelling of perovskite solar cells (PSCs) has improved to a point where it can not only qualitatively reproduce most of the observed responses of these devices to a wide variety of experimental protocols but is also close to being predictive. Nevertheless there remain challenges, not only with the formulation of the charge transport model, but also with its accurate solution in parameter regimes of interest. The first published numerical scheme capable of solving such a model in relevant parameter regimes was presented in [20], while open source 1 software was made available only relatively recently by Courtier et al. [19] and Calado et al. [13]. The aim of this work is to publicise and describe an updated, and improved, version of the first of these two software packages.
The original version of the perovskite solar cell (PSC) simulation tool IonMonger [19] has been successfully used by numerous groups to simulate charge transport and cell performance in these devices. The main goal of the continuum models that IonMonger solves is to improve the understanding of PSCs and to predict device behaviour from its underlying physical properties. For example, IonMonger enables investigation of the effects of altering properties such as charge carrier diffusion lengths and recombination rates on current-voltage (or power) output. IonMonger includes a fully-coupled conservation equation for mobile ion migration in the perovskite layer of a three-layer PSC (see Fig. 1), which has been shown to play a dominant role in device behaviour [67], including steadystate performance [17,18].
Simulation results obtained using IonMonger have enabled researchers to tackle a variety of questions relating to the performance and design of PSCs. Examples of such work include: (1) Cave et al. [15] demonstrated how to determine activation energies for ion vacancy migration in different perovskite compositions using the information contained in current-voltage (JV) measurements; (2) Courtier [17] used IonMonger to verify a novel theory, termed the ectypal diode theory, for the steady-state performance of a PSC, which is notable for its inclusion of mobile ions; (3) Bennett et al. [8] extended this theory to describe the electrochemical impedance response of a PSC and identify a measurable value, analogous to the ideality factor for conventional solar cells, that may be used to diagnose the limiting form of recombination; (4) Riquelme, Castro-Chong, Anta and coworkers [14,54,55] validated this approach by investigating trends in impedance spectra using a combination of IonMonger simulations and investigation of experimental measurements; (5) Diekmann, Le Corre, Stolterfoht and colleagues used IonMonger to study the maximum efficiencies attainable from PSCs incorporating mobile ions in the perovskite layer [23,42], concluding that the presence of mobile ions is detrimental to power conversion efficiency (PCE); (6) Cordoba et al. [21] and Lin [45] evaluated another metric, namely the open-circuit voltage. Lin also simulated the impact of mobile ions on the short-circuit current [46]. (7) Li et al. [43] showed how ion migration can explain the lightsoaking phenomenon observed in experiment [65]; and (8) García-Rodríguez et al. demonstrated that inverted hysteresis can be explained by the drift-diffusion model, comparing IonMonger simulations to experiment [27].
Despite the undoubted success of IonMonger v1.0, and its growing user base, significant issues with its performance and capabilities have been brought to our attention by its early users; these have motivated the improvements to the numerical solver described in this work. Chief amongst the issues identified, and addressed here, are problems with the stability of the code. In addition to performance issues the validity of the Boltzmann model for charge carrier dynamics in the transport layers (particularly where these are composed of organic materials) has been called into question. This has motivated an extension to IonMonger which allows for non-Boltzmann models of the transport layer charge carrier dynamics. From a practical perspective, many experiments are performed with white light and therefore IonMonger v2.0 incorporates a model of charge carrier generation based on multiple wavelength illumination. Furthermore, adapting IonMonger v1.0 to simulate impedance experiments is cumbersome and so version 2.0 includes a module that allows impedance plots to be generated automatically and efficiently.
A development cycle similar to IonMonger's is being carried out to great effect in the field of lithium-ion battery modelling, by the PyBaMM community [34,59]. This builds upon collaborations both within academia and between academia and industry in order to accelerate R &D and serves as a model for the future path of IonMonger. In order to progress towards the ultimate goal of developing stable and commercially viable PSCs, further advances in our understanding of the interplay between ionic and electronic mechanisms, including their effects on degradation, are required. It is envisaged that this goal can be furthered by providing sophisticated and easy-to-use open-source modelling tools [61]. As such, IonMonger 2.0 includes a host of advanced features along with a new user-friendly interface.
Favourable comparisons have been made between Ion-Monger 1.0 and alternative simulation software such as Driftfusion [13], SCAPS (without mobile ions) [23] and a Python implementation of the IonMonger model [6]. As a consensus has not yet been reached on the correct model of charge transport in PSCs, it is vital that simulation software is open source to allow users to investigate the many proposed extensions to the model. The code presented in this paper, therefore is, to our knowledge, the most advanced open-source software for the simulation of charge transport and ion migration in standard three-layer planar PSCs available to-date.
In the next section, Sect. 2, we introduce the impedance spectroscopy module, outlining its features and giving a cursory introduction to its use in practice (a full set of instructions are given in the users' guide). Subsequently, in Sect. 3, we describe how the charge transport model has been extended to incorporate non-Boltzmann statistical models in the description of the transport layers. In Sect. 4, we discuss how steric effects have been accounted for in the equations defining the ion vacancy flux. The purpose of Sect. 5 is to introduce the remaining extensions to the charge transport model, namely; generation from light comprising a spectrum of wavelengths, Auger recombination, immobile ion distributions, parasitic series and shunt resistances (by embedding the charge transport model within an equivalent Fig. 1 Schematic of the basic drift-diffusion model of a planar PSC upon which IonMonger v1.0 was based. The continuum variables modelled in each layer are shown. Carrier generation is assumed to occur only in the perovskite and recombination in the perovskite and at the interfaces between the perovskite and the transport layers circuit), and metal contact band offsets. In Sect. 6, an introduction to IonMonger Lite is given. The penultimate section, Sect. 7, includes details of the improved numerical methods which give rise to better code performance, and how backwards compatibility with IonMonger 1.0 has been ensured. In Sect. 8 we draw our conclusions.
A complete statement of the charge transport model (system of PDEs) being solved in IonMonger 2.0 is given in "Appendix". A schematic diagram showing the new additions to the model is shown in Fig. 2 and these additions are listed in Table 1.

Impedance spectroscopy module
Electrochemical impedance spectroscopy (IS) is a measurement technique that is widely used in the study of electrochemical systems [51]; including in batteries [44,63] and fuel cells [32], where it is used to study corrosion and plating [38] and solar cells [8,53,54,56] where it can be used to identify efficiency and monitor degradation. It is a particularly useful technique as it is non-destructive and relatively easy to perform both in the lab and in the field.  In the context of PSCs, the interpretation of results is still relatively poorly understood, because of the novel complicating feature of ion motion, in addition to charge carrier motion. This motivates the need for accurate impedance simulations from a well-founded drift-diffusion model of the cell. IS is performed by perturbing a device from steady state by applying a low-amplitude voltage oscillation (over a range of frequencies) and measuring the current response. This elicits valuable information due to the different physical processes stimulated at different frequencies, allowing their effects to be disentangled. The applied voltage typically consists of a constant (or 'DC') component and a sinusoidal (or 'AC') component, where V DC is the steady-state voltage, and V p and are the amplitude and the angular frequency of the voltage perturbation, respectively. Keeping the amplitude of the voltage perturbation small (here, typically, < 20 mV) ensures that the current response is an approximately linear function of the applied voltage, taking the general form where J DC is the steady-state current, J p is the amplitude of the sinusoidal component and is the phase relative to the voltage. The complex impedance, Z( ) , is the ratio between the complex representation of the voltage perturbation and that of the current, i.e. (1) A so-called impedance spectrum (i.e. the impedance as a function of frequency, f) is shown in Fig. 3. It is common to decompose the impedance into its real and imaginary components; i.e. the resistance R and the reactance X, respectively, such that Z( ) = R( ) + iX( ) . The projection of a complex impedance curve onto the R-X plane forms a Nyquist plot, as shown in Fig. 4a. Nyquist plots restrict the frequency information that is displayed. Therefore, it is useful to display impedance spectra as their projection into the R-f and X-f planes, referred to here as frequency plots (Fig. 4b). Other representations of impedance such as Bode or capacitance plots can also be used to explicitly show the frequency dependence.
Generally, the impedance spectra of PSCs exhibit two semicircles on a Nyquist plot, one low-and one high-frequency feature. The two features indicate that there are (at least) two species of mobile charge (electronic charges and ionic charges) in a typical PSC [8,53,54]. The highfrequency feature excludes the transient effects of the slow-moving ionic charge and enables information to be obtained about the sources of electron and hole recombination (via the electronic ideality factor as defined in Bennett et al. [8]). The low-frequency feature captures how ions impact the cell response and provides information on the potential barrier to recombination (via the ectypal factor as defined in [17]) as well as the density and mobility of the ionic species [8,54]. In addition to the low-and high-frequency features, other features have been observed in the Nyquist plots of impedance spectra for PSCs. For example, these include a third semicircle [5,28] or a loop between the low-and high-frequency features [24,31]. These less

Simulating impedance spectra
IonMonger allows users to perform virtual IS measurements analogous to those performed by experiment. The voltage protocol is created via the applied_voltage input specified in the parameters file. An impedance simulation consists of the following steps. First, the cell is allowed to equilibrate at the DC voltage, thereby locating steady-state. IonMonger autonomously detects when the steady-state solution has been reached. This solution is then used as initial conditions for the transient stimulation at each different frequency. IonMonger then iterates over all sample frequencies, simulating the response to each voltage perturbation. Thus, an impedance protocol can be specified by six parameters: the minimum and maximum frequencies, the DC applied voltage, the AC voltage amplitude, the number of frequencies to be sampled and the number of complete periods to simulate. Detailed information on how to construct an impedance protocol from these six parameters can be found in the user guide (GUIDE.md).
The resulting structure array of solutions (one for each frequency) is passed to impedance_analysis.m to compute the impedance as a function of frequency. This is achieved by extracting the final two periodic cycles of the current density (at each frequency) and fitting a sinusoid to the data using FourierFit.m.
IS simulations are computed using the same numerical solver as transient simulations, meaning they not only achieve the same speed and accuracy, but are also fully compatible with additions to the charge transport model (detailed in Sects. 3,4,5). If MATLAB's Parallel Computing toolbox is installed, the iterations will be executed on multiple parallel workers, decreasing the compute time. The impedance spectra shown in Fig. 4, are each made up of 100 sample frequencies and took an average time of 48.7s to compute on a desktop computer with a six-core Ryzen 5 5600X CPU. Functions for plotting IS simulations are detailed in Sect. 7.3.

Carrier statistics in the transport layers
In previous versions of IonMonger, electronic carriers in all three layers were assumed to follow Boltzmann statistics, as is common in semiconductor modelling. However, this assumption, in the context of PSCs, has been called into question by recent works [2,61]. It is known that this assumption is not justified in the case of highly doped and/or organic materials, such as the transport layers [4,33,60]. To rectify this, IonMonger now allows carriers to be governed by more general band models, with the option of either a Fermi-Dirac or Boltzmann distribution. In this section we discuss the physical origin of statistical models and the specific options available in IonMonger 2.0 (see Table 2).

General statistical models
The density of free conduction electrons in a semiconductor is the product of the statistical distribution, f(E), and the conduction band density of states (DoS), ĝ c (E) , integrated across all energies, i.e.

Similarly, the density of valence band holes is
These integrals (referred to as statistical integrals) define a relationship between carrier densities and quasi-Fermi levels. As Fermions, electrons adhere to a Fermi-Dirac statistical distribution [48], defined as where E f is the Fermi level. In semiconductor modelling, it is common to approximate f(E) by a Boltzmann distribution [25,30,35,37,40,61], and assume that the DoS is parabolic, resulting in the familiar densities The Fermi-Dirac integral ( F ), the Gauss-Fermi integral ( G s ) and the Blakemore function ( B ) are defined in (24), (26), and (28), respectively Statistical integral Boltzmann approximation where g c,v are the effective densities of states for the conduction and valence bands, respectively, E f n,p are the electron and hole quasi-Fermi levels (QFLs), and E c,v are band edges. This is known as the Boltzmann approximation. Employing a statistical integral that is analytical, differentiable and invertible results in drift-diffusion equations that are significantly easier to solve, leading to widespread adoption of the Boltzmann approximation in semiconductor modelling, including the charge transport model upon which previous versions of IonMonger were based. However, the Boltzmann distribution fails as an approximation to the Fermi-Dirac distribution in the case of materials with high doping levels and/or non-parabolic band shapes, motivating a more general approach. Carrier densities in a general statistical model are given by [29] where S is the statistical integral. The energies E c and E v will henceforth be referred to as reference energies to account for DoS functions that do not possess the same defined edge as parabolic bands.
Electron and hole currents in a semiconductor are driven by gradients in their quasi-Fermi levels, meaning the current densities are defined as where n and p are the electron and hole mobilities, respectively. By rearranging (9), the quasi-Fermi levels can be written as functions of carrier density, where S −1 is the inverse of S such that S S −1 (⋅) ≡ ⋅ . The reference energies, E c,v , are dependent on the local electric potential, leading to the phenomenon known as band-bending. Specifically, we assume that the reference energies take the form E c,v = const. − q . Combined with (11), this allows one to express the electronic current densities as functions of carrier density and electric potential, Note that if carriers are assumed to be Boltzmann distributed in parabolic bands, the statistical integral is simply S(⋅) = exp(⋅) and the familiar expressions for current density, are recovered.

Statistical models in the PSC model
Carrier densities in the transport layers (TLs) of PSCs are typically several orders of magnitude larger than in the perovskite layer due to high effective doping levels and the band offsets between the layers. For this reason, the Boltzmann approximation cannot be guaranteed to be accurate in these layers and we adapt the charge transport model to allow for some general (possibly non-Boltzmann) statistical model. We assume that conduction electrons in the electron transport layer (ETL) are described by some statistical integral S E and valence holes in the hole transport layer (HTL) by S H . Thus the equations for current densities in the transport layers must be updated to reflect the choice of statistical model. We now have in the ETL and in the HTL. The generalised statistical models necessitate that continuity conditions at the edges of the transport layers must be reconsidered. As in previous versions of the model [19], we enforce that the majority carriers' QFLs are continuous across the interfaces between transport layers and the per- expressions for the QFLs in Boltzmann models (7), (8) in the perovskite and general statistical models (11) in the transport layers, we obtain the continuity conditions where a superscript ' + ' denotes a quantity on the right hand side of an interface and '−' the left hand side. The ratios of carrier densities either side of the interfaces when the cell is in equilibrium with no potential difference across it are therefore and we recast the continuity conditions in terms of these ratios as This concludes the changes required to incorporate the option of simulating non-Boltzmann statistics in IonMonger. In the next section we describe the common statistical models that are built-in to IonMonger 2.0.

Available statistical models
In IonMonger 2.0, the user can choose statistical models from three possible band shapes, with either Fermi-Dirac or Boltzmann distributions, to apply to majority carriers in the transport layers. Note that statistical models, determined by the combination of a statistical distribution and a band shape, are commonly given names, summarised in Table 2.
The two transport layers can be assigned different statistical models. Model parameters are used to create reference functions which enable quick and accurate evaluation of the statistical function and its inverse. Note that users who do not wish to make use of alternative statistical models can simply remove the 'stats' structure from the parameter file. In its absence, IonMonger will automatically employ Boltzmann distributions in parabolic bands, as was the case in previous versions of the software. For more information about specifying statistical models in simulations, see the user guide (GUIDE.md). Users also have the option of creating their own statistical models by editing the create_stats_funcs.m file which serves as a template.

The parabolic model
The parabolic band model (shown in Fig. 5a) originates from a Taylor expansion of the dispersion relation near a turning point (see [48] §3.3.5). This Taylor expansion describes the DoS near the edges of the bands (where most carriers exist). This model is usually considered accurate in the case of ordered crystalline semiconductors (typically inorganic).
The parabolic band structure is defined by the DoS functions The conduction and valence bands have defined edges, denoted by E c and E v , respectively. These are the reference energies. Under a Fermi-Dirac distribution, the statistical integral, often (somewhat confusingly) referred to as the Fermi-Dirac integral, is defined to be where and are the dimensionless quasi-Fermi level and state energy, respectively, where both are relative to the band's reference energy and measured in units of the thermal voltage. This function is plotted in Fig. 5b. To clarify the nomenclature, the Fermi-Dirac distribution is the statistical distribution obeyed by all Fermions and the Fermi-Dirac integral is the statistical integral to which only Fermions in parabolic bands adhere. As discussed in the previous section, if carrier densities are approximated by a Boltzmann distribution, the statistical integral becomes S( ) = exp( ) . Both the Fermi-Dirac integral and its Boltzmann approximation are shown in Fig. 5b. It is generally accepted that the Boltzmann approximation to the Fermi-Dirac integral is accurate for quasi-Fermi levels more than three thermal voltages from the band edge (or densities n <

The Gaussian model
Weak molecular bonds in organic materials create a disordered structure in which electron transport takes the form of hopping between discrete molecular sites; this is in contrast to the band transport in ordered crystalline materials. For a system with a large number of molecules, this hopping transport resembles band transport in a Gaussian band [11,22,50,60] where the width of the Gaussian corresponds to the level of disorder in the intermolecular structure. The Gaussian band structure is shown in Fig. 6a. The bands no longer have defined edges, hence the reference energies are now the centres of each Gaussian, denoted by E c and E v , respectively. In organic materials, these are commonly referred to as the lowest unoccupied molecular orbital (LUMO) and the highest occupied molecular orbital (HOMO), and denoted by E L and E H , respectively. This model requires a further parameter, s, to describe the width of the Gaussian (corresponding to the structural disorder). Note that s is dimensionless but is often given in its dimensional form, = sk B T . The Gaussian DoS functions are Under a Fermi-Dirac distribution, the statistical integral, referred to as the Gauss-Fermi integral [52], is defined as If carriers are approximated by a Boltzmann distribution, the statistical integral becomes S( ) = exp( + s 2 ) . Both the Gauss-Fermi statistical integral and its Boltzmann approximation are shown in Fig. 6b for different values of s. Note that the Boltzmann approximation to Gauss-Fermi statistics is far less accurate than the parabolic equivalent, and this is exacerbated for larger values of s. As shown in Fig. 6b, the Boltzmann approximation to G 8 ( ) is still highly inaccurate when the quasi-Fermi level is 50k B T from the reference energy, occurring at a dimensionless carrier density of approximately 10 −9 . Measurements of the Gaussian disorder

The Blakemore model
In the limit of vanishing Gaussian width ( s = 0 ), the Gaussian DoS becomes a Dirac-function, representing a density of g c,v states, all with identical energy E c,v . This is equivalent to hopping transport in materials with no structural disorder. Under a Fermi-Dirac distribution, this DoS yields a statistical integral, referred to as the Blakemore function, and shown in Fig. 7 [9], If carriers are approximated by a Boltzmann distribution, the statistical integral becomes S( ) = exp( ) .

Quasi-Fermi level input
Under Boltzmann distributions, the conversion between a carrier density and a quasi-Fermi level is simple; meaning that users can easily convert experimental measurements of equilibrium quasi-Fermi levels into effective doping densities to use in the parameter set for simulations. To avoid requiring users to perform this conversion for non-Boltzmann statistical models, IonMonger 2.0 can accept either carrier densities or quasi-Fermi levels when setting the TL doping levels. When the user sets the QFLs for the TLs, the relevant doping densities are calculated according to the statistical model in that layer. The two layers are independent, meaning one could set the doping density for one and the QFL for the other. Impact of statistical models on device performance. Finally, we note that the full charge transport model is sufficiently complex that it is impossible to predict how changes to the statistical model in one of the TLs might affect the model's response to any general experimental protocol. While outside the scope of this work, the impact of transport layer statistical models on device-level modelling will be investigated in a forthcoming publication [16].

Steric effects
The standard Poisson-Nernst-Planck (PNP) system, used to govern ionic motion in the previous version of IonMonger has an embedded assumption; namely that there is no shortage of lattice sites for ion vacancies to occupy, i.e. P ≪ P lim where P lim is the density of anion sites (and so is equivalent to the maximum vacancy density). 2 However, in relatively extreme scenarios, e.g. at very large applied voltage and/or in the narrow Debye layers near the edges of the perovskite, it is possible that ion vacancy densities become sufficiently large that steric effects (i.e. ion vacancy crowding, also known as volume exclusion effects) cannot be duly neglected. In such scenarios, the standard PNP system must be amended and replaced with a nonlinear-PNP (nPNP) system with a modified ion flux [1]. Two forms of modified ion flux have been used in the literature, derived using different assumptions. One approach yields a nonlinear (in the vacancy density) drift term in the ion flux whilst the other gives a nonlinear diffusion term. Rather than stipulating one form (or the other) of the ion flux, we allow users to easily switch between the two. We shall now present the two formulations in turn. Bazant [7] presents a general theory for chemical kinetics which describes how an ion flux is driven by the product of the ion density, mobility and gradient of its electrochemical potential = k B T ln  Another appealing justification for this formulation follows from work by Burger et al. [12]. This work considers a hopping model for the diffusion process that retains a non-zero probability that sites adjacent to the ion (or ion vacancy) in question are already occupied. In the derivation of the usual PNP system, this probability is set to zero. When incorporating steric effects, the probability of transition to an occupied site is zero, thereby enforcing that at most one ion can reside on any given lattice site. This constraint results in a densitydependent mobility and hence a modified drift term.

Nonlinear diffusion
Alternatively, steric effects can be modelled via a nonlinearity in the diffusion term, as proposed by Kralj-Iglic and Iglic [41] and Borukhov et al. [10] according to thermodynamic considerations (assuming a constant mobility). In this case, the ion flux is given by Note that diffusion is "enhanced" as P approaches P lim . This formulation has previously been used in the context of PSC modelling [2,13]. Abdel et al. [2] show that the form of diffusion enhancement depends on the choice of statistical function and that the formulation in (30) results from employing a Blakemore model (or, equivalently, a Fermi-Dirac integral of order − 1) to describe the movement of ions.
Both expressions for the ion flux collapse to the form found in the usual PNP system on taking the limit P lim → ∞ . Furthermore, it is important to note that even though the two formulations predict the same steady state ((29) is identical to (30) on setting F p ≡ 0 ), there are marked differences in their dynamics. Figure 8 shows how the dynamics are affected by variations in P lim for the two models, whilst Fig. 9 verifies that IonMonger produces identical steady states for the two different nPNP models as expected.
Steric effects can be activated in IonMonger by setting 'Plim' to a numeric value and 'NonlinearFP' to a string value of either 'Drift' or 'Diffusion'.

Other extensions to the model
In this section we briefly introduce a number of other extensions that have been added during the IonMonger update. See the User Guide for more information.

Absorption spectra
By default, the generation rate G(x, t) follows the Beer-Lambert law of light absorption for a single wavelength and absorption coefficient [48], i.e. Timescales for the current density reaching steady state at an applied voltage of 1 V after a preconditioning step at 0.9 V with different vacancy limitation parameters. Modified drift (29) and diffusion terms (30) show different dynamics but reach the same steady state where F ph is the flux of photons incident on the light-facing perovskite surface (after accounting for reflection) under the equivalent of 1 Sun illumination; is the light absorption coefficient of the perovskite; I s (t) is the intensity of the illumination in Sun equivalents; and, the parameter l determines whether light enters through the ETL ( l = 1 ) or through the HTL ( l = −1 ) for a so-called inverted architecture. In IonMonger 2.0 the user now has the option, via the generation_profile function, to replace (31) with a Beer-Lambert generation profile for light that includes a range of wavelengths, by integrating an absorption coefficient spectrum versus photon energy as follows: In this expression, the incident photon flux and the absorption coefficients are now functions of photon energy. The function can be created from measurement data.

Auger recombination
The total bulk recombination rate R(n, p) can now include contributions from Auger recombination as well as bimolecular and/or SRH recombination. The Auger recombination rate takes the form where A n and A p are the Auger coefficients (with units m 6 s −1 ). This choice ensures that R Auger = 0 at thermal equilibrium, i.e. when np = n 2 i . The total bulk recombination rate is The full expressions for each of these rates are given in the "Appendix".

Immobile ions
IonMonger 2.0 allows for the possibility of immobilising ion vacancies (previously all anion vacancies were assumed to be mobile). The user may set the ion vacancy diffusion coefficient D I equal to zero in order to fix them in their initial distribution. In previous versions, this would cause the characteristic timescale (by which time is rescaled) to become infinite, preventing simulations from running. To correct for this, the numerical scheme now uses a characteristic timescale defined as follows.
The first option is a timescale for electron transport; the second is the ionic timescale used in IonMonger 1.0 [19].

Band offsets at metal contacts
Differences in the band energy levels between the transport layers and the metal contacts are modelled by updating the boundary condition on the majority carrier in each TL as follows.
If the cathode and anode workfunctions are not set, Ion-Monger 2.0 assumes flat-band conditions, as before, which are equivalent to setting

Parasitic resistances
In order to model parasitic resistances, we follow the work of Neukom et al. [49] and embed the PSC drift-diffusion model in an external circuit containing two resistors, as shown in Fig. 2b. Even though this modelling approach (of appending resistors to the PSC element) is essentially ad-hoc, it does allow the model to capture parasitic series and shunt(/ parallel) resistances which may be important in practice. An alternative approach to capturing these resistances would be to upgrade the model to include more spatial dimensions, more of the PSC architecture and a more complex (nonplanar) geometry. However, such complexity would render the model significantly more costly to both solve and parameterise, nullifying some of the benefits of IonMonger.
The effect of the external circuit is taken account of through two modifications to the drift-diffusion model.
Within the PSC, the total current density is independent of x and can be computed at any point in the cell from The first term on the right-hand side accounts for current losses due to shunt resistance in the external circuit. V p denotes the potential difference across the parallel resistor and, equivalently, the PSC element as shown in Fig. 2b. To minimise numerical error, IonMonger automatically calculates the current density at the midpoint of the perovskite layer, where the grid spacing is larger and the solution varies more smoothly. The terms involving time-derivatives give the displacement current density, which can become dominant during very high frequency (IS) measurements. Secondly, we have to take account of the change in V p due to the resistor in series. We choose to incorporate this change at the HTL/metal contact, leading to the boundary condition where R s (Ohm) is the series resistance, R p (Ohm) is the parallel resistance and A (cm −2 ) is the cross-sectional area of the cell perpendicular to the current. We implicitly assume that the displacement current is negligible at the boundary. With these two modifications, the impact of parasitic resistances on cell performance and time-dependent behaviour can be easily investigated alongside the full set of cell parameters.

Easy access via IonMonger Lite
Previously, some knowledge of MATLAB syntax was required to edit the parameters file and run simulations. In contrast, IonMonger Lite (provided as part of Ion-Monger 2.0) is aimed at users with no expertise in the MATLAB language.
IonMonger Lite is a MATLAB live script, i.e. an interactive document that displays formatted text, equations and images [47]. The underlying code can be hidden, thereby providing a more accessible user-experience. The IonMonger Lite interface enables users to run two types of simulation: a single current-voltage ( J-V ) sweep, (38) preceded by a preconditioning step, or impedance spectroscopy (IS). The cell parameters and simulation protocol can be edited in text boxes and drop-down menus. Figure 10 depicts the available protocols, including the adjustable parameters. For the J-V sweep, the six parameters are the steady-state initial voltage, V init ; the voltage at which the sweep begins, V pre ; the target sweep voltage, V tar ; the length of the preconditioning phase, t pre ; the length of the sweep, t sweep ; and the constant light intensity, L. For impedance spectroscopy, the adjustable parameters are the DC voltage, V DC ; the AC voltage amplitude, V p ; the minimum and maximum frequencies, f min and f max ; the number of sample frequencies, n f ; and the constant light intensity, L.
Users are presented with a series of checkboxes to determine which plots will be produced at the end of the simulation. Possible plots for J-V sweeps include the applied voltage and light intensity as functions of time; the current density as a function of either time or applied voltage; and distributions of the anion vacancy density, electron density, hole density or electric potential. For IS simulations, Nyquist and frequency plots such as those in Fig. 4 can be produced. For J-V sweeps, there is also the option to automatically generate an animation of the solution using animate_sections.m (see Sect. 7.3).
While offering a simplified user experience, live scripts lack the flexibility of standard MATLAB code. For this reason, proficient users should continue to run simulations by editing a parameters file and running 'master.m' to make use of IonMonger's full functionality. Useful features of the full version that are not possible in IonMonger Lite include: Fig. 10 The two types of experimental protocol available in Ion-Monger Lite. a Current-voltage sweep protocol with 6 adjustable parameters. b Impedance spectroscopy protocol with two variable voltage parameters, three frequency parameters and one light intensity parameter • advanced simulation protocols involving more sections, asymmetric sweeps, multiple consecutive sweeps, or time-dependent illumination intensity • open-circuit voltage tracking • non-Boltzmann statistical models in the transport layers • running batches of simulations where one or more variables are iterated through a list of values • changing the resolution and error tolerances of the solver • band offsets at the metal contacts • parasitic resistances • resuming saved simulations IonMonger Lite uses the same numerical solver as the full version, meaning it achieves the same performance and accuracy. After the simulation, the solution file is saved in the same format as the full version; thus all of the same postprocessing and analysis/plotting can be performed.

Importing data into Python
Once a simulation is complete, all the data, including the inputs, are saved in a single file. These data can then be imported into any software of the user's choice for further analysis. The transfer of this data from MATLAB to Python is nontrivial because Python has no direct equivalent data type to the structures used by MATLAB. To aid this transfer, we include a Python file named IonMonger_import.py that can import a solution and unpack the data for analysis. This file contains functions to extract all major variables from a saved .mat file as well as example code to generate plots of the data and further instructions. Combined with IonMonger Lite, this function allows users with experience in Python but no expertise in MATLAB to design, run, and analyse simulations with IonMonger, an important step in making IonMonger accessible to a wider scientific community.

Performance, compatibility and accessibility
In this section, we give details of the features of IonMonger 2.0 that enable the improved performance, backwards compatibility and greater accessibility via new functionalities for both running the code and visualising the results.

Performance
In the original version of IonMonger, the Jacobian ( J ) of the DAE system was approximated numerically by ode15s. In IonMonger 2.0, performance is improved by analytic calculation of J in the function AnJac.m. Figure 11 demonstrates the resulting increase in performance for simulations of current-voltage sweeps using the default parameter set. The simulation time for an impedance spectrum with 60 sample frequencies is decreased from 61 to 31 s. Similarly, IonMonger 2.0 benefits from improved stability in response to different parameter sets, due to many small changes that aid the solver, ode15s, in converging on a solution.

Fast and accurate evaluation of statistical integrals and their inverses
The challenge of implementing alternative statistical models is the large number of comparatively costly evaluations of the inverse statistical integral, S −1 (47b,49b). In general, not only is the statistical integral non-invertible, but S also cannot be evaluated exactly. This is the case in both the Fermi-Dirac (24) and Gauss-Fermi (26) integrals. The simplest approach would be to evaluate the integrals numerically and invert the function with a root-finding algorithm. However, even with efficient functions to perform the integration and root-finding, this is computationally expensive and would result in J-V simulations taking minutes to hours instead of seconds to minutes motivating the need for more efficient schemes. Significant efforts have been made to find approximations to the Fermi-Dirac [9,26,64] and Gauss-Fermi integrals [52,58]. However, most approximations are either computationally expensive or lack accuracy. Furthermore, attempts to approximate the inverses of these integrals are rare. In order to maximise both efficiency and accuracy, IonMonger 2.0 utilises interpolation of lookup tables of dimensionless quasi-Fermi levels and their corresponding dimensionless carrier densities, evaluated once at the start of a simulation by numerical integration using the trapezium rule on a high-resolution grid. This approach eliminates the need for a root-finding algorithm as the interpolation can be performed in either direction.
The inverted statistical integrals are evaluated inside an approximated derivative within the DAE system; therefore it is essential that the approximations have continuous first derivatives. For this reason, interpolation is performed using piecewise cubic Hermite interpolating polynomials (PCHIPs), rather than linear interpolation, via MATLAB's interp1 function. This approach results in approximations for F( ) and G s ( ) with a relative error of less than 0.1% in the domains { < 6} and { < 2s , s ≥ 1} , respectively. Figure 12 shows the simulation times for a single current-voltage sweep with and without Boltzmann statistics in the transport layers for a range of spatial grid spacings. The time taken to create the lookup tables at the beginning of each simulation averaged 39 ms for the Fermi-Dirac, and 30 ms for the Gauss-Fermi integral. The average extra time taken for simulations of a 100 mV s −1 J-V sweep of the template parameter set was just 0.54 s. Therefore the effects of alternative statistical models may be investigated using IonMonger without significant negative impact on computation time.

Compatibility
IonMonger 2.0 is fully compatible and consistent with previous versions. Any parameter file from a previous version can be run using version 2.0. If any parameters required for additions to the model are not specified in the parameter file, they will automatically be replaced by their default values, listed in Table 1. An example of the verification that IonMonger 2.0 is indeed consistent with the initial release is shown in Fig. 13, where electric potential distributions calculated on high resolution grids from both versions are shown to overlap exactly. Though not shown, this verification has been performed on all output variables across multiple parameter sets.
Two tests for developers have been written using MAT-LAB's testing framework. The aim of these tests is to verify that future versions of the code are consistent with the initial and subsequent releases. The first is a regression test that compares the output of an example simulation with saved data produced by the original code. The second is an integration test that checks that the current code runs successfully for two different experimental protocols that make use of different subroutines.

New analysis tools
Five new analysis tools are provided as part of IonMonger 2.0 and can be found in the Code\Plotting folder.

Plotting functions
Four of these new tools are plotting functions. For more details, see the README file. Fig. 12 Time taken to compute a 100 mV s −1 J-V sweep of the template parameter set with different statistical models across a range of grid resolutions. Red lines correspond to simulations using Fermi-Dirac statistics in the ETL and Gauss-Fermi in the HTL as the parameter set is based on an inorganic ETL and an organic HTL and blue lines correspond to simulations in which the Boltzmann approximation was employed in both transport layers. Circles were measured on a desktop with a Ryzen 5 5600X CPU and triangles were measured on a laptop with a Ryzen 5 Windows Surface Edition CPU (Color figure online) Fig. 13 Verification of IonMonger 2.0 (dotted lines) against version 1.0 (solid lines) for solutions of electric potential during a voltage sweep on a high resolution grid. Note that the electric potential has been shifted such that = 0 on the HTL contact for clarity and the spatial grid has been stretched such that grid points are linearly spaced in order to examine the behaviour in the Debye layers • plot_recombination.m generates two figures showing the recombination rates as functions of time and space. • plot_dstrbns.m plots spatial distributions of the four model variables ( P, , n, p). • plot_IS.m generates Nyquist and frequency plots as in Fig. 4, as well as a 3D plot of the impedance spectrum as shown in Fig. 3. • plot_bands.m generates a band diagram including the position of references energies and QFLs.

Animating solutions
The function animate_sections.m can be used to animate a solution and save the video as an MP4 file. Details of how to use the function and specify the parameters of the video (length, resolution, frame rate, etc) can be found in the user guide (GUIDE.md). Example videos can be found on the IonMonger Modelling YouTube channel [36] and one is included in the SI.

Conclusions
We have presented the second incarnation of the PSC simulation tool IonMonger, expanding on the success of the original software by: (1) augmenting the charge transport model to account for a diverse list of physical processes, including non-Boltzmann statistics, steric effects for mobile ions, absorption spectra, Auger recombination, and parasitic resistances; (2) adding additional functionality, particularly in regards to simulating impedance responses; (3) improving accessibility and lowering the barrier-to-entry via the introduction of IonMonger Lite (allowing users to leverage many features of the tool with little to no coding or MATLAB experience), and; (4) increased numerical performance and stability. In addition, the performance of the tool has been increased across the board, meaning decreased computation times. This is primarily due to the addition of an analytic Jacobian in the numerical scheme and enables the computation of full, well-resolved, impedance spectra on timescales of seconds to minutes and current-voltage sweeps on timescales of seconds on standard desktop computers. The open-source nature of IonMonger is supported by a dedicated, free-to-join Slack community of users and developers 3 who will further develop both the charge transport model and the code in the future. Thus, the software provides an accessible in-silico laboratory for investigating the physics of PSCs, and it is the authors' hope that it can contribute towards accelerating PSC development.

Appendix: The full charge transport model
In this section we state the full charge transport model solved by IonMonger 2.0. The notation is consistent between versions 1.0 and 2.0; a complete list of symbols and their definitions is provided in the SI. For an overview of the changes between the versions of IonMonger, refer to Fig. 2 and Table 1. In order to make use of model extensions, the parameters listed in Table 1 need to be specified, via the parameters file, otherwise they are set to their default values for backwards compatibility. With or without the model extensions, IonMonger 2.0 outperforms 1.0 due to improvements in the numerical implementation, as described in Sect. 7.1. Equations that have been altered by the model extensions are highlighted in bold.
Perovskite absorber layer ( 0 < x < b ) Following [19], our model comprises statements of conservation of conduction electrons, valence holes, and halide ion vacancies, as well as drift-diffusion equations for the fluxes of these species, Charge carrier photo-generation and recombination within the perovskite are included via the functions G(x, t) and R(n, p), respectively. These are coupled with Poisson's equation for the electric potential, We assume a constant, uniform background density ( N 0 ) of immobile cation vacancies.