Spin-phonon relaxation from a universal ab initio density-matrix approach

Designing new quantum materials with long-lived electron spin states urgently requires a general theoretical formalism and computational technique to reliably predict intrinsic spin relaxation times. We present a new, accurate and universal first-principles methodology based on Lindbladian dynamics of density matrices to calculate spin-phonon relaxation time of solids with arbitrary spin mixing and crystal symmetry. This method describes contributions of Elliott-Yafet and D’yakonov-Perel’ mechanisms to spin relaxation for systems with and without inversion symmetry on an equal footing. We show that intrinsic spin and momentum relaxation times both decrease with increasing temperature; however, for the D’yakonov-Perel’ mechanism, spin relaxation time varies inversely with extrinsic scattering time. We predict large anisotropy of spin lifetime in transition metal dichalcogenides. The excellent agreement with experiments for a broad range of materials underscores the predictive capability of our method for properties critical to quantum information science. First-principles calculations can help design and understand the behaviour of quantum technologies, but this requires the development of accurate methods to predict material properties. Here the authors present a method for calculating the spin-phonon relaxation time of general systems, a key quantity for spintronic devices.

T he manipulation of electron spins is of increasing interest in a wide-range of emerging technologies. The rapidly growing field of spintronics seeks to control spin as the unit of information instead of charge in devices such as spin transistors 1 . Quantum information technologies seek to utilize localized spin states in materials both as single-photon emitters [2][3][4] and as spin-qubits for future integrated quantum computers 5 . Both spintronics and quantum information applications therefore demand a quantitative understanding of spin dynamics and transport in metals and semiconductors. Recent advances in circularly polarized pump-probe spectroscopy 6 , spin injection, and detection techniques 7 have enabled increasingly detailed experimental measurement of spin dynamics in solid-state systems 8,9 . However, a universal first-principles theoretical approach to predict spin dynamics, quantitatively interpret these experiments and design new materials has remained out of reach.
A key metric of useful spin dynamics is the spin relaxation time τ s 1 . For example, spin-based quantum information applications require τ s exceeding milliseconds for reliable qubit operation. Consequently, accurate prediction of τ s in general materials is an important milestone for first-principles design of quantum materials. Spin-spin 10 , spin-phonon 11 , and spin-impurity scatterings, all contribute to spin relaxation, but spin-phonon scattering sets the intrinsic material limitation and is typically the dominant mechanism at room temperature 1 . Further, spin-phonon relaxation arises from a combination of spin-orbit coupling (SOC) and electron-phonon scattering, and is traditionally described by two mechanisms. First, the Elliott-Yafet (EY) mechanism involves spin-flip transitions between pairs of Kramers-degenerate states due to SOC-based spin-mixing of these states 12,13 . Second, the D'yakonov-Perel' (DP) mechanism in systems with broken inversion symmetry involves electron spins precessing between scattering events due to the SOCinduced internal effective magnetic field 14 .
Previous theoretical approaches have extensively investigated these two distinct mechanisms of intrinsic spin-phonon relaxation using model Hamiltonians in various materials 15 . These methods require parametrization for each specific material, which needs extensive prior information about the material and specialized computational techniques, and often only studies one mechanism at a time. Furthermore, most of these approaches require the use of simplified formulae 12,14 and make approximations to the electronic structure (e.g. low spin-mixing) or electron-phonon matrix elements 15 . This limits the generality and reliability of these approaches for complex materials, particularly for the DP mechanism, where various empirical relations are widely employed to estimate τ s 1 . Sophisticated methods based on spin susceptibility 16 and time evolution of density matrix 17 also rely on suitably chosen model Hamiltonians with empirical scattering matrix elements. Therefore, while these methods provide some mechanistic insight, they do not serve as predictive tools of spin relaxation time for the design of new materials.
A general first-principles technique to predict spin-phonon relaxation in arbitrary materials is therefore urgently needed. Previous first-principles studies have addressed the EY mechanism in centrosymmetric semiconductors 18,19 and metals 20 . These methods 18,20 rely on defining a pseudospin that allows the use of Fermi's golden rule (FGR) with only spin-flip transitions 13 . However, this is only well-defined for cases with weak spinmixing such that eigenstates within each Kramers-degenerate pair can be chosen to have small spin-minority components, precluding the study of spin relaxation of states with strong spinmixing, e.g. holes in silicon and noble metals. First-principles calculations have not yet addressed systems with such complex degeneracy structures, where the simple picture of spin-flip matrix elements in a FGR breaks down, or systems without inversion symmetry that do not exhibit Kramers degeneracy. Therefore, a more general first-principles technique without the material-specific simplifying assumptions of these previous approaches is now necessary.
In this work, we establish a new, accurate and unified firstprinciples technique for predicting spin relaxation time based on perturbative treatment of the Lindbladian dynamics of density matrices 21 . Importantly, by covering previously disparate mechanisms (e.g. EY and DP) in a unified framework, this technique is applicable to all materials regardless of dimensionality, symmetry (especially inversion) and strength of spin-mixing, which is critical for new material design. All SOC effects are included self-consistently (and non-perturbatively) in the ground-state eigensystem at the density functional theory (DFT) level, and we predict τ s through a universal rate expression without the need to invoke real-time dynamics. In this article, we first introduce our theoretical framework based on first-principles density-matrix dynamics, and then show prototypical examples of τ s for the broad range of systems, including three with inversion symmetry-silicon, iron, and graphene, and three without inversion symmetry-monolayer MoS 2 , monolayer MoSe 2 , and bulk GaN, in excellent agreement with available experimental data. By doing so, we establish the foundation for quantum dynamics of open systems from first-principles to facilitate the design of quantum materials.

Results
Theory. The key to treating arbitrary state degeneracy and spinmixing for spin relaxation is to switch to an ab initio densitymatrix formalism, which goes beyond specific cases such as Kramers degeneracy or Rashba-split model Hamiltonians. Specifically, we seek to work with density matrices of electrons alone, treating its interactions with an environment consisting of a thermal bath of phonons. In general, tracing out the environmental degrees of freedom in a full quantum Liouville equation of the density-matrix results in a quantum Lindblad equation. Specifically, for electron-phonon coupling 21 based on the standard Born-Markov approximation 22 that neglects memory effects in the environment, the Lindbladian dynamics in interaction picture reduces to where α is a combined index labeling electron wavevector k and band index n, λ is mode index and ± corresponds to q ¼ Ç k À k 0 ð Þ. n ± qλ n qλ þ 0:5 ± 0:5 and n qλ is phonon occupation. G qλ± αα 0 ¼ g qλ± αα 0 δ 1=2 ðε α À ε 0 α ± ω qλ Þ is the electron-phonon matrix element including energy conservation, where ω qλ is the phonon frequency.
This specific form of the Lindbladian dynamics preserves positive definiteness of the density matrix which is critical for numerical stability 21 . In addition, the energy-conserving δfunction above is regularized by a Gaussian with a width γ, which corresponds physically to the collision time. In some cases, the results depend on γ and γ → 0 is not the relevant limit 23 .
Here, the Lindblad master equation with finite smearing parameters corresponding to the collision time can be regarded as the best Markovian approximation to the exact dynamics 23 . In the case of spin relaxation, this is particularly important for systems that exhibit the DP mechanism, as we show below. Consequently, we consistently determine the smearing parameters from ab initio electron-phonon linewidth calculations throughout 24,25 .
The density-matrix formalism allows the computation of any observable such as number and spin density of carriers, and the inclusion of different relaxation mechanisms at time scales spanning femtoseconds to microseconds, which forms the foundation of the general relaxation time approach we discuss below. Given an exponentially relaxing measured quantity O ¼ TrðoρÞ, where o and ρ are the observable operator and the density matrix, respectively, we can define the relaxation rate Γ o and relaxation time where eq corresponds to the final equilibrium state. We note that even when the observables have additional cosðωtÞ oscillation factors, such as due to spin precession with periodicity of ω, the above equation remains an appropriate definition of the overall relaxation rate. For example, for a precessing and relaxing spin system with SðtÞ ¼ S 0 expðÀt=τÞ cosðωtÞ, the initial relaxation rate is _ Sð0Þ ¼ ÀS 0 =τ, which is exactly the same as that of a pure exponential relaxation.
The equilibrium density matrix in band space is ðρ eq Þ k nn 0 ¼ f kn δ nn 0 , where f kn are the Fermi occupation factors of electrons in equilibrium. Writing the initial density matrix ρ = ρ eq + δρ, assuming a small perturbation ||δρ|| ≪ ||ρ eq || and kdiagonal o and δρ, the Lindblad dynamics expression (Eq. (1)) and the definition (Eq. (2)) yield Here, the G is exactly as defined above in Eq. (1), but separating the wavevector indices (k, k 0 ) and writing it as a matrix in the space of band indices (n, n 0 ) alone. Similarly, o and δρ are also matrices in the band space, Tr n and † n are trace and Hermitian conjugate in band space, and ½o; G kk 0 o k G kk 0 À G kk 0 o k 0 , written using matrices in band space. Given an initial perturbation δρ and an observable o, Eq. (3) can now compute the relaxation of expectation value O from its initial value. Even for a specific observable like spin, several choices are possible for the initial perturbation corresponding directly to the experimental measurement scheme. Specifically for spin relaxation rate Γ s,i , the observable is the spin matrix S i labeled by Cartesian directions i = x, y, z, and the initial perturbed state should contain a deviation of spin expectation value from equilibrium. The most general (experiment-agnostic) choice for preparing a spin polarization is to assume that all other degrees of freedom are in thermal equilibrium, which can be implemented using a test magnetic field B i as a Lagrange multiplier for implementing a spin polarization constraint. With a corresponding initial perturbation Hamiltonian of H 1 = −2μ B B i S i /ℏ, where μ B is the Bohr magneton, perturbation theory yields In some cases, S i,k,mn ≈ 0 when ϵ km ≠ ϵ kn . For these cases, In such cases, we can further simplify Eq. (3) to the Fermi Golden rule-like expression, Note that the test field B i etc. drops out of the final expression and only serves to select the direction of the perturbation in the high-dimensional space of density matrices.
where Δs i;knk 0 n 0 s i;kn À s i;k 0 n 0 is the change in (diagonal) spin expectation value for a pair of states. Therefore, in this limit, Eq. (5) reduces to transitions between pairs of states, each contributing proportionally to the square of the corresponding spin change.
See Supplementary Note 1 and 2 for detailed derivations of the above equations. As we show in Supplementary Note 1 and 2, the above equations can be reduced to previous formulae with spin-flip matrix elements in Kramers-degenerate subspaces for systems with inversion symmetry and weak spin-mixing, such as conduction electron spin relaxation in bulk Si, similar to ref. 18 . However, Eq. (3) is much more general, applicable for systems with arbitrary degeneracy and crystal symmetry, and we therefore use it throughout for all results presented below. In addition, the overall framework can also be extended to other observables and can be made to correspond to specific measurement techniques that prepare a different initial density matrix e.g. a circularly polarized pump pulse.
Finally, note that in our first-principles method, all SOCinduced effects (such as the Rashba/Dresselhaus effects) are selfconsistently included in the ground-state eigensystem or the unperturbed Hamiltonian H 0 . This is essential to allow us to simulate τ s by a single rate calculation when there is broken inversion symmetry. On the other hand, if SOC does not enter into H 0 , as in previous work with model Hamiltonians, it must be treated as a separate term that provides an internal effective magnetic field. Consequently, those approaches require a coherent part of the time evolution to describe the fast spin precession induced by this effective magnetic field, which require explicit real-time dynamics simulations even to capture spin relaxation, going beyond a simple exponential decay as in Eq. (2). Using fully self-consistent SOC in a first-principles method is therefore critical to avoid this system-specific complexity and arrive at the universal approach outlined above.
Systems with inversion symmetry: Si and Fe. We first present results for systems with inversion symmetry traditionally described by a Elliot-Yafet spin-flip mechanism. Figure 1a shows that our predictions of electron spin relaxation time (τ s ) of Si as a function of temperature are in excellent agreement with experimental measurements 26,27 . Note that previous first-principles calculations 18 approximated spin-flip electron-phonon matrix elements from pseudospin wavefunction overlap and spinconserving electron-phonon matrix element, effectively assuming that the scattering potential varies slowly on the scale of a unit cell; we make no such approximation in our direct first-principles approach. Importantly, this allows us to go beyond the doubly degenerate Kramers-degenerate case of conduction electrons in Si. In contrast, holes in Si exhibit strong spin-mixing with spin-2/3 character and spin expectation values no longer close to ℏ/2. Figure 1b shows our predictions for the hole-spin relaxation time, which is much shorter than the electron case as a result of the strong mixing (450 fs for holes compared to 7 ns for electrons at 300 K) and is much closer to the momentum relaxation time. In addition, Fig. 1d shows that the change in spin expectation values (Δs) per scattering event (evaluated using Eq. (5)) has a broad distribution for holes in Si, indicating that they cannot be described purely by spin-flip transitions, while conduction electrons in Si predominantly exhibit spin-flip transitions with Δs = 1. We next consider an example of a ferromagnetic metal, iron, which exhibits a complex band structure not amenable for model Hamiltonian approaches. Previous first-principles calculations for ferromagnets employ empirical Elliott relation 28 or FGR formulae with spin-flip matrix elements specifically developed for metals or ferromagnets 20 . Here, we apply exactly the same technique used for the silicon calculations above and predict spin relaxation times in iron in good agreement with experimental measurements (Fig. 1c) [29][30][31] . Our Wannier interpolation also enables systematic and efficient Brillouin zone convergence of these predictions which were not possible previously. Similar to holes in Si, the Δs of Fe also exhibits a broad distribution extending from 0 to ℏ in the contribution to the total spin relaxation rate (Fig. 1d). Therefore, spin relaxation in transition metals are not purely spin-flip transitions, and we expect this effect to be even more pronounced in 4d and 5d metals with stronger SOC than the 3d magnetic metal considered here. Finally, Fig. 1a-c shows that τ s is approximately proportional to momentum relaxation time τ m for both Si and bcc Fe, which is expected for spin relaxation in systems with EY mechanisms 1 .
Systems with inversion symmetry: graphene. Graphene is of significant interest for spin-based technologies, and significant recent work with model Hamiltonians seeks to identify the fundamental limits of spin coherence in graphene 32 . Estimates vary widely from theoretical estimates on the order of microseconds to experiments ranging from picoseconds to nanoseconds [33][34][35][36] , with the discrepancies hypothesized to arise from faster extrinsic relaxation in experiments. However, previous model Hamiltonian studies required parametrization of approximate matrix elements, and focus on specific phonon modes (e.g. flexural modes) for spin-phonon relaxation.
Here we predict intrinsic electron-phonon spin relaxation time for free-standing graphene to firmly establish the intrinsic spin-phonon relaxation limit free of specific model choices or parameters. Figure 2 shows the predicted spin-phonon relaxation times as a function of temperature and Fermi level position. At room temperature, our calculated lifetimes are of the same magnitude (in microseconds) as previous predictions 33 indicating that faster relaxation is likely extrinsic in experiments. However, in addition to the flexural phonon mode 33,35 , in-plane acoustic (A) phonon modes have a strong and non-negligible contribution, while optical modes (O) have an overall smaller effect (Fig. 2b). We also find that the ratio between in-plane and out-of-plane spin relaxation times range from 0.5 to 0.7 (Fig. 2a, c), consistent with experimental measurements 34 . As evident from Fig. 2c, longer spin relaxation time of up to microseconds is achievable at low temperatures in pristine and free-standing graphene. However, at low temperatures, competing effects from substrates and disorder can make overall measured spin relaxation faster than theoretical predictions 35 . Finally note that while the ratio between in-plane and out-ofplane spin relaxation times is nearly 1/2, which is often considered to be a signature of the DP mechanism, freestanding graphene is inversion symmetric and does not exhibit the DP mechanism. Figure 2d shows that the spin relaxation time is mostly insensitive to the extrinsic scattering rates, instead of the linear relation (inverse relation with scattering time) expected for the DP mechanism, as discussed below in further detail. The spin relaxation of graphene may be switched to the DP regime by adding substrates or external electric fields to break inversion symmetry 36,37 , which will be investigated in detail using this theoretical framework in future work.
Systems without inversion symmetry: out-of-plane τ s of MoS 2 and MoSe 2 . The two-dimensional transition metal dichalcogenides (TMDs) exhibit extremely long-lived spin/valley polarization (over nanoseconds) 38 , with long valley-state persistence attributed to spin-valley locking effects. A fundamental understanding of spin/valley relaxation mechanisms is now required to utilize this degree of freedom for valleytronic computing 39 . Next we investigate spin relaxation τ s of systems without inversion symmetry from first-principles, starting with two TMD systems-monolayer In both systems, valence and conduction band edges at K and K′ valleys exhibit relatively large SOC band splitting, with nearly perfect out-of-plane spin polarization. Time-reversal symmetry further enforces opposite spin directions for the band-edge states at K and K′. Previous studies using model Hamiltonians consider the DP mechanism to dominate spin relaxation in these materials 17 , but in our first-principles approach, we do not need to a priori restrict our calculations to EY or DP limits.
In Fig. 3, we show the out-of-plane spin (τ s ) and momentum (τ m ) relaxation time of conduction electrons in two monolayer TMDs as a function of temperature, along with their intervalley/ intravalley contributions and experimental values. First, the overall agreement between our calculations and previous experiments by ultrafast pump-probe spectroscopy is excellent 38,40,41 . Note that ultrafast measurements of TMDs obtain coupled dynamics of spin and valley polarizations according to the selection rules with circularly polarized light, necessitating additional analysis to extract τ s , e.g., a phenomenological model fit to experimental curves in ref. 38 . On the other hand, our firstprinciples method simulates τ s directly without model or input parameters. This provides additional confidence in the experimental procedures of extracting τ s , and lends further insights into different scattering contributions in the dynamical processes as we show below. Moreover, special care is necessary when comparing with certain low temperature measurements with lightly doped samples, which access spin relaxation of excitons rather than individual free carriers, as discussed in refs. 42,43 ; we focus here on spin relaxation of free carriers. Next, comparing the relative contributions of intervalley and intravalley scattering for spin relaxation time, we find that the intravalley process dominates spin relaxation of conduction electrons in both TMDs: the intravalley only spin relaxation time (black squares) in Fig. 3 is nearly identical with the net spin relaxation time (red circles), while the intervalley contribution alone (blue triangles) is consistently more than an order of magnitude higher in relaxation time (lower in rate). Furthermore, with decreasing temperature, the relative contribution of the intervalley process decreases because the minimum phonon energies for wave vectors connecting the two valleys exceed 20 meV, and the corresponding phonon occupations become negligible at temperatures far below 300 K.
Previous theoretical studies of MoS 2 with model Hamiltonians 17 obtained (out-of-plane) τ s two orders of magnitude higher than our predictions, which agree with experimental data 38 . Such significant deviations are possibly because of the approximate treatments of electronic structure and electron-phonon coupling in their theoretical model. In addition, our first-principle calculations treat all phonon modes on an equal footing. Table 1 shows that the relative contributions of each phonon mode to τ s varies strongly with temperature. Full electron and phonon band structure is therefore vital to correctly describe spin-phonon relaxation with varying temperature, while model Hamiltonians that select specific phonon modes have limited range of validity 17 .
Hole-spin relaxation in MoS 2 and MoSe 2 has not been previously investigated in detail theoretically. Figure 4 presents our predictions of hole τ s and τ m in the two TMDs, indicating that hole τ s is much longer than that for electrons at all temperatures, exceeding 1 ns below 100 K. In contrast to the electron case, the intervalley process is relatively much more important and dominates spin relaxation at low temperature in MoS 2 and at all temperatures in MoSe 2 . This is because large SOC splitting at the valence band maximum makes the intravalley transition between two valence bands nearly impossible based on energy conservation in the electron-phonon scattering process. Experimental measurements also observe long spin relaxation times dominated by intervalley scattering in tungsten dichalcogenides 44 , which may facilitate applications in spintronic and valleytronic devices.
External magnetic fields can serve as tools tuning material properties 45 and are an inherent component of spin dynamics measurements 38,44 . Systems with broken inversion symmetry in particular may strongly respond to magnetic fields. We therefore investigate the effects of an external field B on τ s by introducing a Zeeman term (g s μ B /ℏ)B ⋅ S to the electronic Hamiltonian interpolated using Wannier functions (approximating g s ≈ 2), just prior to computing τ s with Eq. (3). Figure 5 shows that the out-of-plane τ s of conduction electrons of MoS 2 decreases with increasing in-plane magnetic field B x , in agreement with experimental work on MoS 2 38 and in general consistency with previous theoretical studies of τ s for systems with broken inversion symmetry 17,46 .
This strong magnetic field response has a simple intuitive explanation: in TMDs, the spin splitting of bands can be considered as the result of the internal effective magnetic field    Supplementary Fig. 4).
B soẑ due to broken inversion symmetry. Applying a finite B x perpendicular to B soẑ will cause additional spin-mixing and increase the spin-flip transition probability, thereby reducing the spin relaxation time. The degree of reduction depends on the detailed electronic structure of MoS 2 and MoSe 2 as shown in Supplementary Figs. 2 and 3: MoSe 2 exhibits a larger spin splitting of conduction bands and a higher internal magnetic field, and is therefore less affected by external B x . Similarly, hole-spin relaxation in both MoS 2 and MoSe 2 (not shown) exhibit very weak dependence on B x because of the large spin splitting and high internal effective magnetic field B so for valence band-edge states compared to those near the conduction band minimum. This insensitivity of hole τ s to magnetic fields is also consistent with experimental studies of hole τ s in WS 2 44 and WSe 2 47 . Finally, out-of-plane magnetic field B z has a negligible effect on spin relaxation for TMDs (not shown), unlike the in-plane magnetic field B x or B y . This is because electronic states around band edges are already polarized along the out-of-plane direction under a strong internal B soẑ . High experimental external magnetic fields~1 Tesla are relatively weak in contrast and only slightly change the spin polarization of the states, rather than introducing a spin-mixing that leads to spin relaxation.
Systems without inversion symmetry: in-plane τ s of MoS 2 . In all cases, the spin-phonon relaxation time decreases with increasing temperature, approximately proportional to the momentum relaxation time τ m . This is expected because both scattering rates, τ À1 s and τ À1 m , are proportional to phonon occupation factors which increase with temperature. The intrinsic inplane spin relaxation time (τ s,xx ) in MoS 2 also shows the same trend with temperature (Fig. 6a), but exhibits a fundamental  Specifically, Fig. 6b shows the dependence of spin relaxation times in MoS 2 conduction electrons as a function of extrinsic scattering rates, which enter Eq. (3) through an additional contribution to the smearing width γ of the energy-conserving δ-functions (in addition to the intrinsic electron-phonon contributions computed from first-principles). This additional smearing physically corresponds to the reduced lifetime and increased broadening of the electronic states in the material due to scattering against defects, impurities etc 37 . Importantly, the inplane spin relaxation time τ s,xx increases linearly with extrinsic scattering rate, or inversely with extrinsic scattering time, which is a hallmark of the DP mechanism of spin relaxation 48 .
Note that this inverse relation competes with the phonon occupation factors in determining the overall temperature dependence of spin relaxation time. At higher temperature, increased phonon occupation factors lower the intrinsic relaxation times of both carrier and spin, as stated above. The lowered intrinsic carrier relaxation time increases the finite smearing width in Eq. (3), which contributes towards increasing the spin relaxation time within the DP mechanism (inverse relation). However, the direct contribution of phonon occupation factors in the spin relaxation rate in Eq. (3) overwhelms this secondary change and results in a net decrease of spin relaxation time, consistent with all calculations above and experiments 49,50 .
In contrast with the in-plane case, the out-of-plane spin relaxation τ s,zz is mostly insensitive to the extrinsic scattering rate (and broadening γ), as all previous spin relaxation results in Kramers-degenerate materials discussed above (e.g. for graphene in Fig. 2d). Note that τ s,xx is also overall much shorter than τ s,zz , because the strong internal magnetic field in TMDs stabilizes spins in the z direction as discussed above. Large anisotropy in spin lifetimes due to a similar spin-valley locking effect has been theoretically predicted 51 and experimentally measured 52 previously in graphene-TMD interfaces as well.
Systems without inversion symmetry: GaN. Finally, we show spin relaxation in GaN as an archetypal example of the DP mechanism. Fig. 7 shows that both in-plane (τ s,xx ) and out-ofplane (τ s,zz ) spin lifetime of GaN are proportional to extrinsic scattering rates, or inversely proportional to extrinsic scattering time. Most importantly, the ratio between τ s,zz and τ s,xx is exactly 1/2 for this material, which is an additional feature of the conventional DP mechanism 1 . Note that, in contrast, the 2D TMDs are more complex due to strong SOC splitting and anisotropy, did not exhibit this 1/2 ratio, and exhibited the extrinsic scattering dependence only for in-plane spin relaxation. Overall, these results indicate that the general density-matrix formalism presented here elegantly captures the characteristic DP and EY mechanism limits, as well as complex cases that do not fit these limits, all on the same footing in a unified framework.

Discussion
In summary, we have demonstrated an accurate and universal first-principles method for predicting spin relaxation time of arbitrary materials, regardless of electronic structure, strength of spin-mixing and crystal symmetry (especially with/without inversion symmetry). Our work goes far beyond previous firstprinciples techniques based on a specialized Fermi's golden rule with spin-flip transitions and provides a pathway to an intuitive understanding of spin relaxation with arbitrary spin-mixing. In TMD monolayer materials, we clarify the roles of intravalley and intervalley processes, which are additionally resolved by phonon modes, in electron and hole-spin relaxation. We predict longlived spin polarization from resident carriers of MoS 2 and MoSe 2 and show their strong sensitivity of electron spin relaxation to inplane magnetic fields.
The predictive power of first-principles calculations is crucial for providing fundamental understanding of spin relaxation in new materials. The same technique can be applied to predict spin relaxation in realistic materials with or without defects useful for quantum technologies, wherever spin relaxation is dominated by electron-phonon scattering. We have already considered the general impact of disorder and electron-impurity scattering on spin-phonon relaxation through carrier broadening, but impurities can contribute an additional channel for spin relaxation, especially in the Kramers-degenerate case and at lower temperatures 18 . The extension of this technique to directly predict electron-impurity scattering for specific defects is relatively straightforward using supercell calculations, but computationally more demanding, while predicting the impact of electron-hole interaction [53][54][55] and electron-electron scattering 56,57 is additionally challenging. Finally, a robust understanding of ultrafast experiments may require simulation of real-time dynamics to capture initial state effects, probe wavelength effects and beyondsingle-exponential decay dynamics, which is a natural next step within the general Lindbladian density-matrix formalism presented here.

Methods
Computational details. All simulations are performed by the open-source planewave code -JDFTx 58 using pseudopotential method, except that the Born effective charges and dielectric constants are obtained from open-source code QuantumESPRESSO 59 . We firstly carry out electron structure, phonon and electron-phonon matrix element calculations in DFT using Perdew-Burke-Ernzerhof exchange-correlation functional 60 with relatively coarse k and q meshes. The phonon calculations are done using the supercell method. We have used supercells of size 7 × 7 × 7, 4 × 4 × 4, 6 × 6 × 1, 6 × 6 × 1, 6 × 6 × 1, 4 × 4 × 4 for silicon, BCC iron, graphene, monolayer MoS 2 , monolayer MoSe 2 and GaN, respectively, which have shown reasonable convergence for each system (<20% error bar in the final spin relaxation estimates). SOC is included through the use of the fully relativistic pseudopotentials 61 . For monolayer MoS 2 and MoSe 2 , the Coulomb truncation technique is employed to accelerate convergence with vacuum sizes 62 .
We then transform all quantities from plane wave to maximally localized Wannier function basis 63 and interpolate them 24,25,64,65 to substantially finer k and q meshes (with >3 × 10 5 total points) for lifetime calculations. Statistical errors of lifetime computed using different random samplings of k-points are found to be negligible (<1%). This Wannier interpolation approach fully accounts for polar terms in the electron-phonon matrix elements and phonon dispersion relations using the approaches of Verdi and Giustino 66 and Sohier et al. 67 for the 3D and 2D systems. In-plane (τ s,xx ) and out-of-plane (τ s,zz ) spin lifetime of GaN as a function of extrinsic scattering rates at 300 K. The experimental τ s,zz at 298 K is from ref. 69 .

Data availability
All relevant data are available from the authors upon request.

Code availability
First-principles methodologies available through open-source software, JDFTx 58 , and post-processing codes available from authors upon request.