How spin relaxes and dephases in bulk halide perovskites

Spintronics in halide perovskites has drawn significant attention in recent years, due to their highly tunable spin-orbit fields and intriguing interplay with lattice symmetry. Here, we perform first-principles calculations to determine the spin relaxation time (T1) and ensemble spin dephasing time (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{2}^{*}$$\end{document}T2*) in a prototype halide perovskite, CsPbBr3. To accurately capture spin dephasing in external magnetic fields we determine the Landé g-factor from first principles and take it into account in our calculations. These allow us to predict intrinsic spin lifetimes as an upper bound for experiments, identify the dominant spin relaxation pathways, and evaluate the dependence on temperature, external fields, carrier density, and impurities. We find that the Fröhlich interaction that dominates carrier relaxation contributes negligibly to spin relaxation, consistent with the spin-conserving nature of this interaction. Our theoretical approach may lead to new strategies to optimize spin and carrier transport properties.


INTRODUCTION
The field of semiconductor spintronics aims to achieve the next generation of low-power electronics by making use of the spin degree of freedom.Several classes of materials for spintronic applications have been discovered, investigated and engineered in the past decade [1][2][3][4][5].Efficient spin generation and manipulation require a large spin-orbit coupling (SOC), with GaAs a prototypical system, whereas long spin lifetimes (τ s ) is mostly found in weak SOC materials, such as graphene and diamond.Materials with large SOC as well as long τ s are ideal for spintronic applications but rare, presenting a unique opportunity for the discovery of new materials.
Halide perovskites, known as prominent photovoltaic [6] and light-emitting materials [7] with remarkable optoelectronic properties, have recently attracted interests also for spin-optoelectronic properties [8][9][10][11][12][13], since these materials exhibit both long lifetimes and large SOC (due to heavy elements).Compared to conventional spintronic materials, the optical accessibility for spin generation and detection of halide perovskites opens up a new avenue for spinoptoelectronics applications.Additionally, with highly tunable symmetry through the organic-inorganic frame-work, large Rashba splitting and high spin polarization have been realized at room temperature, critical for device applications.For example, extremely high spin polarization was produced through charge current in chiral nonmagnetic halide perovskites at room temperature in the absence of external magnetic fields [8] (B ext ), which is a hallmark in semiconductor spintronics.Persistent spin helix states that preserve SU(2) symmetry and that can potentially provide exceptionally long τ s were recently discovered in two-dimensional halide perovskites [12].
Several recent experimental studies have sought to identify the dominant spin relaxation and dephasing mechanisms to further control and elongate τ s of halide perovskites, [8][9][10][11] e.g.via time-resolved Kerr/Faraday rotations.In particular, the bulk halide perovskite such as CsPbBr 3 , which possesses one of the simplest halide perovskite structures, is a good benchmark system to understand the fundamental physical mechanisms but already presents several outstanding questions.First, what is the intrinsic τ s of CsPbBr 3 ?Experimentally this is not possible to isolate due to the unavoidable contributions from defects and nuclear spins.However, the intrinsic τ s are essential as the upper limits to guide the experimental optimization of materials.Next, what scattering processes and phonon modes dominate spin relaxation when varying the temperature, carrier density, etc.? This has been extensively studied for carrier relaxation dynamics, but not yet for spin relaxation dynamics.As we show here, the role of electron-phonon (e-ph) coupling, and especially the Fröhlich interaction known to be impor-tant for carrier relaxation in halide perovskites [14], can be dramatically different in spin relaxation.Lastly, how do electron and hole τ s respond to B ext , and what are the roles of their respective g-factor inhomogeneity?[9,10] To answer these questions, we need theoretical studies of spin relaxation and dephasing due to various scattering processes and SOC, free of experimental or empirical parameters.Previous theoretical work on spin properties of halide perovskites have largely focused on band structure and spin texture [12,15,16], and have not yet addressed spin relaxation and dephasing directly.Here, we apply our recently-developed first-principles real-time densitymatrix dynamics (FPDM) approach [17][18][19][20][21], to simulate spin relaxation and dephasing times of free electrons and holes in bulk CsPbBr 3 .FPDM approach was applied to disparate materials including silicon, (bcc) iron, transition metal dichalcogenides (TMDs), graphene-hBN, GaAs, in good agreement with experiments [17,18,20].We account for ab initio Landé g-factor and magnetic momenta, self-consistent SOC, and quantum descriptions of e-ph, electron-impurity (e-i) and electron-electron (e-e) scatterings.We can therefore reliably predict τ s with and without impurities, as a function of temperature, carrier density, and B ext .

Theory
We simulate spin and carrier dynamics based on the FPDM approach [17,18].We solve the quantum master equation of density matrix ρ (t) in the Schrödinger picture as the following: [18] dρ 12 where the first and second terms on the right side of Eq. 1 relate to Larmor precession and scattering processes respectively.The scattering processes induce spin relaxation via the SOC.H (B ext ) is the electronic Hamiltonian at an external magnetic field B ext .[H, ρ] = Hρ − ρH.H.C. is Hermitian conjugate.The subindex, e.g., "1" is the combined index of k-point and band.P is the generalized scattering-rate matrix considering e-ph, e-i and e-e scattering processes, computed from the corresponding scattering matrix elements and energies of electrons and phonons.
Starting from an initial density matrix ρ (t 0 ) prepared with a net spin, we evolve ρ (t) through Eq. 1 for a long enough time, typically from hundreds of ps to a few µs.
Historically, two types of τ s -spin relaxation time (or longitudinal time) T 1 and ensemble spin dephasing time (or transverse time) T * 2 were used to characterize the decay of spin ensemble or δS tot (t).[22,23] Suppose the spins are initially polarized along Note that without considering nuclear spins, magnetic impurities, and quantum interference effects [24], theoretical τ s (B ext = 0) should be regarded as T 1 .See more discussions about spin relaxation/dephasing in Supporting Information Sec.SI.
Below we first show theoretical results of T 1 and T * 2 of bulk (itinerant or delocalized) carriers.For bulk carriers of halide perovskites, T 1 are mainly limited by Elliott-Yafet (EY) and D'yakonov-Perel' (DP) mechanisms [25,26].EY represents the spin relaxation pathway due to mostly spin-flip scattering (activated by SOC).DP is caused by randomized spin precession between adjacent scattering events and is activated by the fluctuation of the SOC fields induced by inversion symmetry broken (ISB).Different from T 1 , T * 2 is additionally affected by the Landé-g-factor fluctuation at transverse B ext .We later generalize our results for other halide perovskites by considering the ISB and composition effects.We at the end discuss T * 2 of localized carriers due to interacting with nuclear spins.By simulating T 1 and T * 2 , and determining the dominant spin relaxation/dephasing mechanism, we provide answers of critical questions raised earlier and pave the way for optimizing and controlling spin relaxation and dephasing in halide perovskites.

Spin lifetimes at zero magnetic field
Intrinsic τ s , free from crystal imperfections and nuclear spin fluctuation, is investigated first, which sets up the ideal limit for experiments.At B ext = 0, bulk CsPbBr 3 possesses both time-reversal (nonmagnetic) and spatial inversion symmetries, resulting in Kramers degeneracy of a pair of bands between (pseudo-) up and down spins.Spin relaxation in such systems is conventionally characterized by EY mechanism [26].To confirm if such mechanism dominates in CsPbBr 3 , the proportionality between τ s and carrier lifetime (τ p , τ s ∝ τ p ) is a characteristic signature, as is discussed below.Even for intrinsic τ s , varying temperature (T ) and carrier density (n c ) would lead to large change, and its trend is informative for mechanistic understanding.
Fig. 1a and 1b show theoretical τ s at B ext = 0, including e-ph and e-e scatterings, as a function of T and n c respectively, for free electrons and holes (SI Fig. S7).Note that although bulk CsPbBr 3 crystal symmetry is orthorhombic, the spin lifetime anisotropy along three principle directions is weak (see SI Fig. S8).Therefore only τ s along the [001] direction is presented here.We have several major observations as summarized below.) FIG. 1. Spin lifetime τs of electrons of CsPbBr3.We compare electron and hole τs in Supplementary information (SI) Fig. S7 and they have the same order of magnitude at all conditions we investigated.(a) τs due to both the electron-phonon (e-ph) and electron-electron (e-e) scatterings calculated as a function of T at different electron densities ne compared with experimental data.In Fig. S6, we show τs versus T using log-scale for both y-and x-axes to highlight low-T region.Exp.A are our experimental data of T * 2 of free electrons in bulk CsPbBr3 at a small external transverse magnetic field.For Exp. A, the density of photo-excited carriers is estimated to be about 10 18 cm −3 .Exp. B are experimental data of exciton τs of CsPbBr3 films from Ref. 11. Exp.C and Exp.D are experimental data of spin relaxation time T1 of bulk CsPbBr3 and CsPbBr3 nanocrystals measured by the spin inertia method from Ref. 9 and 27 respectively.In Ref. 27, it was declared that quantum confinement effects do not modify the spin relaxation/dephasing significantly (see its Table 1), so that their T1 data can be compared with our theoretical results.For Exp. C and D, the measured lifetimes cannot be unambiguously ascribed to electrons or holes and can be considered as values between electron and hole T1.The carrier densities are not reported for Exp.First, a clear decay of τ s as increasing T is observed.As τ s with and without e-e scattering (SI Fig. S7) has little difference, this indicates e-ph scattering is the dominant spin relaxation mechanism (without impurities and B ext ).Note that with increasing T , phonon occupations increase, which enhances the e-ph scattering and thus lowers both carrier (τ p ) and spin (τ s ) lifetime.
Next, τ s steeply decreases with increasing n c at low T but is less sensitive to n c at high T , as shown in Fig. 1b.The trend of T 1 decreasing with n c is consistent with the experimental observation of T 1 decreasing with pump power/fluence in halide perovskites [28][29][30][31].At 4 K, τ s decreases steeply by three orders of magnitude with n c increasing from 10 16 cm −3 to 10 19 cm −3 .Such phenomenon was reported previously for monolayer WSe 2 [18,32], where spin relaxation is dominated by EY mechanism, the same as in CsPbBr 3 .The cause of such strong n c -dependence at low T is discussed below in more details, attributing to n c effects on (averaged) spin-flip matrix elements.As a result, at low T and low n c , τ s of CsPbBr 3 can be rather long, e.g., ∼200 ns at 10 K and ∼8 µs at 4 K.This is in fact comparable to the ultralong hole τ s of TMDs and their heterostructures [18,33,34], ≥2 µs at ∼5 K, again suggesting the advantageous character of halide perovskite in spintronic applications.
Importantly, good agreement between theoretical results and several independent experimental measure-ments is observed.Our theoretical results agree well with experimental T 1 of bulk CsPbBr 3 [9] (Exp.C) assuming n c ≈ 10 18 cm −3 , and CsPbBr 3 nanocrystal [27] (Exp.D) assuming n c ≈ 10 16 cm −3 , respectively.We further compare theoretical results with our own measured T * 2 (at B ext =100 mT; Exp.A).Excellent agreement is observed at T ⩾10 K with n c around 10 18 cm −3 (estimated from the experimental averaged pump power).The agreement however becomes worse at T <10 K.The discrepancy is possibly due to nuclear-spin-induced spin dephasing of carriers, as will be discussed in the last subsection.
We then study the effects of the e-i scattering on τ s for various point defects.We find that at low T , e.g., T <20 K, the e-i scattering reduces τ s , consistent with EY mechanism (which states increasing extrinsic scatterings reduces spin lifetime).With a high impurity density n i , e.g., 10 18 cm −3 , the e-i scattering may significantly reduce τ s below 10 K, seemingly leading to better agreement between theoretical τ s and experimental data from Exp. A, as shown in SI Fig. S9.However, as will be discussed below in the subsection of magnetic-field effects, a relatively high n i predicts incorrect values of T * 2 and worse agreement with experimental data (Exp.A) on B ext -dependence.Therefore, the discrepancy between our theoretical τ s and our measured T * 2 below 10 K is probably not explained by the impurity scattering effects.
In addition, the electron and hole τ s have the same order of magnitude (Fig. S7), consistent with experiments, but in sharp contrast to conventional semiconductors (e.g., silicon and GaAs [35]), which have longer electron τ s than hole owing to band structure difference between valence and conduction band edges.
Finally, we also predict the spin diffusion length (l s ) of pristine CsPbBr 3 in the low-density limit, which sets the upper bound of l s at different T .We use the relation l s = √ Dτ s , where D is diffusion coefficient obtained using the Einstein relation, with carrier mobility µ from first-principles calculations [19] (more details in Sec.SVII).Excellent agreement between theoretical and experimental carrier mobility is found for CsPbBr 3 (SI Fig. S12a).We find l s is longer than 10 nm at 300 K, and possibly reach tens of µm at T ≤ 10 K (see details in Sec.SVII and Fig. S12 in SI).

Analysis of spin-phonon relaxation
To gain deep mechanistic insights, we next analyze different phonon modes and carrier density effects on spin relaxation through examining spin-resolved e-ph matrix elements.
In Fig. 2, we compare the contribution of different phonon modes to τ s and τ p .First, we find that at a very low T -4 K, only acoustic modes (A1-A3) contribute to spin and carrier relaxation.This is simply because the optical phonons are not excited at such low T (corresponding k B T ∼0.34 meV much lower than optical energy ≳ 2 meV).At T ⩾10 K, optical modes are more important for both spin and carrier relaxation (green and blue dashed lines closer to black line (all phonons) in Fig. 2).
In particular, from Fig. 2b, we find that two special optical modes -57th and 58th modes (O57-O58, modes ordered by phonon energy with their phonon vector plots in SI Fig. S3) dominate carrier relaxation at T ⩾50 K, because τ p due to O57-O58 (blue dashed line) nearly overlaps with τ p due to all phonon modes (black line).These two optical modes are mixture of longitudinal and transverse vibration as shown in SI Fig. S3.In contrast, for spin relaxation in Fig. 2a, at T ⩾10 K O57-O58 are less important than other optical modes (green dashed line).More specifically, in this temperature range, there are tens of phonon modes (with energies ranging from 2 meV to 18 meV), contributing similarly to spin relaxation.This is contradictory to the simple assumption frequently employed in previous experimental studies [9,36,37] that a single longitudinal optical (LO) phonon with a relatively high energy (e.g.∼18 meV for CsPbBr 3 in Ref. 9) dominates spin relaxation over a wide T range, e.g., from 50 K to 300 K, through a Fröhlich type e-ph interaction.
In the simplified picture of Fermi's golden rule (FGR), τ −1  Here carrier density nc is set to be 10 18 cm −3 .We note that special optical phonon modes O57 and O58 are dominant in carrier relaxation above 50 K (panel b), consistent with the usual Fröhlich interaction picture, but are not important in spin relaxation (panel a).
From Fig. 3a, we find that spin-flip ME is dominated by "other optical modes" (blue line), opposite to the spinconserving ME in Fig. 3b (i.e.instead, dominated by special optical phonon modes O57 and O58 (red line)).This well explains the different roles of optical O57-O58 modes in carrier and spin relaxation.Moreover, spin-conserving ME for O57-O58 in Fig. 3b diverges at q → 0, which indicates its dominant long-range nature, consistent with the common long-range Fröhlich interaction picture [38], mostly driving carrier relaxation in polar materials at high T (e.g., 300 K).On the contrary, the small magnitude of spin-flip ME for O57-O58 modes indicates that Fröhlich interaction is unimportant for spin relaxation.This is because all spin-dependent parts of the e-ph interaction are short-ranged, while Fröhlich interaction is the only long-range part of the e-ph interaction but is spin-independent.This important conclusion again emphasizes the sharp difference between spin and carrier relaxations in polar materials.
To explain the strong n c dependence of τ s at low T , we

. The analysis of the e-ph matrix elements (ME). (a)
The q-resolved modulus square of spin-flip e-ph ME | g ↑↓ | 2 (q) at a high temperature -300 K with a part of or all phonon modes.(b) The same as panel (a) but for spin-conserving e-ph ME g Ω = ± g Ω if the excess/excited spin δS tot (t) precesses along ±δS tot (t) × B ext .g Ω is close to the averaged g-factor g defined in Eq. 10. (c) The effective amplitude of the fluctuation of g factors -∆ g defined in Eq. 11 as a function of carrier density at 10 K.
edge.In "Methods" section, we have proven that at the low density limit, since carrier occupation satisfies Boltzmann distribution, both | g ↑↓ | 2 and D S are µ F,c and n c independent.
Landé g-factor and transverse-magnetic-field effects At B ext , the electronic Hamiltonian reads where µ B is Bohr magneton; g 0 is the free-electron gfactor; S and L are the spin and orbital angular momentum respectively.The simulation of L is nontrivial for periodic systems and the details are given in Method section and Ref. ? .Having H (B ext ) at a transverse B ext perpendicular to spin direction, T * 2 is obtained by solving the density-matrix master equation in Eq. 1.
The key parameters for the description of the magnetic-field effects are the Landé g-factors.Their values relate to B ext -induced energy splitting (Zeeman effect) ∆E k (B ext ) and Larmor precession frequency Ω k , satisfying Ω k ≈ ∆E k = µ B B ext g k with g k the k-resolved Landé g-factor.More importantly, the g-factor fluctuation (near Fermi surface or µ F,c ) leads to spin dephasing at transverse B ext , corresponding to T * 2 .Fig. 4a shows g k of electrons and holes at k-points around the band edges.
g k are computed using ∆E k (B ext ) (Eq. 8 and 9) obtained from H k (B ext ).Our calculated electron g k are larger than hole g k , and the sum of electron and hole g k range from 1.85 to 2.4, in agreement with experiments [9,35].Furthermore, g k are found sensitive to state energies and wavevectors k, and the fluctuation of g k is enhanced with increasing the state energy.In Figs.4b and 4c, we show the global g-factor g Ω and the amplitude of the g-factor fluctuation (near the Fermi surface) ∆ g (Eq.11) as a function of n c .Both g Ω and ∆ g are insensitive to n c at n c <10 16 cm −3 , but sensitive to n c at n c ⩾10 16 cm −3 .
In Fig. 4, we show ab initio g-factors computed with the PBE functional [39].We further compare g-factors computed using different exchange-correlation functionals (V xc ) in SI Sec.SV.It is found that the magnitude of ∆ g and the trend of g-factor change with the state energy are both insensitive to V xc .Since T * 2 only depends on ∆ g, our predictions of T * 2 should be reliable.Next, we discuss magnetic-field effects on τ s in Fig. 5, calculated from our FPDM approach, and analyze them with phenomenological models.At transverse B ext , the total spin decay rate is approximately expressed by (i) Free induction decay (FID) mechanism if τ p ∆Ω ≳ 1 (weak scattering limit).We have where C ∆g is a constant and often taken as 1 or 1/ √ 2 ≈ 0.71 [9,22,36,[40][41][42][43]. The latter assumes a Gaussian distribution of g-factors and the scattering being completely absent [22,41,43].
As a result, we only discuss τ s at B ext ̸ = 0 below 20 K, specifically at 4K afterwards.From Fig. 5b, we can see that magnetic-field effects on electron and hole τ s are quite similar, which is a result of their similar band curvatures, e-ph scattering, and ∆ g (Fig. 4c), although their absolute g-factors are quite different, as shown in Fig. 4a and 4b.
We further examine magnetic-field effects on τ s at 4 K in Fig. 5c.As discussed above, τ −1 s (B ext ) increases with B ext .More specifically, we find that the calculated τ −1 s (B ext ) is proportional to (B ext ) 2 at low B ext (details in SI Fig. S13) following the DP mechanism (Eq.5), but linear to B at higher B following the free induction decay mechanism(Eq.4).
The comparison of calculated τ −1 s (B ext ) with experimental data (orange diamond in Fig. 5b) is reasonable with n e around 10 18 cm −3 (the experimental estimated average n c ).However, their B ext -dependence is not the same in the small B ext range, e.g. at (as shown in SI Fig. S13), whereas the experimental τ −1 s (B ext ) is more likely linear to B ext .In principles, extremely small B ext will lead to ∆Ω small enough falling in the DP regime ((τ ∆Ω s ) −1 proportional to (B ext ) 2 ).However, experimental results still keep in the FID regime ((τ ∆Ω s ) −1 linear dependent on B ext ) at small B ext .This inconsistency implies additional magnetic field fluctuation contributes to ∆Ω and/or other faster spin dephasing processes exist at small external B ext .It may originate from nuclear spin fluctuation, magnetic impurities, carrier localization, chemical potential fluctuation, etc. [9,35] in samples, which are however rather complicated for a fully first-principles description.In this work, we focus on spin dephasing of bulk carriers due to Zeeman effects and various scattering processes.
Moreover, in Fig. 5c, we find that at B ext ̸ = 0, τ s decreases with n c , similar to the case at B ext = 0.But the origin of the strong n c dependence at high B ext is very different from τ s at B ext = 0.When B ext ≥0.4 Tesla, τ s is dominated by the FID mechanism (Eq.4), thus its n c dependence is mostly from ∆ g's strong n c dependence shown in Fig. 4c.
Finally, we show τ −1 s (B ext ) as a function of B ext at 4 K with the e-i scattering in Fig. S14.We find that with relatively strong impurity scattering (e.g, with 10 17 cm −3 V Pb neutral impurities), the B ext -dependence of τ s becomes quite weak, in disagreement with experiments, indicating that impurity scattering is probably weaker in those experiments.See more discussions in Sec.SVIII.

Inversion symmetry broken (ISB), composition effects and hyperfine coupling
For halide perovskites, ISB may present due to ferroelectric polarization, strain, surface, applying electric fields, etc.One of the most important effects from ISB is inducing k-dependent SOC fields (called B in ).B in can change the electronic energies and spin textures, which may significantly modify the spin relaxation/dephasing.The linear relation between ensemble spin dephasing rate and B ext was frequently found and used in previous experimental measurements [9,36,40,41].
To understand the ISB effects, we simulate τ s with two important types of B in -Rashba and PSH (persistent spin helix) ones.Rashba SOC presents in many 2D and 3D materials, e.g., wurtzite GaN and graphene on SiO 2 /hBN.PSH exhibits SU(2) symmetry [44,45] which is robust against spin-conserving scattering, and was recently realized in 2D hybrid perovskites [12].Their effects are considered by introducing an additional SOC term to the electronic Hamiltonian perturbatively.The specific forms of Rashba and PSH SOC Hamiltonians are given in Eq. 18-22 in "Methods" section.
From Fig. 6a, we find that τ s is reduced by Rashba SOC and the reduction is significant when the SOC coefficient α ≥ 0.5 eV Å.This is because Rashba SOC induces a nonzero ∆Ω ∝ α and then induces an DP/FID spin decay channel additional to the EY one.Similar to Eq. 3, the total rate τ −1 ∆Ω becomes large, so that τ s is significantly reduced from τ −1 s | α=0 .τ s keeps decreasing with α but its low limit is bound by τ p .On the other hand, with PSH SOC, τ s (along PSH B in -B PSH , which is along y direction here) is unchanged at α ≤ 2 eV Å, and increases at a larger α.The reason is: with PSH SOC, spins are highly polarized along B PSH , so that τ s along B PSH is still dominated by EY mechanism (no spin precession).One critical effect of B PSH is then modifying the energies (spin split energies).At small α, the energy changes are not significant compared with k B T , so that the e-ph scattering contribution to spin relaxation is not modified much; as a result, τ s is close to τ s | α=0 .From Fig. 6b, we can see that at large α (e.g., 7 eV Å) the band structure is however significantly changed.The valence band maxima are now at two opposite k-points away from Γ and with opposite spins.Therefore, at large α, spin relaxation is dominated by the spin-flip scattering processes between two opposite valleys away from Γ.This can lead to longer τ s since the spin-flip processes within one valley (intravalley scattering) are suppressed, essentially a spin-valley locking condition is realized [12,19].Our FPDM simulations with model SOC suggest that Rashba SOC likely reduces τ s while PSH SOC can enhance τ s as anticipated in previous experimental study [45].Note that in practical materials, the ISB effects may not be completely captured by model SOC fields as introduced here.Although in general, we include self-consistent SOC in our FPDM calculations instead of perturbatively, but since the studied equilibrium bulk structure has inversion symmetry, we therefore have to include model ISB SOC perturbatively to simulate such effects induced by various causes.Therefore, further FPDM simulations of materials with ISB structures are important for comprehensive understanding of the ISB effects.
Furthermore, it is crucial to understand the chemical composition effects to improve our understandings of spin dynamics and transport in many other kinds of halide perovskites beside CsPbBr 3 .As an initial study, we performed FPDM simulations of τ s of holes of pristine bulk CsPbCl 3 , CsPbI 3 , MAPbBr 3 and CsSnBr 3 as a function of temperature, at the same carrier density.We consider the inversion-symmetric orthorhombic phase for all systems, the same as CsPbBr 3 here, in order to study chemical composition effect alone.From Fig. 7, our FPDM simulations show that the differences of τ s of CsPbBr 3 , CsPbCl 3 , CsPbI 3 , MAPbBr 3 and CsSnBr 3 are mostly tens of percent or a few times in the wide temperature range from 4 K to 300 K. Specifically, τ s of MAPbBr 3 is found always shorter than CsPbBr 3 .τ s of CsSnBr 3 is found slightly longer than CsPbBr 3 at 300 K but becomes shorter than CsPbBr 3 at T <100 K.A trend of hole τ s is found for CsPbX 3 : τ s (CsPbCl 3 ) > τ s (CsPbBr 3 ) > τ s (CsPbI 3 ), indicating that the lighter the halogen atom, the longer the spin lifetime.This trend may be partly due to two reasons: (i) For the band gap, we have CsPbCl 3 > CsPbBr 3 > CsPbI 3 (1.40,1.03 and 0.75 eV respectively at PBE), so that spin mixing due to the conduction-valence band mixing is reduced at lighter halogen compound, which usually weakens the spin-phonon interaction; (ii) The lighter halogen atom reduces the SOC strength of the material (weaker SOC reduces the spin mixing between up and down states).Additionally, we find that for all these inversion-symmetric orthorhombic materials, the anisotropy of τ s along different crystalline directions is rather weak (see SI Fig. S8).Overall, our results indicate that the chemical composition effects on τ s are not very strong when comparing with the effects of the symmetry change (e.g.broken inversion symmetry resulting in Rashba or PSH discussed in Fig. 6).A more systematic study of the composition, symmetry, and dimensionality effects is of great importance and will be our future work.
Above we focus on spin relaxation/dephasing of bulk (or itinerant or delocalized) electrons, for which hyperfine coupling is usually unimportant [22,46].In actual samples, due to polarons, ionized impurities, etc., a considerable portion of electron carriers are localized.It is known that hyperfine coupling can induce spin dephasing of localized electrons through spin precessions about randomly-distributed nuclear-spin (magnetic) fields B Nuclear .When nuclear spins are weakly polarized (because of weak B ext ), T * 2 of localized electrons -T * 2,loc is often estimated based on FID mechanism 1/T * 2,loc ∼ σ Ω N [47][48][49], where Ω N is Larmor frequency of a localized electron due to B Nuclear and σ Ω N is the parameter describing its fluctuation or determining its distribution (Eq.28 for B ext =0).According to Refs. 9, 47-49, σ 2 Ω N ∼ C loc /V loc (Eq.30), where V loc is the localization volume.At B ext =0, C loc is determined by hyperfine constant A, nuclear spin I, isotope abundance and unitcell volume (Eq.31).See detailed formulae and our estimates of the above quantities in "Methods" section.Our estimated C loc is ∼180 and ∼530 nm 3 ns −2 for electrons and holes respectively.The estimated localization radii halide perovskites are 2.5-14 nm [50][51][52][53], which lead to T * 2,loc (B ext = 0) ∼0.6-8.0 ns for electrons and ∼0.35-4.6 ns for holes.Since bulk and localized carriers coexist in materials, T * 2,loc roughly gives the lower bound of the effective carrier T * 2 .In addition to the hyperfine coupling for spin dephasing of localized carriers above, the fluctuation of hyperfine coupling for bulk (delocalized) carriers at different kpoints may lead to spin dephasing when nuclear spins are polarized along a non-zero transverse B ext .This effect is however rather complicated (probably requiring the difficult L contribution [54] to hyperfine coupling), beyond the scope of this work.
In summary, through a combined ab initio theory and experimental study, we reveal the spin relaxation and dephasing mechanism of carriers in halide perovskites.Using our FPDM approach and implementing ab initio magnetic momenta and g-factor, we simulate free-carrier τ s as a function of T and B ext , in excellent agreement with experiments.The transverse magnetic-field effects are found only significant at T <20 K.We predict ultralong T 1 of pristine CsPbBr 3 at low T , e.g., ∼200 ns at 10 K and ∼8 µs at 4 K.We find strong n c dependence of both T 1 and T 2 at low T , e.g.τ s can be tuned by three order of magnitude with n c from the low density limit to 10 19 cm −3 .The reasons are attributed to the strong electronic-energy-dependences of spin-flip e-ph matrix elements and ∆ g for T 1 and T * 2 respectively.From the analysis of e-ph matrix elements, we find that contrary to common belief, Fröhlich interaction is unimportant to spin relaxation, although critical for carrier relaxation.We further study ISB and composition effects on τ s of halide perovskites.We find that ISB effects can significantly change τ s , i.e. spin lifetime can increase with PSH SOC, but not with Rashba SOC.The composition effect is found not very strong and only changes τ s by tens of percent or a few times in a wide temperature.Our work provides fundamental insights on how to control and manipulate spin relaxation in halide perovskites, which are vital for their spintronics and quantum information applications.

METHODS
Spin dynamics and transport.Spin dynamics and spin lifetime τ s are simulated by our recently developed first-principles density-matrix dynamics (FPDM) method [17][18][19][20][21]. Starting from an initial state with a spin imbalance, we evolve the time-dependent density matrix ρ (t) through the quantum master equation with Lindblad dynamics for a long enough simulation time, typically from ns to µs, varying with systems.After obtaining the excess spin observable δS tot (t) from ρ (t) and fitting δS tot (t) to an exponentially oscillating decay curve, the decay constant τ s and the precession frequency Ω are then obtained (Eq.S3 and Fig. S1 in SI).All required quantities of FPDM simulations, including electron energies, phonon eigensystems, e-ph and e-i scattering matrix elements, are calculated on coarse k and q meshes using the DFT open source software JDFTx [55], and then interpolated to fine meshes in the basis of maximally localized Wannier functions [56][57][58].The e-e scattering matrix is computed using the same method given in Ref. 18.More theoretical background and technical details are given in Ref. 19 and 18, as well as the Supporting Information.
Using the same first-principles electron and phonon energies and matrix elements on fine meshes, we calculate the carrier mobility by solving the linearized Boltzmann equation using a full-band relaxation-time approximation [59] and further estimate spin diffusion length based on the drift-diffusion model (SI Sec.SVII).
Orbital angular momentum.With the Blöch basis, the orbital angular momentum reads where m and n are band indices; ϵ and u are electronic energy and the periodic part of the wavefunction, respectively; H 0 is the zero-field Hamiltonian operator.Eq. 7 can be proven equivalent to L = 0.5 * (r × p − p × r) with r and p the position and momentum operator respectively.The detailed implementation of Eq. 7 is described in Ref. ? .Our implementation of L has been benchmarked against previous theoretical and experimental data for monolayer MoS 2 (Table S1).g-factor of free carriers.In experimental and model Hamiltonian theory studies [9,35], g-factor is defined from the ratio between either B ext -induced energy splitting ∆E (B ext ) or Larmor precession frequency Ω (B ext ) to µ B B. Therefore, we define g-factor of an electron or a hole at state k, where g S k is g-factor defined based on spin expectation values.B ext is the unit vector along B ext , where S exp k,h B ext and S exp k,l B ext are the spin expectation value (exp) of the higher (h) and lower (l) energy band at k projected to the direction of B ext respectively.However, in previous theoretical studies of perovskites [35,60], g-factors were defined based on pseudo-spins related to the total magnetic momenta J at , which are determined from the atomic-orbital models.The pseudo-spins can have opposite directions to the actual spins.Most previous experimental studies adopted the same convention for the signs of carrier g-factors.Therefore, to compare with g-factors obtained in previous theoretical and experimental studies, we introduce a correction factor C S→J and define a new g-factor: C S→J = m at S /m at J with m at J and m at S the total and spin magnetic momenta respectively, obtained from the atomic-orbital model [35].C S→J is independent from kpoint, and is ∓1 for electrons and holes respectively for CsPbBr 3 .
g k is different at different k; therefore we define its statistically averaged value (depending on temperature T and chemical potential µ F,c ) as and its fluctuation amplitude as where f ′ k is the derivative of the Fermi-Dirac distribution function.Here for simplicity the band index of f ′ k is dropped considering both valence and conduction bands are two-fold degenerate.
We have further defined a more general g-factor as a tensor and its fluctuation amplitude in SI Sec.SV.For CsPbBr 3 , we find different definitions predict quite similar values (differences are not greater than 10%).
Analysis of e-ph matrix elements.For EY spin relaxation, in the simplified picture of Fermi's golden rule (FGR), τ −1 s is proportional to the modulus square of the spin-flip scattering matrix element.As the e-ph scattering plays a crucial role in spin relaxation in CsPbBr 3 , it is helpful to analyze the spin-flip e-ph matrix elements.
Note that most matrix elements are irrelevant to spin relaxation and we need to pick the "more relevant" ones, by defining a weight function related to occupation and energy conservation.Therefore we propose a T and µ F,c dependent effective modulus square of the spin-flip e-ph matrix element | g ↑↓ | 2 as where g ↑↓,qλ k,k−q is e-ph matrix element, related to a scattering event between two electronic states of opposite spins at k and k − q through phonon mode λ at wavevector q; n qλ is phonon occupation; f k is Fermi-Dirac function; ω c is the characteristic phonon energy specified below, and w k,k−q is the weight function.Here we drop band indices for simplicity, as CsPbBr 3 bands are twofold Kramers degenerate and only two bands are relevant to electron and hole spin/carrier dynamics.
The matrix element modulus square is weighted by n qλ since τ −1 s is approximately proportional to n qλ according to Eq. 5 of Ref. 17.This rules out high-frequency phonons at low T which are not excited.ω c is chosen as 4 meV at 10 K based on our analysis of phonon-moderesolved contribution to spin relaxation.The trends of | g ↑↓ | 2 are found not sensitive to ω c as checked.w k,k−q selects transitions between states separated by ω c and around the band edge or µ F,c , which are "more relevant" transitions to spin relaxation.
We also define a q-resolved modulus square of the spinflip e-ph matrix element | g ↑↓ | 2 (q) as Note that for spin relaxation, only states around the band edges are relevant.Thus we restrict |ϵ k −ϵ edge | <180 meV for the calculation of Eq. 14, which is about 7k B T at 300 K relative to the band edge energy (ϵ edge ).
Analysis of the EY spin relaxation rate.According to FGR, the EY spin relaxation rate of an electronic state should be also proportional to the density of pair states allowing spin-flip scattering between them.Therefore, we propose a scattering density of states D S (which is T and µ F,c dependent), D S can be regarded as an effective density of spinflip or spin-conserving e-ph transitions satisfying energy conservation between one state and its pairs (considering that the number of spin-flip and spin-conserving transitions are the same for Kramers degenerate bands).
When ω c = 0 (i.e.elastic scattering), we have We then discuss µ F,c dependence of τ −1 s at low n c limit.For simplicity, we only consider conduction electrons.At low n c limit, we have exp Therefore, according to Eq. 12, 13 and 15, both | g ↑↓ | 2 and D S are independent from µ F,c (as exp and n c at low n c region, e.g.much lower than 10 16 cm −3 for CsPbBr 3 .We can similarly define spin conserving matrix elements | g ↑↑ | 2 and | g ↑↑ | 2 (q) by replacing g ↑↓,qλ k,k−q to g ↑↑,qλ k,k−q in Eq. 12 and 14.Then we have the approximate relation for carrier relaxation rate due to e-ph scattering, The Hamiltonian for model SOC.In general, the Hamiltonian for model SOC reads For the Rashba field, − → Ω model k in Eq. 18 is defined in the plane (xy plane here) surrounding Γ point, where α R is the Rashba SOC strength coefficient.f cut (k/k cut ) is 1 at small k but vanishes at large k.It is introduced to truncate the SOC fields at k > k cut smoothly in order to avoid unphysical band structures around first Brillouin zone boundaries.It is taken as k cut is taken 0.12 bohr −1 for CsPbBr 3 .This value is about half of the length of the shortest reciprocal lattice vector, about 0.28 bohr −1 for orthorhombic CsPbBr 3 .We can see that f cut is almost 1 at k = Γ but almost vanishes at first Brillouin zone boundaries.
Persistent Spin Helix (PSH) was first proposed by Bernevig et al. [44].PSH has SU(2) symmetry which is robust against spin-conserving scattering.In general, for PSH SOC, where directions i and j are orthogonal.PSH fields are all along the same axis (y axis here).We take where α PSH is the PSH SOC strength coefficient.T * 2,loc due to nuclear spin fluctuation.The Hamiltonian of hyperfine coupling between an electron spin and nuclear spins approximately reads [9,48] where − → Ω N is Larmor precession vector, related to the effective hyperfine field (called Overhauser field) generated by all nuclei and acting on electron spin.s is the spin operator of the electron.Eq. 24 specifically refers to the hyperfine Fermi contact interaction between an electron and nuclear spins.The sum in Eq. 24 goes over all nuclei.I j is the spin operator of nucleus j.V u is the unit cell volume.A j is the hyperfine coupling constant considering only the Fermi contact contribution, which was assumed to be the dominant contribution in Refs.9, 48, 49 for CsPbBr 3 and GaAs.µ j and I j are the magnetic moment and spin of nucleus j, respectively.µ B is the Bohr magneton.ψ (R j ) and u c (R j ) are the electron envelope wave function and the electron Bloch function at the j-th nucleus respectively, whose product gives the electronic wavefunction With this definition, |u c (R j )| 2 ∝ 1/V u , therefore, from Eq. 25, The value of A j depends on the isotope of the nucleus.For CsPbBr 3 , it was found that the relevant isotopes are 207 Pb with natural abundance of about 22% for holes, and 79 Br and 81 Br for electrons [9].Since the total abundance of 79 Br and 81 Br is almost 100% and their nuclear spins are both 3/2, 79 Br and 81 Br can be treated together.
According to the proportional relation in Eq. 27, A j of orthorhombic CsPbBr 3 is approximately 1/4 of A j of cubic CsPbBr 3 , considering that their Bloch functions at the band edges are similar [61] (e.g., their hole Bloch functions both have significant Pb-s-orbital contribution), and V u of orthorhombic CsPbBr 3 is about 4 times of that of cubic CsPbBr 3 .Therefore, using estimated A j of cubic CsPbBr 3 in Ref. 9, we obtain that A j of 207 Pb for holes is about 25 µeV and A j of Br for electrons is about 1.75 µeV.
When nuclear spins are not polarized (due to B ext =0), the nuclear field is zero on average.However, due to the finite number of nuclei interacting with the localized electron, there are stochastic nuclear spin fluctuations, which are characterized by the probability distribution function [47] where σ Ω N determines the dispersion of hyperfine field, and the angular brackets denotes the statistical averaging: Ω 2 N = 3σ 2 Ω N /2.For the independent and randomly oriented nuclear spins, we have (at B ext =0) where j u is nucleus index in the unit cell, ξ is the isotope, and c is the unit cell index in the whole system.α ξ is the abundance of isotope ξ.Since ψ (R juc ) usually varies slowly on the scale of a unit cell, V u c |ψ (R juc )| 4 can be replaced by an integral in the whole system -|ψ (r)| 4 dr.Define V loc = 1/ |ψ (r)| 4 dr.V loc is the localization volume.Therefore (at B ext =0), C loc = 2V u 3 juξ α ξ I juξ (I juξ + 1) A 2 juξ .
With σ Ω N , T * 2,loc is often estimated based on FID mechanism [47][48][49] (Eq.4) T * 2,loc ∼ σ −1 Ω N .As α ξ , I juξ and V u can be easily obtained and with A juξ estimated above, we obtain C loc ∼180 and ∼530 nm 3 ns −2 for electrons and holes respectively.V loc can be estimated from the localization radii r loc of localized carriers, In Table S2, we listed values of the parameters used to calculate T * 2,loc .Experimental synthesis.Growth of CsPbBr 3 single crystal: Small CsPbBr 3 seeds were first prepared with fresh supersaturated precursor solution at 85 • C. Small and transparent seeds were then picked and put on the bottom of the vials for large crystal growth.The temperature of the vials was set at 80 • C initially with an increasing rate of 1 • C/ h, and was eventually maintained at 85 • C. Vials were covered with glass slides to avoid fast evaporation of the DMSO.So the growth driving force is supersaturation achieved by slow evaporation of DMSO solvent.After 120-170 hours, a centimeter-sized single crystal was picked from the solution, followed by wiping the residue solution on the surface.
Experimental Spin Lifetime Measurement.For measuring the spin lifetime in CsPbBr 3 single crystals, we have used the ultrafast circularly-polarized photoinduced reflectivity (PPR) method at liquid He temperature under the influence of a magnetic field.The experimental setup was described elsewhere [62? ].It is a derivative of the well-known 'pump-probe' technique, where the polarization of the pump beam is modulated by a photoelastic modulator between left (δ + ) and right (δ − ) circular polarization, namely LCP and RCP, respectively.Whereas the probe beam is circularly polarized (either LCP or RCP) by a quarter-wave plate.The transient change in the probe reflection, namely c-PPR(t), was recorded.The 405 nm pump beam, having 150 femtoseconds pulse duration at 80 MHZ repetition rate, was generated by frequency doubling the fundamental at 810 nm from the Ti:Sapphire laser (Spectra Physics) using a SHG BBO crystal.The 533 nm probe beam was generated by combining the 810 nm fundamental beam with the 1560 nm infrared beam from an OPA (optical parametric amplifier) onto a BBO type 2 SFG (Sum Frequency Generation) crystal.The pump/probe beams having average intensity of 12 Wcm −2 and 3 Wcm −2 , respectively were aligned onto the CsPbBr 3 crystal that was placed inside a cryostat with a built-in electro-magnet that delivered a field strength, B up to 700 mT at temperatures down to 4 K.Using this technique we measured both t-PPR responses at both zero and finite B to extract the Bdependent electron and hole spin lifetimes.From the c-PPR(B,t) dynamics measured on (001) facet with B directed along [010] [62] (see example c-PPR(B,t) dynamics in SI Fig. S15), we could obtain the electron and hole T * 2 by fitting the transient quantum beating response with two damped oscillation functions: 2,e cos(2πf 1 t + ϕ 1 ) + A 2 e −t T * 2,h cos(2πf 2 t + ϕ 2 ), (33) where T * 2,e and T * 2,h are the spin dephasing times of the electrons and holes; f1 and f2 are the two QB frequencies that can be obtained directly from the fast Fourier transform of the c-PPR dynamics.

DATA AVAILABILITY
The input files of all simulations (including the groundstate DFT simulations, Wannier fitting and interpo-lation, and the real-time density-matrix simulations), python post-processing scripts, example output files and necessary source data files (for plotting) generated for this study are available in the SI repository (LINK GITHUB).
FIG.1.Spin lifetime τs of electrons of CsPbBr3.We compare electron and hole τs in Supplementary information (SI) Fig.S7and they have the same order of magnitude at all conditions we investigated.(a) τs due to both the electron-phonon (e-ph) and electron-electron (e-e) scatterings calculated as a function of T at different electron densities ne compared with experimental data.In Fig.S6, we show τs versus T using log-scale for both y-and x-axes to highlight low-T region.Exp.A are our experimental data of T * 2 of free electrons in bulk CsPbBr3 at a small external transverse magnetic field.For Exp. A, the density of photo-excited carriers is estimated to be about 10 18 cm −3 .Exp. B are experimental data of exciton τs of CsPbBr3 films from Ref. 11. Exp.C and Exp.D are experimental data of spin relaxation time T1 of bulk CsPbBr3 and CsPbBr3 nanocrystals measured by the spin inertia method from Ref. 9 and 27 respectively.In Ref.27, it was declared that quantum confinement effects do not modify the spin relaxation/dephasing significantly (see its Table1), so that their T1 data can be compared with our theoretical results.For Exp. C and D, the measured lifetimes cannot be unambiguously ascribed to electrons or holes and can be considered as values between electron and hole T1.The carrier densities are not reported for Exp.C and D. (b) τs due to both the e-ph and e-e scatterings as a function of ne at different T .The vertical dashed line in panel (b) corresponding to ne with chemical potential µF,c at the conduction band minimum (CBM).

FIG. 2 .
FIG.2.Phonon-mode contribution analysis.(a) Spin lifetime τs and (b) carrier lifetime τp due to different phonon modes."A" and "O" denote acoustic and optical modes respectively.The number index is ordered by increasing phonon energies.The phonon dispersion is given in SI Fig.S2.Here carrier density nc is set to be 10 18 cm −3 .We note that special optical phonon modes O57 and O58 are dominant in carrier relaxation above 50 K (panel b), consistent with the usual Fröhlich interaction picture, but are not important in spin relaxation (panel a).

3 0FIG. 4 .
FIG. 3. The analysis of the e-ph matrix elements (ME).(a)The q-resolved modulus square of spin-flip e-ph ME | g ↑↓ | 2 (q) at a high temperature -300 K with a part of or all phonon modes.(b) The same as panel (a) but for spin-conserving e-ph ME| g ↑↑ | 2 (q).(c) | g ↑↓ | 2 , | g ↑↑ | 2 and | g ↑↓ | 2 D S ofconduction electrons as a function of carrier density at a low T -10 K compared with the spin relaxation rates 1/τs.| g ↑↓ | 2 and | g ↑↑ | 2 are the T and µF,c dependent effective (averaged around the band edge or µF,c) modulus square of spin-flip and spin-conserving e-ph ME, respectively (see Eq. 12).D S is the scattering density of states (Eq.15).The vertical dashed line corresponding to µF,c at CBM.

where τ 0 s− 1
is the zero-field spin relaxation rate due to EY mechanism; τ ∆Ω s −1 is induced by the Larmorprecession-frequency fluctuation (∆Ω = µ B B ext ∆ g), and can be described by different mechanisms depending on the magnitude of τ p ∆Ω[22,26]:

FIG. 5 . 1 s
FIG. 5.The effects of transverse B ext (perpendicular to spin direction) on calculated τs of free carriers of CsPbBr3.(a) The ratio of τs at B ext =1 T and τs at B ext =0 as a function of T .Here electron carrier density ne=10 18 cm −3 .(b) Spin decay rates (τ −1 s ) of electrons and holes as a function of B ext at 4 K with nc = 10 18 cm −3 .(c) τ −1 s as a function of B ext at 4 K at different ne."Exp."(orange open diamond) represent our experimental data (with B ext along [010] direction), where the density of photo-excited carriers is estimated about 10 18 cm −3 .The orange dashed line is the linear fit of experimental data.The linear relation between ensemble spin dephasing rate and B ext was frequently found and used in previous experimental measurements[9,36,40,41].

FIG. 6 .
FIG. 6.The effects of model SOC fields.(a) textures in the kx − ky plane of the CsPbBr3 system with model Rashba SOC. S exp ≡ S exp x , S exp y , S exp z

are
Larmor precession vectors induced by k-dependent B in .s k is spin operator.With the total electronic Hamiltonian H k = H 0,k + H model k , τ s considering the effects of model SOC is obtained by solving the density-matrix master equation in Eq. 1.
ϵ) / dϵ − df dϵ D (ϵ) with D (ϵ) density of electronic states (DOS).So D S can be roughly regarded as an weighted averaged DOS with weight − df dϵ D (ϵ).With | g ↑↓ | 2 and D S , we have the approximate relation for spin relaxation rate,