Device Model for Methylammonium Lead Iodide Perovskite With Experimentally Validated Ion Dynamics

Being based on mixed ionic‐electronic semiconductors, the operation of perovskite solar cells depends on many parameters. To develop an experimentally validated numerical device model, it is therefore necessary to isolate individual physical phenomena. To this end, the dynamics of ion motion in lead halide perovskites is investigated by measuring impedance spectra and the electric displacement as a function of frequency in dark. The displacement response is fully reproduced by a numerical device model that combines electronic and ionic conduction. For a quantitative description of the displacement, it is critical to consider the frequency‐dependent apparent dielectric constant, the ion concentration and the ion diffusion coefficient. The numerical simulations enable to quantify the effect of ion motion and voltage scan speed on the electric field distribution in MAPbI3 based devices, laying the foundations for an experimentally validated perovskite device model.


Device Model for Methylammonium Lead Iodide Perovskite With Experimentally Validated Ion Dynamics
Mohammad Sajedi Alvar,* Paul W.M. Blom, and Gert-Jan A. H. Wetzelaer* DOI: 10.1002/aelm.201900935 functional theory vacancy assisted I − diffusion was predicted to be the fastest diffusion, where positively charged I − vacancies are the mobile species instead of I − ions. [7] For these vacancies, a room temperature diffusion coefficient of 10 −16 m 2 s −1 was calculated. Using transient ion-drift measurements on MAPbI 3 PSCs, Futscher et al. observed that both MA + and I − are mobile and reported an ion density of 10 21−22 m −3 and diffusion coefficient of 10 −16 and 10 −13 m 2 s −1 , respectively. [18] The hysteresis in the J-V characteristics are in accordance with such slow moving positive vacancies. [19] To improve the understanding of the device physics of perovskite solar cells, a numerical device model, including electronic and ionic properties, is indispensable. However, for such a model to be accurate, the large amount of input parameters should be carefully quantified. So far, several device models have been reported for describing the operation of PSCs. [9,14,[20][21][22][23][24][25][26] Reenen et al. explained the hysteretic I-V of PSCs by combining ion motion and charge trapping at the interfaces of perovskite layer. [14] However, their choice of ion density (10 24 m −3 ) and relative dielectric constant (6.5) is not consistent with other reports. [7,9,21] Calado et al. simulated optoelectronic transient measurements on PSCs and reported the same mechanisms for hysteretic I-V, but with higher ion density (10 25 m −3 ) and much lower ion diffusivity (2.5 × 10 −18 m 2 s −1 ). [9] Richardson et al. reported the motion of iodide vacancies as the origin of I-V hysteresis and estimated the ion diffusion coefficient to be around 10 −16 m 2 s −1 . [21] Another device model was developed by Sherkar et al, where they studied the influence of ferroelectric polarization and grain boundaries on the performance of PSCs. [25,26] However, ion dynamics were not included in their model. [25] Recently, Neukom et al. applied an electronic-ionic device model to simulate a series of experiments, but a validation of the multiple fit parameters was missing. [22] In all the aforementioned models the question remains what are the proper input parameters for the ionic and electronic properties. It is clear that the output of a device model strongly depends on the input parameters. Due to the significant impact of ion motion on the electric field distribution in PSCs, all the physical phenomena in PSCs can be influenced by the ionic properties of perovskites. As an example, the movement of ions towards a biased electrode has a strong effect on the electronic properties of a PSC. [27] In a PSC with electron-and hole selective contacts the built-in electric field at zero bias will for example move positive ions towards the negatively charged electrode, as schematically indicated in Figure 1a. Consequently, the electric Being based on mixed ionic-electronic semiconductors, the operation of perovskite solar cells depends on many parameters. To develop an experimentally validated numerical device model, it is therefore necessary to isolate individual physical phenomena. To this end, the dynamics of ion motion in lead halide perovskites is investigated by measuring impedance spectra and the electric displacement as a function of frequency in dark. The displacement response is fully reproduced by a numerical device model that combines electronic and ionic conduction. For a quantitative description of the displacement, it is critical to consider the frequency-dependent apparent dielectric constant, the ion concentration and the ion diffusion coefficient. The numerical simulations enable to quantify the effect of ion motion and voltage scan speed on the electric field distribution in MAPbI 3 based devices, laying the foundations for an experimentally validated perovskite device model.

Introduction
Solar cells based on lead halide perovskites have demonstrated a rapid increase in power conversion efficiency from 4% to more than 25% in less than ten years. [1,2] However, understanding of the operation mechanism of perovskite solar cells (PSCs) has been hampered by the coexistence of electronic and ionic conduction in lead halide perovskites. [3,4] The presence of mobile ions and their accumulation at the interfaces has a large effect on the electric-field distribution inside the perovskite layer, [5][6][7][8][9] which affects the charge extraction in a solar cell, [6,8] In addition, knowledge of the electric-field distribution is also vital for quantifying the steady-state charge-carrier mobility. Furthermore, the occurrence of hysteresis in the current-voltage characteristics, which affects the accuracy of the solar-cell efficiency measurement, [10][11][12] has been attributed to slow-moving ions, [6,[13][14][15] Several candidates for ion motion in lead halide perovskites have been addressed, [4,7] Experimentally, motion of I − ions with a diffusion coefficient of 2.4 × 10 −12 m 2 s −1 was proposed as well as movement of positive MA + with a diffusion coefficient of only ≈10 −16 m 2 s −1 . [16,17] Theoretically, from density field is not homogeneous but mainly drops at the electrode, such that the electric field in the bulk of the active layer is (partially) screened. The amount of screening depends on the dielectric constant, the ion concentration and their ability to follow the applied field, which is governed by the ion diffusion coefficient. In case of weak screening, the photogenerated charge carriers at short-circuit conditions are swept out by the built-in field, the photocurrent is drift dominated. In contrast, for complete screening the photogenerated carriers are driven towards the contacts by diffusion. Since metal-halide perovskites are affected by both electronic charge-carrier and ion movement, many parameters are simultaneously influencing the solar-cell currentvoltage characteristics, such as the ion concentration and diffusion coefficient, the dielectric constant, charge-carrier mobilities, charge trapping, contact barriers, and recombination rates.
Here we present a combined experimental and device modelling study on MAPbI 3 , providing an independent quantification and validation of the ion dynamics, governed by the dielectric constant, ion diffusion coefficient and ion concentration. Experimentally, we use two different methods to quantify the ion diffusion, ion density, and the apparent dielectric constant of MAPbI 3 , independently of the electronic properties. We demonstrate that these parameters can be directly obtained from impedance spectroscopy measurements on Au/MAPbI3/ Au capacitors. The ion dynamics were further investigated by measuring and modeling the electric displacement loops of Au/ MAPbI 3 /Au parallel-plate capacitors. All measurements were carried out under dark conditions. With a numerical device model combining electronic and ionic conduction, both the magnitude and frequency dependence of the electric displacement as a function of electric field are excellently described using anion diffusion coefficient and concentration similar to the impedance data. The mixed electronic-ionic model, with the ionic parameters independently validated by impedance and displacement measurements, is employed to simulate the time-dependent electric field distribution in MAPbI 3 based devices and paves the way for a device model based on experimentally validated input parameters.

Impedance Spectroscopy
To obtain the ion concentration, ion-diffusion coefficient, and frequency-dependent apparent dielectric constant, we first use impedance spectroscopy. In order to avoid parasitic effects on the measurements from surrounding electron and hole transport layers that are typically existing in solar cells, we fabricated bare MAPbI 3 parallel-plate capacitors based on an Au/ MAPbI 3 /Au. As studied by Courtier et al., [28] charge-transport layers significantly affect the electrical characteristics of perovskite-based devices. We have therefore purposely avoided such layers to isolate the ionic parameters of MAPbI 3 thin films. A schematic depiction of the bare perovskite parallel-plate capacitors, along with the employed equivalent circuit for the analysis, is displayed in Figure 1. In Figure 2a the measured real and  imaginary part of the complex impedance Z are shown at zero bias as a function of frequency, together with the resulting dielectric loss ε r ″ and constant ε r ′ (Figure 2b,c), defined as ε r ′-jε r ″ = 1/jωZ, and dielectric loss tangent ε r ″/ε r ′ (Figure 2c) with ω the angular frequency ω = 2πf. The amplitude of the ac voltage used for the impedance measurement was 100 mV. The real part Z′ strongly decreases for frequencies higher than 70 Hz, whereas the imaginary part Z″ is characterized by a maximum at f 1 = 34.4 Hz and a minimum at f 2 = 5.3 Hz. At f = 1 Hz we obtain an apparent dielectric constant ε r ′ of 990, which is identical to earlier reported values. [29][30][31] With increasing frequency ε r ′ strongly decreases to a value of 67 for frequencies higher than 1 kHz. Furthermore, the dielectric loss ε r ″ exhibits a 1/f dependence in the frequency range 1 Hz to 1 kHz, in agreement with Jonscher's law for ionic conductors. [32] Jonscher's law behavior has been reported previously for MAPbI 3 . [33] For further analysis of the data and extraction of the ion mobility we use a model for ionic liquids and electrolytes developed by Bandara and Mellander. [34] We assume that positively charged defects (i.e., iodide vacancies) are mobile, whereas negative defects are fixed. The electrical behavior of our MAPbI 3 thin films between two (ion) blocking contacts can be described by an equivalent circuit as shown in Figure 1b.
We assume that the positively charged ionic species are mobile, and the negatively charged species are immobile. Upon application of an electrical signal with frequency f (or angular frequency ω = 2πf) positive ions will move to the negatively biased electrode, leaving fixed negatively charged ions behind. As schematically indicated in Figure 1, the accumulated positive ions at the electrode form a Helmholz double layer and behave as a parallel RC circuit , where C I and R I represent the capacitance and resistance of the Helmholz layer. This parallel RC element is placed in series with another parallel RC-circuit, where C and R represent the capacitance and resistance of the bulk MAPbI 3 , respectively (C I >> C). We note that the equivalent circuit of Figure 1b can only be applied to analyze the impedance data under the condition that C I >> C, meaning that accumulation of ions at the electrode should take place. [35] When ions are absent or immobile the corre-sponding equivalent circuit would further simplify to a parallel RC circuit. For C I >> C the calculated Z′ and Z″ of such a circuit for a hypothetical case of R I = 5 MΩ C I = 0.5 µF, R = 1.8 MΩ and C = 2.5 nF are shown in Figure 3a. For such a circuit Z′ and Z′′ are given by In the ideal case, the Nyquist plot of a double parallel RC circuit would give two semicircles in the whole range of frequency. The time constants of the two parallel RC elements are given by For the frequency range of ω > 1 R C I I , τ 1 can be calculated using Z′ = R/2. For the low frequency regime, where ω < 1 RC , the characteristic time constant of τ 2 occurs at ′ = + 2 Z R RI . Figure 3 shows the total Z″ as a function of frequency. From the figure two peaks at ω 1 and ω 2 can be distinguished, corresponding to each of parallel RC elements. Between the two peaks, Z" exhibits a minimum at a frequency of ω 3 , which obviously appears at the frequencies higher than ω 2 . The time constant for this characteristic minimum is given by The maximum at ω 1 and minimum at ω 3 thus provides us with two characteristic time constants, representing relaxation Adv. Electron. Mater. 2020, 6,1900935  of the space charge in the bulk and interface, τ 1 and τ 3 respectively, which are related via with δ defined as L/λ, where L is the sample thickness and λ is the Debye length. The Debye length is equal to √(Dτ 3 ). As a result, with τ 1 and τ 3 known, the ion diffusion coefficient D ion is then directly given by We note that for very thick samples [17] ω 3 can become lower than 0.1 Hz and moves out of the accessible measurement range. In that case the impedance data resemble a typical parallel RC circuit ( Figure S1, Supporting Information).
Furthermore, as displayed in Figure 3b, the frequency dependence of the apparent dielectric constant ε r ′ in the model of Bandara and Mellander [34] is given by The experimental impedance data shown in Figure 2 qualitatively exhibit the same functional dependence on frequency as compared to the proposed equivalent circuit. Not only Z″ exhibits a minimum (ω 3 ) and a maximum (ω 1 ), but also Z′ is equal to Z″ at ω 1 . The second maximum occurs in a very low frequency regime, which is out of experimentally accessible frequency range. In Figure 2, the experimental impedance data are fitted with an equivalent circuit model with values of C = 2.7 nF, C I = 0.115 μF, R = 1.7 MΩ, R I = 4.5 MΩ, showing that this simple equivalent circuit captures most of the essential physics. We note that more complex circuits have been used for describing the interfaces in solar cells including charge transport layers. [36] The agreement between theory and experiment allows determination of the ion diffusion coefficient directly from the characteristic frequencies. From the observed frequencies f 1 = 34.4 Hz and f 2 = 5.3 Hz it follows that τ 1 = 4.6 ms and τ 2 = 30 ms, resulting in δ = 42. With a sample thickness L = 240 nm, an ion diffusion coefficient D ion of 1 × 10 −15 m 2 s −1 is obtained via Equation (7). This value is in the same range as the diffusion coefficient obtained from fitting a circuit including a Warburg element [37] (3 × 10 −16 m 2 s −1 ), as well as the coefficient needed to simulate the hysteresis in the J-V characteristics, [21] and furthermore also in agreement with the value obtained from ab initio theory. [7] Using the obtained values for δ and τ 1 we also compare the expected frequency dependence of ε r ′ according to Equation (8) (solid line) with the experimentally obtained results (symbols), as shown in Figure 2c. The high frequency value ε∞′ amounts to 67, in agreement with earlier reported results. [30,33] The calculated increase of ε r ′ at lower frequencies is in reasonable agreement with experiment, given the fact that we approximate the ion conduction in MAPbI 3 by a relatively simple equivalent circuit. Furthermore, the model predicts a static (DC) ε r ′ of about 2900.
As next step, we show that also the ion concentration N ion can be obtained from the impedance data. With the ion diffusion coefficient D ion and thus ion mobility µ known via the Nernst-Einstein relation (µ = eD/kT), knowledge of the ion conductivity σ would suffice to calculate N ion = σ/eµ. The conductivity can be obtained from Jonscher's law for ionic conductors, via the equation ε r ″ = σ/ε 0 ω. A fit using this relation is shown by the black solid line in Figure 2b. Considering a slight deviation from 1/ω behavior in the experimental data, the fit yields an ionic conductivity of σ = (1.4 ± 0.3) × 10 −7 S m −1 , accounting for the fitting error. The ion diffusion coefficient D ion of 1 × 10 −15 m 2 s −1 yields an ion mobility µ of 4 × 10 −14 m 2 V −1 s −1 . Using σ = eN ion µ, an ion concentration N ion of 2 × 10 25 m −3 is obtained, remarkably similar to the value predicted by ab initio theory. [38] From the impedance data we obtain R I = 4.5 MΩ and R = 1.7 MΩ, showing that there is a large resistance present both in the bulk and at the interfaces, most likely due to the presence of an injection barrier. [27] As a result, the electronic current at steady state is very low, having negligible influence on the extracted ionic parameters.
To verify that the deduced ion diffusion coefficient is not influenced by the intensity of the perturbation signal, the impedance spectra of the MAPbI 3 capacitor were measured for different amplitudes of the AC voltages. We note that in order to extract reliable results from impedance measurements, the amplitude of the perturbation signal has to be sufficiently small. In this way, the system is still close to its steady-state condition. Figure S2 in the Supporting Information shows the spectra of the imaginary (a) and real (b) part of the impedance and the dielectric loss (c) as a function of frequency, as well as the extracted ion diffusion coefficient (d) for five different values of V AC . By varying V AC from 10 to 100 mV, the positions of the maximum and minimum in the spectra are approximately constant, yielding similar time constants and, as a result, similar values for the extracted ion diffusion coefficient. However, at a higher V AC of 200 mV, the height of the low frequency peak decreases and the position of the minimum between peaks shifts to lower frequencies. Consequently, the extracted diffusion coefficient at high V AC deviates from the value of diffusion coefficient determined from measurements with lower V AC . These measurements demonstrate that the diffusion coefficient can be extracted reliably for voltage amplitudes of 100 mV and lower. The extracted conductivity was not affected by changing V AC .

Electric Displacement
An independent confirmation of these ionic properties by directly modelling the current-voltage characteristics of solar cells is not trivial since a device model contains many additional parameters, such as contact barriers, charge-carrier mobilities and concentrations, trap concentrations and trap depths, and even an apparent dielectric constant that depends on the frequency and thus the scan rate. Highly desirable is a method that characterizes the ionic properties independently, such that the ionic part of the device model can be validated. Therefore, we measure the electric displacement with a Sawyer-Tower circuit in perovskite based capacitors www.advancedsciencenews.com www.advelectronicmat.de under dark conditions. A schematic diagram of the Sawyer-Tower setup is shown in Figure S3 in the Supporting Information. By measuring the voltage over a reference capacitor with known capacitance, which is placed in parallel with the Au/MAPbI 3 /Au device, the charge density on the plates of the perovskite capacitor can be directly obtained. This charge density equals the electric displacement field D, which can be recorded as a function of the applied bias. As mentioned above, upon application of a bias voltage positive ions drift towards the negatively biased electrode and form an accumulation layer. This accumulation layer screens the electric field in device generated by the amount of electronic charge on the capacitor plates. For large ion concentrations the electric field is almost completely screened. In that case, the amount of ionic charge at the perovskite surface is nearly equal to the amount of (free) electronic charge on the plate, which is represented by the electric displacement D. Consequently, the amount of free charges on the electrodes of the Au/MAPbI 3 /Au capacitors is directly linked to the build-up of ionic charges at the negatively biased electrode. The measured amount of free charges then depends on the ion concentration N, whereas the frequency dependence of the D-E loops depends on how fast ions can follow the applied field, which is governed by the ion diffusion coefficient D I . We note that the use of two gold electrodes prevents the results from being affected by a large built-in electric field as well as chemical reaction between the electrode and perovskite.
The experimental D-V characteristics of the MAPbI 3 -based capacitor as a function of frequency are presented in Figure 4a. The measurements have been carried out in dark. A large hysteresis loop can be observed at 1 Hz, which shows a different shape from a typical loop of a ferroelectric material, excluding ferroelectricity as the origin of the hysteresis. [39] Depending on the applied voltage, at 1 Hz the magnitude of the electric displacement varies between +150 to −200 mC m −2 . At a frequency of 10 Hz, the hysteresis in the electric displacement reduces, with the loop only showing negative displacement. In addition, the magnitude of the electric displacement at 10 Hz compared to 1 Hz decreases by almost one order of magnitude. Upon increasing the frequency to 100 Hz, the hysteresis almost disappears. At higher frequencies of 1 kHz and 10 kHz the hysteresis is absent and the electric displacement of the device exhibits a linear dependence on the electric field with a slope of ε 0 ε r (see Equation (10)).

Numerical Device Model
For a quantitative analysis of the displacement loops, numerical device simulations are required that include the presence and movement of ions. In order to calculate the D-V characteristics corresponding to the experimental results, an ionic-electronic drift-diffusion model was developed. In the device model, it is assumed that the perovskite layer includes a certain density of positively and negatively charged ions. The positively charged ionic species are assumed to be mobile and the negatively charged ionic species are immobile. Furthermore, the ion flow through the metal-perovskite interface is considered to be zero. A detailed description of the model is presented in the Supporting Information.
Due to the fast dynamics of electronic charges and considering the frequency regime of our study, namely 1 Hz to 10 kHz, the transport of electrical charges does not play a significant role in the electric displacement. Therefore, the choice of electronic properties does not significantly influence the simulation results. Importantly, in this frequency regime, the results are a direct consequence of the ionic properties of the perovskite. This enables us to isolate the effect of ion diffusion on the ion distribution and electric field inside the device. In a system in which both ions and charge carriers are present, the electric displacement has two components. One component is related to the polarization charges and the second component is the ionic part. The total electric displacement is the summation of the polarization and ionic part and can be described by following equation where D is the total electric displacement and D polarization and D Ionic are the polarization and ionic portion of the electric displacement. The polarization part of the displacement is defined as follows Adv. Electron. Mater. 2020, 6, 1900935 Here, ε ∞ is the high frequency dielectric constant of MAPbI 3 , ε 0 is vacuum permittivity, and E is the electric field. In addition, the ionic electric displacement can be defined by the following equation where ε r is the apparent dielectric constant and ∆E is difference between the electric field at the interface, where the ions accumulate, respect to the bulk electric field For a quantitative simulation of the perovskite capacitors, the key parameters are the ion concentration, ion diffusion coefficient, and the apparent dielectric constant ε r . The experimentally measured values for ε r , shown in Figure 2c, were used as input parameters in the drift-diffusion model to simulate the D-V characteristics. Using an ion diffusion coefficient of 8 × 10 −16 m 2 s −1 and an ion concentration of 1.9 × 10 25 m −3 the calculated D-V characteristics are shown in Figure 4b for a range of frequencies from 1 Hz to 10 kHz. Similar to the experimental results, the simulated electric displacement shows strong frequency dependence, with the simulations being in excellent agreement with the experimentally obtained characteristics. The simulations reveal that the large hysteresis observed at 1 Hz is due to the slow dynamics of ionic charges. At such low frequency, the ionic component of the electric displacement is dominant, while the polarization part does not play a significant role. In addition, the magnitude of the electric displacement is voltage dependent and similar to the experimental results it varies between +150 to −200 mC m −2 . Moreover, in accordance with experiments, the hysteresis and also the magnitude of the electric displacement drastically reduce at a frequency of 10 Hz. In fact, at a frequency of 10 Hz, the slow ions can marginally follow the variations of the applied electric field. Therefore, the ionic portion of the displacement has a small voltage-dependent variation which leads to a small hysteresis in the D-V characteristics of the device. For frequencies above 100 Hz, no significant hysteresis can be observed. In this frequency regime, the ionic charges are not sufficiently mobile to respond to the applied electric field. Thus, the ionic portion of the electric displacement becomes constant and the polarization part of the electric displacement is the only voltage-dependent term.
To explain the observed displacement loops, the rearrangement of ionic and electronic charges inside MAPbI 3 under the applied bias is visualized in Figure 5. Figure 5a shows the distribution of ions as a function of position at 1 Hz for different bias voltages. After an initial D-V sweep (0 V → +2.5 V → −2.5 V → 0 V) at 0 V a population of ions is accumulated at the interface of MAPbI 3 and the right electrode (blue line), which has a slightly higher work function. This results in a large and negative electric displacement of ≈−200 mC m −2 . As the bias again increases to +2.5 V, the ions start moving toward the opposite (left) electrode. Therefore, the population of the ions at the right electrode gradually decreases, while the ionic population starts increasing on the left electrode. As a result, the electric displacement also reduces (see Figure 4). At a bias voltage of +2.5 V, a certain number of ions accumulate at the left electrode (red line), which leads to a positive electric displacement of ≈+20 mC m −2 . In the subsequent downward scan, from +2.5 V to 0 V, the bias remains positive and consequently the number of accumulated ions at the left electrode increases (green line). As a result, at 0 V the electric displacement has further increased to approximately +130 mC m −2 . Upon further deacreasing the applied bias, the ions start moving toward the right electrode again and concomitantly the electric displacement gradually deacreases. At the maximum bias voltage of −2.5 V, the presented distribution of charges (black line) gives a negative electric displacement of approximately −100 mC m −2 . This displacement continues to decrease while sweeping back to 0 V. Figure 5d shows the position-dependent distribution of mobile ions at various applied biases at frequency of 10 kHz, for which the displacement loops do not show hysteresis. Due to a small built potential, some ions are accumulated at the higher work function (right) electrode. In contrast to what is observed at 1 Hz, in the high frequency regime the distribution of ions approximately persists at all voltages and the hysteresis in the D-V loops disappears. This is a result of the ions not being able to follow the changing electric field at high frequency. Therefore, the ionic part of the electric displacement remains constant, which in combination with the polarization part gives rise to a linear dependence of the electric displacement on the electric field with a slope of ε ∞ ε 0 . From the measurements a slope of 67ε 0 is obtained, which is the high frequency permittivity of MAPbI 3 ( Figure S4a, Supporting Information). Figure 5b,e presents the position-dependent distribution of holes at different applied biases at frequencies of 1 Hz and 10 kHz. The density of holes inside the device is several orders of magnitude lower than the density of mobile ions. Furthermore, the dynamics of electronic charge carriers is much faster than the frequency regime of our study. In other words, the response time of the electronic charges is much faster than the variation of applied bias. The low density and the fast dynamics of holes means that the displacement response is dominated by ionic charges.
The position-dependent distribution of electrons is presented in Figure S5 in the Supporting Information. The injection barrier for electrons is assumed to be about 1 eV, meaning the electron density in the MAPbI 3 layer is negligible compare to the density of holes and ions. Therefore, it has no effect on the observed displacement. Furthermore, the sensitivity of the calculated displacement on the hole-injection barrier was examined. The electric displacement at a frequency of 1 Hz was simulated for a range of injection barriers from 0.2 to 0.5 eV. The results are shown in Figure S7c in the Supporting Information. It can be seen that the magnitude of the hole injection barrier does not have a major effect on the electric displacement loop. This demonstrates that the displacement response of the device is dominated by the ionic contribution.
In order to investigate the contribution of the electronic current experimentally, the J-V characteristics of the MAPbI 3 -based capacitor was measured at different frequencies. For the J-V measurements, the reference capacitor in the Sawyer-Tower circuit ( Figure S3c, Supporting Information) was replaced by a reference resistor. Figure 6 shows the J-V characteristics of the device at different frequencies ranging from 1 Hz to 10 kHz, shown in more detail in Figure S8 in the Supporting Information. Here the frequency f refers to the frequency of the applied voltage which is related to the scan rate via = × 1 ( ) 4 f t scan rate V V s m (where t and V m are the period and the amplitude of the applied voltage, respectively). At a low frequency of 1 Hz, the current is dominated by the low electronic current and the displacement current is negligible. However, as the frequency increases the contribution of the displacement current increases. At a high frequency of 10 kHz, the current is fully dominated by the displacement contribution, while the electronic current has a negligible contribution to the total current. This is evident from the large rectangular hysteretic loops appearing at higher frequencies.
Using the device model with the same set of parameters that were used for calculating the D-V loops, the J-V characteristics of the MAPbI 3 -based capacitor were simulated at different frequencies. The simulations reproduce the experiments, showing the appearance of large hysteresis with increasing frequency.
The simulated distribution of the electric field inside the MAPbI 3 capacitor at 1 Hz and at 10 kHz is shown in Figure 5c,f, respectively. The results for frequencies of 10, 1000, 10 000 Hz are provided in the video in the Supporting Information ( Figure S6, Supporting Information). At 10 kHz, the ions are not sufficiently mobile to respond the applied voltage, upon variation of the applied voltage the difference between electric field at the interface with respect to the bulk value, resulting from ion movement as a result of the small built-in voltage, persists. This effect results in a constant (negative) ionic component in the electric displacement, as shown in Figure 4. At 1 Hz, when no bias is applied at steady state, the ions are mostly accumulated at the higher work function (right) electrode interface and consequently a narrower Debye layer is formed at this interface ( λ = εε β 0 2 k T Ne ) (blue line). From 0 to +2.5 V, the ions start to move toward the opposite electrode. Therefore, the number of ions at high work function electrode decreases, while it enhances on the opposite electrode. Consequently, the sign of the electric field at the electrodes changes and the screening effect is strongly reduced when an applied bias of +2.5 V is reached, the field is nearly constant (red line). From +2.5 V to 0 V, more ions reach the left electrode and the screening effect again increases, while the sign of the electric field in the bulk changes from negative to positive value (green line). As the applied voltage increases from 0 to −2.5 V, the direction of ion motion changes toward the right electrode, meaning the ion density at the low work function (left) electrode decreases and the (positive) electric field in the bulk is further enhanced due to reduced screening (black line).
Adv. Electron. Mater. 2020, 6,1900935  We also verified the need for taking into account a frequency dependent apparent dielectric constant [35] to model the displacement loops. As an example, we have re-calculated the electric displacement as a function of frequency with a fixed apparent dielectric constant of 67 ( Figure S4b, Supporting Information). Keeping the rest of the input parameters similar we observe that for an apparent dielectric constant of 67 the magnitude of D at 1 Hz has strongly decreased to a range of +80 to −90 mC m −2 , which is much lower than the experimental observations. When using the frequency dependent apparent dielectric constant as measured with impedance spectroscopy, the correct value for the electric displacement is obtained at all frequencies. As a result, a frequency dependent apparent dielectric constant is required to reproduce the experimental electric displacements.
In order to demonstrate the sensitivity of the calculated displacement loops on ion concentration and diffusion coefficient, the electric displacement at 1 Hz was calculated for a range of ion densities and diffusion coefficients ( Figure S7, Supporting Information). It is observed that changing the ion concentration within an order of magnitude drastically modifies the magnitude of the electric displacement. For example, an ion density of 5 × 10 24 m −3 gives a variation in electric displacement of +20 to −50 mC m −2 which is much smaller than what is observed in experiment. Conversely, a large ion density of 5 × 10 25 m −3 results in a much larger electric displacement ranging from +350 to −400 mC m −2 , much larger than the experimental values. The best fit of the experimental results is obtained with an intermediate ion concentration of 1.9 × 10 25 m −3 .
As a next step the electric displacement was calculated for a range of ion diffusion coefficients, the ion density and the apparent dielectric constant were kept at the optimum value ( Figure S7b, Supporting Information). For a low ion diffusion coefficient to 8 × 10 −17 m 2 s −1 , there is only small hysteresis in the electric displacement, since the ions are not able to follow the variations in the applied electric field. An ion diffusivity of 8 × 10 −16 m 2 s −1 gives the maximum hysteresis, which quantitatively matches the corresponding experimental D-V loop. In this case, the ionic charges can partially follow the alteration of the electric field with some delay, giving rise to the characteristic hysteresis observed in experiment. By further increasing the ion diffusivity to than 8 × 10 −14 m 2 s −1 . the interface charge density can immediately follow the applied field and consequently negligible hysteresis in the electric displacement is observed. Most importantly, the obtained ion density of 1.9 × 10 25 m −3 and diffusion coefficient of 8 × 10 −16 m 2 s −1 from modelling of the electric displacement are in excellent agreement with the values of 2 × 10 25 m −3 and 1 × 10 −15 m 2 s −1 obtained from impedance analysis performed on the same device layout. As a result, the electric displacement measurements and device modelling independently confirm the ionic parameters directly obtained from impedance measurements.
We note that for ion diffusion coefficients > 10 −13 m 2 s −1 hysteresis is fully absent at a measurement frequency of 1 Hz, since the ions are sufficiently fast to follow the modulation completely. We note that from theoretical and experimental studies for I − ions (or iodide vacancy) a broad range of activation energies (≈0.1-0.6 eV) [7,15,[40][41][42][43][44][45][46][47][48] and diffusion coefficients ranging from 10 −16 to 10 −7 m 2 s −1 , have been repo rted. [7,[16][17][18][19]42,46,49,50] For MA + ions, activation energies of ≈0.4-1.1 eV and diffusion coefficients in the range of 10 −15 -10 −20 m 2 s −1 have been reported. [7,[15][16][17][18]40,41,48,49] Due to the much higher activation energy and lower diffusivity of Pb 2+ ions, they are approximately immobile at room temperature. [7] Here we observed that for high ion diffusion coefficients (>10 −12 m 2 s −1 ) the maximum hysteresis in the electric displacement loops would occur around 10 kHz or higher, which is not observed in our experimental data. Although we cannot accurately distinguish the type of ions with our techniques, our experimental and simulation data show only ionic species with a diffusion coefficient of about 10 −15 m 2 s −1 are present in the perovskite layer. Such slow moving ions are consistent with the observed hysteresis in both electric displacement and J-V characteristics [21] at low frequencies.
Regarding ion concentrations reported values from simulations and theoretical models are rather high in the range of 10 24 m −3 -10 25 (m −3 ), [7,9,21,24] whereas the experimentally measured values are often a few order of magnitude lower (in the range of 10 21 m −3 -10 23 m −3 ). [5,8,18,19,35,43,51] As stated above, measurement of the electric displacement is a direct method to quantify ion concentrations, since the number of free carriers on the capacitor plate is nearly equal to the amount of accumulated ions at the interface. The obtained ion concentration of 2 × 10 25 m −3 is furthermore verified by another independent experimental method. Our experiments combined with simulations show that the observed high electric displacement cannot be generated by ion densities in the range of only 10 21 -10 23 m −3 . Only the higher range of ion density is consistent with both experiment and simulation.
Having demonstrated that the displacement characteristics can be accurately reproduced by simulations when using the Adv. Electron. Mater. 2020, 6, 1900935 Figure 6. Experimental (circles) and simulated (solid lines) current density-voltage characteristics of an MAPbI 3 capacitor at different frequencies ranging from 1 Hz to 10 kHz. At low frequencies, the total current is small and dominated by the electronic current. At 10 kHz, the current is fully dominated by the displacement current.
correct apparent dielectric constant, ion concentration, and ion diffusion coefficient, we now turn to the implications of ionic movement for the potential distribution in perovskite solar cells. The potential profile of PSCs has been experimentally measured and numerically simulated. [8,9,14,21,23] Richardson et al. simulated the potential profile of a PSC and they reported an extremely narrow Debye layer with the width of 1.5 nm. [21,23] Weber et al. successfully measured the potential distribution across a PSC based on mixed-cation mixed halide and recorded the time dependence of electric field distribution after a voltage step. [8] According to their observations, the relaxation time for the ions is slightly below one second. In addition, they reported a much larger value of Debye length. [8] Using our device model with experimentally validated ion dynamics, we have simulated the "time-dependent" potential profile in a perovskite solar cell after applying a voltage step. We assumed a solar cell device with built-in potential of 1 V. Initially, the device was kept under a bias of 1 V to compensate the effect of the built-in potential on the distribution of ionic charges. Then, a voltage step was applied to reduce the applied voltage to 0.5 V across the device. The potential profile across the device in response to the voltage step was calculated as a function of time, as displayed in Figure 7. At short timescales (10 ms) the ions are not sufficiently mobile to redistribute across the film. As a result, the electric potential linearly decreases across the device, which is equivalent to a uniformly distributed electric field across the perovskite layer. After a longer period of time the positive ions gradually reorganize and partially move toward the biased contact. For a relaxation time of 1 s, the system reaches to a steady state and the ions are approximately settled down in a new arrangement with a large population of ions accumulated at the interface. Although Weber et al. used a different composition of perovskite in their solar cell, this relaxation time is in the same range as their experimental value of 0.7 s. [8] As a result, we here demonstrate that the dynamical potential profile across a perovskite solar cell can be simulated with our device model. Qualitative comparison of the simulated potential profile by our device model and other reports of the experimental potential profile, shows that our validated ionic parameters are realistic.

Conclusion
In conclusion, we have developed a numerical device model with experimentally validated ion dynamics. The apparent dielectric constant, diffusion coefficient and concentration of ions in MAPbI 3 were quantified using impedance spectroscopy. An ion diffusion coefficient D ion = 1 × 10 −15 m 2 s −1 and an ion concentration N ion = 2 × 10 25 m −3 were obtained. These numbers were confirmed independently by frequency-dependent displacement measurements and simulations. Using the developed electronic-ionic drift-diffusion model, the experimental D-V loops of perovskite based capacitors were accurately reproduced with identical ion concentration and diffusion coefficient as extracted from the impedance measurements. The frequency dependence of the impedance and displacement measurements allows us to isolate the ionic contribution from the electronic contribution by the charge carriers.
The accurate description of the experimental electric displacement loops by our numerical simulations is an important step in the development of an experimentally validated perovskite device model. As such a device model consists of many parameters, it is important to first isolate the ionic contribution to the device physics, prior to modeling the charge transport and full current-voltage characteristics of a perovskite solar cell. As first application of the model we have simulated the timedependent potential profile in a MAPbI 3 solar cell, based on experimentally validated parameters, which is dominated by the movement of ions in the layer.

Experimental Section
Device fabrication: Cr/Au (1 nm/50 nm) was deposited as the bottom contact on a glass substrate. A mixture of lead acetate trihydrate (Pb(Ac) 2 ) and methylammonium iodide (MAI) in a 1:3 molar ratio in N,N-dimethylformamide (DMF) was prepared and spin coated on the substrate inside a nitrogen-filled glovebox. The films were annealed at 100 °C for 30 min. Finally, the capacitors were finished by the deposition of 50 nm top Au electrodes, giving rise to a device area of 1 mm 2 .
Device Characterization: Impedance measurements were conducted under controlled N 2 atmosphere, using a computer-controlled Solartron 1260 impedance analyzer. Electric displacement measurements were carried out an output of Sawyer-Tower method ( Figure S3, Supporting Information). As function generator a Tektronix AFG 3022B has been used, the input voltage and output voltage were recorded with a Waverunner LT374 LeCroy Digital Oscilloscope.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.