Influences of dielectric constant and scan rate on hysteresis effect in perovskite solar cell with simulation and experimental analyses

In this work, perovskite solar cells (PSCs) with different transport layers were fabricated to understand the hysteresis phenomenon under a series of scan rates. The experimental results show that the hysteresis phenomenon would be affected by the dielectric constant of transport layers and scan rate significantly. To explain this, a modified Poisson and drift-diffusion solver coupled with a fully time-dependent ion migration model is developed to analyze how the ion migration affects the performance and hysteresis of PSCs. The modeling results show that the most crucial factor in the hysteresis behavior is the built-in electric field of the perovskite. The non-linear hysteresis curves are demonstrated under different scan rates, and the mechanism of the hysteresis behavior is explained. Additionally, other factors contributing to the degree of hysteresis are determined to be the degree of degradation in the perovskite material, the quality of the perovskite crystal, and the materials of the transport layer, which corresponds to the total ion density, carrier lifetime of perovskite, and the dielectric constant of the transport layer, respectively. Finally, it was found that the dielectric constant of the transport layer is a key factor affecting hysteresis in perovskite solar cells.

esis phenomenon. In this section, the structure is Au/Spiro-OMeDA/MAPbI 3 /SnO 2 /ITO, and modeling parameters are shown in Table 1. To evaluate the effect of ion accumulation, the default scan rate was set to 10 5 mV/s . The complete modeling sequence involved the following steps: (1) pre-condition at 0 V for 10 s to demonstrate ion accumulation at the interface; (2) calculate the J-V curve under forward scanning conditions from 0 V to V oc at a scan rate of 10 5 mV/s (to demonstrate the situation when the ion accumulates at the interface, a faster scan rate is needed in this case); (3) pre-condition at V oc for 10 s to evaluate ion accumulation at the interface; and (4) calculate the J-V curve under reverse scanning conditions at a scan rate of 10 5 mV/s . Figure 1a,b show the net www.nature.com/scientificreports/ ion distribution and the conduction band after pre-conditioning at 0 V and V oc for 10 s, respectively. To achieve a steady-state, the ions redistribute to screen the built-in electric field, generating an ion-induced electric field with a direction opposite the built-in electric field. Figure 1a,b show that the built-in electric field is screened by the ion-induced electric field; thus, ion accumulation at the interface can be observed. Because the direction of the built-in electric field is opposite that for the pre-condition at 0 V and V oc , the net ion distribution is also opposite. Bias was then applied from 0 V to V oc (forward scanning) and from V oc to 0 V (reverse scanning). Figure 1e shows the J-V characteristics under the forward and reverse scanning conditions, and apparent hysteresis behavior can be observed. There is a cross between V oc under forward and reverse scanning. This reason is that the ion distribution after the forward scanning and before the reverse scanning is like effective enhanced p-type doping near the HTL interface and enhanced n-type doping near the ETL layer (as shown in Fig. 1b). Hence, we can notice that there is a difference in the turn-on voltage between forward and reverse scanning in Fig. S9a in the dark current condition. The reverse scanning condition in the dark current condition has a relatively smaller contact resistance, where carriers are more easily injected and recombined. This effect would be screened if the photo current is higher, so we only can observe this phenomenon when the current density is around 0 (voltage is around Voc) under illumination or under dark. Also, Fig. S9b reflects this phenomenon and we can notice that the Voc is almost the same for case of net ion density being 10 16 cm −3 . The hysteresis phenomenon can be explained by the conduction bands in the two cases. Figure 1c,d show the conduction bands for the forward and reverse scanning, respectively. For the forward scanning, the direction of the ion-induced electric field is from the HTL to the ETL, the same direction as the applied electric field. Figure 1c shows that electrons are accumulated at the side of the HTL, leading to a decrease in carrier extraction. It is not an ideal condition for carrier extraction when the direction of the ion-induced electric field is the same as that of the applied electric field. For the reverse scanning, the direction of the ion-induced electric field is from the ETL to the HTL, the opposite direction as the applied electric field. Under this condition, carriers can be easily extracted because the applied electric field is screened by the ion-induced electric field. Unlike the forward scanning, the ion-induced electric field can improve carrier extraction in this situation. The difference in the non-radiative recombination rate distribution between the forward and reverse scanning is shown in Fig. 1f. Hence, the hysteresis phenomenon can be attributed to the difference in non-radiative recombination between the forward and reverse scanning, which is caused by the ion-induced electric field.
In summary, interfacial ion accumulation is caused by the built-in electric field and leads to an opposite ioninduced electric field in the active layer. Because carrier extraction is affected by the applied electric field, the ion-induced electric field has the opposite direction under the forward and reverse scanning. Therefore, carrier extraction differs between these two cases, leading to different non-radiative recombination at the same voltage. This results in different current densities when measured at the same voltage. Hence, the hysteresis phenomenon is observed. In the next section, the effect of the voltage scan rate is discussed.  To discuss this issue clearly, in this work, both experiments and modeling are demonstrated, the structures are the same as Au/Spiro-OMeDA/MAPbI 3 /SnO 2 /ITO, and modeling parameters are shown as Table 1. Furthermore, to express whole process of forward and reverse scanning clearly, there are animations including scan rates of 10, 700, and 10,000 mV/s in supporting information.
Before discussing the effect of scan rate, it is necessary to define a factor to quantify the degree of hysteresis. Several such factors have been defined in previous studies 41,60 . In this work, the hysteresis index (HI) is defined as: where PCE FS and PCE RS are power conversion efficiency under the forward and reverse scanning conditions, respectively. Figure 2a schematically shows the definition of HI. Figure 2b shows that both of the experimental and modeling results show a similar trend. The experiment shows there are different hysteresis degrees under the different scan rates, and the hysteresis pattern like a Gaussian distribution, which means a lower scan rate (below 100 mV/s) shows lower hysteresis, then hysteresis phenomenon would get more serious with scan rate increasing. But hysteresis decrease as the scan rate is higher than 700 mV/s. Figure 2c-f show that the degree of hysteresis changes with scan rate (see the Supporting Information Section Experimental Data to find all experimental data in J-V curve). In the following section, the simulation program was utilized to explain the mechanism of these results. The modeling results aim to demonstrate the effect of scan rate on hysteresis along with the associated mechanism.
The protocol for the simulation described in the next section involved the following steps (different from the previous section, there is no pre-bias in this section): (1) calculate the current density from a bias of 0 V to V oc under a specific scan rate; and (2) after forward scanning, immediately calculate the current density from a bias of V oc to 0 V using the same scan rate. The scan rate in this simulation ranged from 1 to 10 5 mV/s. The modeling results of HI at scanning rates of 1-10 5 mV/s are shown in Fig. 2b. The maximum HI is observed at the scan rate of 700 mV/s, and the degree of hysteresis gradually decreases as the scan rate increases or decreases from this value. As mentioned before, Tress et al. reported this phenomenon in a previous work 59 . Levine et al. also reported similar behavior; they termed the appropriate scan rate as the "system response region" and found that hysteresis also occurs in the p-i-n structure but with a different system response region than for the n-i-p structure 61 . To understand the mechanism of this effect, three scan rates (10, 700, and 10 4 mV/s , corresponding to slow, appropriate, and fast scan rates, respectively) were chosen for further evaluation. The conduction bands, net ion densities, and electric fields for these cases are shown in Fig. 3a-f. www.nature.com/scientificreports/ Figure 3a shows that under the forward and reverse scanning conditions at the scan rate of 10 mV/s, the net ion densities are similar, and the conduction bands are flat. Because the ion movement responds to the scan rate, the net ion distribution depends on the specific bias when the scan rate is slow. Figure 3d shows that the electric field of the active layer is quite small under both forward and reverse scanning conditions. Because the ions have sufficient time to accumulate at the interface, the built-in electric field is completely screened by the ion-induced electric field. Hence, the J-V characteristics are nearly the same under the forward and reverse scanning conditions, and the hysteresis is weak. For the case of the fast scan rate (Fig. 3b), the velocity of ion redistribution cannot respond to the high scan rate (10 4 mV/s ), and less ion accumulation is observed compared to the other two scan rates. Therefore, the effect of the screening of the built-in electric field is not apparent. Figure 3e shows that the electric field is almost the same under the forward and reverse scanning conditions, and the hysteresis is weak. Figure 3c shows the net ion density and conduction band distribution under forward and reverse scanning conditions at 0.6 V. The difference in ion density between the forward and reverse scanning is quite large; thus, the conduction band also differs between the forward and reverse scanning. Figure 3f shows that the difference in the electric field (green line) between the forward and reverse scanning at the same voltage is larger at the scan rate of 700 mV/s than at the other two scan rates. Hence, the scan rate of 700 mV/s is the system response scan rate in this case.
In summary, hysteresis cannot occur when the scan rate is larger or smaller than the system response scan rate. At a scan rate that is lower than the system response rate, the net ion distribution depends on the bias, and the ions have sufficient time to accumulate at the interface. Therefore, the measurement results are the same under the forward and reverse scanning. For scan rates larger than the system response rate, the ions cannot respond to the voltage scan, and the ion effect is thus quite weak. Thus, the measurement results are again the same under the forward and reverse scanning. As mentioned earlier, some studies have shown that the system response region is related to the degree of perovskite dissociation, the quality of perovskite, and the properties of the transport layer. These are the possible reasons for the observed difference in HIs between the n-i-p and p-i-n structures. The following section presents the simulation results for different net ion densities, carrier lifetimes, and properties of the transport layer.
Effect of transport layer's dielectric constant on hysteresis. As mentioned above, PSCs with n-i-p structures exhibit significant hysteresis, while hysteresis is weak in devices with p-i-n structures. Some studies have indicated that the properties of the transport layer are important in the degree of hysteresis 34,38 . Fang et al. reported that PCBM, a fullerene material, plays an important role in the p-i-n structure. They attempted to remove PCBM and adjust the thickness of PCBM and discovered that using PCBM can reduce the hysteresis due to the trapping of ions by PCBM. However, some studies have shown that hysteresis can be further reduced by replacing TiO 2 with SnO 2 as the ETL in the n-i-p structure [62][63][64] . Previous studies have attributed the reduced hysteresis to the permittivity of the transport layer [39][40][41] . In this section, both experiments and modeling were developed to support this viewpoint and analyze the mechanism. We replaced the materials of ETL from SnO 2 ( ǫ r = 9) to TiO 2 ( ǫ r = 31) in the experiment, a significant increase of hysteresis index can be observed (see Supporting Information Section Experimental Data for the complete experimental data under different scan rates). The experimental data of TiO 2 shows a serious hysteresis. The reason is surface traps at interface between TiO 2 and perovskite layer besides a higher dielectric constant 65 . By investigating and comparing reported cases in which the hysteresis was reduced and the experiments in this work, we can note that the dielectric constant is an important factor. The dielectric constants of PCBM, SnO 2 , and TiO 2 are 3, 9, and 31, respectively. The dielectric constant reflects the degree of polarization in the material; thus, a higher dielectric constant corresponds to a lower electric field and vice versa. The dielectric constants of metallic oxides and organic materials are significantly different. Therefore, the difference in dielectric constant is a possible explanation for the difference in hysteresis degree.
In this section, in order to discuss the effects of different dielectric constants purely, the dielectric constant of the ETL was changed from 3 (organic material) to 31 (metallic oxide), while the other parameters including conduction and valence bands of ETL remained the same as Table 1. Figure 4b shows the changes in HI with scan rate for the dielectric constants of 9 and 31 in this simulation. Hysteresis is reduced when the dielectric constant changes from 31 to 9. The scan rate of 700 mV/s was chosen to analyze the effect of dielectric constant on HI. As shown in Fig. 4b, the HI is lower when the dielectric constant is 3 (corresponding to PCBM as the ETL) than when the dielectric constant is 31 (TiO 2 as the ETL). The HI for the dielectric constant of 9 (SnO 2 as the ETL) is lower than that for the dielectric constant of 31 (TiO 2 ). The cases with the dielectric constants of 3 and 31 were chosen to explain the observed effect, as shown in Fig. 4c. Figure 4d,e show the conduction band and net ion density for the two cases at the voltage of 0.6 V, respectively. For the conduction band, the band bending of the ETL is stronger when the dielectric constant is 3 than when it is 31 (Fig. 4d). The lower dielectric constant results in a stronger degree of polarization. Thus, the electric field of the ETL is higher when the dielectric constant is 3 than when it is 31. The entire device can be simplified as a series circuit; hence, the built-in electric field of perovskite increases when the electric field of the ETL decreases. Therefore, more ion accumulation also occurs, as observed in Fig. 4d,e. Figure 4f shows the electric field in the middle of the active layer. Compared to when the dielectric constant is 3, the electric field is stronger when the dielectric constant is 31, and the degree of hysteresis is thus higher.
In conclusion, the ETL plays an important role in the degree of hysteresis. The effect of the ETL can be attributed to various properties of the ETL, especially the dielectric constant. Due to the polarization of the ETL, the electric field of the ETL is weak when the ETL material has a high dielectric constant (e.g., metallic oxide). This causes the built-in electric field of perovskite to be strong, leading to serious hysteresis because the strong builtin electric field results in serious ion accumulation. This explains why studies have found reduced hysteresis when the ETL is changed from TiO 2 to another metallic oxide with a lower dielectric constant, such as SnO 2 .  66 has shown that the light intensity would affect the ion density in perovskite. Hence, it is important to investigate the influences of these two parameters. However, discussing these two factors is quite hard in experiments. The reason is that degradation and quality of perovskite depend on the type of perovskite, doping other materials (such as K, Rb, and Li), or crystallizing of perovskite [67][68][69] . Hence, modeling is utilized to discuss these two factors without experiments in this section. The structure is Au/Spiro-OMeDA/MAPbI 3 /SnO 2 /ITO, and modeling parameters are shown in Table 1 except ion density in perovskite.
A previous study reported a total ion density of nearly 10 18 cm −3 in perovskite 26 . However, recent studies have shown that the total ion density is determined by multiple factors, including the process, material, and structure. As we mentioned before, some materials have been doped into the perovskite layer to capture ions ( I − ), including K, Rb, and Li, thereby reducing hysteresis. In this section, the total ion density ( N ion ) is changed from 10 16 to 5 × 10 18 cm −3 while the other parameters remain the same. Figure 5a shows the changes in HI with scan rate for different N ion . The results can be divided into two behaviors: (1) N ion = 10 16 to 10 17 cm −3 ); (2) N ion = 10 18 to 10 19 cm −3 ). For the lower total ion density [ N ion of 10 16 and 10 17 in Fig. 5a, the hysteresis gradually increases as the total ion density increases, as shown in Fig. 5b. The hysteresis is quite weak when N ion is only 10 16 cm −3 . Figure 5c shows that the difference in the conduction band of the perovskite layer (black line) is not obvious between the forward and reverse scanning at the same voltage of 0.6 V. Compared to the other cases, the ion accumulation (blue line in Fig. 5c) is relatively low in this case ( ≈ 10 16 cm −3 ). Because the total ion density is not sufficient to screen the built-in electric field, hysteresis is quite weak. When N ion is increased to 10 17 cm −3 , the hysteresis becomes obvious. Figure 5d shows that the change in the conduction band of the perovskite layer (black line) is greater than in the case of lower total ion density, while the ion accumulation (blue line) is higher in this case ( ≈ 10 17 cm −3 ). Because the extra ion density results in a greater degree of screening of the built-in electric field compared to in the case of lower total ion density, the hysteresis is obvious.
For the higher total ion density ( N ion of 10 18 and 10 19 in Fig. 5a), the patterns in HI are similar at the different scan rates. The system response region shifts with total ion density and the point of maximum HI changes; the maximum HI values are 700 and 5000 mV/s for N ion = 1 × 10 18 and 10 19 , respectively. As shown in Fig. 5e, the J-V characteristics are the same for the two cases; thus, the HI values are also the same. However, Fig. 5f shows that the net ion density distribution is not the same for the three cases. The net ion density at the interface is higher when the total ion density is higher. Although the distributions are different, the summations are almost the same, indicating that the ion-induced electric fields are the same. This phenomenon can be explained by the different timescales of ion accumulation. Because the net ion distribution is decided by the ion drift-diffusion equation [Eq. (14)], the time required to achieve sufficient ion accumulation to screen the built-in electric field is shorter than in the other cases when the total ion density is quite high. Since the time is shorter and the voltage www.nature.com/scientificreports/ is the same, the system response scan rate is higher than in the other cases. In other words, it takes less time to accumulate sufficient ions to screen the built-in electric field. Figure 6a,b show that hysteresis is more serious at the carrier lifetime of 20 ns than in cases with longer carrier lifetimes. To conveniently examine this situation, we chose the scan rate with the maximum HI (700 mV/s ). Some studies have attributed strong hysteresis in PSC devices to high non-radiative recombination, and poor crystalline quality 70 . Indeed, the quality of perovskite can affect the performance of the PSC device; however, the reason for the high degree of hysteresis is less well studied. Figure 6c,d show the electron/hole distributions and conduction bands for carrier lifetimes of 20 and 200 ns under forward and reverse scanning, respectively. Under the forward scanning, the carrier distribution of both electrons and holes is less for the carrier lifetime of 20 ns than for the lifetime of 200 ns. Furthermore, at the voltage of 0.6 V, the conduction bands are similar for the forward and reverse scanning. Figure 6e shows that the net ion distribution is almost the same; thus, the difference in the degree of hysteresis has nothing to do with ion accumulation. Figure 6f shows that the non-radiative recombination is stronger for the carrier lifetime of 20 ns than for the lifetime of 200 ns. This is the primary reason that perovskite quality affects the performance of the PSC device. The different degrees of hysteresis under different carrier lifetimes are related to the proportions of non-radiative recombination. The non-radiative recombination of the active layer for each case is shown in Table 2. While ion accumulation can improve the performance under reverse scanning conditions in both cases, the performance in the case of a 20 ns carrier lifetime is quite low under forward scanning conditions. The improvement in non-radiative recombination from ion accumulation is relatively greater than the case of 200 ns. Hence, the improvement of performance is greater than the case of 200 ns, and it causes a stronger hysteresis.

Conclusion
In this work, perovskite solar cells (PSCs) with different ETLs were fabricated to understand the hysteresis phenomenon under a series of scan rates. The experimental results show that hysteresis phenomenon would be affected by the dielectric constant of ETL and scan rate significantly. Also, a modified Poisson-DD solver coupled with a fully time-dependent ion migration model was used to analyze ion migration in different situations and reveal the mechanisms of the observed effects. Hysteresis was shown to be caused by ion accumulation. Due to the built-in electric field of perovskite, ions accumulate at both sides of the interface and form an ion-induced (e,f) Net ion density distribution and non-radiative recombination rate, respectively, for the two cases at a voltage of 0.6 V for the forward and reverse scanning. Table 2. Non-radiative recombination rates at 0.6 V.  www.nature.com/scientificreports/ electric field. Therefore, the J-V curves are not the same under the forward and reverse scanning. The different dielectric constants will lead to different electric fields in the active layer and help us to verify this concept. Generally, the performance of the PSCs device is better under the reverse scanning than under the forward scanning because carrier extraction is enhanced by the ion-induced electric field under reverse scanning, while it decreases under forward scanning. This is the primary factor that leads to hysteresis. Generally, the ion would cause detriment to performance because it would cause other effects in the real case. For example, it would generate more defects or reduce the lifetime of the device. These factors might worsen the device performance and lead to a stronger hysteresis effect. The change in HI with scan rate was demonstrated, and the curve of HI vs. scan rate was shown to be Gaussian-shaped rather than monotonic. Hysteresis is weak when the scan rate is lower or higher than the system response scan rate. At low scan rates, the net ion distribution depends on the bias. At fast scan rates, the ions cannot respond to the voltage scan, and the ion effect is thus quite weak. Hence, the HI is low under both slow and fast scan rates, although the mechanisms are different. Hysteresis only occurs when the scan rate lies in the system response region. Possible reasons leading to different degrees of hysteresis in previous works include the quality of the perovskite, degree of perovskite degradation, and properties of the ETL; these factors correspond to the carrier lifetime of perovskite, total ion density, and dielectric constant, respectively. The mechanisms of the effects of the above-mentioned factors were demonstrated and explained. The effect of total ion density is explained in a relatively intuitive fashion. Hysteresis is weak when the total ion density is low; thus, doping with K, Rb, and Li can significantly reduce the hysteresis. The effect of carrier lifetime is related to the different proportions of non-radiative recombination under the forward and reverse scanning. Hence, hysteresis is stronger when the lifetime of perovskite is lower (i.e., the perovskite quality is poor). Finally, the properties of the ETL, especially the dielectric constant, play an important role in hysteresis in PSC devices. According to the modeling results, using an ETL with a low dielectric constant can reduce the degree of hysteresis. This is because the electric field of the ETL increases as the dielectric constant of the ETL decreases; thus, decreasing the dielectric constant can reduce the built-in electric field of perovskite and hence reduce hysteresis. Device characterization. The device performances were measured using a source measurement unit instrument (Model 2420 Source Meter, Keithley Instruments) and a solar simulator (Sun 2000 Class A, ABET technologies) under the standard AM 1.5G (100 mW/cm 2 ). The measurement method of cases with different scan rates adjusts the measurement voltage interval to change the scanning speed. In this work, the measurement voltage range is from + 1.5 to − 0.5 V for all cases, and the range of scanning rate is from 10 to 1200 mV/s. The number of measurement data points will decrease, leading to the J−V curve's characteristic not being clear if the scanning rate is more significant than 1200 mV/s.

Modeling section. Simulation of optical properties.
To model the light absorption in the device, the FD-TD method was implemented to model wave propagation and absorption. The FD-TD program was developed by our lab. The first FD-TD algorithm, which was proposed by Prof. Yee, is based on solving Maxwell equations in the time and space domains 71 . The details of the FD-TD modeling, including the simulation domain and parameters, and the modeling results are summarized in the Supporting Information (Section FD-TD modeling).
Poisson and drift-diffusion solver. After the optical properties were obtained using the FD-TD solver, the Poisson-DD solver was applied to calculate the electrical characteristics of the device 72,73 . Briefly, this solver is based on a self-consistent solution of the Poisson equation: the drift-diffusion equations: and (2) ∇ r · (ε∇ r V (r)) = q(n − p + N anion − N cation ), (3) J n = −qµ n n(r)∇ r V (r) + qD n ∇ r n(r) where q is the elementary charge; ε is the dielectric constant; V is the potential of the device; n and p are the electron density and hole density, respectively; and N cation and N anion are the cation and anion density, which are calculated by the ion migration model. The drift-diffusion equations are shown in Eqs. (3) and (4), where J is the current density, µ is the mobility, and D is the diffusion coefficient. The subscripts e and h indicate electrons and holes, respectively. The continuity equation is shown in Eq. (5), where R is the recombination rate, which is determined by Shockley-Read-Hall (SRH) non-radiative recombination, and the radiative recombination coefficient (B) is shown in Eqs. (6) and (7): where τ n and τ p are the non-radiative carrier lifetimes of electron and hole, respectively; n i is the intrinsic carrier density; and G opt is the generation rate, which is obtained from Eq. (8): where ω is the photon energy; P is the Poynting vector; n is the refractive index; k is the extinction coefficient; and E is the electric field.
Gaussian density of state. The band edges of organic material are called the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO); these are similar to the valance band and conduction band. However, different from semiconductors, the band edges are not abrupt. Hence, some tail states exist between the LUMO and HOMO, meaning that the density of state is a disordered distribution. Bässler et al. proposed that the disordered density of state in an organic material has a Gaussian-like shape 46 . As we had mentioned before, using Gaussian density of state can simulate the transport layer of organic material more correctly, which shows in Fig. 7b,c. Figure 7b shows that carriers can not be extracted efficiently as Gaussian density of state is unimplemented. Using heavily doping in transport layers is a common and feasible method to solve this problem, but it can not present the real case. Figure 7b shows using Gaussian density of state can achieve a similar result which implementing heavily doping, and the advantage of this method can be observed in Fig. 7c. The conduction and valence band bending can be demonstrated by using Gaussian density of state, which is a key factor to discuss the mechanism of hysteresis. Hence, Gaussian density of state [Eq. (9)] is implemented to describe the tail states in this work: R =SRH + Bnp, www.nature.com/scientificreports/ where N t is the total density of state; σ is the broadening factor of the Gaussian shape; and E t is the difference between the center of the Gaussian density of state and the LUMO or HOMO. And the Gaussian density of state that have successfully utilized in our previous works 35,36 . Details of the Gaussian density of state, including a schematic and the parameters used in this work, are shown in the Supporting Information (Section Gaussian Density of State). When a Gaussian density of state is applied, the carrier density can be expressed as otherwise, the carrier density is expressed as and where m e,h is the effective mass of the electron or hole, respectively; E c,v is the conduction or valance band, respectively; and f e,h is the Fermi-Dirac distribution function for the electron or hole, respectively.
Time-dependent ion migration model. To evaluate the effect of the voltage scan rate on hysteresis, the timedependent ion migration model was coupled with the Poisson-DD solver. The ion migration model considers both the anion (iodide, I − ) and cation (methyl ammonium, MA + ). However, some studies have shown that the mobility of I − is 10 4 times faster than the mobility of MA +58,74 . Due to the low mobility of MA + , MA + cannot respond to typical scan rates, which range from 10 to 10 4 mV/s. Hence, MA + is almost immobile in this model. Therefore, although both types of ions are considered in the model, the net ion density is defined as the sum of I − and MA + density. In this work mainly reflects the movement of I − density, which is shown in Fig. 7d and the net ion density is zero before the ion movement. To evaluate cation and anion migration, the following cation and anion drift-diffusion equations were used: and where N is the ion density distribution; µ is the ion mobility; E is the electric field; and D is the diffusion coefficient, which is determined by the Einstein relation D = µk B T , where k B is the Boltzmann constant, and T is room temperature. The subscripts cation and anion indicate MA + and I − , respectively. The flowchart of the simulation process and the details of the time-dependent model are shown in the Supporting Information (Section. Timedependent Ion Migration Model).

Data availability
The supporting information is available on the website including figures and animations and the modeling tool was developed by our laboratory. More detailed information about the modeling tool can be found on our website (http:// yrwu-wk. ee. ntu. edu. tw/).