Driftfusion: an open source code for simulating ordered semiconductor devices with mixed ionic-electronic conducting materials in one dimension

The recent emergence of lead-halide perovskites as active layer materials for thin film semiconductor devices including solar cells, light emitting diodes, and memristors has motivated the development of several new drift-diffusion models that include the effects of both electronic and mobile ionic charge carriers. In this work we introduce Driftfusion, a versatile simulation tool built for modelling one-dimensional ordered semiconductor devices with mixed ionic-electronic conducting layers. Driftfusion enables users to model devices with multiple, distinct, material layers using up to four charge carrier species: electrons and holes plus up to two ionic species. The time-dependent carrier continuity equations are coupled to Poisson’s equation enabling transient optoelectronic device measurement protocols to be simulated. In addition to material and device-wide properties, users have direct access to adapt the physical models for carrier transport, generation and recombination. Furthermore, a discrete interlayer interface approach circumvents the requirement for boundary conditions at material interfaces and enables interface-specific properties to be introduced.


Introduction
Accurate models of semiconductor devices are essential to further our understanding of the key physical processes governing these systems and hence rationally optimise them. One approach to modelling devices on the mesoscopic scale is to use continuum mechanics, whereby charge carriers are treated as continuous media as opposed to discrete particles. Typically, electronic carriers are modelled at discrete energy levels with a transport model describing the dynamics of carriers in response to an electric field (drift) and carrier density gradients (diffusion). This drift-diffusion (Poisson-Nernst-Planck) treatment leads to a system of coupled partial differential equations (the van Roosbroeck system [1]): a set of continuity equations, defining how the density of each charge carrier changes with time at each spatial location, are coupled with Poisson's equation (Gauss' Law), which relates the space-charge density to the electrostatic potential. For many architectures of thin-film semiconductor device (with the notable exception of transistors), provided that the materials are homogeneous and isotropic, it is sufficient to model devices with properties that vary in a single spatial dimension. In all but the most elementary of cases the resulting system of equations must be solved numerically.

The emergence of lead-halide perovskites and recent progress in mixed electronic-ionic conductor device models
The recent emergence of lead-halide perovskites (referred to herein as perovskites) as active layer materials for thin film semiconductor devices including solar cells, light emitting diodes (LEDs), and memristors has motivated the development of several new drift-diffusion models that include mobile ionic species in addition to electronic carriers. [2,3,4,5,6,7,8,9,10] Ab initio calculations and experimental evidence has shown that the charge density distribution, and consequently the electric field, in perovskite materials is dominated by high densities of relatively slow-moving mobile ionic defects. [11,12,13,14] This has a profound impact on the optoelectronic response of devices with perovskite active layers, leading to strong hysteresis effects in experimental measurements on timescales from microseconds to hundreds-of-seconds. [15,16] To date, both experimental and theoretical research into perovskites has primarily focussed on their application as a photovoltaic absorber material for solar cells and we now review the recent advances in device-level modelling in this field of application. Van Reenen, Kemerink & Snaith were the first to publish perovskite solar cell (PSC) simulations using a coupled model that included continuity equations for three charge carriers: electrons, holes and a single mobile ionic species. [3] They found that current-voltage (J-V ) hysteresis in PSCs could only be reproduced by including a density of trap states close to one of the interfaces acting as a recombination centre. [3] Later calculations of a 1.5 nm Debye length 1 [4] suggested, however, that the choice of a 4 nm mesh spacing in the simulations was too coarse to properly resolve the ionic charge profiles at the perovskite active layer-transport layer interfaces (described herein simply as interfaces). Richardson and co-workers overcame the numerical challenge of high ionic carrier and potential gradients at the interfaces by using an asymptotic analytical model to calculate the potential drop in the Debye layers of a single mixed electronic-ionic conducting material layer. [2,4] While this approach enabled the reproduction of hysteresis effects using high rates of bulk recombination, the inability to accurately model interfacial recombination limited the degree to which the simulation could represent real-world devices. [4] In a later publication by the same group modelling dark current transients, interfacial recombination was implemented, but only at the inner boundary of the Debye layer. [17] Furthermore, since these models were limited to a single layer, unrealistically large ionic charge densities were calculated at the interfaces as compared to three-layer models with discrete electron and hole transport layers (ETL and HTL respectively).
Our own work simulating PSCs began with a threelayer p-i-n dual homojunction model in which the p-and n-type regions simulated the HTL and ETL and where interfacial recombination was approximated by including high rates of recombination throughout these layers. Our results supported van Reenen et al.'s conclusion that both mobile ions and high rates of interfacial recombination are required to reproduce J-V hysteresis effects and other comparatively slow transient optoelectronic phenomena in p-i-n solar cells. [18] Shortly after Neukom et al. published a modelling study with similar conclusions. [5] They used the commercial package SETFOS [19] to solve for 1 Based on an ion density of 1.6 × 10 19 cm −3 . electronic carriers in combination with a separate MATLAB code that solved for the ionic carrier distributions. More recently, Courtier et al. published results from IonMonger, a freely-available, fully-coupled, three-layer device model that included a single ionic charge carrying species and boundary conditions at the interfaces such that surface recombination of electronic carriers at interfaces between the different layers could be explicitly included. [8,20] There remained some limitations with the model however; only the majority carriers were calculated in the ETL and HTL, excluding the possibility of simulating single carrier devices, and intrinsic or low-doped transport layers such as organic semiconductors; ions were confined to the perovskite layer and; users only have the possibility to simulate three-layer devices. Jacobs et al. also published results from a three-layer coupled electronic-ionic carrier simulation implemented using COMSOL Multiphyics® [21] and MATLAB Livelink™. [22,23] Most recently, Tessler & Vaynzof published impressive results from a similar three-layer PSC device model that included the option to use either Boltzmann or Fermi-Dirac statistics.
[24] Notwithstanding, the methodological details from both Jacobs et al. [23] and Tessler & Vaynzof [24] are sparse and at the time of writing neither code is publicly available.

Driftfusion:
An open source code for simulating ordered semiconductor devices with mixed ionic-electronic conducting materials in one-dimension Here we present a comprehensive guide to Driftfusion, our open source simulation tool designed for simulating semiconductor devices with mixed ionic-electronic conducting layers in one dimension. The software (based in MATLAB) enables users to simulate devices with any number of distinct material layers and up to four charge carrying species: electrons and holes by default plus up to two ionic species. The time-dependent continuity equations are fully-coupled to Poisson's equation enabling transient optoelectronic measurements to be accurately simulated. In addition to common material parameters, users have direct access to adapt the carrier transport, recombination and generation models as well as the system's initial and boundary conditions.[25] Driftfusion uses a discrete interlayer interface approach for junctions between material layers (heterojunctions) such that energetic and carrier density properties are graded between adjacent layers using a range of grading options. This method has the added benefits that it both circumvents the requirement for boundary conditions at heterojunctions, and enables interface-specific properties to be defined within the interface regions. While the example architectures and outputs given in this work use PSCs as a model system, Driftfusion can, in principle, be used to model any ordered one-dimensional mixed ionic-electronic semiconductor or redox system for which the drift-diffusion approach is valid. This work is divided into four main sections; we begin with a general overview of the simulation tool in Section 2; in Section 3 the default physical models for charge carrier transport, generation and recombination are outlined; Section 4 provides a detailed description of the system architecture and a step-by-step guide of the important commands and functions that will enable readers to get started with using Driftfusion; we conclude in Section 5 by comparing calculations from Driftfusion to two analytical and two numerical models to validate the simulation solutions and the discrete interface approach.

Workflow
A flow diagram summarising Driftfusion's general workflow is given in Figure 1. The system is designed such that the user performs a linear series of steps to obtain a solution; (1) the process begins with the creation of a semiconductor device object for which both device-wide properties, such as the carrier extraction coefficients, and layer-specific properties, such carrier mobilities are defined. A user-definable physical model comprised of one-dimensional generation, recombination and transport models determines the continuity equation for each charge carrier (see Section 3 for the default expressions); (2) the continuity equations are solved simultaneously with Poisson's equation (Equation 17) to obtain a solution for the electron density, hole density, cation density (optional), anion density (optional), and electrostatic potential distributions at equilibrium; (3) an experimental protocol such as a current-voltage scan is defined using the appropriate input parameters e.g. scan rate and voltage limits. The protocol generates time-dependent voltage and light conditions that are subsequently applied to the device, typically using the equilibrium solution as the initial conditions. In more sophisticated protocols the solution is broken into a series of steps whereby intermediate solutions are fed back into the solver. Likewise, protocols can be cascaded such that the solution from one protocol supplies the initial conditions for the next; (4) the desired solution is output as a MATLAB data structure (see Solution structures highlighted box); (5) the solution structure can be analysed to obtain calculated outputs such as the charge carrier currents, quasi-Fermi levels etc; (6) a multitude of plotting tools are available to visualise the simulation outputs. Instructions on how to run each step programmatically and further details of the system architecture, protocols, solutions, and analysis functions are given in Section 4.

Licensing information
The front end code of Driftfusion has been made opensource under the GNU Affero General Public License v3.0 in order to accelerate the rate of development and expand our collective knowledge of mixed ionic-electronic conducting devices. [26] It is important to note, however, that Driftfusion currently uses MATLAB's Partial Differential Equation solver for Parabolic and Elliptic equations (pdepe), licensed under the MathWorks, Inc. Software License Agreement, which strictly prohibits modification and distribution. If you use Driftfusion please consider giving back to the project by providing feedback and/or contributing to its continued development and dissemination.
We now proceed to describe the default physical models underlying this release of Driftfusion.
[26] Herein relevant Driftfusion functions and commands are highlighted using boxes and referred to using command line typeface.
Solution structures Following successful completion of the steps given in Figure 1, Driftfusion outputs a MATLAB structure sol (known herein as a solution structure) containing the following elements: -The solution matrix u: a three-dimensional matrix for which the dimensions are [time, space, variable]. The order of the variables are as follows: 1. Electron density 2. Hole density 3. Cation density (where 1 or 2 mobile ionic carriers are stipulated) 4. Anion density (where 2 mobile ionic carriers are stipulated) -The spatial mesh x.
-The time mesh t.
-The parameters object par.
As illustrated in Section 2.2, sol can be used as the input argument for analysis functions contained within dfana or plotting functions within dfplot. See Section 4 for further details. 1. Define layer and device-wide properties e.g. semiconductor band energies, layer thickness etc.
2. Solve drift diffusion equations to obtain equilibrium (dark) state electron, hole, cation, anion and electrostatic potential profiles 3. Apply experimental protocol defining light and voltage conditions to device e.g. current voltage sweep. See Table 2 for further examples 4. Solve drift diffusion equations to obtain electron, hole, cation, anion and electrostatic potential profiles as a function of time for the conditions defined in (3) 5. Analyse the solution to obtain calculated quantities e.g. charge carrier currents, quasi Fermi levels etc. 6. Plot the outputs e.g. energy level diagrams, current-voltage curves etc. The device equilibrium state is solved for using the given recombination and transport models; 3. An experimental protocol, which defines time-dependent voltage and optical generation conditions, is applied to the equilibrium solution; 4. A solution is obtained that may be fed back into the protocol until the desired solution is reached; 5. The solution is analysed to obtain calculated outputs; 6. The outputs are visualised using plotting tools.

Implementation of established semiconductor theoretical principles in Driftfusion
The device physics implemented in Driftfusion is principally based on established semi-classical semiconductor transport and continuity principles, which are well described in Sze & Kwok [27] and Nelson[28]. Elements of this section are adapted from Ref. [29] and are provided here as a direct reference for the reader. The equations described herein are written in terms of a single spatial dimension and can only be applied to devices with onedimensional architectures and for which the material layers are homogeneous.
Driftfusion evolved from a diffusion-only code written to simulate transient processes in dye sensitised solar cells [30] and uses MATLAB's [22] built-in pdepe solver. [31] The code solves the continuity equations and Poisson's equation for electron density n, hole density p, cation density c (optional), anion density a (optional), and the electrostatic potential V as a function of position x, and time t.
The full details of the numerical methods employed by the pdepe solver for discretising the equations are given in Skeel & Berlizns 1990.[32] 3.1 Semiconductor energy levels Figure 2a shows the energy levels associated with an idealised intrinsic semiconductor. The electron affinity Φ EA and ionisation potential Φ IP are the energies required to add an electron to the conduction band (CB) from the vacuum level E vac and to remove an electron from the valence band (VB) to E vac respectively. Note that in contrast to the established convention and the description given here, in Driftfusion Φ EA and Φ IP are input as negative values for consistency with other energetic properties referenced using the electron energy scale.
The electronic band gap E g of the material can be defined as:

The vacuum energy
The vacuum energy, E vac is defined as the energy at which an electron is free of all forces from a solid including atomic and external potentials.
[28] Spatial changes in the electrostatic potential, V are therefore reflected in E vac such that, at any point in space, where q is the elementary charge.

Conduction and valence band energies
The conduction and valence band energies E CB and E VB are defined as the difference between the vacuum energy and Φ EA and Φ IP respectively.
The band energies then, include both the potential energy associated with the specific molecular orbitals of the solid and the electrostatic potential arising from the existence of charge both within and external to the material.

The occupation probability distribution function and electronic equilibrium carrier densities
At equilibrium the net exchange of mass and energy into and out of, as well as between different locations within a system is zero. Under these conditions the average probability, f that an electron will occupy a particular state of energy, E at equilibrium in a semiconductor at temperature T is given by the Fermi-Dirac (F-D) distribution function: where k B is Boltzmann's constant. The equilibrium Fermi energy E F0 defines the energy at which a hypothetical electronic state has a 50% probability of occupation. At equilibrium, and for V = 0, the Fermi energy is identical to the chemical potential of the material. For an intrinsic semiconductor, E F0 lies close to the middle of the gap.
Where the semiconductor is p-type, dopant atoms accept electrons from the bands, shifting E F0 towards the valence band ( Figure 2b). Similarly, where the semiconductor is n-type, dopants donate electrons to the bands, shifting E F0 towards the conduction band ( Figure 2c). Note that herein the subscript '0' denotes the value or expression of properties at equilibrium e.g. E CB,0 is the conduction band energy at equilibrium. 2 To obtain the density of free electrons in the conduction band, n, the product of the probability distribution function, f and the conduction band density of states (DOS) function, g CB is integrated across energies above the conduction band edge, E CB : Similarly, to obtain the density of free holes in the valence band p, the product of the average probability that an electron is not at energy E (i.e. (1 − f )) with the valence band DOS function, g VB is integrated across all energies up to the valence band edge E VB : For semiconductor materials, g CB and g VB are typically modelled as parabolic functions with respect to electron energy at energies close to the band edges. Specifcally, , where m * e and m * h are the effective electron and hole masses and h is Planck's constant. Unfortunately, closedform solutions cannot be found to Equations 6 and 7 using the parabolic band approximation and the F-D distribution function (Equation 5). Hence, we use Blakemore's approximation[40] to the above integrals to obtain closedform expressions for the equilibrium carrier densities: Intrinsic equilibrium n-type equilibrium p-type equilibrium where N CB and N VB are the temperature-dependent 3 effective density of states (eDOS) 4 of the conduction and valence bands respectively. γ is a constant defining how close the approximation is to the Boltzmann regime (note that Equations 8 and 9 reduce to the Boltzmann approximation for γ = 0). Following Farrell et al., we set γ = 0.27 by default for ordered materials.
[41] This results in a close agreement to F-D statistics for E VB − 1.3k B T < E F0 < E CB + 1.3k B T , such that degenerate semiconductor states are permissible within the scope of the model.

Equilibrium Fermi levels in doped materials
In charge-neutral n-type materials the equilibrium electron density is approximately equal to the density of donor dopant atoms such that n 0 ≈ N D . Similarly in p-type materials, p 0 ≈ N A , where N A is the density of acceptor dopants. In Driftfusion users input values for E F0 for each material layer and the corresponding equilibrium carrier and doping densities, n 0 , p 0 , N D , and N A are calculated during creation of the device parameters object (see Section 4.2) according to Equations 8 and 9.
The equilibrium carrier densities n0 and p0 and Fermi levels EF0 for individual material layers are calculated and stored as a function of position in the device structures dev and dev sub of the device parameters object par. See Subsection 4.2.5 for further details. Note that when the material layers with different equilibrium Fermi levels are brought into contact, n 0 , p 0 , and E F0 become position-dependent owing to the creation of a space charge regions and associated electric fields. 3 For simplicity, the temperature-dependence of N CB and N VB has been omitted from the equations herein and it should further be noted that this temperature dependence is not explicitly dealt with in this release of Driftfusion.

Quasi-Fermi levels
A key approximation in semiconductor physics is the assumption that, under external optical or electrical bias, the electron and hole populations at any particular location can be treated separately, with individual distribution functions and associated quasi-Fermi levels (QFLs), E Fn and E Fp (Figure 2d). This is permitted because thermal relaxation of carriers to the band edges is typically significantly faster than interband relaxation, resulting in quasi-equilibrium states for each population.
[28] Under these circumstances a similar approach to that taken for the true equilibrium state can be used to derive expressions for the QFLs: It can be helpful to conceptualise the QFLs as the sum of the electrostatic (V ) and average chemical potential energy (k B T ln(n/N CB − γ) − Φ EA for electrons, −k B T ln(p/N VB − γ) − Φ IP for holes) components of the carriers at each location. It follows that the gradient of the QFLs provides a convenient way to determine the direction of the current since, from the perspective of the electron energy scale, electrons move 'downhill', and holes move 'uphill' in response to electrochemical potential gradients. Moreover, the electron and hole currents, J n and J p , can be expressed in terms of the product of the electron and hole QFL gradients with the corresponding carrier conductivities, σ n and σ p : Here, the conductivities are the product of the electronic carrier mobilities µ n and µ p with their corresponding concentrations and the elementary charge:

Open circuit voltage
The open circuit voltage, V OC is the maximum energy per unit charge that can be extracted from an electrochemical cell for a given charge state at open circuit. The V OC can be calculated using the difference in the electron QFL at the location of the cathode (x cathode ) and the hole QFL at the location of the anode (x anode ) with the cell at open circuit.
The open circuit voltage can be output using the command: Here, sol OC is an open circuit solution obtained either by applying V app = V OC or approximated by setting the external series resistance, R S to a high value (e.g. 1 MΩ cm 2 ) using the lightonRs protocol (see Section 4.4 for further information on protocols).

Poisson's equation
Poisson's equation (Gauss's Law) relates the electrostatic potential to the space charge density ρ and the material dielectric constant ε r via the Divergence Theorem. The space charge density is the sum of the mobile carrier and static charge densities at each spatial location. Doping is simulated via the inclusion of fixed charge density terms for ionising donor N D and acceptor N A atoms. In the default version of Driftfusion mobile ionic carriers are modelled as Schottky defects [33] for which every ion has an oppositely charged counterpart, maintaining overall ionic defect charge neutrality within the device. 5 The mobile cation density c is initially balanced by a uniform static counter-ion density N cat and the mobile anion density a is similarly balanced by a static density N ani . For the one-dimensional system described, Poisson's equation can be explicitly stated as: where ε 0 is permittivity of free space. We emphasise that p, n, c, and a represent mobile species, while N A , N D , N cat and N ani are static ion densities. z c , and z a are the integer charge states for the ionic species (by default z c = 1, and z a = −1).
Terms can easily be added or removed from Poisson's equation by editing the S potential term in the Equation Editor in dfpde subfunction of the core df code. See Subsection 4.5 and Listing 1 for further details.
The space charge density rho can be output from a Driftfusion solution structure sol using the command: rho is output as a two dimensional matrix for which the dimensions are [time, space].

Charge transport: Drift and diffusion
As the name suggests, the drift-diffusion (Poisson-Nernst-Planck) model assumes that charge transport within semiconductors is driven by two processes: 1. Drift arising from the Lorentz force on charges due to an electric field F, where F = −dV /dx. 2. Diffusion arising from the entropic drive for carriers to move from regions of high to low concentration.

Bulk transport
Within the bulk of material layers the expressions for the flux density of electrons j n , holes j p , anions j a , and cations j c with mobility µ y and diffusion coefficient D y (where y denotes a generic charge carrier) are given by: Figure S.1 illustrates how the direction of electron and hole flux densities is determined from gradients in the electric potential and charge carrier densities. An analogous diagram can be drawn for mobile ionic species by substituting cations for holes and anions for electrons. The carrier (particle) currents are calculated as the product of the flux densities with the specific carrier charge qz y such that J y = qz y j y .
The electric field calculated from the gradient of the potential (FV) and by integrating the space-charge density a (Frho) can be obtained from a Driftfusion solution structure sol using the syntax:

Diffusion enhancement
The implementation of electronic carrier statistics beyond the Boltzmann approximation (Section 3.2.1) necessitates the inclusion of a generalised Einstein relation to define the relationship between the carrier mobilities and diffusion coefficients as a function of band state occupancy.
[41] The result is a non-linear diffusion enhancement as the QFLs approach and move into the bands. Under Blakemore's approximation,[40] the diffusion coefficient-mobility relationships for electrons and holes can be expressed using the closed-forms: .
We use similar expressions to those above for the ionic carriers from a model proposed by Kilic et al.
[34] to account for steric effects at high ion densities: Here a max and c max denote the limiting anion and cation densities. In the first instance these are set to the lattice cite density for the corresponding ions.

Transport across heterojunctions
At the interface between two different semiconductor materials there is a change in the band energies and electronic density of states. In Driftfusion we choose to model the mixing of states at the interface using a smooth transition in material properties over a discrete interlayer region, in contrast to the commonly employed abrupt interface model 6 (see Figure 4 below for a schematic illustrating the difference between the two models). To accommodate this approach, Equations 18 and 19 are modified to include additional gradient terms for spatial changes in Φ EA , Φ IP , N CB , and N VB . This leads to an adapted set of flux equations for electrons and holes within the interfaces: Values of between 1 − 2 nm have been extensively tested for the interfacial region thickness and are used in the example parameter files accompanying Driftfusion. By default, Φ EA and Φ IP are graded linearly, while N CB and N VB are graded exponentially within the interfacial regions.
The transport equations of Driftfusion can be edited using the carrier flux terms F n, F p, F c, and F a of the Equation Editor in the dfpde subfunction of the core df code. See Subsection 4.5 for further details.

Displacement current
The displacement current J disp , as established in the Maxwell-Ampere law, is the rate of change of the electric displacement field, ∂ D/∂t. In terms of the electric field the displacement current can be expressed as:

Total current
The total current, J is the sum of the individual current components at each point in space and time: Fluxes and currents are calculated from the Driftfusion solution structure sol using the command: [J , j , xout ] = dfana . calcJ ( sol , " sub ") J is a structure containing the individual carrier particle currents J.n, J.p, J.c, and J.a, the displacement current J.disp, and the total current J.tot at each spatial location and time calculated by integrating the continuity equations. j is a structure containing the corresponding carrier and total fluxes. The second input argument "sub" indicates that the currents and fluxes are requested on the subinterval spatial mesh, which is output as xout (see Subsection 4.2.3).

Validity criteria for the drift-diffusion equations
The drift-diffusion approach set out above is valid for semiconductor materials that satisfy the following criteria: [28] 1. The electron and hole populations are at quasi-thermal equilibrium. 2. The electron and hole population temperatures are the same as that of the lattice. 3. Changes in state occupancy are more likely to be due to scattering collisions within a band than generation and recombination events between bands or trapping events. 4. The electron and hole states can be described by a quantum number, k. 5. The mean free path length of carriers,L is significantly shorter than the layer thickness, d (L << d).

Charge continuity
The continuity equations are a set of 'book-keeping' equations, based on the conservation of charge, describing how charge carrier densities change in time at each location. In one-dimension, the continuity equation for a generic carrier density y with flux density j y , and source/sink term S y can be expressed as: For electronic carriers S is composed of two components; 1. Generation, g of carriers by both thermal and photo excitation and; 2. Recombination, r of carriers through radiative (photon emission) and non-radiative pathways. determined by the generation, recombination, and difference between the incoming and outgoing flux density of carriers.
Where chemical reactions take place within devices, additional generation and recombination terms for carriers may also contribute to S. In the current version of Driftfusion, mobile ionic charge carriers are treated as inert such that g c = g a = r c = r a = 0. Users are, however, free to edit the default source terms using the Equation Editor (Section 4.5). A guide describing how to do this is included in the Supplemental Information Section S. 6.
In one-dimension the continuity equations for electrons, holes, cations and anions are given by: Equation 17 and Equations 31 -34 then form the complete set of equations to be solved.

Steady-state approximation to electronic carrier densities and fluxes within the interfacial regions
To better understand the discrete interface model employed by Driftfusion we solve the electron and hole continuity equations (Equations 26, 27, 31 and 32) to obtain analytical expressions for the electronic carrier densities within the discrete interfacial regions using the following approximations and assumptions: 1. Carriers within an interface are at steady-state with respect to the surroundings layers (dn/dt = 0, d p/dt = 0). 2. There is no optical generation within the interface (g = 0). 3. The electric field can be treated as approximately constant throughout the interfacial region (dF/dx = 0). 4. The recombination rate, r within the interfacial region is constant and distributed uniformly.

QFLs remain within the Boltzmann regime (E CB
As detailed in the Supplemental Information Section S.3.1, using the boundary conditions n(x n = 0) = n s , p(x p = 0) = p s , j n (x n = 0) = j n,s , and j p (x p = 0) = j p,s (see Figure 4b), the following expressions can be obtained for the carrier densities within the interfacial regions: Carrier density (log scale) n(x n ) = n s e αx n + j n,s where, The corresponding fluxes are given by: As illustrated in Figure 4b, the translated co-ordinates x n and x p are taken to be in the direction for which α and β are negative and typically the direction for which carrier densities decay.
Example solutions comparing the analytical approximations to numerical solutions calculated using Driftfusion under different transport and recombination regimes are given in the Supplemental Information, Section S.3.2. Where transport is a limiting factor within the interfaces the solutions become strongly dependent on the boundary flux and recombination rates. It is noteworthy however that in the limiting case of infinitely fast transport (µ n,p → ∞) Equations 35 and 36 converge towards purely exponential forms for which the carrier densities change by a Boltzmann factor (∆ n = N CB e αd int and ∆ p = N VB e β d int ) across the width of the interface, d int . For the special case where F = 0, the result is a change in carrier densities equivalent to that expected from an abrupt interface model using Boltzmann statistics. The results presented in this section are applied below in Section 3.7.3 to the interfacial volumetric surface recombination model.

Electronic carrier generation
Two optical models for electronic carrier generation are currently available for use in Driftfusion; uniform generation for which a uniform volumetric generation rate, g 0 is defined for each layer (excluding interfacial regions) and; Beer-Lambert law generation as described below. Irrespective of the choice of optical model the generation rate is zeroed within the interfacial regions to avoid potential stability issues.

Beer-Lambert law generation
The Beer-Lambert law models the photon flux density as falling exponentially within a material with a characteristic photon energy-dependent absorption coefficient α abs . The volumetric generation rate g, over a range of photon energies E γ with incident photon flux density ϕ 0 , is given by the integral across the spectrum: where κ is the reflectance. For simplicity, we assume that a single electron-hole pair is generated by a single photon.

Arbitrary generation profiles
An arbitrary generation profile can be inserted following creation of the parameters object for users who wish to use profiles calculated from different models using an external software package. Details on how to do this are given in Section 4.2.7.

Recombination
By default, two established models for recombination are included in Driftfusion: band-to-band recombination and trapmediated Shockley-Read-Hall (SRH) recombination. Figure  5 is a simplified energy level schematic illustrating these mechanisms. The recombination expressions can be modified in the source terms of the Equation Editor (Section 4.5).

Band-to-band recombination
The rate of band-to-band recombination r btb (also commonly termed direct, radiative or bimolecular recombination) is proportional to the product of the electron and hole densities at a given location such that: where B is the band-to-band recombination rate coefficient. The n 2 i term is equivalent to including an expression for thermal generation and ensures that np n 2 i at steady-state.

Shockley-Read-Hall (SRH) recombination
Recombination via trap states is modelled using a simplified Shockley-Read-Hall (SRH) recombination[37] expression r SRH for which the capture cross section, mean thermal velocity of carriers, and trap density are collected into SRH time constants, τ n,SRH and τ p,SRH for electrons and holes respectively: Here, n t and p t are parameters that define the dependence of the recombination rate to the trap level and are given by the electron and hole densities when their respective QFLs are at the position of the trap energy, E t : It should be noted that Equation 43 is valid only when trapped carriers are in thermal equilibrium with those in the bands. It follows that the rate of trapping and de-trapping of carriers is assumed to be fast compared to the timescale being simulated, such that the approximation is reasonable. In the current version of Driftfusion we also assume that the quantity of trapped carriers is negligible compared to that of the free carriers such that trapped carriers can be neglected in Poisson's equation. To more completely simulate the dynamics of capture and emission of carriers, and the associated contribution to the chemical capacitance of devices, one or more additional variables could be included.
The volumetric recombination rate can be obtained from a Driftfusion solution structure using the command: r is a structure containing three matrices r.btb, r. srh, r.vsr, and r.tot in which the band-to-band, SRH, volumetric surface recombination (VSR, see below) and the total recombination rates are stored as two dimensional matrices with dimensions [ time, space].
The recombination models used in the simulation can be edited using the carrier source terms S electron, S hole, S cation, and S anion in the Equation Editor in dfpde subfunction of the core df code. Note that the models used in dfana must also be updated in accordance with any changes to dfpde as these functions are not coupled. See Subsection 4.5 for further details.

Surface recombination at interfaces
Abrupt interface models typically use a SRH surface recombination model to determine the recombination flux, R int between majority carriers n s and p s at the interface between two materials (see Figure 4a): [20] R int (t) = n s (t)p s (t) − n 2 i 1 s n (p s (t) + p t ) + 1 s p (n s (t) + n t ) .
Here, s n and s p are the surface recombination velocities for electrons and holes at the interface. This model implies that the electron and hole populations in the two materials have delocalised wave functions that overlap significantly such that recombination events are probable.
Since Driftfusion uses discrete interfacial regions, in order to obtain an equivalent recombination flux to the abrupt interface model, we convert Equation 46 into a volumetric surface recombination rate, r vsr by distributing the recombination uniformly across a zone of thickness, d vsr within the interface (see Figure S.4), such that r vsr = R int /d vsr . By default the recombination zone is automatically located next to the interface with the highest minority carrier density at equilibrium. To obtain an expression for r vsr , Equations 35 and 36 can be rearranged to express n s and p s in terms of n(x n ) and p(x p ): For sufficiently high values of µ n , and µ p the n(x n ) and p(x p ) terms dominate Equations 47 and 48 and the carrier density profiles within the interfaces tend towards purely exponential functions, such that n s ≈ n(x n )e −αx n and p s ≈ p(x p )e −β x p . In many instances it is then sufficient to approximate the volumetric surface recombination rate within the recombination zone as: where α and β are given in Equations 37 and 38, and the d vsr term is subsumed into the volumetric surface recombination time constants, τ n,vsr and τ p,vsr such that: We stress here that this is not a physically motivated model in the sense that we do not anticipate recombination to be distributed uniformly throughout a recombination zone in reality. This approach does however result in a good approximation to the established abrupt interface surface recombination model for a wide variety of devices and conditions (see Section 5.4).
The volumetric surface recombination model can be toggled on and off by using the property par. vsr mode. Where par.vsr mode = 1, E t and the n i , n t and p t terms are set to constant and calculated from the energy levels defined for the interfacial region. This ensures that r vsr remains approximately constant throughout the recombination zone. Where par.vsr mode = 0, the standard SRH expression in Equation 43 is assumed. In this case E t is graded linearly and n i , n t and p t are graded exponentially.
The assumptions used in the derivation of Equation 49 breakdown when the transport within the interface is limited or where recombination fluxes are particularly high (see Section S.3). Since both the flux and recombination terms in Equations 47 and 48 are unknowns without welldefined limits, Driftfusion performs a check for selfconsistency directly following calculation of the solution when VSR mode is switched on: the function compare rec flux calculates the sum of the interfacial recombination fluxes using the values of n s and p s from the solution and the SRH model given in Equation 46. This sum is compared to that of the integrated recombination rate calculated using the VSR model (Equation 49) across all interfaces. If the fractional difference in the two calculations is greater than par.RelTol vsr, for fluxes above par.AbsTol vsr, a warning is displayed. In such cases users could consider increasing the electronic carrier mobilities within the interfacial regions or reducing the recombination coefficients.

Initial conditions
At present two sets of initial conditions are used in Driftfusion, dependent on the number of layers. These conditions are designed to be consistent with the boundary conditions and to minimise the error in the space charge density at junctions which can lead to large electric fields and convergence failure.

Single layer device
A linearly varying electrostatic potential and exponentially varying electronic carrier densities over the layer thickness d are used as the initial conditions (Equations 53, 54, and 52) when simulating a single layer. Uniform ionic carrier density profiles are used throughout the layer to guarantee ionic defect charge neutrality (Equations 56 and 55).
n(x) = n 0,l exp ln n 0,r n 0,l Here, the built-in potential V bi of the device is determined by the difference in boundary electrode workfunctions Φ l and Φ r :

Multilayer device
For multilayer devices the electrostatic potential is set to fall uniformly throughout the device (Equation 58), while the electronic carrier densities are chosen to be the equilibrium densities for the individual layers (n 0 and p 0 ). As with the single layers, the ionic carriers are given a uniform density (Equations 59 -61), thus guaranteeing local electroneutrality.
Here, the device thickness d dev is the sum of the individual layer thicknesses d i (d dev = ∑ i d i ). Driftfusion auto-detects the number of layers in the device and uses the appropriate set of initial conditions when running the equilibrate protocol to obtain the equilibrium solutions for the device (Section 4.3).
The initial conditions of the simulation can be edited in the dfic subfunction of the core df code. See Subsection 4.5 for further details.

Boundary conditions
Solving Equation 17 and Equations 31 -34 requires two constants of integration for each variable, which are provided by the system boundary conditions. For the charge carriers, Neumann (defined-flux value) conditions are used to set the flux density into and out of the system. The electrostatic potential uses Dirichlet conditions (defined-variable value) such that the potential is fixed at both boundaries at each point in time as detailed in Section 3.9.1. The details of these boundary conditions are discussed in the following subsections.

Electrostatic potential boundary conditions
In Driftfusion the electrostatic potential at the left-hand boundary is set to zero (Equation 63) and used as the reference potential. The applied electrical bias, V app and an effective potential arising from series resistance V Rs are applied to the right-hand boundary as described in Equation 64 .
Here, Ohm's law is used to calculate V Rs from the electron and hole flux densities: where R s is the area-normalised series resistance, given by the product of the external series resistance and the device active area. Setting R s to a high value (e.g. R s = 10 6 Ω cm 2 ) approximates an open circuit condition for devices with metal electrodes. Technically this can be achieved using the lighton Rs protocol (see Section 4.4 for a description of protocols). 7 The boundary conditions of the simulation can be edited in the dfbc subfunction of the core df code. See Subsection 4.5 for further details.
The function generator fun gen defines the applied potential Vapp as a function of time t, which can be recalculated from a Driftfusion solution structure sol using the command:

Carrier selectivity and surface recombination at the system boundaries
Many architectures of semiconductor device, including solar cells and LEDs employ selective contact layers that block minority carriers from being extracted (or injected) via energetic barriers. These are known variously as transport layers, blocking layers, blocking contacts, or selective contacts. For solar cells semiconductor layers are typically sandwiched between two metallic electrodes constituted of metals or highly-doped semiconductors. Such materials can be numerically challenging to simulate owing to their high charge carrier densities and thin depletion widths. Consequently, a common approach is to use boundary conditions defining charge carrier extraction and recombination flux densities to simulate the properties of either the contact or electrode material. It should be noted, however, that the employment of fixed electrostatic potential boundary conditions (as defined in Section 3.9.1) implies that the potential falls only within the discrete system and not within the electrodes. This approximation is only realistic for contact materials with vanishingly small depletion widths (infinite interfacial capacitances) i.e. metals and highly doped semiconductors. It follows that to accurately simulate semiconductor contacts layers with finite depletion regions these layers must also be included within the discretised system. For electronic carriers the surface recombination velocity coefficients s n and s p determine the carrier extraction/recombination rate at the boundaries of the system. For majority carriers in solar cells, high values of s n and s p (e.g. > 10 7 cm s −1 [38]) are advantageous for carrier extraction, while low values imply poor contact extraction properties. For minority carriers, high values of s n and s p are typically undesirable as they imply high rates of surface recombination at the electrode. In Driftfusion the expressions for electronic boundary carrier flux densities, j n and j p are given by the typical first-order expressions: j n,r (t) = s n,r (n r (t) − n 0,r ), where n 0,l , n 0,r , p 0,l , and p 0,r are the equilibrium carrier densities at the left (x = 0) and right-hand (x = d) boundaries, calculated using Equations 8 and 9 under the assumption that the semiconductor QFLs are at the same energy as the electrode Fermi energy, which is assumed to remain constant. For the left-hand boundary n 0,l and p 0,l are given by Equations 70 and 71.
where Φ l is the left-hand electrode work function. Analogous expressions are used for n r and p r at the right-hand boundary. Extraction barriers can also be modelled with this approach by including a term for the barrier energy in the exponent of Equations 70 and 71. At present, however, quantum mechanical tunnelling and image charge density models for energetic barriers at the system boundaries are not accounted for in Driftfusion.

Ionic carrier boundary conditions
In the simplest case, ionic carriers are confined to the device and do not react at the electrode boundaries. This leads to a set of zero flux density boundary conditions for mobile anions and cations: Where an infinite reservoir of ions exists at a system boundary (such as an electrolyte), a Dirichlet boundary condition defining a constant ion density could alternatively be imposed.
This concludes our description of the physical models employed in Driftfusion. In the following section the system architecture and key commands are introduced as well as a guide on how to get started with using Driftfusion.

System architecture and how to use Driftfusion
Driftfusion is designed such that the user performs a linear sequence of simple procedures to obtain a solution. The key steps are summarised in Figure 6; Following initialisation of the system, the user defines a device by creating a parameters object containing all the individual layer and device-wide properties; The equilibrium solutions (soleq. el and soleq.ion) for the device are then obtained before applying a voltage and light protocol, which may involve intermediate solutions; Once a desired solution (sol) has been obtained, analysis and plotting functions can be called to calculate outputs and visualise the solutions. Below, the principal functions are discussed in further detail.

Initialising the system: initialise df
At the start of each MATLAB session, initialise df needs to be called from within the Driftfusion parent folder (not one of its subfolders) to add the program folders to the file path and set plotting defaults. This action must be completed before any saved data objects are loaded into the MATLAB workspace to ensure that objects are associated with their corresponding class definitions.

Defining device properties and creating a parameters object: pc(file path)
The parameters class, pc contains the default device properties and functions required to build a parameters object, which we shall denote herein as par. The parameters object defines both layer-specific and device-wide properties.
Layer-specific properties must be a cell or numerical array containing the same number of elements as there are layers, including interface layers. For example, a three-layer device with two heterojunctions requires layer-specific property arrays to have five elements. Examples can be found in the /Input files folder (see Section 4.2.1). Figure 7 shows the processes through which the main components of the parameters object are built; The user is required to define a set of material properties for each layer. The comments in pc describe each of the parameters in detail and give their units (also see Supplemental Information, Table S.2 for quick reference); Several subfunctions (methods) of pc then calculate dependent properties, such as the equilibrium carrier densities for example, from the choice of probability distribution function (prob distro function) and other user-defined properties. The treatment of properties in the device interfaces is dealt with, and can be changed, using the device builder build device (see below).

Importing properties
Typically, the most important user-definable properties are stored in a .csv file, which is easily editable with a spreadsheet editor such as LibreOffice or Microsoft Excel. The file path to the .csv file can then be used as an input argument for pc, for example: This allows the user to easily create and store sets of key device parameters without editing pc. Default values for properties set in the parameters class pc will be overwritten by the values in the .csv file during creation of the parameters object par. New properties defined in pc can easily be  added to the .csv file provided that they are also included in import properties, which tests to see which properties are present in the text file and reads them into the parameters object where present. Following the properties read-in step performed by import properties, the number of rows in the layer type column (see Section 4.2.2) is used for error checking all other entries to confirm that properties have been defined for each discrete layer of the system (this does not include the electrode rows, which are pseudo-layers).
To avoid potential incompatibility, new user-defined material properties defined in pc, which have distinct values for each layer, should be included in the .csv file and added to the list of importable properties in import properties, as well as build device. An example of how to do this is given in the Supplemental Information, Section S.6.

Layer types
Layer types, set using the layer type property, flag how each layer should be treated. Driftfusion currently uses three layer types: 1. 'electrode': A pseudo-layer which defines the boundary properties of the system. These are not discrete layers and do not appear in visualised outputs. 2. 'layer': A slab of semiconductor for which all properties are spatially constant.
3. 'active': As layer but flags the active layer of the device. The number of the first layer designated 'active' is stored in the active layer property and is used for calculating further properties such as the active layer thickness, d active. Flagging the active layer proves particularly useful when automating parameter explorations as the active layer can easily be indexed. 4. 'interface': An interfacial region between two different material layers. The properties of the interface are varied according to the specific choice of grading method as defined in build device (see below). It is critical that interfacial layers are included between material layers with different energy levels and eDOS values i.e. at heterojunctions. See the included default input files for examples of how to set up devices with heterojunctions.

Spatial mesh
The computational grid is divided into N intervals with N − 1 subintervals, where the position of the subintervals is defined by pdepe solves for the variable values u i+1/2 on the subintervals (x i+1/2 ) and their associated flux densities j i on the integer intervals (x i ) as illustrated in Figure 8. Table 1 Key to properties contained in external .csv parameters files and their constraints. 1 Note that contrary to convention, Phi EA and Phi IP take negative values for consistency with other energies referenced to the electron energy scale. 2 These properties are required for electrode pseudo-layers, but all other properties are ignored in these rows. 3 For electrode layers, entries for EF0 are stored in the parameters object (par) as the distinct properties Phi left and Phi right rather than as part of the EF0 array. Phi left and Phi right take negative values for consistency with other energies referenced to the electron energy scale. 4 For electrode layers, entries for sn and sp are stored in the parameters object (par) as the distinct properties sn l, sn r, sp l, and sp r rather than as part of the sn and sp arrays.
Uniform Owing to the use of a finite element discretisation scheme the details of the spatial mesh in Driftfusion are of critical importance to ensure fast and reliable convergence. In this release of Driftfusion two types of spatial mesh are available: 1. 'linear': Linear piece-wise spacing 2. 'erf-linear': Mixed error function (bulk regions)linear (interfacial regions) piece-wise spacing meshgen x generates integer and subinterval spatial meshes, x and x sub respectively (see Figure 8), based on the layer thickness, the number of layer points defined in the device properties, and the xmesh type. The solution is interpolated for the integer grid points when generating the output solution matrix sol.u (see Section 4.6). The property xmesh coeff controls the spread of points for regions where an error function is used for point spacing i.e. layer and active layer types: higher values determine higher point densities close to the layer boundaries. In general we recommend using 'erf-linear' for devices with relatively high ionic defect densities as, where ionic carriers are confined, high point densities are required at the layer boundaries to resolve the carrier distributions. Where the depletion of ionic carriers extends into the bulk xmesh ceoff can be reduced to increase the bulk point density.

Time mesh
pdepe uses an adaptive time step for forward time integration and solution output is interpolated for the user-defined time mesh. Convergence of the solver is weakly-dependent of the user-defined time mesh interval spacing and stronglydependent on the maximum time and the maximum time step. These can be adjusted by changing the tmax and MaxStepFactor properties of the parameters object (see Section 4.2). The options for different time mesh types (tmesh type) are: 1. 'linear' or 1 2. 'log 10' or 2 Typically, the time mesh is changed frequently for intermediate solutions within protocols to accommodate the different timescales on which carriers move. For example, in some cases the ionic carriers may be frozen to obtain a stable short timescale solution for the electronic carriers, before a second, longer timescale, solution is calculated with the ionic carriers mobile. Similarly the tmesh type is adjusted dependent on the voltage and light conditions. For example during a J − V scan a linear time mesh is used in keeping with the linear change in the applied voltage with time. In other instances a logarithmic mesh is more appropriate in order to resolve time periods for which the carrier time derivatives are larger.

The device structures
build device and build property are called during creation of the parameters object to build two important data structures, which we call the 'device structures': dev, defined on the integer grid intervals (x i ), is used to determine the initial conditions only, while dev sub, defined on the subintervals (x i+ 1 2 ), is used by the pdepe solver function.
dev and dev sub contain arrays defining all spatially varying properties at every location within the device including the interfacial regions. build property enables the user to specify different types of interface grading for each property listed in build device. At present there are four generic grading option types: 1. 'zeroed': The value of the property is set to zero throughout the interfacial region. 2. 'constant': The value of the property is set to be constant throughout the interfacial region. Note that the user must define this value in the .csv file. 3. 'lin graded': The property is linearly graded using the property values of the adjoining layers. 4. 'exp graded': The property is exponentially graded using the property values of the adjoining layers.
For the volumetric surface recombination scheme described in Section 3.7.3, we also introduce parameterspecific grading schemes for the VSR time constants: The interface grading type for each property is set in the build device function e.g. the default grading type for the electron affinity, Φ EA is 'lin graded'.
Dependent on the choice of grading option, input values may, or may not, be required for a given material property in the interfacial layers. For example, when using the 'lin graded' option a value is not needed because the interface property values are calculated from those of the adjacent layers. By contrast, when using the 'constant' option, a property value does need to be specified for the interfacial layer to avoid an error. Since property values that are not required are ignored, we recommend that users specify all property values for all layers to future-proof against problems arising when experimenting with grading options.

The electronic carrier probability distribution function
The class distro fun defines the electronic carrier probability distribution function for the model and calculates equilibrium boundary, and initial carrier densities when building the parameters object. At present there are two available options for the choice of probability distribution function: While the Boltzmann approximation results in marginally faster calculations owing to the absence of a diffusion enhancement we recommend using Blakemore statistics for their extended domain of validity.

The electronic carrier generation profile function
generation calculates the two generation profiles gx1 and gx2, which can be used for a constant bias light and a pulse source for example, at each spatial location in the device for the chosen optical model and light sources (see Section 3.6 and the flow diagram in the Supplemental Information Figure S.5). The light sources can be set using the light source1 and light source2 properties. There are two options for the optical model: 1. 'uniform': a uniform volumetric generation rate defined by the property g0 is applied layers with layer and active layer types. The multiplier properties int1 and int2 define the intensities for light source 1 and 2 respectively. 2. 'Beer Lambert': The generation profile follows the Beer-Lambert law as detailed in Section 3.6.1 with material layer optical properties taken from ./ Libraries/Index of Refraction library.xls for the corresponding materials defined in the par. material cell array. As with uniform generation, the generation rate profile is multiplied by the intensity properties int1 and int2 before being applied.
A arbitrary generation profile calculated using an external program (such as Solcore [42] or the McGeeHee Group's Transfer Matrix code[43]) can be inserted into the parameters object par by overwriting the generation profile properties gx1 or gx2 following creation of the object. The profile must be interpolated for the subinterval (x i+1/2 ) grid points defined by x sub (see Section 4.2.3). We also recommend that the generation rate is set to zero within the interfacial regions to avoid stability issues.
The light source time-dependencies are controlled using the g1 fun type and g2 fun type properties, which define the function type (e.g. sine wave) and the g1 fun arg and g2 fun arg properties, which are coefficient arrays for the function generator (see Section 4.5.5) e.g. the frequency, amplitude, etc.

Protocols: equilibrate
Once a device has been created and stored in the MATLAB workspace as a parameters object the next step is to find the equilibrium solution for the device. The function equilibrate starts with the initial conditions described in Section 3.8 and runs through a number of steps to find the equilibrium solutions, with and without mobile ionic carriers, for the device described by par. The output structure soleq contains two solutions (see Section 4.6);  2. soleq.ion: electronic and ionic carriers are mobile and are at equilibrium.
Storing both solutions in this way allows devices with and without mobile ionic charge to be compared easily.

Protocols: General
A Driftfusion protocol is defined as a function that contains a series of instructions that takes an input solution (initial conditions) and produces an output solution. For many of the existing protocols, listed in Table 2, with the exception of equilibrate, the input is one of the equilibrium solutions, soleq.el or soleq.ion. Figure 9 is a flow diagram illustrating the key functions called during execution of a protocol. Protocols typically start by creating a temporary parameters object that is a duplicate of the input solution parameters object. This temporary object can then be used to write new voltage and light parameters that will be used by the function generator to define the generation rate at each point in space and time and the potential at the boundary at each point in time. Additional parameters, for example those defining the output time mesh or carrier transport are also frequently adjusted. The approach is to split a complex experimental protocol into a series of intermediate steps that facilitate convergence of the solver. For example, where ionic mobilities are separated by many orders of magnitude from electronic mobilities, a steady-state solution is easiest found by accelerating the ions to similar values to the electronic carriers using the K a and K c properties. The simulation can then be run with an appropriate time step and checked to confirm that a steady-state has been reached. The possibilities are too numerous to list here and users are encouraged to investigate the existing protocols listed in Table 2 in preparation of writing their own.

The Driftfusion master function df
df is the core of Driftfusion. The function takes the device parameters, and voltage and light conditions, calls the solver and outputs the solution (Figure 9). The underlying equations in Driftfusion can be customised for specific models and applications using the Equation Editor in the dfpde subfunction. df contains three important subfunctions: dfpde, dfic, and dfbc which we now proceed to describe.

Driftfusion Partial Differential Equation function dfpde and the Equation Editor
dfpde defines the equations to be solved by MATLAB's Partial Differential Equation: Parabolic and Elliptic solver toolbox (pdepe).
[31] For one-dimensional Cartesian co-ordinates, pdepe solves equations of the form: Here, u is a vector containing the variables V , n, p, c, and a at each position in space x and time t. C is a vector defining the prefactor for the time derivative, F is a vector determining the flux density terms (note: F does not denote the electric field in this section), and S is a vector containing the source/sink terms for the components of u. By default C = 1 for charge carriers but could, for example, be used to change the active volume fraction of a layer in a mesoporous structure. The equations can be easily reviewed and edited in the dfpde Equation Editor as shown in Listing 1. Here, device properties that have a spatial dependence are indexed with the variable i to obtain the corresponding value at the location given by x sub(i). A step-by-step example of how to change the physical model using the Equation Editor is given in the Supplemental Information, Section S.6.

Driftfusion Initial Conditions dfic
dfic defines the initial conditions to be used by pedpe. If running equilibrate to obtain the equilibrium solution or running df with an empty input solution, the first set of initial conditions are as described in Section 3.8. Otherwise, the final time point of the input solution is used.

Driftfusion Boundary Conditions dfbc
dfbc defines the system boundary conditions. The boundary condition expressions are passed to pdepe using two coefficients P and Q, with N u elements, where N u is the number of independent variables being solved for. The boundary conditions are expressed in the form: For Dirichlet conditions P must be non-zero to define the variable values, whereas for Neumann conditions Q must be non-zero to define the flux density. Listing 2 shows how the default Driftfusion boundary conditions (described in Section 3.9) are implemented in the current version of the code.
In addition to these subfunctions, df also calls two important external functions: the time mesh generator, meshgen t and the function generator, fun gen.

Calculating outputs: dfana
dfana is a collection of functions (methods) that enable the user to calculate outputs such as carrier currents, quasi-Fermi levels, recombination rates etc. from the solution matrix u, the parameters object par and the specified physical models. The use of a class in this instance enables the package syntax dfana.my calculation(sol) to be used. For example the command: Due to the computational cost of calling functions external to pdepe, the physical model described in the Equation Editor is not coupled to that used in the analysis functions. Users should, therefore, take great care when adapting the physics of the simulation to make certain that the models defined in dfana and df are consistent.  Figure 15a). If no second argument is given then only the final time point is plotted. dfplot also includes the generic property plotting function dfplot.x2d to allow users to easily create new two-dimensional plots.
For variables plotted as a function of time the second argument defines the position. For example the command: dfplot . Jt ( sol , 1e -5) ; plots the current density for each carrier as a function of position at x = 10 −5 cm. For plots where variables are integrated over a region of space, the second argument is a vector containing the limits [x 1 , x 2 ]. Further details can be found in the comments of dfplot.

Getting started and the example scripts
While the underlying system may appear complex, Driftfusion has been designed such that with a few simple commands, users can simulate complex devices and transient optoelectronic experiment protocols. Table 2 is a complete list of protocols available at the time of writing. In addition to the brief guide below, a quick start with up-todate instructions can be found in the README.md file in the Driftfusion GitHub repository,[26] and a series of example scripts for running specific protocols are also presented in the Scripts folder. New users are advised to study these scripts and adapt them to their own purposes. In addition, we have written an introductory workshop to guide students through the process of building basic semiconductor devices and applying optical and voltage biases to them. This can be found in the Getting-started-workshop-2021 branch of the Driftfusion GitHub repository.

How to build a device object, find the equilibrium solution and run a cyclic voltammogram
In this section some commonly-used commands are put together to show new users how to create a device object, obtain the device equilibrium solutions, and run a protocol, which in this example simulates a cyclic voltammogram (CV). The doCV protocol applies a triangular wave voltage function to the device, with optional constant illumination, for a set number of cycles enabling the device current-voltage characteristics at a given scan rate to be calculated.
At the start of each session, the system must be initialised by typing the command: initialise_df ; To create a parameters object using the default parameters for a Spiro-OMeTAD/perovskite/TiO 2 perovskite solar cell the parameters class, pc is called with the file path to the appropriate .csv file as the input argument: par = pc ( ' Input_files / spiro_mapi_tio2 . csv ') ; The equilibrium solutions with and without mobile ionic carriers for the device can now be obtained by calling the equilibrate protocol: soleq = equilibrate ( par ) ; As discussed in Section 4.3, the output structure soleq contains two solutions: soleq.el and soleq.ion. In this example we are interested in seeing how mobile ionic carriers influence the device currents so we will use the solution including mobile ionic charge carriers, soleq.ion.
To perform a cyclic voltammogram simulation from 0 to 1.2 to −0.2 to 0 V at 50 mVs −1 , under 1 sun illumination we call the doCV protocol with the appropriate argument values as detailed in the protocol comments shown in Listing 4. The second argument of dfplot.JVapp is the precalculated dependent property d midactive, the value of which is equal to the position at the midpoint of the active layer of the device. The resulting plot is given in Figure 10 for reference.

How to change the physical model
A detailed step-by-step example of how to modify the physical model to account for the possible effects of photogenerated mobile ionic charge carriers is described in Section S.6 of the Supplemental Information. In some situations where device properties, such as layer widths, are changed in a user-defined script or function users may need to rebuild the device structures dev and dev sub, and spatial meshes x and x sub. To maintain code performance this is not performed automatically since non-devicerelated parameters are frequently changed within protocols. meshgen x and build device can be rerun and stored using the new parameter set ('refreshed') using the syntax: par = r efresh_device ( par ) ; To further illustrate how to use refresh device, refresh device script, a script describing the nec-essary steps to change the interfacial recombination parameters, is provided in the Scripts folder of the Driftfusion repository.

Parallel computing and parameter exploration: explore
Driftfusion's parameter exploration class explore, takes advantage of MATLAB's parallel computing toolbox, to enable multiple simulations to be calculated using a parallel pool. explore script is an example script demonstrating how to use explore to run an active layer thickness versus light intensity parameter exploration and plot the outputs using explore's embedded plotting tools.

Validation against existing models
To verify the numerical accuracy of the simulation, results from Driftfusion were compared against those from two analytical and two numerical models. In Section 5.1 current-voltage characteristics obtained using analytical and numerical solutions for a p-n junction solar cell are compared. In Section 5.2 the simulation's time integration is verified by calculating the transient photovoltage response of a single, field-free layer and comparing it to the solution obtained using a zero-dimensional kinetic model. In Section 5.3, numerical solutions for three-layer, dual heterojunction devices obtained using Driftfusion are compared with those from the Advanced Semiconductor Analysis (ASA) simulation tool, an established, commercially available package. [35]. Finally, in Section 5.4, J-V characteristics calculated using Driftfusion for devices dominated by bulk and interfacial recombination processes are compared with those of IonMonger [20], a recently published, free-touse, mixed ionic-electronic carrier semiconductor device simulator.
The location of the MATLAB scripts and the parameter sets used to obtain the results in this section can be found in the Supplemental Information, Section S.9.

The depletion approximation for a p-n junction
The p-n junction depletion approximation The Depletion Approximation (DA) allows the continuity equations and Poisson's Equation (Equations 31, 32 and 17) to be solved analytically for a p-n homojunction.
[45] The depletion region at the junction of the device is assumed to have zero free carriers such that the space charge density ρ can be described using a step function with magnitude equal to the background doping density (see Figure 11, top panel). Transport and recombination of free carriers in the depletion region is also ignored. Poisson's equation can then be solved by applying fixed carrier density (p(x = −∞) = p 0 and n(x = ∞) = n 0 ) and zero-field boundary conditions to obtain the depletion widths for n-and p-type regions, w n and w p respectively:[27] Solving the DA for the current flowing across the junction under the assumption that the diffusion length of both carriers is significantly greater than the device thickness (L n,p >> d) yields the Shockley diode equation: where J SC and J 0 are the short circuit and dark saturation current densities respectively. Here we use the convention that a positive applied (forward) bias generates a positive current flowing across the junction.
To make the comparison between numerical and analytical current-voltage (J-V ) characteristics, values for J 0 and J SC need to be related to input parameters for the simulation. J 0 embodies the recombination characteristics of the device. For the contribution to recombination in the quasi-neutral region, J 0 can be related to the electron and hole diffusion lengths L n and L p , minority carrier lifetimes τ n and τ p , and diffusion coefficients D n and D p , according to: [27] J 0 = qD p p 0,n−type L p + qD n n 0,p−type L n where L p = τ p D p and L n = √ τ n D n . The material band gap and AM1.5 solar spectrum were used to calculate the theoretical maximum current density J SC,max and corresponding uniform generation rate throughout the depletion region g 0 . Figure  The limiting short circuit photocurrent for a perfectly absorbing semiconductor of band gap E g is given by: where η is the external quantum efficiency. If η = 1 for photon energies E γ ≥ E g , and η = 0 for E γ < E g , the maximum theoretically achievable short circuit current for a single junction J SC,max is given by the integral of φ 0 from the bandgap energy to infinity. Figure S.8 of the Supplemental Information shows the maximum achievable current density J SC,max over a range of band gap energies. For the comparison we use the bandgap of silicon E g = 1.12 eV resulting in J SC,max = 42.7 mA cm −2 .
Simulation methods A p-n junction was devised in Driftfusion with very thick n-and p-type layers (≈ 100 µm, p 0 ≈ N A = 9.47 × 10 15 cm −3 , n 0 ≈ N D = 2.01 × 10 15 cm −3 ) to approximate the assumptions and boundary conditions used in the DA. Furthermore, ionic carriers and trapped electronic charges are not included in the simulations in this section. The surface recombination velocity was set to zero for minority carriers at both boundaries. A special recombination scheme was implemented using the following simplified first order expressions: where τ n and τ p are the electron and hole lifetimes respectively. For simplicity, τ n and τ p were set equal to one another and, for consistency with the DA, recombination was switched off in the depletion region.
To convert the value obtained for J SC,max into a uniform carrier generation rate, the short circuit flux density ( j SC,max = J SC,max /q) was divided by the depletion region thickness d DR , yielding g 0 = j SC,max /d DR . The complete parameter sets for the simulations in this section are given in Tables S.5 and S.6.
Results Both the analytical and numerical solutions for the space charge density, electric field, and electric potential are shown in Figure 11: The space charge widths, field strength and potential profiles all show good agreement. Figure 12 shows the light and dark current-voltage curves for the analytical solution obtained using equation 80  for three values of τ n,p . For τ n,p = 10 −6 and τ n,p = 10 −7 the agreement is very good, even at current densities as low as 10 −14 A cm −2 . The solutions begin to diverge at τ n,p = 10 −8 s, for which L n,p = 2.3 µm such that L n,p << d and the underlying assumptions of the DA are no longer valid. The deviation from the ideal model in Driftfusion at higher current densities (Figure 12, inset) is expected due to the absence of series resistance in the DA.  produce a small additional photovoltage ∆V . As detailed in Supplemental Information Section S.9.2.1, it can be shown using a kinetic model that ∆V is given by:

Transient photovoltage response of a single layer field-free device
for − t pulse < t 0, (85) where t pulse is the length of the laser pulse, ∆ g is the additional volumetric generation rate due to the excitation pulse, n OC is the steady-state open circuit carrier density, and k TPV is the decay rate constant of the TPV signal.
Simulation methods A 100 nm field-free single layer of semiconductor with E g = 1.6 eV was simulated in Driftfusion with the zero flux density boundary conditions for electronic carriers representing perfect blocking contacts. The constant volumetric generation rate was set to g 0 = 1.89 × 10 21 cm −3 s −1 = 1 sun equivalent based on the integrated photon flux density for the AM1.5G solar spectrum and a step function absorption. In this special case the splitting of the quasi-Fermi levels in the simulation is solely attributable to changes in chemical potential, such that no series resistance term is needed for the potential boundary conditions. It should be noted that for instances where an electric field is present in the device these boundary conditions will not represent an open circuit condition; a high value for R s or a mirrored cell approach (see Ref. [18] for details) should be used instead.
The second order band-to-band recombination coefficient was set to B = 10 −10 cm 3 s −1 and SRH recombination was switched off. Since the underlying kinetic theory demands that the TPV perturbation must be small such that the additional carrier density ∆ n << n OC , we used t pulse = 1 µs and set the pulse intensity equivalent to be 20% of the bias light intensity. The complete parameter sets for the simulations are given in Tables S.7 and S.8.
Results A comparison of the results from the analytical and simulation models for bias light intensities of 0.1, 1, and 10 sun equivalent is shown in Figure 13. The steady-state charge carrier densities and open circuit voltage obtained using Driftfusion agreed to within 10 decimal places with the analytical values calculated using Equations S.15 and S.14 (n = 4.35 × 10 15 cm −3 , V OC = 1.08 V at 1 sun equivalent). The transient photovoltage perturbations also behaved as predicted, with the rate constants extracted from fitting the TPV decays correct to within 3 significant figures of the values calculated using the analytical expression in Equation 86.

Numerical solution for three-layer devices with electronic carriers
Simulation methods To verify that the graded treatment of the interfaces in Driftfusion produces similar results to models that use abrupt interfaces, a HTL/absorber/ETL device was simulated using both Driftfusion and the Advanced Semiconductor Analysis (ASA) simulation tool. [35] To maintain consistent layer dimensions the linear grid spacing for the ASA simulations and the interface thickness in Driftfusion were set to be 1 nm. The base material parameters for the active layer were based loosely on those for a perovskite material excluding mobile ionic charge. The parameters for the contact layers were not chosen to simulate real materials but rather to produce a large built-in potential and to vary as many properties as possible including the layer thickness, dielectric constants, recombination coefficients, mobilities etc. For optical generation the Beer-Lambert option (without back contact reflection) was chosen and the same optical constant and incident photon flux density spectrum data were used in both simulators. Four different parameter sets (PS) were compared, the key differences for which are summarised in the first four columns of Table 3.
The complete parameter sets for the simulated devices are given in Tables S.9 -S.12.
Results: Beer-Lambert optical model The results for the integrated generation rate profiles are shown in Figure S.10. Despite the wavelength-dependent generation rates appearing to be very closely matched between the two simulators ( Figure S.11), the integration across photon energies resulted in marginally different generation profiles in the two simulations corresponding to a total difference in generation current of 1.24 mA cm −2 . As a means to ensure that the input generation profiles were identical, the generation profile from ASA was inserted into the Driftfusion parameters objects using the method described in Section 4.2.7.
Results: Current-voltage characteristics To compare the current outputs from the different simulation tools currentvoltage scans were performed from V app = 0 to 1.3 V. Since ASA solves for the steady-state current, the J-V scan rate was set to be low (k scan = 10 −10 V s −1 ) in Driftfusion to minimise contributions from the displacement current. Figure 14a shows a comparison of the J-V characteristics obtained from the two simulation tools for Parameter Sets (PS) 1a and 2a. While the different parameter sets result in distinctly different characteristics for the two devices, the data show excellent agreement between the two simulators despite the significant difference in discretisation schemes and interface treatment. Closer examination of the results from the two simulators ( Figure 14b) reveals that, the percentage difference (calculated as 100 × (J ASA − J DF )/J ASA ) for PS 1a is on the order of 1 % for current densities beyond J = 10 −12 A cm −2 .
Halving the active layer thickness has little impact on this difference (PS 1b). The percentage difference calculated for PS 2a, however, is much greater, with a maximum of ≈ 5 % for current densities J > 10 −12 A cm −2 . The larger difference between the two simulators in this instance can be attributed to a change of over 7 orders of magnitude in the electron density at the absorber-ETL interface (Figure S.14) due to both a transition in the conduction band eDOS from N CB = 10 18 to 10 20 cm −3 and a change in the conduction band energy of 0.3 eV. Under these circumstances the difference in discretisation schemes between the two tools becomes apparent: the linear discretisation method used in the PDEPE and Driftfusion cannot calculate the change in carrier densities within the interfaces to as high a degree of accuracy as the internal boundary conditions used in ASA (see Section S.8.1 for further details). The deviation can be reduced significantly, however, by using a uniform eDOS of N CB = 10 18 cm −3 across all layers as the results for PS 2b show. With respect to device characteristics, despite the percentage difference in results for PS 2a, the key metrics of the J-V curve such as the V OC , ideality factor and fill-factor are all preserved. Figure 15 shows a comparison of the electrostatic potential and hole density profiles calculated using Driftfusion and ASA for PS 1a under illumination at increasing applied bias. The electron density is given in the Supplemental Information, Figure S  ) they have no effect in these regions. For the comparison we ran J-V scan simulations at a variety of scan speeds using similar parameter sets to those published in Ref. [8]: a device dominated by bulk recombination and a device dominated by interfacial recombination using the volumetric surface recombination scheme described in Section 3.7.3.   is somewhat sacrificed in Driftfusion in favour of greater flexibility, which enables the physical models to be easily edited, devices with any number of material layers to be simulated and a range of interface-specific properties and grading functions to be specified.

Results
We have verified Driftfusion against two analytical and two existing numerical models and found that in all cases the calculated results are in reasonably good agreement. The results show that the discrete interface approach produces results comparable to models with abrupt interfaces, albeit with marginal errors introduced owing to the combination of the interface treatment and the linear discretisation scheme used in Driftfusion.

Conclusions
We have developed an efficient and powerful onedimensional drift-diffusion simulation tool for modelling semiconductor devices with mixed ionic-electronic conducting layers based on MATLAB's pdepe toolbox. Distinct from existing codes, in addition to electronic carriers, Driftfusion can include up to two ionic carrier species and virtually any number of material layers. This flexibility is made possible by a discrete interlayer interface approach whereby the material properties can be graded between two adjoining semiconductor layers. This method has the added advantage that interface-specific properties can easily be specified.
The default physical models underlying the simulation were described and analytical approximations to the carrier densities and fluxes within the interfacial regions were given. Using these solutions, a method for approximating surface recombination within interfacial regions using a volumetric recombination scheme was derived.
The system architecture was presented and the processes by which users can define device properties and change the simulation's physical model were outlined. Protocol functions, determining time-dependent voltage and light conditions, were described as well as the simulation solution structures and how to use built-in analysis and plotting functions to calculate and visualise outputs. A step-by-step guide was presented detailing how to create a device, find an equilibrium solution, and calculate current-voltage characteristics using a cyclic voltammogram protocol. Lastly, advanced features, including how to calculate multiple solutions in parallel, were briefly introduced.
Driftfusion was verified by comparing calculated solutions with two analytical and two existing numerical models which tested different aspects of the simulation. In all cases the agreement was good and the general device behaviour was reproduced. The discrete treatment of the interfaces resulted in small variations in the calculated currents as compared to other simulators using abrupt layer boundaries.
The ease with which the underlying models can be changed and new material layers can be introduced places Driftfusion in a unique space compared to other free-to-use simulation codes for which an intricate knowledge of the numerical mechanics is required to adapt the physical models and device architecture. By making Driftfusion both accessible and free-to-use our hope is that this work will advance our collective understanding of mixed ionic-electronic conducting materials and devices.
Acknowledgements The authors would like to thank the UK Engineering and Physical Sciences Research Council for funding this work (Grant No. EP/J002305/1, No. EP/M025020/1, No. EP/M014797/1, No. EP/N020863/1, No. EP/R020574/1, No. EP/R023581/ 1, and No. EP/L016702/1). J.N. thanks the European Research Council for support under the European Union's Horizon 2020 research and innovation program (grant agreement No. 742708). We would also like to thank: Emilio Palomares and Davide Moia for lending the expertise of their students to work on this project, and Diego AlonsoÁlvarez for his invaluable advice on preparing and improving the code.

Funding
The authors would like to thank the UK Engineering and Physical Sciences Research Council (grant No. EP/J002305/1, EP/M025020/1, EP/M014797/1, EP/R020574/1, EP/R023581/1, and EP/L016702/1) and the European Union's Horizon 2020 research and innovation program grant agreement (grant No. 742708) for their support in funding this work.

Conflicts of interest
The authors declare that they have no conflicts of interest.

Availability of data and material
The data presented in this work is available at reasonable request from the authors.

Supplemental information for Driftfusion
An open source code for simulating ordered semiconductor devices with mixed ionic-electronic conducting materials in one dimension Philip Calado 1 · Ilario Gelmetti 2 · Benjamin Hilton 1 · Mohammed Azzouzi 1 · Jenny Nelson 1 · Piers R. F. Barnes 1  Table S.2 Table of symbols. †The absorption and photon flux data are loaded directly from libraries according to the layer names defined in stack and the choice of light source using the properties light source1 and light source2. * Solely for the Methods-depletion-approximation branch. * * For the .csv file the name thickness can also be used. For a complete description of all user-defined and dependent properties please see the comments in the parameters class pc. For information on calculated outputs see the analysis class dfana. VSR denotes auto-generated coefficients used in the volumetric surface recombination scheme.

Driftfusion
Variable Name Variable type Name Unit Vres Solver variable Voltage drop across external circuit V

S.3.1 Solution
To examine how the carrier densities change within the interfacial regions we begin with the continuity equation for electrons and assume no generation such that: For simplicity, we will assume that electrons within the interfaces are in a steady-state condition such that ∂ n(x,t)/∂t = 0 and where: where, To solve Equation S.2 we make the following approximations and assumptions; 1. Transport within the interface is fast with respect to the surrounding materials such that the carriers within the interface can be treated as begin at steady-state. 2. Given that the interfacial regions are largely depleted of carriers, the electric field can be treated as constant throughout the interfaces. 3. The recombination rate r is uniform throughout the thickness of the interface such that it can be treated as constant for a given majority carrier boundary density. 4. The non-generalised Einstein relation (D = k B T µ) can be applied (i.e. Boltzmann statistics are valid).
Under these conditions the following general solution to Equation S.2 can be obtained: C I , C II are constants, and x is the translated spatial co-ordinate such that x = x − x 1 , where x 1 is the position of the left-hand interface boundary as shown in Figure 4 of the Main Text. Using the boundary conditions n(0) = n s , and flux j n (0) = j n,s the following expressions for the electron density and flux within the interfaces can be obtained: j n (x ) = j n,s − rx .

(S.8)
A similar method can be used to obtain analogous expressions for the hole carrier and flux densities: where,

S.3.2 Validation under different transport and recombination regimes
Repository branch: Methods-interface-solutions Parameters file: ./Input files/3 layer test vary.csv Script: ./Scripts/interface numerical analytical To validate the analytical solutions derived in Section S.3.1, here we present results comparing analytical and numerical solutions of the carrier densities within the interfacial regions for a three-layer device including ionic carriers in the active layer. Solely for testing purposes a constant and uniform recombination rate, r con was set throughout the interfacial regions. Four different mobility and recombination regimes were tested; 1. High mobility, no recombination; 2. Low mobility, no recombination; 3. High mobility, high recombination; and 4. Low mobility, high recombination. Current-voltage scans from 0 to 1 to 0 V at a scan rate of 100 mV s −1 were executed in each case. To calculate the analytical values the boundary values n s , p s , j n,s and j p,s were taken from the numerical solution (the nearest interior point on the subinterval mesh) for each interface. The solutions in Figure S The results from the analytical and numerical solutions show close agreement. The small differences between the two can be accounted for by the fact that the fluxes are calculated on the subinterval mesh. In the low mobility, high recombination case, negative carrier densities result from the unphysical constant recombination rate. In practice this can be avoided by using relatively high mobilities within the interfacial regions such that the r term in Equations S.7 and S.9 is minimised. As discussed in Section 3.7.3, the use of a recombination zone further minimises this issue by including recombination only in the region with the highest minority carrier density. Position, x (nm)

S.6.2 Adding new device properties
In order to be able to adjust our generation term we will introduce two new device properties k iongen and k ionrec that characterise the ionic generation and recombination rates. For simplicity we will assume that these terms are material specific and equal to zero in the transport layers of our three-layer device. k iongen and k ionrec need to be added as arrays to the parameters class pc and the device builder before we can create a parameters object for the device. First we create the properties in pc with default values as shown in Listing S1.
Listing S2 shows how new properties can then be added to the device builder. The new coefficient arrays k iongen and k ionrec are used as the first input arguments. The device building code build device then uses these values with our chosen interface grading option to define the values of these properties at every point in the device and stores the resulting arrays in the device structures par.dev and par.dev sub. Here the 'lin graded' option is chosen to linearly grade the new properties in the interface regions.
Once the new properties have been added to the pc and build device, they can additionally be added to the .csv parameters file and import properties to enable them to be easily imported. Figure  New properties can be inserted anywhere the user chooses into the table as import properties checks for matches between the column headers and listed property names. Lastly, the new properties are added to import properties using try statements (Listing 3). If a match to a column heading is found in the .csv file the the properties are imported, otherwise the default values in pc are used.
This completes the steps necessary to add new properties to the parameters object and make them easily accessible from the .csv file. We will now look at how to adapt the physical model using the Equation Editor to include our new ionic carrier generation and recombination terms.

S.6.3 Using the equation editor
For simplicity, we will assume that ionic carriers are generated at a rate proportional to the optical intensity and hence the electron and hole generation rate g. We will further assume that the ionic carriers recombine at a rate proportional to their respective densities. We stress here that there is no physical basis for this choice of models and they are only been used for illustrative purposes. To adapt the Driftfusion master code df to our purposes we first need to unpack the new variables from the parameters object by adding the to the list of variables at the start of the code as shown in Listing S4.
While this approach is laborious, creating individual variables in the workspace in this These terms are then simply added to the appropriate lines in the Equation Editor as shown in Listing S5.

S.6.4 Adapting analysis functions in dfana to the new physical model
It is particularly important to bear in mind that the physical models in the master code and the analysis functions are not coupled. Users therefore need to adapt any of the analysis functions in dfana to account for changes to the physical model. For example, if one were to change the recombination model, dfana.calcr would need to be updated in accordance with the changes made to the expressions in the Equation Editor.

S.7 Troubleshooting
We hope that your experience using Driftfusion will be a relatively trouble-free experience. Along with the possibility of minor bugs that come with experimental research software, you may encounter two common error messages: spatial discretization and time integration failures.

S.7.1 Spatial discretization has failed
This error typically occurs for one of two reasons: 1. The point density of spatial mesh is too high: Try reducing the number of points in individual layers or using a different spatial mesh type (see meshgen x for possible options). 2. The boundary conditions are not consistent with initial conditions: The initial value of a variable must be consistent with those set in the boundary conditions. If the initial or boundary conditions are changed from the default, ensure that they are consistent with one another.

S.7.3 Unexpected values calculated using dfana
As discussed in the Main Text, owing to the computational cost of using functions external to df for the solving the equations, the physical model described in the Equation Editor is not coupled to that used in the analysis functions. For greater consistency an obvious change to the programming architecture would be to use centralised flux and source functions to describe the physical model, which could then be called during both the solving and analysis stages. Multiple tests have shown, however, that this method comes with a significant increase in computation cost owing to the large number of calls made by the solver to such functions. Hence the present method of adapting the model in both df and dfana separately, while more prone to user error, is considerably more efficient.

S.7.4 Bug reporting
If you find a bug with Driftfusion please raise an issue using the Issues tab on the Driftfusion GitHub repository page. [1] This is the best way to ensure that other users can see which bugs have been addressed and understand how they have been fixed.

S.8.1 Linear discretisation
Many existing drift-diffusion models use Scharfetter-Gummel finite volume discretisation scheme [2][3] to account for the exponential change in carrier density in a constant electric field. In order to take advantage of MATLAB's Partial Differential Equation Parabolic and Elliptic (PDEPE) toolbox and the ease with which transport models can be altered with it, a major trade off is the use a simplified finite element discretisation scheme for which carrier densities are assumed to change linearly between neighbouring grid points. This can be compensated for somewhat with the use of a logarithmically-spaced spatial grid and a higher density of grid points. There are however limits to the allowed point density and while the comparison simulations that use the Scharfetter-Gummel scheme (Section 5.3 of the Main Text) is generally good, the user should be aware that the use of lower point densities will increase this error. One consequence of this simplified scheme is that the currents may not always be close to zero at equilibrium. Following exploratory work at lower point densities, users are encouraged to use as high a point density as possible for obtaining the final results, especially since the solver is highly efficient. Communication with Mathworks has suggested the maximum point density is a function of various factors and therefore cannot be distilled to a single value. A degree of trial and error is therefore required on the part of the user to find the maximum allowed point density for a given device configuration.

S.8.2 Numerical errors with interfacial currents
Following from the brief discussion in Subsection S.8.1, and as noted in the comparison of Driftfusion and ASA results in Section 5.3, the linear discretisation method and large gradients in the band energies and eDOS within the interfacial interlayers of Driftfusion makes these regions particularly prone to numerical errors. As shown in Section S.9.4, errors can be minimised by reducing band energy and eDOS gradients as well as increasing the number of spatial mesh points within the interfaces. However, where accurate calculation of very low current densities (< 10 −12 A cm −2 for example, based on the ASA comparison calculations) across one or more heterojunction interfaces is required we presently do not recommend the use of Driftfusion for this purpose.

S.8.3 Normalisation of the dielectric constant
As may be noted from the Equation Editor (Listing 1, Main Text) the dielectric constant must be normalised for the electrostatic potential flux and source terms to avoid a spatial discretisation error. The origin of this problem is currently unknown as normalisation has little effect on the magnitude of the input values.

S.8.4 High extraction coefficients do not result in constant charge density at the boundaries
Under some uncommon operating conditions e.g. large preconditioning voltages, the carrier densities at the boundaries do not tend to their limiting value when using high extraction coefficients. Under most circumstances this can be solved by switching to fixed carrier densities at the boundaries. This does however present a problem when calculating currents for a single carrier device as the minority carrier fluxes can no longer be used as the boundary value for the fluxes. This issue is currently ongoing and under investigation. S.8.5 Drift and diffusion currents do not sum to give the total current Drift and diffusion currents calculated using dfana.Jddxt do not sum to give the correct current. This is related to the way in which fluxes are calculated in the solver but a solution has yet to present itself. Total carrier currents are correctly calculated using the continuity equations in dfana.calcJ and users are recommended to use this method instead, reserving dfana.Jddxt only for individual analysis of the approximate drift and diffusion currents.

S.8.6 lightOnRs protocol instability
At the time of writing lightOnRs protocol does not converge for all devices parameter sets at high values of series resistance.

S.9 Validation against existing models
The following sections provide references to the relevant Github Driftfusion repository branches, parameters files and scripts required to reproduce the results in Section 5 of the Main Text.
S.9.1 The depletion approximation for a p-n junction Repository branch: Methods-depletion-approximation-comparison Parameters file: ./Input files/TPV test.csv Script: df methods depletion approx pn junction.m This section has been largely reproduced from reference [5]. Under the assumption that the parabolic band and Boltzmann approximations are valid, 2 the charge densities of electrons n and holes p for an intrinsic semiconductor with intrinsic carrier density n i and Fermi energy E Fi are given by: [6] Using Equations S.12 and S. 13, the open circuit (OC) voltage V OC in a field free, zerodimensional material, where electron and hole charge carrier concentrations are equal (n = p), can be expressed as: Here, n OC is carrier density at OC and is dependent on the generation rate g and the recombination model. For example, in an idealised device with band-to-band recombination only, the recombination rate r is given by r = B(n 2 − n 2 i ), where B is the band-to-band recombination coefficient. At OC steady-state, generation is equal to recombination (r = g) leading to: During a small perturbation measurement a small additional charge density ∆ n << n OC is injected into the device such that the state of the system is not significantly altered. In a transient photovoltage measurement this charge is generated by an excitation light pulse imposed onto a background bias light ( Figure S.9a). In the most general case, following a small perturbation ending at t = 0, the rate of change of addition charge d∆ n dt can be expressed using a small perturbation rate constant, k TPV (see reference [7] , pp. 101 − 103 for full derivation): where i and j are summation indices and the exponents γ i and γ j can take any value (including non-integers). k i j is the associated rate coefficient for the reaction order. The solution to Equation S.16 is an exponential: for t > 0 (S.18) Using Equation S.14, the change in open circuit voltage ∆V OC produced by a light pulse introducing an additional charge ∆ n can be expressed as: Substituting for ∆ n using Equation S.18 leads to: A similar method can be used to find the change in charge carrier density and voltage rise during a pulse of duration t pulse : where ∆ g is the additional generation rate from the pulse light.
In experimental measurements, a single empirical reaction order γ, with the corresponding rate constant k γ , is typically assumed to dominate recombination such that: Accordingly, the slope of the log(k TPV ) vs. log(n OC ) plot can be used to determine γ:     Table S.14 Interface-specific parameters for simulation comparisons with IonMonger. (a) and (b) denote the bulk and interfacial recombination dominated devices respectively.* The relative dielectric constant is set artificially high in this parameter set to reduce the electric field strength within the interfacial regions: this modification facilitates a more direct comparison with the abrupt interface model.         The actual recombination fluxes calculated in Driftfusion, using the volumetric scheme described in Section 3.7.3 are labelled as DF volumetric (dashed coloured curves). The recombination fluxes calculated using IonMonger [8] are also back-calculated using the expression in Ref [9], where additional factors are included with the n 2 i , which have been neglected here. The volumetric interfacial recombination scheme in Driftfusion shows a high level of self-consistency with the 2-D model. Some variation is observed between Driftfusion and IonMonger, which can be attributed to the difference in electron density at the AL/HTL interface owing to small differences in the ionic carrier distributions shown in Figure Figure 16 of the Main Text. The accuracy of the volumetric recombination scheme can be improved by doubling the interface thickness to 2 nm (thus reducing α and β ) and increasing the number of points in the interfaces from 60 to 100 (Coloured, dotted curves).