Bayesian parameter estimation for characterising mobile ion vacancies in perovskite solar cells

To overcome the challenges associated with poor temporal stability of perovskite solar cells, methods are required that allow for fast iteration of fabrication and characterisation, such that optimal device performance and stability may be actively pursued. Currently, establishing the causes of underperformance is both complex and time-consuming, and optimisation of device fabrication thus inherently slow. Here, we present a means of computational device characterisation of mobile halide ion parameters from room temperature current-voltage (J-V) measurements only, requiring $\sim 2$ hours of computation on basic computing resources. With our approach, the physical parameters of the device may be reverse modelled from experimental J-V measurements. In a drift-diffusion model, the set of coupled drift-diffusion partial differential equations cannot be inverted explicitly, so a method for inverting the drift-diffusion simulation is required. We show how Bayesian Parameter Estimation (BPE) coupled with a drift-diffusion perovskite solar cell model can determine the extent to which device parameters affect performance measured by J-V characteristics. Our method is demonstrated by investigating the extent to which device performance is influenced by mobile halide ions for a specific fabricated device. The ion vacancy density $N_0$ and diffusion coefficient $D_I$ were found to be precisely characterised for both simulated and fabricated devices. This result opens up the possibility of pinpointing origins of degradation by finding which parameters most influence device J-V curves as the cell degrades.


Introduction
Lead-halide perovskites (LHPs) have attracted intense research activity as light harvesting layers in solar cell devices due to power conversion efficiency (PCE) increases of 9-25.8% since 2009 [1], rivalling the record efficiency reported for silicon-based solar cells of 26.7% [2].Additionally, perovskite solar cells (PSCs) may be fabricated using low-energy, low-cost manufacturing techniques such as solution processing [3].Despite the progress to large power conversion efficiencies, problems with the reproducibility and stability of PSCs remains [4,5].The decrease of PCE over time originates in extrinsic factors such as exposure to heat [6], light [7], humidity [8] and oxygen [9] and intrinsic factors such as the structural stability of the perovskite [10] and migration of mobile ion vacancies across the perovskite layer [11].
To find strategies and fabrication techniques in addressing these problems, a large number of parameters is required to specify layer widths and materials parameters for the perovskite layer and the electron and hole transport layers (ETL, HTL respectively).These parameters are often not accurately known and may vary from device to device.Identifying which of these parameters is sensitive to the fabrication method and extent of degradation of a specific fabricated device gives important insights on the origins of subpar device performance.The PSCs can be characterised by experimental measurements of, for example, current-voltage (J-V) characteristics at varying voltage scan rates.These characteristics do not however give directly the device and materials parameter values.Here we present a method for reverse modelling, namely mapping from J-V device measurements to device and materials parameters.We do so by solving the inverse parameter problem where a drift-diffusion (DD) simulation is combined with Bayesian parameter estimation, BPE.It is not possible to invert explicitly the coupled partial differential equations in the DD model.
Despite strong evidence for the migration of halide ions, and an understanding of its effect on device performance, it is not clear how to determine whether performance losses in a specific device may be attributed to ion migration.For example, a recent study has been made of degradation pathways linked to iodide diffusion and iodide reactions.[12] The electric field in the perovskite layer results from the built in voltage of the device V bi and the externally applied voltage V app , such that the electric field across the perovskite E = (V bi − V app )/b p where b p is the width of the perovskite layer, once potential drops across the charge transport layers are allowed for.The electric field causes the mobile ion vacancies to migrate across the perovskite layer, resulting in a build up of negatively charged ions at the boundary between the perovskite and ETL and positively charged ion vacancies at the HTL (see Figure 1).In turn, the resulting non-uniform charge distribution partially screens the total electric field in the perovskite layer -severely affecting charge transport, reducing the power output of the device and causing current-voltage hysteresis [13][14][15].The effect of halide ion migration can be detrimental to cell performance, causing significant reduction in PCE over time on the order of seconds to days [16].
As shown by density functional theory (DFT) calculations [11,17], halide ions (e.g.iodide I − ) are the most likely candidate to undergo migration through vacancy pathways due to the ions' large diffusion coefficients in comparison to other ionic species.Migration of halide ions within the perovskite layer has also been confirmed experimentally [18,19].J-V characteristics such as hysteresis are strongly dependent on mobile ion vacancy parameters, the mobile ion vacancy density N 0 and mobile ion diffusion coefficient D I [20,21].Here, we demonstrate our reverse modelling approach by mapping from J-V device measurement to ion vacancy parameters.Our method characterises mobile ion vacancy parameters (N 0 and D I ) from room temperature J-V device measurements only, without further experiment.Knowledge of N 0 and D I for a specific device can then be used to minimise the effect of ion migration on the performance and stability of PSCs, such as optimising experimental processing [22,23], tailoring perovskite composition [24,25], and tuning transport layer properties [20,26], can be tested.Without appropriate characterisation of mobile ion vacancy parameters for a specific device, it is difficult to analyse how the reduction in performance, potentially caused by mobile ion vacancies, depends on processing and fabrication methods.This limits the rate at which experimental optimisation can be achieved and restricts progress towards improving the temporal stability of PSCs.
In our approach, BPE [27][28][29] was used to iteratively query a drift-diffusion PSC model to obtain a posterior distribution over the ion vacancy parameters most likely to reproduce experimental J-V curves.BPE has previously been used in solar cell applications to characterise the charge carrier mobility, minority carrier lifetime and conduction and valence band offsets between the light-harvesting and transport layers [30].Kirchartz et al. demonstrated that BPE may be applied to PSCs and also utilised a neural network model to approximate the drift-diffusion simulations [31].However, the effect of ion migration present in PSCs is not included in reference [31].Our choice of BPE method requires only single temperature J-V measurements at two voltage scan-rates, in comparison to the multiple temperature [30] and illumination [31] measurements used in previous studies.We model both charge carriers and mobile ion vacancies, and allow all parameters associated with the device (perovskite and transport layer parameters, see Table 1) to influence performance.Additionally, we determine posterior distributions over all fundamental device parameters, allowing for a greater understanding of the underlying physical processes affecting performance.

Methods
Solving the inverse problem of mapping from J-V device measurement to device parameters requires a model of the forward relationship (device parameters to J-V measurement) which may then be inverted.Here, this forward relationship has been approximated by utilising a coupled electron-hole-ion drift diffusion device model, namely IonMonger [32,33].IonMonger is capable of simulating the current across the device as a function of the experimental conditions (illumination, applied bias, voltage scan-rate) and fundamental physical parameters (see Table 1).We propose a Bayesian Parameter Estimation (BPE) approach to solve the inverse problem; the advantages and disadvantages of this method are discussed in section 3. BPE applies Bayes' Rule to inversely derive a posterior distribution p(θ|y) over a set of parameters θ associated with a measured outcome y, p(y|θ) is the likelihood of obtaining the measured outcome y from θ, p(θ) is the prior distribution over θ, and θ p(y|θ)p(θ) is a normalising factor.In this application, θ is the set of fundamental physical device parameters associated with a PSC device (see Table 1) and y is the set of characteristics obtained from a measured J-V curve; specifically y = [J sc , V r oc , V f oc , P r max , P f max ], where J sc is the short-circuit current density, V oc is the open-circuit voltage, and the superscripts r and f indicate quantities measured on the reverse and forward scans, respectively.
The prior distribution over perovskite parameters p(θ) was assumed to be a uniform distribution over each parameter, representing the inherent epistemic uncertainty prior to any measurement of the device.Table 1 lists an estimated parameter range within each prior distribution for the TiO 2 -MAPbI 3 -Spiro PSC architecture.The likelihood p(y|θ) of obtaining a measured set of J-V characteristics y from a set of parameters θ was determined by querying IonMonger.Let the prediction of y by IonMonger be ỹ for given θ.The likelihood that θ produces the observed characteristics y can therefore be modelled as a multivariate Gaussian distribution with mean y, where the covariance is determined by the set of experimental uncertainties in each measured J-V characteristic σ, assumed to be ±0.1%.The dominating uncertainty originates from measurement of the current density, as shown by Brandt et al. [30] to be of the order ±0.1%, and σ was therefore set to this value.Note that the covariance matrix is diagonal, quantifying an assumption that the measured J-V characteristics for a specific device are independently distributed.The posterior distribution p(θ|y) was obtained by sampling from the product of the likelihood and the prior p(y|θ)p(θ) using the Metropolis-Hastings algorithm [34,35], a Markov-Chain Monte-Carlo (MCMC) sampling method [36,37].Thus, a distribution over the parameters of the device likely to produce the experimental J-V measurement may be determined, without an explicit model of this inverse dependence.
It was also found possible to decrease the variance of final posterior distributions over mobile ion vacancy parameters by including J-V characteristics measured from multiple voltage scan-rates.The set of characteristics y used in the method above was therefore the concatenation of two sets of characteristics for scan-rates of 0.1V/s and 1.0V/s.Previous studies [21] have shown that hysteresis is maximised at scanrates of 0.1V/s and 1.0V/s and were therefore selected for this method to maximise the information on mobile vacancy dynamics within the likelihood.An analysis of this observation is provided in section 3 and describes how the method may be augmented to probe all regions of the mobile ion vacancy parameter space, enabling characterisation of an arbitrary device.

Results and Discussion
To evaluate the method, a TiO 2 -MAPbI 3 -Spiro architecture was simulated by IonMonger and J-V measurements taken at 0.1V/s and 1.0V/s scan-rates.Using the BPE method presented, posterior distributions over mobile ion vacancy parameters, vacancy density N 0 and diffusion coefficient D I , were determined (see Figure 2).Specifically, 6 Markov Chains were iterated and allowed to search the device parameter space for 500 Metropolis-Hasting steps, resulting in 3000 total device parameter samples.This process required ∼ 2.5 hours of computation on a 32 core Intel Xeon 2.40GHz processor (2012), with IonMonger parallelised across all cores.
It can be seen from Figure 2 that the resulting posterior distribution obtained over N 0 and D I is a precise (i.e.non-zero probability density varying over ∼ 1 order of  magnitude) characterisation of the true mobile ion vacancy parameters associated with the observed J-V measurements, as simulated by IonMonger.The resulting distribution allows an understanding of internal mobile ion vacancy dynamics; specifically, whether poor/optimal device performance may be attributed to the number density of mobile ion vacancies (large N 0 ) and/or the vacancy diffusion coefficient D I is at a value likely to affect device performance at a given scan-rate.Additionally, posterior distributions were obtained from J-V characteristics taken at each scan-rate (0.1V/s and 1.0V/s) individually at a temperature of 295K.The variance associated with the individual posteriors is significantly increased in comparison to the final posterior distribution obtained using J-V characteristics at both scan rates to specify the likelihood.This may be understood by considering that the posterior distribution obtained from both scanrates must be determined from the overlap of the posterior distributions from individual scan-rates.The probability that a set of mobile ion parameters is associated with a set of J-V characteristics must now be determined from the probability the parameters are associated with J-V characteristics measured at 0.1V/s and 1.0V/s.This overlap of the two likelihoods at each scan-rate therefore greatly reduces the probable mobile ion parameter space and reduces the uncertainty in characterisation.
Further, observed J-V characteristics exhibit a strong dependence on the scan-rate at which the device was measured [20,21].Importantly, the matching of time-scales between ion diffusive motion and the rate at which the voltage is scanned determines whether J-V hysteresis is observed.Cave et al. [21] showed that by varying the scanrate and locating the rate that causes the largest J-V hysteresis, a prediction of the mobile ion diffusion constant may be obtained.Similarly, here the scan-rate may be varied to 'probe' an arbitrary region of the (N 0 , D I ) parameter space, with resulting J-V characteristics specifying the likelihood of the BPE method.As J-V hysteresis is strongly dependent on the dynamics of mobile ion vacancies, maximising observed J-V hysteresis was found to provide the highest characterisation precision within the BPE method (see Figure 2(c) and (d)).This freedom to increase characterisation precision in any region of mobile ion parameter space by varying a single parameter (scan-rate), allows for characterisation of an arbitrary device with arbitrary values of N 0 and D I .
As shown in Table 1, all physical device parameters were allowed to vary within a uniform prior range.Consequently, characterisation of all device parameters was obtained via BPE and a visualisation of the posterior distributions across all parameters is shown in Figure 3. Figure 3 shows that parameters N 0 and D I are characterised with the least variance; the layer width b, permittivity of perovskite layer ϵ p and absorption coefficient α are also consistently well characterised but with larger uncertainties.The plausible parameter ranges over other physical device parameters have been considerably reduced, with varying degrees of characterisation precision.The method presented has been optimised to characterise N 0 and D I by specifying the likelihood with measured J-V characteristics that are strongly dependent on these parameters (e.g.P f max ).However, it is known that device performance, as characterised by J-V measurements, is also dependent on other physical parameters such as the doping densities of the transport layers [20].It is therefore important that all device parameters were allowed to vary and affect performance (J-V characteristics) to ensure that any under-performance is not wrongly attributed to the effect of ion motion.In principle, the precision of characterisation over any device parameter may be improved by incorporating further device measurements within the likelihood that are known to be strongly dependent on that device parameter.For example, pulse J-V measurements have been shown to allow for characterisation of PSCs while minimising the effect of ion migration [43].By specifying the likelihood of the BPE method with characteristics measured by the pulse J-V technique, the influence of mobile ion vacancy parameters on the measurement can be minimised.Other device parameters may therefore be more precisely characterised as their effect on performance is amplified.Here, the J-V characteristics selected to specify the likelihood lead to high characterisation precision for mobile ion parameters N 0 and D I and lower characterisation precision for transport layer parameters.This indicates that the selected J-V characteristics are less dependent on changes in the transport layer parameters than mobile ion parameters.Indeed, high characterisation precision of N 0 and D I has been obtained despite poorer characterisation of the transport layer parameters.This suggests the potential to de-couple the influence of mobile ion Figure 3: Boxplot over posterior parameter samples obtained via BPE method for each normalised physical device parameter (see Table 1).Lower, median and upper quartiles shown by vertical lines of central box.Whiskers are inclusive of extremal parameter samples and therefore indicate the extent of the Markov-Chain search over prior uniform range for each parameter.Black dots show true device parameters associated with the simulated J-V curves given in Figure 2. parameters on performance by selection of a relevant set of performance characteristics (here, J-V characteristics) strongly dependent on the mobile ion parameters and weakly dependent on other device parameters.Again this reasoning may be extended to the characterisation any device parameter via a BPE method.
The method was tested experimentally with J-V measurements on a TiO 2 -MAPbI 3 -Spiro architecture from Cave et al. [21] where the device was measured at scan-rates of 0.178V/s and 1.0V/s at 295K.A posterior distribution over mobile ion parameters N 0 and D I was obtained using the BPE method and the result can be seen in Figure 4 as well as the two J-V measurements taken by Cave et al. [21] Figure 4 shows how N 0 and D I are deduced to have median values of 4.9 × 10 23 m −3 and 3.5 × 10 −14 m 2 s −1 , respectively.Assuming the ion diffusion coefficient D I has an Arrhenius dependence on temperature T , where the high temperature diffusion The vacancy hop activation energy E a may be predicted from this relationship and was calculated to be 0.35 eV.This prediction lies within experimentally obtained activation energies of 0.36-0.43eV.[60] Cave et al [21] predicted 0.41 eV for E a in a PSC with HTL, ETL of spiro, TiO 2 respectively and the DFT prediction of 0.37eV (0.51eV) for tetragonal (cubic) MAPbI 3 , allowing for a reaction enthalpy of 0.07 eV [61], and noting that the cubic-to-tetragonal phase transition in MAPI occurs around 54 °C, 327K [62].However, individual device characterisation, and this method, assume that mobile ion vacancy parameters vary between devices as a result of fabrication and device history and DFT predictions for E a are approximate.The single DFT prediction for all TiO 2 -MAPbI 3 -Spiro architectures assumes a level of uniformity across devices that oversimplifies the resulting mobile ion dynamics in PSCs.Using the BPE method, the prediction for the activation energy of mobile ion vacancy diffusion is device specific.This therefore allows for an analysis to be constructed, correlating fabrication methods with resulting mobile ion vacancy dynamics.
Additionally, an analysis of the assumptions of the numerical simulation model, IonMonger, leads to further potential insight.Specifically, IonMonger assumes that all ion vacancies are mobile -an assumption which is known not to hold in fabricated PSCs due to, for example, grain boundary dominated migration [63].The resulting posterior over vacancy density N 0 predicted by this BPE method is therefore a quantification of the density of mobile vacancies.The number of total ion vacancies may be larger due to immobile ion vacancies.For example, DFT calculations predict an iodide vacancy density for cubic MAPbI 3 to be 1.6 × 10 25 m −3 at room temperature [11,44]; whereas Figure 4 predicts a iodide vacancy density of 4.9 × 10 23 m −3 .This discrepancy may be due to either a lower number of total ion vacancies resulting from fabrication or a lower number of mobile ion vacancies, with migration heavily limited by the density of grain boundaries.
The efficiency of the BPE method, and most importantly the efficiency of the Metropolis-Hasting algorithm, was observed to be strongly dependent on jumping distribution variance (see Figure 5).Within the Metropolis-Hastings algorithm, each Markov chain proposes a new position in the parameter space by sampling a point in this space from a jumping distribution centered on the current point.Here, the jumping distribution was specified by a multivariate Gaussian distribution with the standard distribution in each parameter given by a fixed percentage, q, of the prior distribution range (see Table 1).For large values of q, the algorithm takes large steps around the parameter space, potentially finding a region of the space with high-likelihood more quickly.Yet, due to the large jumping distribution variance, the acceptance rate of proposed states will be low and the algorithm is unlikely to remain precisely in this region.Conversely, if the value of q is too low the acceptance rate will increase but the algorithm may take too long to find a region of the space with high-likelihood.Optimal scaling theory provides a prediction of the asymptotically optimal acceptance rate of the Metropolis-Hasting algorithm for a multi-dimensional parameter space to be 23.4% [64].The variance of the jumping distribution was therefore tailored to achieve this acceptance rate on average, which was found to be for q ∼ 1%.
Future work may look to improve the efficiency of this method via superior MCMC algorithms such as Hamiltonian Monte Carlo (HMC) [65] or No-U-Turn Samplers (NUTS) [66].These methods utilise the derivative of the likelihood to efficiently search the parameter space, requiring fewer iterations/evaluations of the posterior and therefore, in this application, fewer device model queries.Active learning methods [67], may also be investigated.Active learning methods look to specify the set of device parameters with which to query the drift-diffusion PSC model such that the decrease in entropy of the posterior distribution is maximized, potentially resulting in high query/data efficiency.
Application specific efficiency improvements may also be obtained by constraining the device parameter space with prior device knowledge.In the examples presented here, prior distributions over device parameters were uninformative uniform distributions over the full range of plausible values for MAPbI 3 .Reducing the parameter space over which the method is performed would decrease the number of MCMC iterations required and may be achieved by application-specific informative prior distributions.For example, the prior distribution over perovskite layer width ranged from 400nm-900nm; however, the layer width is commonly characterized previously to J-V measurement.The prior distribution over this parameter may therefore be specified by a normal distribution centered on the measured value with variance given by the measurement uncertainty.This may be extended to any number of device parameters that have previously been characterized, greatly decreasing the searched parameter space.

Conclusions
A method for fast (∼ hours) reverse modelling, mapping from J-V device measurements to device and materials parameters, has been presented.We have demonstrated that mobile halide ion parameters can be deduced from room temperature J-V measurements only, requiring no further experimental protocols.In our method, Bayesian Parameter Estimation is coupled with a drift-diffusion PSC model, to determine posterior distributions over mobile ion parameters which were found to precisely characterise the Figure 5: Example Markov Chains for jumping distribution standard deviations given by q = 1% (left) and q = 4% (right).Red points indicate initial randomly sampled starting point of the Markov Chain.The resulting posterior (q = 1%) is given in Figure 2(e).ion vacancy density N 0 and diffusion coefficient D I for both simulated and fabricated devices.With this method, an understanding of the internal dynamics of mobile ions, a phenomenon known to cause J-V hysteresis and significantly hinder optimal chargetransport, may be obtained for a specific device.This understanding allows fabrication techniques and other potential device optimisations to be correlated with the resulting internal dynamics of mobile ion vacancies.Further, by performing characterisation over all physical device parameters, the influence of mobile ion vacancies on performance may be decoupled from the effect of other physical parameters.The relative influence of physical processes within the device may therefore be established.
The method presented here can also be applied in a digital twin-style setting, where a virtual model of a laboratory device is created and continuously updated to reflect the current output of the device.The digital twin would use our drift diffusion code combined with Bayesian Parameter Estimation to find values for the drift-diffusion model input parameters that reproduce measurements of the PSC being studied, such as its current-voltage (J-V) characteristics.If the device changes over time, e.g.due to degradation, this method would be able to identify the causes of changes in device performance and track these over time.This approach could be of value, not only in current efforts to improve PSC by helping to determine underlying physical processes, but also in the future by employing it in operational settings, where methods for monitoring and optimising devices would be needed.
High-performance computing (HPC) may be leveraged to further reduce the characterisation time, potentially allowing for full device characterisation in minutes, requiring only J-V measurements as input.Further methodological improvements may also be investigated such as improved MCMC sampling and active learning methods.

Figure 1 :
Figure 1: (a) Illustration of halide ion vacancy migration in perovskite solar cells.Mobile halide ions (purple) migrate through ion vacancies (green) to the electron transport layer (ETL) boundary, producing an electric field E ion .This field partially screens the electric field due to V bi and V app .Electrons and holes are indicated by e − and h + , respectively.(b) Simulated example of J-V hysteresis where the current density on the forward voltage-scan (dashed line), differs from the reverse voltage-scan (solid line); arrows indicate voltage-scan direction.The maximum power point (P max ) for both the reverse (superscript r) and forward (superscript f) scan are shown, indicating the operating power losses as a result of hysteresis.

Figure 2 :
Figure 2: Simulated J-V measurements (a) 0.1mV/s, (b) 1.0mV/s and associated posterior distributions shown in (c) and (d), respectively, where the likelihood p(y|θ) has been specified by single scan-rate J-V characteristics.(e) Posterior distribution obtained with the likelihood specified by J-V characteristics obtained from both scan-

Figure 4 :
Figure 4: Experimental J-V measurements at scan rates of (a) 0.178V/s and (b) 1.0V/s.(c) Associated posterior distribution over mobile ion vacancy parameters, N 0 and D I .High posterior density indicated by yellow, fading to dark purple for zero density.