The dominant role of surfaces in the hysteretic behavior of hybrid perovskites

Journal Pre-proof The dominant role of surfaces in the hysteretic behavior of hybrid perovskites Claudia Caddeo, Alessio Filippetti, Alessandro Mattoni PII:

S2211-2855(19)30869-9

DOI:

https://doi.org/10.1016/j.nanoen.2019.104162

Reference:

NANOEN 104162

To appear in:

Nano Energy

Received Date: 17 August 2019 Revised Date:

27 September 2019

Accepted Date: 1 October 2019

Please cite this article as: C. Caddeo, A. Filippetti, A. Mattoni, The dominant role of surfaces in the hysteretic behavior of hybrid perovskites, Nano Energy (2019), doi: https://doi.org/10.1016/ j.nanoen.2019.104162. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

The Dominant Role of Surfaces in the Hysteretic Behavior of Hybrid Perovskites Claudia Caddeoa, Alessio Filippettia, b, Alessandro Mattoni * a a

Istituto Officina dei Materiali (CNR - IOM) Cagliari, Cittadella Universitaria, I-09042 Monserrato (Ca), Italy

b

Dipartimento di Fisica, Università degli Studi di Cagliari, Cittadella Universitaria, I-09042 Monserrato (Ca), Italy *

[email protected] Abstract Hysteresis and memory-effects are an issue for the development of a reliable solar technology based on hybrid perovksites. In this context, we study the effect of surfaces on the electric mobility of charged point-defects at room conditions by molecular dynamics simulations. At variance with the bulk, for which the ionic current is linear and reversible upon inversion of the electric field, we find that the presence of surfaces gives rise to accumulation of charged point-defects and memory effects. Surfaces act as trapping regions for iodine defects with fast accumulation and much slower release kinetics (by at least two orders of magnitude). The calculated release times are of hundreds of millisecond for both vacancies and interstitials at MAI- and PbI-terminated surfaces, respectively and are consistent with the decays of ionic currents observed in time-resolved Kelvin Probe Force Microscopy experiments. As we do not find ferroelectric polarization and memory effects on ionic currents in perfect crystalline bulks at room conditions, we conclude that surfaces (and likely interfaces or boundaries) are key factors giving rise to hysteresis. Present results are complemented by an analytic model of the effect of surfaces on point defects and of the kinetics of ionic polarization for the design of functional films with better control of hysteretic ionic currents.

Keywords: defects, simulations, model potential, polarization, memory effects

1. Introduction Perovskite-based solar cells promise to be a leading technology in photovoltaics due to their power conversion efficiencies well above 20% combined to simple and inexpensive solution processes[1]. Despite hybrid perovskites such as CH 3 NH 3 PbI 3 (MAPI) show elevated microstructural disorder their electronic properties are highly resilient to defects. This has been attributed to the fact that most common defects only create shallow levels inside the band gap which do not enhance carrier recombination[2, 3, 4, 5, 6]. Material disorder is however responsible for undesired memory effects such as current-voltage hysteresis, causing strong dependence of the performances on the history of the sample. There is a general agreement among researchers on the fact that such phenomena are connected to the migration of charged ionic species under the effect of electric fields. The high concentration of point defects in MAPI (especially solution processed samples)[7, 8, 9, 10], which can easily migrate throughout the device[11, 12, 13, 14] has been in fact associated to several peculiar phenomena, such as switchable photovoltaic effect[15, 13], hysteresis in the current-voltage curves[16, 17, 18], and degradation[19, 20, 21]. It is also worth noting that defect-mediated ion migration can also have beneficial effects[22] and enable new applications of MAPI (other than the solar cells) exploiting ion diffusion. Examples include the use of hybrid perovskites for resistive switching and memory devices[23, 24]. Concerning memory effects and current density-voltage (J-V) hysteresis they are ultimately connected phenomena, and there is a general consensus that they depend on migration of various ionic species[25]. Hysteresis was originally attributed to the reorientation of organic cations with permanent electric dipole and it was proposed to substitute MA ions with zero dipole molecules[26, 27]. Mixed organic-inorganic halide perovskites featuring organic cations with reduced dipole moments (e.g., formamidinium and/or guanidinium ions) have been realized [28, 29] partly eliminating hysteresis. However, current-voltage hysteresis occurs over timescales (ms-s) much longer then molecular reorientation (1-100 ps)[25, 30, 31]. The beneficial effect of zero dipole molecules on hysteresis is likely indirect. For example, the reduction of the electric capacity at the

interface between the perovskite active layer and the electron transport layer, which in turn reduces the ion accumulation at that interface[32]. Different explanations of J-V hysteresis have been proposed, and they all describe the formation of charged layers of mobile species which migrate with different timescales, screening the electric field inside the perovskite[33, 12, 16, 34]. Due to this timeevolving screening, the charge extraction efficiency and, consequently, the photocurrent depend on the history of the sample. The effect of surfaces on hysteresis is debated, and it is not clear if charge trapping at grain boundaries or at the interfaces between the perovskite and the electron/hole transport layers (ETL/HTL) is crucial[35, 36]. Recently, Weber et al. have observed the fast (few ms) formation of a localized positive charge close to the ETL upon application of an electric field. This localized charge is stable for over 500 ms after switching off the bias, generating an opposite electric field inside the sample which is time-dependent. It is believed that the localized charge is due to trapped iodine vacancies, but a direct experimental demonstration of the nature of the trapping is missing. Atomistic simulations are suitable to study specific point-defects. Ab initio methods have been applied to study electronic effects of point-defects either in the bulk[38, 39, 40, 41] or at the surfaces[42, 41], but, due to the associated computational cost, they cannot easily afford large-scale models necessary to study diffusion of point-defects in presence of interfaces. Model potential molecular dynamics (MPMD)[30], on the other hand, is suitable for the study of dynamical properties and the microstructure evolution of large atomistic models [43, 44, 45, 46, 47, 48, 49]. MPMD has been applied to study the thermally activated diffusivity of point-defects in crystals showing that iodine vacancies and interstitials are the most mobile defects[11, 50, 51] in absence of electric fields in periodic crystals. In this work, with the aim of understanding the origin of hysteresis, we study by MD the electric mobility of iodine point defects at room conditions both in the bulk and in crystalline samples of finite thickness. A reversible kinetics of native point-defects is found in the bulk where we exclude ferroelectric polarization and global molecular ordering at room conditions. Conversely, we show that surfaces (either polar or compensated) perturb the lattice and, independently on their termination, act as trapping layers with fast accumulation and slow release kinetics of both vacancy and interstitial point defects. Accordingly, history-dependent ionic currents are expected in the material in presence of surfaces. Considering that charged defects are able to screen externally applied voltages and contribute to ionic currents, we conclude that the asymmetry in accumulation and release of defects at surfaces (or interfaces) with long release times is at the origin of the observed hysteresis, as also suggested by recent experimental findings[37]. Present results are complemented by an analytic model of the effect of surfaces on point defects and of the kinetics of ionic polarization for the design of functional films with better control of hysteretic ionic currents.

2. Results and discussion 2.1. Response of the crystalline bulk to electric fields We first consider a periodic crystal and its response to electric field. The objective of present section is to establish whether transient currents and memory effects are possible in the perfect bulk. We consider in this section two possible sources of hysteresis: (i) the occurrence of bulk macroscopic polarization; (ii) the mobility of ionic defects. Concerning the former point, it has been extensively anticipated in literature[45, 31] that MA dipoles undergo fast rotations in the picosecond time scale excluding the possibility of a ferroelectric state with complete ordering of molecular dipoles. Here we directly simulate the crystal under electric fields, up to very large values of 2 MV/cm. In response to the electric field a finite polarization arises (i.e. a dielectric response) consisting of ionic displacements and partial orientation of dipoles. In order to study whether such a linear polarization can be associated to hysteresis we considered alternating electric field periodically switched between +E and −E , with periods in the range 0.1 − 10 ns. At each inversion, the total bulk polarization is inverted with transient times as short as a few picoseconds and this is consistent with the angular correlation time (i.e. reorientational time) of methylammonium molecules at room temperature (

1 − 10 ps)[30]. Accordingly, polarization transients are nine orders of magnitude faster than currentvoltage hysteresis (picoseconds vs. milliseconds) Figure 1: Defect drift velocity against applied electric field. Symbols: calculated velocities; lines: model of Eq.1, with Γ of Tab.1. Having clarified that bulk molecular polarization cannot explain slow current transients and memory effects, we focus on the response to electric field of bulk crystals containing ionic point-defects. The ionic mobility μ of a specific defect in MAPI can be directly calculated by applying an electric field E to a periodic crystal containing the defect and calculating the average defect velocity v (i.e. drift) along the field direction. In linear regime and for a single charge species it is expected

v = μE as derived from the fact that ionic current density j satisfies j = σ E = nq μ E on one hand and j = nqv on the other, where σ is the ionic conductivity, n is the defect concentration and q its charge. The fields E have magnitude ranging from 0.5 to 2 MV/cm. Smaller values of E can be applied but the larger statistical errors at small fields would require longer simulations. Figure 1 shows the iodine vacancy and interstitial velocities calculated by MD as a function of E . The drift velocities can be described over the whole range of E by using a nonlinear model based on rate theory[52]

⎡ ⎛⎜ d0 qE ⎞⎟ ⎛⎜ − d0 qE ⎞⎟ ⎤ 2k T 2k T v = d 0 f + − d 0 f − = d 0 ⋅ Γ ⋅ ⎢ e⎝ B ⎠ − e⎝ B ⎠ ⎥ = ⎢ ⎥ ⎣ ⎦ (1) ⎛ d qE ⎞ 2d 0 ⋅ Γ ⋅ sinh ⎜ 0 ⎟ ⎝ 2 k BT ⎠

Γ is the frequency of jumps between neighboring sites at distance d 0 in absence of electric fields; f + and f − are the probability to jump in a direction parallel and antiparallel to the applied electric field, respectively; qEd0 is the energy difference between neighboring sites. At zero electric field the

velocity is zero as f + and f − are equal. At small electric fields sinh ( E )  E , the velocity is linearly

dependent on the bias and the mobility becomes a constant proportional to diffusivity

v (E )

qd 02 Γ qD μ= ≈ = E E → 0 k BT k BT

(2)

At large electric fields the velocity and the mobility grow exponentially with E , i.e.

v ≈ d0 ⋅ Γ ⋅ e

⎛ d 0 qE ⎞ ⎜ ⎟ ⎝ 2 k BT ⎠

. The fit of atomistic data by the model is reported in Fig.1; fitted Γ and the mobility at small fields are reported in Table 1. Present calculations provide very high mobility of ionic defects consistent with previous results based on ionic diffusivity[11]. The very high mobility of defects in the bulk demonstrates that ionic defects are important in the overall conductivity of the material but it does not explain alone the occurrence of current transients and hysteresis in the bulk due to ionic polarization. As a matter of fact, by switching the electric field between +E and −E the ionic current is inverted almost suddenly (within tens of picoseconds) without showing long relaxation times comparable with timescales of observed hysteresis (at least milliseconds). We conclude that in the ideal bulk (without surfaces) neither molecular polarization nor ionic currents due to point-defects can explain hysteresis. Table 1: Jump frequencies Γ and defects mobilities μ for vacancies and interstitials.

Γ (s −1 ) μ (cm 2 /V ⋅ s) V 1.6 ⋅1010

4.7 ⋅10−4

I

2.6 ⋅10−4

1.0 ⋅1010

2.2. Effect of free surfaces Figure 2: Left: system under study (top) and energy profiles (bottom) after that the macroscopic polarization field has been removed. Right: schematic representation of the three different zones in the samples and the energy profile under external bias Vext . The above scenario for ionic currents is dramatically altered by the presence of surfaces. Surfaces perturb the bonding structure of the crystal giving rise to electronic defects and surface states; furthermore, surfaces can induce long-range effects such as lattice strain or electric fields. The latter are expected when there is a net polarization charge at the surfaces as in case of polarized surfaces of ionic crystals (such as MAPI). Electric and strain effects are long-ranged and can have important implications on the energetics and diffusion of point-defects. To study the above effects we consider a crystal slab of size ∼ 3.5x3.5x25 nm 3 , with (001) surfaces (PbI 2 -terminated on one end and MAI-terminated on the other one) orthogonal to the z direction (see top left panel of Fig.2). We choose the (001) surface that is considered the most likely cut in real samples as indicated by powder X-ray diffraction measurements[53] and by theoretical calculations based on present classical force-field[46] and first-principles[54, 55]. The effect of surfaces on point defects can be understood by calculating the energy profile of an iodine vacancy (or interstitial) placed at different distance from the surfaces (see bottom left panel of Fig.2). At each distance the atomic forces and positions are fully optimized and the potential energy of the system is calculated. The energy profile for the vacancy (interstitial) increases (decreases) linearly with the z coordinate. Considering that the vacancy (interstitial) is positively (negatively) charged, the linear behavior reveals a constant electric field within the slabs. Though the nominal charge of both surfaces (PbI 2 and MAI) is zero, there is a charge unbalance between negative PbI 2 and positive MAI planes. This difference is found to be 0.25e per surface formula unit (e is the electronic charge) by ab initio calculations[56] and it is well reproduced by the present force field providing 0.23e . These charge values correspond to a residual polarization in the film of the order of

2.65μ C/cm 2 . Notice that such a polarization is related to surface polarity and it cannot be switched by the application of external electric fields. Polarized surfaces can in principle contribute to hysteresis in MAPI since the velocity of charge carriers is enhanced (reduced) when the current is parallel (antiparallel) to the internal electric field. However, polarized surfaces are energetically unfavored in real samples as the potential energy increases with the slab thickness. Already for nanometer-thick layers several phenomena tend to compensate the surface charge (screening due to mobile charges, metallization, reconstructions or Zener breakdown[57]). Considering that hybrid perovskites contain large concentrations of mobile defects (either ionic or electronic) the surface polarity is likely compensated in real layers. The concentration of point defects necessary to compensate the surfaces is calculated in Appendix ( ∼ 1013 cm −2 ) and it is comparable with the density of defects expected in real conditions. When the superficial charge is fully compensated, the surfaces have nevertheless sizable interaction with point-defects. The surface-defect interaction can be extracted from the atomistic energy profile once the linear dependence on z (i.e. the potential of the internal electric field) is subtracted. This corresponds to cancel the long-range electric field due to the slab terminations. The result, labeled as U 0 ( z ) in Fig.2 for vacancies (red) and interstitials (blue), shows a non monotonic dependence of the energy on the defect-surface distance that is roughly symmetric with respect to the slab center. Two

maxima occur close to the surfaces and two deep wells at the boundaries where point-defects can be trapped. Notice that both point-defects are attracted at each surface independently on their charge state (negative interstitials and positive vacancies). This is a clear signature that the surface-defect interaction is not due to electrostatic interactions but a more complex combination of strain, dangling bonds, molecular and ionic rearrangements. For example, on the one hand, the formation energy of defects is smaller on the surface since atoms are undercoordinated and strain is higher; on the other hand the electrostatic energy is lower in the bulk where the defect charge is more efficiently screened. These opposite factors give rise to the energy maxima observed. Starting from each surface we can identify three regions in the energy profiles (see Fig.2, left, for the interstitial case): the outmost superficial region S, (consisting of a well of width LS and height E A ); a

j ) characterized by a smooth decrease of the potential subsequent superficial region S (of length L S energy; a bulk-like region B (of length LB ) with approximately flat energy EB . We consider to be in the B region when the value of the potential energy has dropped by ∼ 90% of its maximum value. The overall energy profile U 0 ( z ) of the slab can be fitted by the convolution of two curves (one at each of the surfaces) of functional form ⎛ z − z0 ⎞ λ ⎟⎠

⎛ z − z0 ⎞ ⎜− U 0 ( z ) = qV0 ( z ) = E0 ⋅ ⎜ + 1 ⎟ ⋅ e⎝ ⎝ λ ⎠

+ EB

(3)

corresponding to the universal relation (UER)[58, 59]. z0 is the position of the maximum. The energy barrier E A is related to E0 through the relation E A = E0 − U ( z S ) , where z S is the coordinate of the

surface. The UER reproduces the cohesive energy dependence on interatomic distance in most materials but it is also useful for surface-related properties (e.g. surface-surface energy[60] or molecule-molecule interactions on surfaces[61]).

j are in nm. Table 2: Fitting parameters for eq.3. E0 , E A and EB are in eV, LS , z0 and L S MAI surface

PbI 2 surface

E0

EA

j L S

z0

LS

E0

EA

j L S

z0

LS

EB

V

0.23

0.37

5.92

1.53

1.58

0.51

0.56

8.09

22.37

2.02

0.13

I

0.44

0.52

6.38

1.83

1.70

0.27

0.49

5.23

22.27

1.62

0.22

The parameter λ controls the curvature close to the maximum and the width of LS : moving from z0

⎛5⎞ by LS = 4λ , the energy decreases by a factor ⎜ 4 ⎟  90% and approaches the asymptotic bulk ⎝e ⎠

value. The fitting provides estimate of all parameters (see Table 2) of the MD energy profile and in j are ∼ 2 nm and ∼ 6 nm, respectively, for both surfaces (MAI and PbI ) and particular LS and L 2 S

j are superficial properties of MAPI almost defect type (vacancies and interstitials). LS and L S independent on the sample size. Conversely, the extension of the bulk-like region, LB , scales with the

(

)

j . We have verified that the defect-surface sample length L, being equal to L − 2 ⋅ LS + L S interaction calculated above is independent on the terminations of the slab (either MAI-PbI 2 or PbI 2 PbI 2 or MAI-MAI), once the long-range potential due to the corresponding surface charges are removed.

The importance of the above energetic analysis is to provide information on accumulation and release of defects on/from the surfaces. In Fig.2, right, the effect of external bias Vext is schematically depicted during accumulation (top) and release (bottom) of point-defects at one surface. The represented quantity is the total potential, i.e. the sum of internal U 0 ( z ) / q and external potential

Vext ( z ) . In the bulk-like region, where the internal energy is flat, the defects are driven by the external potential V ( z ) = V ( z ) = −E ⋅ z . In the superficial regions, i.e. S and S , the internal ext

contributions are dominant with respect to external ones. Typical external bias are of the order of 1 V/μm = 1 mV/nm which are smaller than the internal ones being ∼ 50 mV/nm. Experiments show that when an external bias is applied (removed) the potential of the slab gives rise to current (and potential) transients driving the system to stationarity. We can explain the origin of the transient by considering the bulk region of the slab as a capacitor of capacitance C = ε

A and LB

−1

⎛ A⎞ resistance R = ⎜ σ ⎟ , where A is the surface area: when a bias is applied to the slab the defects ⎝ LB ⎠ accumulate from the bulk towards the boundaries (regions S and S) and the capacitor undergoes V di i charging. The resistance-capacitance model satisfying equation + = ext δ ( t − t0 ) (where dt RC R δ ( t ) is the delta-function and t0 is the switch-on time) predicts an exponential transient current

Vext −(t −t0 ) / RC (where θ ( t ) is the Heaviside step function) and voltage V ( t ) = Ri ( t ) e R with time of charging/discharging RC corresponding to the resistance-capacitance product; for the i ( t ) = θ ( t − t0 )

MAPI slab the equivalent time is

tB =

ε

σ

(4)

that is, in principle, size independent. This is expressed in terms of the dielectric constant and defects conductivity in MAPI. By using ε  62 ⋅ ε0 [62] and σ = nqμ , with μ in the range 10 −5 − 10−4 cm 2 V −1

s −1 and n in the range 1014 − 1015 cm

−3

[50], the timescale for charging the slab turns out to be of

the order of 0.1 − 10 ms, consistent with experimental observations. When the bias is turned off, the defects are released from the surfaces and the capacitor undergoes discharging with a time constant that is still RC . If the material and its ionic conductivity are reversible after the charging, then the accumulation and release times are expected to be the same. The calculated potential at the superficial regions suggests however a different scenario from the simple resistance-capacitance model, in particular for the release kinetics as it will be shown below. During the charging stage, as soon as the defects accumulate close to the superficial layers, they can

diffuse through S (see Fig.2) and reach the surface. The time for diffusing to the surface tS turns out to be as short as few nanoseconds, that is much smaller than charging time t B . This time can be estimated as

j2 L tS = S D

(5)

where D is the local diffusion coefficient. Though the energy profile at the superficial region reduces the effective diffusivity D with respect to the equilibrium bulk diffusivity D,

D = D ⋅ e

−

d U 0 d0 dz k BT

qVext d 0 k BT L

e

 D⋅e

−

d U0 d0 dz

(6)

a very fast diffusion to surface is actually found. This due to the bias and the nanometric extension of the S region. Once defects have overcome the maximum they will spontaneously fall into the superficial well in a time tS ∼

q ⋅ L2S < 1 ns. μ EA

Let us consider now the release of defects (see Fig2 right bottom) from the surface. An inverse field E > 0 now points towards the bulk zone (see lower panel of Fig. 2). For the discharging to occur it is necessary first that the trapped defects in region S overcome a large barrier E A ( ∼ 0.5 eV) and drift through S to the bulk. According to the transition state theory the release time per defect is ⎛ EA ⎞

1 ⎜k T ⎟ t'S = ⋅ e⎝ B ⎠ (7) Γ where Γ is the attempt frequency calculated in the previous section. Defects that have overcome the

barrier will drift through S region towards the bulk driven by both the external and internal potentials ( Vext + V0 ) in a short timescale of

j2 L S t' S = μ (Vext + V0 )

(8)

Once that defects have detrapped from the surface, the discharge can occur within the transient time t' B = t B . In conclusion, the overall accumulation τ a and release τ r times for an idealized crystalline slab can be estimated as

τ a = tS + tS + tB τ r = t'S + t'S + tB

(9)

where each term in the sums can be calculated through Eq.s 4-8 using data obtained from atomistic j , the energy barrier E , the bulk diffusivity D and mobility μ) and simulations (i.e. lengths LS , L A S from experiments (i.e. the external bias Vext and the layer thickness L). As discussed above, the leading term in accumulation is the capacitor charging time t B , while for defects release the transition state time is larger or comparable to discharging:

τ a ∼ εσ −1 1 Γ

τ r ∼ εσ −1 + ⋅ e E

A / k BT

(10)

The times τ r and τ a calculated for the four defect-surface combinations are reported in Table 3. It is found that accumulation is always faster than release. The slowest kinetics is that of vacancies at MAI surfaces and of interstitials at PbI 2 surfaces, for which release times are of the order of hundred milliseconds, while the remaining kinetics are faster. We point out that, according to the model, we expect exponential kinetics for both accumulation and release. In fact, the accumulation is controlled by the RC time and charge density transients are exponential

na ( t ) ∼ na ( 0 ) ⋅ e

−

t

τa

On the other hand the release time τ r is obtained within a transition state theory and the defect concentration n ( t ) (per unit volume) leaving the surface in the interval dt is dn = −

1

τr

n ( t ) dt

giving

nr ( t ) ∼ nr ( 0 ) ⋅ e

−

t

τr

i.e. an exponential transient, too. The above analysis based on static atomistic calculations is supported by finite-temperature dynamics of defects in the MAPI slab. To this aim we have first compensated the slab surfaces by introducing an equal number of negative interstitials and positive vacancies at the surfaces with the aim to minimize the internal electric field. An additional single point defect (interstitial or vacancy) has been introduced in the center of the slab and its drift-diffusion simulated at room temperature under the effect of external electric fields. Consistent with the static analysis we found that the defects easily migrate to the surfaces where they can be trapped in a time as short as few nanoseconds. In Fig.3 we report three snapshots of an iodine interstitial moving from the bulk (left), close to the MAI surface (middle) and trapped on the surface (right). The short time for accumulating the defect at the surface is consistent with the result of the static model tiS + tS if we consider that the time t B due to capacitive effects requires interaction between positive and negative ions and is not included when simulating a single defect. Once the defects were trapped on the surfaces, it was not possible to observe their release from the surface into the slab in the time affordable by MD ( ∼ 100 ns), even by inverting the electric field (from E to −E ). This is consistent with the static analysis predicting much longer release times ( ∼ 100 milliseconds). Similar results were obtained for slabs with different compensated terminations (e.g. PbI2-PbI2 or MAI-MAI). Figure 3: Snapshots taken during a MPMD simulation of an interstitial point defect, which migrates from the bulk (left), close to the MAI surface (center) and trapped on the same surface (right) in a few ns. Once that the defect reaches the surface, it remains there during the remaining simulation time. It is notable that this simple slab model (without interfaces) already catches the physical origin of the long current transients. As a matter of fact, it is the asymmetry between accumulation and release of defects, together with the long release times, the physical key explaining current hysteresis observed experimentally. At the origin of the phenomenon there is the energy gain of defects forming or trapped at the surfaces; in presence of residual surface polarity the effects could be also enhanced; conversely the bulk-like regions do not trap defects and do not contribute to slow and asymmetric current transients. Table 3: Estimated time needed for point defects to accumulate from bulk to surface ( τ a ) and released from surface to bulk ( τ r ); conductivity σ = 4.7 ⋅106 ( Ωm )

−1

has been estimated by using

−3

n = 10 cm . 15

MAI

PbI 2

τa

τr

τa

τr

V

0.07 ms

0.55 ms

0.07 ms

319 ms

I

0.12 ms

110 ms

0.12 ms

35 ms

Figure 4: Top: schematic representation of charged defects migration in response to switching on and off an external field; arrows represent the ionic current. Bottom: comparison between experimental and calculated variation of the CPD signal after switching on ( ton ) and off ( toff ) the

external voltage; the shaded area covers the decay times between the highest values of τ r for interstitial (blue) and vacancy (red). Zones A, B, C, D are explained in detail in the main text. The calculated release time, τ r , can be compared with time-resolved Kelvin Probe Force Microscopy (tr-KPFM)[63, 64, 65, 37, 66] measurements. In particular, very recently, Weber et al.[37] have recorded the temporal evolution of the local contact potential difference (CPD) at a position close to interface between perovskite and oxide upon the switching on ( ton ) and off ( toff ) of a voltage or light pulse. The experimental CPD data[37] are reported in bottom panel of Fig.4 (green line) where we have labeled the different time stages (A equilibrium, B bias, C unbias, D equilibrium). Top panels (A-D) represent schematically the expected behavior of ionic defects during the corresponding stage. At equilibrium (panel A) point defects are distributed inside the MAPI bulk-like regions. The corresponding CPD potential is zero. When the external bias V1 is turned on (panel B), the defects are separated and trapped at the surfaces within time τ a . The surface charge screens the external bias and it is observed a CPD sharp transient from V 1 to V 2 (stage B, t > ton ). Upon switching off the bias (C), the surface charge generates a response CPD potential opposite to the initial bias (V 2 -V 1 ); the CPD potential exponentially decays slowly back to zero ( t > toff ) as the charge is released from surfaces. Consistent with theoretical prediction the voltage overshoot is a single exponential. The experimental decay time τ r =125.7 ms is in the range of calculated release time τ r = 110 -319 ms (see shadowed region). Considering that the atomistic calculations focus on ideal surfaces and not on specific interfaces, the agreement is remarkable. Theory predicts the correct ordering τ r > τ a with release much slower then accumulation; furthermore only one defect species per surface type (vacancies at MAI and interstitials at PbI 2 surfaces) is expected to control the decay times and this is consistent with indications from recent experimental works[67, 37]. We conclude present study by showing explicitly how the slow release of defects from the surfaces can explain hysteresis. To this aim we observe that an ionic charge Q localized at the surface of area A corresponds to a (ionic) polarization P =

Q . The dynamics of P ( t ) can be derived from the A

following analytical model

dP P0 ( E ) − P ( t ) = dt τ

(11)

Q0 is the polarization driven by the external bias (i.e. the charge density induced at A surfaces by the field E when stationary conditions are achieved, dP / dt = 0 ), P ( t ) is the

where: P0 =

instantaneous polarization, and τ is the relaxation time that controls the time evolution of P towards target P0 . The differential equation states that the current (i.e.

dP ) is inversely proportional to the dt

relaxation time and directly proportional to the unbalance between instantaneous and equilibrium values of charge at the surfaces. To solve the model we assume that P0 has linear dependence on the electric field P0 = K ⋅ E . The actual value of the constant K is unnecessary here though it can be related to the ionic susceptivity χion through the relation K = ε0

χ ion

χ ion + 1

valid for dielectrics. We

finally observe that the relaxation time τ is different depending whether we consider accumulation ( P ( t ) < P0 ( E ) ) or release ( P ( t ) > P0 ( E ) ):

⎪⎧τ if P0 ( E ) − P ( t ) > 0 τ =⎨ a ⎪⎩τ r if P0 ( E ) − P ( t ) < 0 The two-times model is solved numerically for the typical electric field E ( t ) used in experiments, i.e. a periodic signal that varies linearly from zero to the maximum value and vice versa during each period τ s (scan time). Figure 5: Calculated instantaneous polarization P ( t ) as a function of the applied electric field for different values of accumulation τ a and release τ r times. The times are normalized to the scan time

τ s , while the electric field is normalized to its maximum value. In Figure 5 we report the calculated polarization P as a function of E for different values of accumulation τ a and release τ r times. When both τ a , τ r are small with respect to the scan time (bottom panel τ a = τ r = 0.001τ s ) no hysteresis is observed and polarization is a single-valued function of E . This is the behavior expected in the ideal bulk where defects cannot be trapped. As soon as the release time becomes longer (e.g. τ r = 0.05τ s ) there appear closed loops in the P- E figure. The loop area further increases for larger release time τ r = 0.5τ s . This is actually the typical experimental conditions for which scan times are of the order of seconds and release times from surfaces are hundreds of milliseconds. This model is a direct theoretical evidence that the slow release of defects from surfaces (or interfaces) can explain the onset of hysteretic figures of the polarization. The effect of accumulation time τ a is also reported in the top panels of Figure 5. When also accumulation is long and comparable to τ s (for example this can be the case of specific interfaces), the formation of open loops in the P- E plane can be observed.

3. Conclusions By model potential molecular dynamics, we were able to describe the effect of surfaces on the electric mobility of iodine point-defects in MAPI, showing that experimentally observed memory effects cannot be ascribed to bulk and can originate from surface or interface effects. Mobile point-defects are easily trapped on the surfaces due to the spontaneous occurrence of superficial energy wells. In particular, the calculated times for release of vacancies from MAI-terminated and of interstitials from PbI 2 -terminated surfaces are higher than for accumulation and in the hundreds of millisecond timescale, in good agreement with recent tr-KPFM experiments. We propose that the observed hysteresis in perovskite solar cells originates from trapping of iodine defects at surfaces and/or interfaces, which give rise to an electric field able to screen externally applied voltages. We further confirm that no ferroelectric polarization is present in the bulk as the ionic current can be reversed in the picosecond time scale by reversing the applied electric field. By clarifying and modeling the asymmetric interaction of defects with surfaces, present study can help the design of novel functional interfaces and nanodevices based on ionic transport.

4. Computational details Model potential molecular dynamics simulations were performed by using the LAMMPS code [68]. MAPI interatomic forces are described by the MYP potential[30]. The long-range electrostatic forces were calculated by the particle-particle particle-mesh Ewald algorithm [69], and van der Waals interactions were cutoff at 10 Å. The velocity-Verlet algorithm [70] with a time step as small as 0.5 fs was used to solve the equations of motion. Trajectories were analyzed by using the VMD 1.9 code [71]. After a conjugate gradient minimization and a short (5 ps) low temperature (1K) relax, the systems were heated up to the temperatures of interest. Temperature was controlled by the NoseHoover thermostat[72, 73].

The velocity vz of the iodine defects under constant electric field E is calculated from all the iodine

velocities vi , z ( t ) and total current

∑v ( t ) q as follows i,z

i

vz = ∫dt ∑vi , z ( t )

(12)

i

The average velocity vz also corresponds to the time derivative of the polarization of the defected slab vz =

A ⋅ L dPz ( t ) under bias. q dt

5. Acknowledgements We acknowledge MIUR for funding through project PON04a2_00490 M2M Netergit and PRACE for awarding access to Marconi KNL at CINECA, Italy, through projects UNWRAP (2016153664) and DECONVOLVES (2018184466)

Appendix A. Charge compensation We have previously observed that a macroscopic polarization field, resulting from the cut of the surfaces, will be present in MAPI samples. This is true at zero temperature, but also at room temperature. This macroscopic field can be calculated by placing a charged defect ( q = ±1.13 e − ) at varying distance from the surface, and observing the energy variation ΔU along the sample length L after relaxing. Given ΔU ∼ 8.9 eV, the polarization field EP can be estimated as

EP =

ΔU V = 3 ⋅10−2 = 3 MV / cm (A1) q⋅L Å

The resulting interaction basin is reported for vacancies in Fig.A1 (black filled circles). Figure A1: Residual field upon displacement of iodine ions from the PbI 2 to the MAI surface (legend is in 10 13 defects/cm 2 ); the curves have been transalted for better visualization. The slab has been relaxed and equilibrated at room temperature and pressure for 1 ns before calculating the interaction basins. Due to this charge unbalance, charged defects are expected to migrate towards the oppositely charged surfaces. Furthermore, in MAPI solar cells there are not free surfaces but rather interfaces with hole/electron transport layers, which can in turn be a source of mobile charges, which will migrate until they compensate the electric fields. Fig. A1 shows also the residual field upon displacement of up to four iodine ions from the PbI 2 to the MAI surface. The intrinsic field decreases upon displacement of up to three iodine ions, while displacement of a fourth ion has the effect to invert the field. Displacement of atoms from one surface to another corresponds to the introduction of a surface defect density, which in our case goes from zero to 3.18 ⋅1013 cm −2 , and a surface charge going from zero to 5.76 μC/cm −2 . By fitting the data at different defect concentration, we estimate that the intrinsic electric field is expected to be fully compensated with a surface defect density of ∼ 2.6 ⋅1013 cm −2 , corresponding to a polarization of σ c = ± 4.7 μC/cm −2 . For typical grain sizes of the order of 100-200 nm, this corresponds to a volumetric density of defects around 1018 cm −3 . From the same analysis at T=0 K we obtain σ c = ± 2.65 μC/cm 2 . Being equal to the surface charge density needed to compensate the intrinsic field, this value is, with opposite sign, the surface charge present on MAPI free (001) surfaces, and it is in very good agreement with the values calculated ab initio at T=0 K for the same surface terminations[56].

References

[1] A. K. Jena, A. Kulkarni, T. Miyasaka, Halide perovskite photovoltaics: background, status, and future prospects, Chemical reviews 119 (2019) 3036–3103. doi:10.1021/acs.chemrev.8b00539. [2] W.-J. Yin, T. Shi, Y. Yan, Unusual defect physics in ch3nh3pbi3 perovskite solar cell absorber, Appl. Phys. Lett. 104 (2014) 063903. URL: https://doi.org/10.1063/1.4864778. doi:10.1063/1.4864778. arXiv:https://doi.org/10.1063/1.4864778. [3] J. Kim, S.-H. Lee, J. H. Lee, K.-H. Hong, The role of intrinsic defects in methylammonium lead iodide perovskite, J. Phys. Chem. Lett. 5 (2014) 1312–1317. URL: https://doi.org/10.1021/jz500370k. doi:10.1021/jz500370k. arXiv:https://doi.org/10.1021/jz500370k, pMID: 26269973. [4] M. H. Du, Efficient carrier transport in halide perovskites: theoretical perspectives, J. Mater. Chem. A 2 (2014) 9091–9098. URL: http://dx.doi.org/10.1039/C4TA01198H. doi:10.1039/C4TA01198H. [5] M.-H. Du, Density functional calculations of native defects in ch3nh3pbi3: Effects of spin-orbit coupling and self-interaction error, J. Phys. Chem. Lett. 6 (2015) 1461–1466. URL: https://doi.org/10.1021/acs.jpclett.5b00199. doi:10.1021/acs.jpclett.5b00199. arXiv:https://doi.org/10.1021/acs.jpclett.5b00199, pMID: 26263152. [6] D. Meggiolaro, S. G. Motti, E. Mosconi, A. J. Barker, J. Ball, C. A. R. Perini, F. Deschler, A. Petrozza, F. De Angelis, Iodine chemistry determines the defect tolerance of lead-halide perovskites, Energy Environ. Sci. 11 (2018) 702–713. [7] A. Filippetti, A. Mattoni, Hybrid perovskites for photovoltaics: Insights from first principles, Phys. Rev. B 89 (2014) 125203. URL: https://link.aps.org/doi/10.1103/PhysRevB.89.125203. doi:10.1103/PhysRevB.89.125203. [8] J. Burschka, N. Pellet, S.-J. Moon, R. Humphry-Baker, P. Gao, M. K. Nazeeruddin, M. Gratzel, Sequential deposition as a route to high-performance perovskite-sensitized solar cells, Nature 499 (2013) 316–319. URL: http://dx.doi.org/10.1038/nature12340. doi:10.1038/nature12340. [9] A. Polman, M. Knight, E. C. Garnett, B. Ehrler, W. C. Sinke, Photovoltaic Materials: Present Efficiencies and Future Challenges, Science (80-. ). 352 (2016) aad4424–aad4424. URL: http://www.sciencemag.org/cgi/doi/10.1126/science.aad4424. doi:10.1126/science.aad4424.

[10] M. I. Saidaminov, V. Adinolfi, R. Comin, A. L. Abdelhady, W. Peng, I. Dursun, M. Yuan, S. Hoogland, E. H. Sargent, O. M. Bakr, Planar-integrated single-crystalline perovskite photodetectors, Nat. Commun. 6 (2015) 8724. doi:https://doi.org/10.1038/ncomms9724. [11] P. Delugas, C. Caddeo, A. Filippetti, A. Mattoni, Thermally activated point defect diffusion in methylammonium lead trihalide: Anisotropic and ultrahigh mobility of iodine, J. Phys. Chem. Lett. 7 (2016) 2356–2361. URL: https://doi.org/10.1021/acs.jpclett.6b00963. doi:10.1021/acs.jpclett.6b00963. arXiv:https://doi.org/10.1021/acs.jpclett.6b00963, pMID: 27237630. [12] T.-Y. Yang, G. Gregori, N. Pellet, M. Grätzel, J. Maier, The significance of ion conduction in a hybrid organicÄìinorganic lead-iodide-based perovskite photosensitizer, Angew. Chem., Int. Ed. 54 (2015) 7905–7910. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/anie.201500014. doi:10.1002/anie.201500014. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/anie.201500014 . [13] Y. Yuan, J. Chae, Y. Shao, Q. Wang, Z. Xiao, A. Centrone, J. Huang, Photovoltaic switching mechanism in lateral structure hybrid perovskite solar cells, Adv. Energy Mat. 5 (2015) 1500615. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/aenm.201500615. doi:10.1002/aenm.201500615. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/aenm.201500615 . [14] J. Bisquert, G. Garcia-Belmonte, A. Guerrero, Impedance characteristics of hybrid organometal halide perovskite solar cells, in: Organic-Inorganic Halide Perovskite Photovoltaics, Springer, 2016, pp. 163–199. [15] Z. Xiao, Y. Yuan, Y. Shao, Q. Wang, Q. Dong, C. Bi, A. Sharma, P.and Gruverman, J. Huang, Giant switchable photovoltaic effect in organometal trihalide perovskite devices, Nat. Mater. 14 (2015) 193–198. doi:10.1038/nmat4150. [16] H. J. Snaith, A. Abate, J. M. Ball, G. E. Eperon, T. Leijtens, N. K. Noel, S. D. Stranks, J. T.-W. Wang, K. Wojciechowski, W. Zhang, Anomalous hysteresis in perovskite solar cells, J. Phys. Chem. Lett. 5 (2014) 1511–1515. URL: https://doi.org/10.1021/jz500113x. doi:10.1021/jz500113x. arXiv:https://doi.org/10.1021/jz500113x, pMID: 26270088.

[17] Y. Shao, Z. Xiao, C. Bi, Y. Yuan, J. Huang, Origin and elimination of photocurrent hysteresis by fullerene passivation in ch3nh3pbi3 planar heterojunction solar cells, Nat. Commun. 5 (2014) 5784. doi:10.1038/ncomms6784. [18] B. Chen, M. Yang, S. Priya, K. Zhu, Origin of j-v hysteresis in perovskite solar cells, J. Phys. Chem. Lett. 7 (2016) 905–917. URL: https://doi.org/10.1021/acs.jpclett.6b00215. doi:10.1021/acs.jpclett.6b00215. arXiv:https://doi.org/10.1021/acs.jpclett.6b00215, pMID: 26886052. [19] T. Supasai, N. Rujisamphan, K. Ullrich, A. Chemseddine, T. Dittrich, Formation of a passivating ch3nh3pbi3/pbi2 interface during moderate heating of ch3nh3pbi3 layers, Appl. Phys. Lett. 103 (2013) 183906. URL: https://doi.org/10.1063/1.4826116. doi:10.1063/1.4826116. arXiv:https://doi.org/10.1063/1.4826116. [20] A. Dualeh, N. Tétreault, T. Moehl, P. Gao, M. K. Nazeeruddin, M. Grätzel, Effect of annealing temperature on film morphology of organicÄìinorganic hybrid pervoskite solid-state solar cells, Adv. Funct. Mater. 24 (2014) 3250–3258. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/adfm.201304022. doi:10.1002/adfm.201304022. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/adfm.201304022 . [21] J. Xing, Q. Wang, Q. Dong, Y. Yuan, Y. Fang, J. Huang, Ultrafast ion migration in hybrid perovskite polycrystalline thin films under light and suppression in single crystals, Phys. Chem. Chem. Phys. 18 (2016) 30484–30490. URL: http://dx.doi.org/10.1039/C6CP06496E. doi:10.1039/C6CP06496E. [22] C. Ran, J. Xu, W. Gao, C. Huang, S. Dou, Defects in metal triiodide perovskite materials towards high-performance solar cells: origin, impact, characterization, and engineering, Chem. Soc. Rev. 47 (2018) 4581–4610. URL: http://dx.doi.org/10.1039/C7CS00868F. doi:10.1039/C7CS00868F. [23] E. J. Yoo, M. Lyu, J.-H. Yun, C. J. Kang, Y. J. Choi, L. Wang, Resistive switching behavior in organicÄìinorganic hybrid ch3nh3pbi3?xclx perovskite for resistive random access memory devices, Adv. Mater. 27 (2015) 6170–6175. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/adma.201502889. doi:10.1002/adma.201502889. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.201502889 .

[24] C. Gu, J.-S. Lee, Flexible hybrid organic-inorganic perovskite memory, ACS Nano 10 (2016) 5413–5418. URL: https://doi.org/10.1021/acsnano.6b01643. doi:10.1021/acsnano.6b01643. arXiv:https://doi.org/10.1021/acsnano.6b01643, pMID: 27093096. [25] S. Meloni, T. Moehl, W. Tress, M. Franckevičius, M. Saliba, Y. H. Lee, P. Gao, M. K. Nazeeruddin, S. M. Zakeeruddin, U. Rothlisberger, et al., Ionic polarization-induced current–voltage hysteresis in ch 3 nh 3 pbx 3 perovskite solar cells, Nat. Commun. 7 (2016) 10334. [26] G. Giorgi, J.-I. Fujisawa, H. Segawa, K. Yamashita, Organic–inorganic hybrid lead iodide perovskite featuring zero dipole moment guanidinium cations: a theoretical analysis, The Journal of Physical Chemistry C 119 (2015) 4694–4701. [27] G. Giorgi, K. Yamashita, Zero-dipole molecular organic cations in mixed organic–inorganic halide perovskites: possible chemical solution for the reported anomalous hysteresis in the current– voltage curve measurements, Nanotechnology 26 (2015) 442001. [28] N. De Marco, H. Zhou, Q. Chen, P. Sun, Z. Liu, L. Meng, E.-P. Yao, Y. Liu, A. Schiffer, Y. Yang, Guanidinium: a route to enhanced carrier lifetime and open-circuit voltage in hybrid perovskite solar cells, Nano letters 16 (2016) 1009–1016. [29] D. J. Kubicki, D. Prochowicz, A. Hofstetter, M. Saski, P. Yadav, D. Bi, N. Pellet, J. Lewiński, S. M. Zakeeruddin, M. Grtzel, et al., Formation of stable mixed guanidinium–methylammonium phases with exceptionally long carrier lifetimes for high-efficiency lead iodide-based perovskite photovoltaics, Journal of the American Chemical Society 140 (2018) 3345–3351. [30] A. Mattoni, A. Filippetti, M. I. Saba, P. Delugas, Methylammonium rotational dynamics in lead halide perovskite by classical molecular dynamics: The role of temperature, J. Phys. Chem. C 119 (2015) 17421–17428. URL: http://dx.doi.org/10.1021/acs.jpcc.5b04283. doi:10.1021/acs.jpcc.5b04283. [31] T. Chen, B. J. Foley, B. Ipek, M. Tyagi, J. R. D. Copley, C. M. Brown, J. J. Choi, S.-H. Lee, Rotational dynamics of organic cations in the CH 3 NH 3 PbI 3 perovskite, Phys. Chem. Chem. Phys. 17 (2015) 31278–31286. URL: http://arxiv.org/abs/1506.02205␣http://xlink.rsc.org/?DOI=C5CP05348J . doi:10.1039/C5CP05348J. arXiv:1506.02205. [32] D. Yao, C. Zhang, N. D. Pham, Y. Zhang, V. T. Tiong, A. Du, Q. Shen, G. J. Wilson, H. Wang, Hindered formation of photoinactive δ-fapbi3 phase and hysteresis-free mixed-cation planar heterojunction perovskite solar cells with enhanced efficiency via potassium incorporation, The journal of physical chemistry letters 9 (2018) 2113–2120.

[33] W. Tress, Metal halide perovskites as mixed electronic[start]â[end][start] [end][start] [end][start]Ã[end][start] [end][start]Ã[end][start]¬[end]io nic conductors: Challenges and opportunities[start]â[end][start] [end][start] [end][start]Ã[end][start] [end][start]Ã[end][start]®[en d]from hysteresis to memristivity, J. Phys. Chem. Lett. 8 (2017) 3106–3114. URL: https://doi.org/10.1021/acs.jpclett.7b00975. doi:10.1021/acs.jpclett.7b00975. arXiv:https://doi.org/10.1021/acs.jpclett.7b00975, pMID: 28641009. [34] P. Calado, A. Telford, D. Bryant, X. Li, J. Nelson, B. O’Regan, P. Barnes, Evidence for ion migration in hybrid perovskite solar cells with minimal hysteresis, Nat. Commun. 7 (2016) 13831 EP –. URL: https://doi.org/10.1038/ncomms13831. doi:Article. [35] E. L. Unger, E. T. Hoke, C. D. Bailie, W. H. Nguyen, A. R. Bowring, T. Heumüller, M. G. Christoforo, M. D. McGehee, Hysteresis and transient behavior in current–voltage measurements of hybrid-perovskite absorber solar cells, Energy Environ. Sci. 7 (2014) 3690–3698. [36] P. Yadav, S.-H. Turren-Cruz, D. Prochowicz, M. M. Tavakoli, K. Pandey, S. M. Zakeeruddin, M. Grätzel, A. Hagfeldt, M. Saliba, Elucidation of charge recombination and accumulation mechanism in mixed perovskite solar cells, J. Phys. Chem. C 122 (2018) 15149–15154. URL: https://doi.org/10.1021/acs.jpcc.8b03948. doi:10.1021/acs.jpcc.8b03948. arXiv:https://doi.org/10.1021/acs.jpcc.8b03948. [37] S. A. L. Weber, I. M. Hermes, S.-H. Turren-Cruz, C. Gort, V. W. Bergmann, L. Gilson, A. Hagfeldt, M. Graetzel, W. Tress, R. Berger, How the formation of interfacial charge causes hysteresis in perovskite solar cells, Energy Environ. Sci. 11 (2018) 2404–2413. URL: http://dx.doi.org/10.1039/C8EE01447G. doi:10.1039/C8EE01447G. [38] J. M. Azpiroz, E. Mosconi, J. Bisquert, F. De Angelis, Defect migration in methylammonium lead iodide and its role in perovskite solar cell operation, Energy Environ. Sci. 8 (2015) 2118–2127. [39] D. Meggiolaro, E. Mosconi, F. De Angelis, Modeling the interaction of molecular iodine with mapbi3: a probe of lead-halide perovskites defect chemistry, ACS Energy Lett. 3 (2018) 447–451. [40] D. Meggiolaro, F. De Angelis, First-principles modeling of defects in lead halide perovskites: Best practices and open issues, ACS Energy Lett. 3 (2018) 2206–2222. URL: https://doi.org/10.1021/acsenergylett.8b01212. doi:10.1021/acsenergylett.8b01212. arXiv:https://doi.org/10.1021/acsenergylett.8b01212.

[41] L. Zhang, P. H.-L. Sit, Ab initio study of the dynamics of electron trapping and detrapping processes in the ch3nh3pbi3 perovskite, J. Mater. Chem. A 7 (2019) 2135–2147. URL: http://dx.doi.org/10.1039/C8TA09512D. doi:10.1039/C8TA09512D. [42] H. Uratani, K. Yamashita, Charge carrier trapping at surface defects of perovskite solar cell absorbers: A first-principles study, J. Phys. Chem. Lett. 8 (2017) 742–746. URL: https://doi.org/10.1021/acs.jpclett.7b00055. doi:10.1021/acs.jpclett.7b00055. arXiv:https://doi.org/10.1021/acs.jpclett.7b00055, pMID: 28129504. [43] J. Yu, M. Wang, S. Lin, Probing the soft and nanoductile mechanical nature of single and polycrystalline organic–inorganic hybrid perovskites for flexible functional devices, ACS nano 10 (2016) 11044–11057. [44] M. Wang, S. Lin, Anisotropic and ultralow phonon thermal transport in organic–inorganic hybrid perovskites: Atomistic insights into solar cell thermal management and thermoelectric energy conversion efficiency, Adv. Funct. Mater. 26 (2016) 5297–5306. [45] A. Mattoni, A. Filippetti, C. Caddeo, Modeling hybrid perovskites by molecular dynamics, J. Phys. Condens. Mat. 29 (2017) 043001. URL: http://stacks.iop.org/09538984/29/i=4/a=043001. [46] C. Caddeo, M. I. Saba, S. Meloni, A. Filippetti, A. Mattoni, Collective molecular mechanisms in the ch3nh3pbi3 dissolution by liquid water, ACS Nano 11 (2017) 9183–9190. URL: https://doi.org/10.1021/acsnano.7b04116. doi:10.1021/acsnano.7b04116, pMID: 28783296. [47] C. Caddeo, C. Melis, M. I. Saba, A. Filippetti, L. Colombo, A. Mattoni, Tuning the thermal conductivity of methylammonium lead halide by the molecular substructure, Phys. Chem. Chem. Phys. 18 (2016) 24318–24324. URL: http://dx.doi.org/10.1039/C6CP04246E. doi:10.1039/C6CP04246E. [48] C. Caddeo, D. Marongiu, S. Meloni, A. Filippetti, F. Quochi, M. Saba, A. Mattoni, Hydrophilicity and water contact angle on methylammonium lead iodide, Adv. Mater. Interfaces 6 (2019) 1801173. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/admi.201801173. doi:10.1002/admi.201801173. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/admi.201801173 .

[49] E. Mosconi, P. Salvatori, M. I. Saba, A. Mattoni, S. Bellani, F. Bruni, B. Santiago Gonzalez, M. R. Antognazza, S. Brovelli, G. Lanzani, et al., Surface polarization drives photoinduced charge separation at the p3ht/water interface, ACS Energy Letters 1 (2016) 454–463. [50] D. Barboni, R. A. De Souza, The thermodynamics and kinetics of iodine vacancies in the hybrid perovskite methylammonium lead iodide, Energy Environ. Sci. 11 (2018) 3266–3274. URL: http://dx.doi.org/10.1039/C8EE01697F. doi:10.1039/C8EE01697F. [51] R. A. De Souza, D. Barboni, Iodide-ion conduction in methylammonium lead iodide perovskite: some extraordinary aspects, Chem. Commun. 55 (2019) 1108–1111. [52] B. Bagley, The field dependent mobility of localized electronic carriers, Solid State Commun. 8 (1970) 345 – 348. URL: http://www.sciencedirect.com/science/article/pii/0038109870904643. doi:https://doi.org/10.1016/0038-1098(70)90464-3. [53] T. Baikie, Y. Fang, J. M. Kadro, M. Schreyer, F. Wei, S. G. Mhaisalkar, M. Graetzel, T. J. White, Synthesis and crystal chemistry of the hybrid perovskite (ch3nh3)pbi3 for solid-state sensitised solar cell applications, J. Mater. Chem. A 1 (2013) 5628–5641. URL: http://dx.doi.org/10.1039/C3TA10518K. doi:10.1039/C3TA10518K. [54] J. Haruyama, K. Sodeyama, L. Han, Y. Tateyama, Termination dependence of tetragonal ch3nh3pbi3 surfaces for perovskite solar cells, J. Phys. Chem. Lett. 5 (2014) 2903–2909. URL: xhttp://dx.doi.org/10.1021/jz501510v. doi:10.1021/jz501510v. [55] J. Haruyama, K. Sodeyama, L. Han, Y. Tateyama, Surface properties of ch3nh3pbi3 for perovskite solar cells, Acc. Chem. Res. 49 (2016) 554–561. URL: http://dx.doi.org/10.1021/acs.accounts.5b00452. doi:10.1021/acs.accounts.5b00452. [56] L. She, M. Liu, D. Zhong, Atomic structures of ch3nh3pbi3 (001) surfaces, ACS Nano 10 (2016) 1126–1131. URL: https://doi.org/10.1021/acsnano.5b06420. doi:10.1021/acsnano.5b06420. arXiv:https://doi.org/10.1021/acsnano.5b06420, pMID: 26643387. [57] C. Noguera, Polar oxide surfaces, J. Phys. Condens. Mat. 12 (2000) R367. [58] J. H. Rose, J. R. Smith, F. Guinea, J. Ferrante, Universal features of the equation of state of metals, Phys. Rev. B 29 (1984) 2963–2969. URL: https://link.aps.org/doi/10.1103/PhysRevB.29.2963. doi:10.1103/PhysRevB.29.2963.

[59] A. Mattoni, M. Ippolito, L. Colombo, Atomistic modeling of brittleness in covalent materials, Phys. Rev. B 76 (2007) 224103. URL: https://link.aps.org/doi/10.1103/PhysRevB.76.224103. doi:10.1103/PhysRevB.76.224103. [60] J. H. Rose, J. Ferrante, J. R. Smith, Universal binding energy curves for metals and bimetallic interfaces, Phys. Rev. Lett. 47 (1981) 675. [61] G. Malloci, A. Petrozza, A. Mattoni, Atomistic simulations of thiol-terminated modifiers for hybrid photovoltaic interfaces, Thin Solid Films 560 (2014) 34 – 38. URL: http://www.sciencedirect.com/science/article/pii/S0040609013015150. doi:https://doi.org/10.1016/j.tsf.2013.09.043, european Materials Research Society (E-MRS) Spring Meeting 2013 Symposium B: Organic and hybrid interfaces in excitonic solar cells: from fundamental science to applications. [62] I. Anusca, S. Balčiūnas, P. Gemeiner, Š. Svirskas, M. Sanlialp, G. Lackner, C. Fettkenhauer, J. Belovickis, V. Samulionis, M. Ivanov, et al., Dielectric response: Answer to many questions in the methylammonium lead halide solar cell absorbers, Adv. Energy Mat. 7 (2017) 1700600. [63] R. Berger, A. L. Domanski, S. A. Weber, Electrical characterization of organic solar cell materials based on scanning force microscopy, Eur. Polym. J. 49 (2013) 1907 – 1915. URL: http://www.sciencedirect.com/science/article/pii/S0014305713001213. doi:https://doi.org/10.1016/j.eurpolymj.2013.03.005. [64] V. W. Bergmann, S. A. L. Weber, F. Javier Ramos, M. K. Nazeeruddin, M. Grätzel, D. Li, A. L. Domanski, I. Lieberwirth, S. Ahmad, R. Berger, Real-space observation of unbalanced charge distribution inside a perovskite-sensitized solar cell, Nat. Comm. 5 (2014) 5001. URL: https://doi.org/10.1038/ncomms6001. doi:10.1038/ncomms6001. [65] L. Collins, M. Ahmadi, T. Wu, B. Hu, S. V. Kalinin, S. Jesse, Breaking the time barrier in kelvin probe force microscopy: Fast free force reconstruction using the g-mode platform, ACS Nano 11 (2017) 8717–8729. URL: https://doi.org/10.1021/acsnano.7b02114. doi:10.1021/acsnano.7b02114. arXiv:https://doi.org/10.1021/acsnano.7b02114, pMID: 28780850. [66] Z. Kang, H. Si, M. Shi, C. Xu, W. Fan, S. Ma, A. Kausar, Q. Liao, Z. Zhang, Y. Zhang, Kelvin probe force microscopy for perovskite solar cells, Sci. China Mater. 62 (2019) 776–789. URL: https://doi.org/10.1007/s40843-018-9395-y. doi:10.1007/s40843-0189395-y. [67] L. Bertoluzzi, R. A. Belisle, K. A. Bush, R. Cheacharoen, M. D. McGehee, B. C. O’Regan, In situ measurement of electric-field screening in hysteresis-free ptaa/fa0.83cs0.17pb(i0.83br0.17)3/c60

perovskite solar cells gives an ion mobility of ∼ 3x10 −7 cm 2 /(v s), 2 orders of magnitude faster than reported for metal-oxide-contacted perovskite cells with hysteresis, J. Am. Chem. Soc. 140 (2018) 12775–12784. URL: https://doi.org/10.1021/jacs.8b04405. doi:10.1021/jacs.8b04405. arXiv:https://doi.org/10.1021/jacs.8b04405, pMID: 30189142. [68] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comp. Phys. 117 (1995) 1–19. URL: http://lammps.sandia.gov. [69] W. C. Swope, H. C. Andersen, P. H. Berens, K. R. Wilson, A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters, J. Chem. Phys. 76 (1982) 637–649. [70] L. Verlet, Computer experiments on classical fluids. i. thermodynamical properties of lennardjones molecules, Phys. Rev. 159 (1967) 98–103. doi:10.1103/PhysRev.159.98. [71] W. Humphrey, A. Dalke, K. Schulten, Vmd: Visual molecular dynamics, J. Mol. Graphics 14 (1996) 33 – 38. doi:http://dx.doi.org/10.1016/0263-7855(96)00018-5. [72] S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys. 81 (1984) 511–519. doi:10.1063/1.447334. [73] W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A 31 (1985) 1695–1697. doi:10.1103/PhysRevA.31.1695.

• Hysteresis and memory effects in hybrid perovskites are not bulk phenomena but determined by the effect of surfaces (or interfaces) on ionic dynamics; • Surfaces act as trapping regions for highly mobile iodine point-defects; • the slow release of point-defects from surfaces is at the origin of hysteresis; • Surface-defects interaction is calculated by classical molecular dynamics and fitted by analytical models providing accumulation and release times in agreement with experiments; • Present results can help the design of hysteresis-free perovskites layers;