# Towards exact predictions of spin-phonon relaxation times: an ab initio
implementation of open quantum systems theory

Alessandro Lunghi<EMAIL_ADDRESS>School of Physics, AMBER and CRANN
Institute, Trinity College, Dublin 2, Ireland

###### Abstract

Spin-phonon coupling is the main drive of spin relaxation and decoherence in
solid-state semiconductors at finite temperature. Controlling this interaction
is a central problem for many disciplines, ranging from magnetic resonance to
quantum technologies. Spin relaxation theories have been developed for almost
a century but often employ a phenomenological description of phonons and their
coupling to spin, resulting in a non-predictive tool and hindering our
detailed understanding of spin dynamics. Here we combine fourth-order time-
local quantum master equations with advanced electronic structure methods and
perform predictions of spin-phonon relaxation time for a series of solid-state
coordination compounds based on both transition metals and lanthanide Kramers
ions. The agreement between experiments and simulations demonstrates that an
accurate, universal and fully ab initio implementation of spin relaxation
theory is possible, thus paving the way to a systematic study of spin-phonon
relaxation in solid-state materials.

Relaxation of localized magnetic moments in condensed-phase is mainly driven
by three interactions: spin-spin, spin-conduction electron and spin-phonon
couplingWu  _et al._ (2010). In the case of semiconductors at temperatures
above a few K, the latter is the most important one. In these circumstances,
the oscillatory motion of atoms acts as a time-dependent perturbation on the
spin degrees of freedom, due to the presence of spin-orbit and spin-spin
interactions, and leads to spin thermalization and loss of spin coherence. On
the one hand a deep understanding of spin-phonon relaxation and spin-phonon
coupling has the potential to open up new ways to address central problems for
modern technologies such as magnetic information storageMaehrlein  _et al._
(2018), spintronicsŽutić  _et al._ (2004), quantum information
scienceWolfowicz  _et al._ (2021); Atzori and Sessoli (2019), magnetic
resonanceEaton and Eaton (2002), and more. On the other hand, the study of
spin relaxation in condensed-phase represents a fundamental benchmark for our
understanding of the theory of open quantum systems. Although a formally exact
description of the quantum dynamics of open systems is achievable through
time-convolutionlessBreuer  _et al._ (2001); Timm (2011) or Nakajima-
ZwanzigKoller  _et al._ (2010) equations, their application is often done
either in a parametric fashion or for simplified models, therefore escaping
the ultimate benchmark of a full ab initio implementation for realistic
systems. Spin-phonon coupling represents a paradigmatic example of open-
quantum-system dynamics and offers a unique test-bed for open quantum system
theory and ab initio simulations at the same time.

Spin-phonon relaxation has now been debated for almost a century starting with
the seminal works of WallerWaller (1932), Van VleckVan Vleck (1940),
RedfieldRedfield (1957) and OrbachOrbach (1961). Redfield relaxation theory
and its adaptations form the fundamental building block of our understanding
of spin relaxation and have served as the basis to model both nuclear and
electronic spin relaxation experiments. Redfield theory is generally adapted
to the case of solid-state systems by accounting for a Debye-like phonon
density of state and by including all the details about spin-phonon coupling
under a single average phenomenological coupling constantOrbach (1961).
Despite the incredible usefulness of such a formalism for interpreting
experiments, these approximations have the major drawback of making it
impossible to use the theory as a real predictive tool.

Figure 1: Molecular structures. Left, central and right panels report the
structures of 1, 2 and 3, respectively. Color code: pink for Vanadium, Purple
for Cobalt, light green for Dysprosium, green for Carbon, blue for Nitrogen,
red for Oxygen, yellow for Sulphur and white for Hydrogen.

Only recently, fully-ab initio methods to predict spin-phonon relaxation in
solid-state materials has been proposedEscalera-Moreno  _et al._ (2017);
Lunghi  _et al._ (2017a); Goodwin  _et al._ (2017); Lunghi and Sanvito
(2019a); Park  _et al._ (2020); Xu  _et al._ (2020). At the heart of these
methods there is a combination of Redfield theory of relaxation and electronic
structure methods. The latter are used to define all the coefficients needed
to set up the former on a material-specific case and thus enable a non-
parametric treatment of all the relevant interactions, such as spin and phonon
spectra, and all the spin-phonon coupling coefficients. Applications of this
method to magnetic molecules have recently generated large interest, with
examples of studies in both molecular qubitsEscalera-Moreno  _et al._ (2017);
Albino  _et al._ (2019); Lunghi and Sanvito (2019a); Mirzoyan and Hadt (2020);
Lunghi and Sanvito (2020a); Kazmierczak  _et al._ (2021) and single-ion
magnetsLunghi  _et al._ (2017a, b); Goodwin  _et al._ (2017); Moseley  _et
al._ (2018); Ullah  _et al._ (2019); Lunghi and Sanvito (2020b); Reta  _et
al._ (2021); Blockmon  _et al._ (2021); Briganti  _et al._ (2021). Indeed,
coordination compounds offer a versatile playground to test relaxation
theories as their chemical structure and properties can be finely controlled
and characterized experimentallyZabala-Lekuona  _et al._ (2021); Moreno-Pineda
and Wernsdorfer (2021). At the same time they allow for the use of advanced
electronic structure methods to describe their lattice and magnetic
propertiesAtanasov  _et al._ (2015); Ungur and Chibotaru (2017); Neese  _et
al._ (2019).

However, until now, a fully quantitative prediction of relaxation times have
been elusive and deviations of orders of magnitude with respect to
experimental data have been observed. Such a situation casts shadows on our
understanding of spin-phonon relaxation and on the reliability of ab initio
methods for blind prediction of spin relaxation time in absence of
experimental validation.

In this contribution we advance the field of ab initio spin dynamics by
demonstrating that quantitative predictions of spin relaxation rate can be
obtained for realistic solid-state compounds. We review the derivation of the
dynamical equations for the spin reduced density matrix in solid-state
materialsLunghi and Sanvito (2019a, 2020a, 2020b) and individuate the key
importance of including all the terms up to fourth-order perturbation theory
and accounting for the dynamics of the entire density matrix, i.e. including
the explicit description of coherence terms’ dynamics. We perform simulations
for three systems: (1) an $S=1/2$ molecular qubit [VO(dmit)2]2-
(dmit=1,3-dithiole-2-thione-4,5-dithiolate)Atzori  _et al._ (2016), (2) an
$S=3/2$ mononuclear complex [CoL2]2- (H2L =
1,2-bis(methanesulfonamido)benzene)Rechkemmer  _et al._ (2016); Rechkemmer
(2016), and (3) a $J=15/2$ mononuclear complex [DyCp${}^{\textrm{ttt}}_{2}$]+
(Cp${}^{\textrm{ttt}}$=[C5H${}_{5}^{\;\;\textrm{t}}$Bu3-1,2,4])Goodwin  _et
al._ (2017). The molecular structures of (1)-(3) are reported in Fig. 1. These
compounds represent the state-of-the-art in molecular magnetism for what
concerns molecular qubits and single-molecule magnets with long spin lifetime.
Moreover, their properties span a large range of magnetic splitting and
relaxation times, providing an extremely robust proof-of-concept of the
accuracy of ab initio spin dynamics.

## Results

Ab initio theory of spin-phonon relaxation.

The total system of spin plus phonon is described by the Hamiltonian

$\hat{H}=\hat{H}_{\mathrm{s}}+\hat{H}_{\mathrm{ph}}+\hat{H}_{\mathrm{s-ph}}\,$
(1)

where the first two terms correspond to the spin and phonons Hamiltonians,
respectively, and the third one represents the coupling between the two
subsystems. Under a weak-coupling assumption, the interaction between a spin
system and an ensemble of phonons, $q_{\alpha}$, can be modelled with a spin-
phonon coupling Hamiltonian of the form,

$\displaystyle\hat{H}_{\mathrm{s-ph}}=$
$\displaystyle\sum_{\alpha}\left(\frac{\partial\hat{H}_{\mathrm{s}}}{\partial
q_{\alpha}}\right)_{0}q_{\alpha}(t)$ (2) $\displaystyle+$
$\displaystyle\sum_{\alpha\geq\beta}\left(\frac{\partial^{2}\hat{H}_{\mathrm{s}}}{\partial
q_{\alpha}\partial q_{\beta}}\right)_{0}q_{\alpha}(t)q_{\beta}(t)\>,$ (3)

where the summations are understood to run over both phonon band indexes and
reciprocal lattice $\mathbf{q}$-points. The first- and second-order
derivatives of $\hat{H}_{\mathrm{s}}$ are generally called spin-phonon
coupling coefficients. The use of the spin Hamiltonian formalism in Eqs. 1, 2,
and 3 does not constitute an approximation as long as an exact mapping between
the low-lying electronic states of the system and the eigenstates of the spin
Hamiltonian exists, which is generally the caseAtanasov  _et al._ (2015);
Ungur and Chibotaru (2017).

The formalism of time-convolutionless master equations provides an exact
description of the time-evolution of the spin reduced density matrix in the
interaction pictureBreuer  _et al._ (2001); Timm (2011), $\hat{\rho}_{s}$,
i.e. once the explicit dynamics of the phonon degrees of freedom have been
integrated out from the total Liouville equation

$\frac{d\hat{\rho}_{\mathrm{s}}}{dt}=\hat{\hat{\mathbf{R}}}(t)\hat{\rho}_{\mathrm{s}}(t)\>.$
(4)

The relaxation super-operator $\hat{\hat{\mathbf{R}}}(t)$ can be
systematically expanded in a perturbative wayBreuer  _et al._ (2001); Timm
(2011),
$\hat{\hat{\mathbf{R}}}(t)=(\hat{\hat{\mathbf{R}}}\mathbf{2}(t)+\hat{\hat{\mathbf{R}}}\mathbf{4}(t)+...)$,
and in principles it can be implemented exactly if all the terms of Eq. 1 are
known. Assuming a separation of timescales between spin relaxation and the
intrinsic relaxation rates of the phonons, mostly due to anharmonic phonon-
phonon scattering, it is possible to assume that the latter are always at
thermal equilibrium and therefore apply the Born-Markov approximation. Namely,
we assume that phonons have a much larger specific heat than the spin system
(no phonon bottleneckTesi  _et al._ (2016)) and that they are able to rapidly
equilibrate while exchanging energy with the spin. Under these assumptions and
considering only the second-order contribution
$\hat{\hat{\mathbf{R}}}\mathbf{2}(t)$, Eq. 4 becomes identical to the well-
known Redfield equationsLunghi and Sanvito (2019a)

$\frac{d\rho^{\mathrm{s}}_{ab}(t)}{dt}=\sum_{cd}e^{i(\omega_{ac}+\omega_{db})t}R2^{n-\mathrm{ph}}_{ab,cd}\rho^{\mathrm{s}}_{cd}(t)\>.$
(5)

Considering only the linear term of Eq. 2,
$\hat{\hat{\mathbf{R}}}\mathbf{2}^{n-\mathrm{ph}}$ becomes

$\displaystyle R2^{1-\mathrm{ph}}_{ab,cd}=-\frac{\pi}{2\hbar^{2}}$
$\displaystyle\sum_{\alpha}\Big{\\{}\sum_{j}\delta_{bd}V^{\alpha}_{aj}V^{\alpha}_{jc}G^{1-\mathrm{ph}}(\omega_{jc},\omega_{\alpha})-V^{\alpha}_{ac}V^{\alpha}_{db}G^{1-\mathrm{ph}}(\omega_{bd},\omega_{\alpha})$
$\displaystyle-V^{\alpha}_{ac}V^{\alpha}_{db}G^{1-\mathrm{ph}}(\omega_{ac},\omega_{\alpha})+\sum_{j}\delta_{ca}V^{\alpha}_{dj}V^{\alpha}_{jb}G^{1-\mathrm{ph}}(\omega_{jd},\omega_{\alpha})\Big{\\}}\>,$
(6)

where the terms $V^{\alpha}_{ab}$ are a short-hand notation for $\langle
a|(\partial\hat{H}_{\mathrm{s}}/\partial q_{\alpha})|b\rangle$, and
$|a\rangle$ and $|b\rangle$ are eigenstates of $\hat{H}_{\mathrm{s}}$.
Finally, $\omega_{ab}=(E_{a}-E_{b})/\hbar$, where $E_{a}$ and $E_{b}$ are
eigenvalues of $\hat{H}_{\mathrm{s}}$. The function $G^{\mathrm{1-ph}}$ that
appears in Eq. 6 is the Fourier transform of the single-phonon correlation
function and reads

$G^{\mathrm{1-\mathrm{ph}}}(\omega,\omega_{\alpha})=\delta(\omega-\omega_{\alpha})\bar{n}_{\alpha}+\delta(\omega+\omega_{\alpha})(\bar{n}_{\alpha}+1)\>,$
(7)

where
$\bar{n}_{\alpha}=[\mathrm{exp}(\hbar\omega_{\alpha}/\mathrm{k_{\mathrm{B}}}T)-1]^{-1}$
is the Bose-Einstein distribution of thermal population,
$\hbar\omega_{\alpha}$ is the $\alpha$-phonon energy and $\mathrm{k_{B}}$ is
the Boltzmann constant. The function $G^{1-\mathrm{ph}}$ accounts for the
spectral density and population of phonons. Eqs. 5 and 6, describe one-phonon
resonant spin transitions, such as those responsible for Direct and Orbach
relaxation. Considering now the quadratic term of Eq. 3,
$\hat{\hat{\mathbf{R}}}\mathbf{2}^{n-\mathrm{ph}}$ becomesLunghi and Sanvito
(2020a)

$\displaystyle R2^{2-\mathrm{ph}}_{ab,cd}=-\frac{\pi}{4\hbar^{2}}$
$\displaystyle\sum_{\alpha\geq\beta}\Big{\\{}\sum_{j}\delta_{bd}V^{\alpha\beta}_{aj}V^{\alpha\beta}_{jc}G^{2-\mathrm{ph}}(\omega_{jc},\omega_{\alpha},\omega_{\beta})-V^{\alpha\beta}_{ac}V^{\alpha\beta}_{db}G^{2-\mathrm{ph}}(\omega_{bd},\omega_{\alpha},\omega_{\beta})$
$\displaystyle-V^{\alpha\beta}_{ac}V^{\alpha\beta}_{db}G^{2-\mathrm{ph}}(\omega_{ac},\omega_{\alpha},\omega_{\beta})+\sum_{j}\delta_{ca}V^{\alpha\beta}_{dj}V^{\alpha\beta}_{jb}G^{2-\mathrm{ph}}(\omega_{jd},\omega_{\alpha},\omega_{\beta})\Big{\\}}\>,$
(8)

where $V^{\alpha\beta}_{ab}$ now stands for $\langle
a|(\partial^{2}\hat{H}_{\mathrm{s}}/\partial q_{\alpha}\partial
q_{\beta})|b\rangle$. It should be noted that the term $\alpha=\beta$ is only
included for double absorption/emission terms (vide infra). The function
$G^{2-\mathrm{ph}}$ accounts for three possible processes involving two
phonons: absorption of two phonons, emission of two phonons and simultaneous
emission of one phonon and absorption of another one

$\displaystyle G^{2-\mathrm{ph}}(\omega_{ba},\omega_{\alpha},\omega_{\beta})=$
$\displaystyle\delta(\omega_{ba}-\omega_{\beta}+\omega_{\alpha})\bar{n}_{\beta}(\bar{n}_{\alpha}+1)+\delta(\omega_{ba}+\omega_{\beta}-\omega_{\alpha})(\bar{n}_{\beta}+1)\bar{n}_{\alpha}+$
(9)
$\displaystyle\delta(\omega_{ba}-\omega_{\beta}-\omega_{\alpha})\bar{n}_{\beta}\bar{n}_{\alpha}+\delta(\omega_{ba}+\omega_{\beta}+\omega_{\alpha})(\bar{n}_{\beta}+1)(\bar{n}_{\alpha}+1)\>.$
(10)

As evident from Eq. 10, in two-phonon transitions, only the sum or difference
of the two phonons’ energies must be resonant with the spin transition, while
the energy of the single vibrations could be anything. Spin relaxation
mechanisms due to two-phonon processes generally falls under the definition of
Raman relaxation mechanism.

At this point of the derivation of the Redfield equations, the secular
approximation is generally performed. This approximation involves setting to
zero the matrix elements of $R2^{n-\mathrm{ph}}_{ab,cd}$ for which
$(\omega_{ac}+\omega_{db})\neq 0$. This is justified by the fact that the
terms $(\omega_{ac}+\omega_{db})$ appear as arguments of an oscillatory term
in Eq. 5. These oscillating terms are averaged out for times commensurable to
the spin relaxation time, $\tau$, provided that the period of the natural spin
oscillations ($\omega_{ab}$) is much shorter than $\tau$. Bearing in mind that
spin-phonon relaxation in magnetic molecules often occurs on timescales of ns
to hundreds of seconds, depending on the value of $T$ and the specific system,
the secular approximation is generally fulfilled for systems with zero-field
splitting or in the presence of small external or dipolar magnetic fields.
Under these conditions, only the terms $R2^{n-\mathrm{ph}}_{aa,bb}$ and
$R2^{n-\mathrm{ph}}_{ab,ab}$, namely population transfer and coherence
relaxation terms, contribute to the dynamics of $\hat{\rho}_{s}$. However,
when dealing with Kramers systems in zero external magnetic field, additional
terms of $R2^{n-\mathrm{ph}}_{ab,cd}$ also survive the secular approximation
in virtue of the spectrum degeneracy that makes additional terms
$(\omega_{ac}+\omega_{db})$ vanish. For instance, the population-coherence
transfer process $R2^{n-\mathrm{ph}}_{aa,bc}$, with $E_{b}=E_{c}$, and
coherence transfer terms $R2^{n-\mathrm{ph}}_{ad,bc}$, with $E_{b}=E_{c}$ and
$E_{a}=E_{d}$, would not be washed out by the secular approximation and could
in principles contribute to the spin-phonon dynamicsEastham  _et al._ (2016);
Gösweiner  _et al._ (2020).

Extending the expansion of $\hat{\hat{\mathbf{R}}}(t)$ up to the fourth-order
and considering only the linear term of spin-phonon coupling, we obtain
another source of two-phonon spin transitions. A set of equations describing
the diagonal population terms of $\hat{\rho}_{s}$ due to this contribution
readsTimm (2011)

$\frac{d\rho^{\mathrm{s}}_{aa}(t)}{dt}=\sum_{bb}R4^{2-\mathrm{ph}}_{aa,bb}\rho^{\mathrm{s}}_{bb}(t)\>,$
(11)

where

$\displaystyle
R4^{2-\mathrm{ph}}_{aa,bb}=\frac{\pi}{2\hbar^{2}}\sum_{\alpha\geq\beta}\left[A_{\alpha\beta}W^{--}_{ab}(\alpha\beta)+A_{\alpha\beta}W^{++}_{ab}(\alpha\beta)+B_{\alpha\beta}W^{+-}_{ab}(\alpha\beta)+B_{\alpha\beta}W^{-+}_{ab}(\alpha\beta)\right]\>,\>\text{and}$
(12) $\displaystyle W^{--}_{ab}(\alpha\beta)=\Big{|}\sum_{c}\frac{\langle
a|\hat{V}_{\alpha}|c\rangle\langle
c|\hat{V}_{\beta}|b\rangle}{E_{c}-E_{b}-\hbar\omega_{\beta}}+\frac{\langle
a|\hat{V}_{\beta}|c\rangle\langle
c|\hat{V}_{\alpha}|b\rangle}{E_{c}-E_{b}-\hbar\omega_{\alpha}}\Big{|}^{2}\bar{n}_{\alpha}\bar{n}_{\beta}\delta(\omega_{ab}-\omega_{\alpha}-\omega_{\beta})$
(13) $\displaystyle W^{++}_{ab}(\alpha\beta)=\Big{|}\sum_{c}\frac{\langle
a|\hat{V}_{\alpha}|c\rangle\langle
c|\hat{V}_{\beta}|b\rangle}{E_{c}-E_{b}+\hbar\omega_{\beta}}+\frac{\langle
a|\hat{V}_{\beta}|c\rangle\langle
c|\hat{V}_{\alpha}|b\rangle}{E_{c}-E_{b}+\hbar\omega_{\alpha}}\Big{|}^{2}(\bar{n}_{\alpha}+1)(\bar{n}_{\beta}+1)\delta(\omega_{ab}+\omega_{\alpha}+\omega_{\beta})$
(14) $\displaystyle W^{+-}_{ab}(\alpha\beta)=\Big{|}\sum_{c}\frac{\langle
a|\hat{V}_{\alpha}|c\rangle\langle
c|\hat{V}_{\beta}|b\rangle}{E_{c}-E_{b}+\hbar\omega_{\beta}}+\frac{\langle
a|\hat{V}_{\beta}|c\rangle\langle
c|\hat{V}_{\alpha}|b\rangle}{E_{c}-E_{b}-\hbar\omega_{\alpha}}\Big{|}^{2}\bar{n}_{\alpha}(\bar{n}_{\beta}+1)\delta(\omega_{ab}-\omega_{\alpha}+\omega_{\beta})$
(15) $\displaystyle W^{-+}_{ab}(\alpha\beta)=\Big{|}\sum_{c}\frac{\langle
a|\hat{V}_{\alpha}|c\rangle\langle
c|\hat{V}_{\beta}|b\rangle}{E_{c}-E_{b}-\hbar\omega_{\beta}}+\frac{\langle
a|\hat{V}_{\beta}|c\rangle\langle
c|\hat{V}_{\alpha}|b\rangle}{E_{c}-E_{b}+\hbar\omega_{\alpha}}\Big{|}^{2}(\bar{n}_{\alpha}+1)\bar{n}_{\beta}\delta(\omega_{ab}+\omega_{\alpha}-\omega_{\beta})\>,$
(16)

where $A_{\alpha\beta}=(1-3/4\delta_{\alpha\beta})$ and
$B_{\alpha\beta}=(1-1/2\delta_{\alpha\beta})$. One can recognize in Eq. 16 the
same contributions of $G^{2-ph}$, except for different pre-factors, which are
now accounting for the presence of excited spin states, i.e. the virtual
states often invoked in perturbation theory. The processes describe by Eq. 16
also contributes to the so-called Raman relaxation. It is important to remark
that Eq. 16 should also include additional terms that cancel the divergence of
the denominatorsKoller  _et al._ (2010); Timm (2011). We neglect those terms
in light of the fact that the numerical sparsity of the phonon spectrum never
leads to such divergences. Importantly, Eq. 16 differs from a previously
presented one, where not all the terms had been includedLunghi and Sanvito
(2020b); Briganti  _et al._ (2021). Finally we note that Eq. 16 does not
account for the dynamics of coherence terms and is therefore only strictly
valid for non-degenerate spin spectra.

As long as the Born-Markov and secular approximations hold, the equations just
detailed provide an exact description of spin relaxation up to two-phonon
processes. The only remaining challenge to their implementation lies in the
definition of the many coefficients that enter these equations. This challenge
is tackled with electronic structure methods. Methods such as Complete Active
Space Self Consistent Field (CASSCF) and Density Functional Theory (DFT) have
now reached a high degree of sophistication in the prediction of spin
Hamiltonian parameters of magnetic moleculesAtanasov  _et al._ (2015); Ungur
and Chibotaru (2017) as well to predict phonons’ frequencies and normal modes
of vibrationsGarlatti  _et al._ (2020).

Spin-phonon relaxation in $S=1/2$ systems.

Compound (1) is a V4+ coordination complex with an $S=1/2$ ground state and a
nuclear spin $I=7/2$Atzori  _et al._ (2016). The spin Hamiltonian can be
modelled as

$\hat{H}_{\mathrm{s}}=\mu_{\mathrm{B}}\vec{\mathbf{B}}\cdot\mathbf{g}\cdot\vec{\mathbf{S}}+\gamma_{\mathrm{N}}\vec{\mathbf{B}}\cdot\mathbf{\vec{I}}+\vec{\mathbf{S}}\cdot\mathbf{A}\cdot\vec{\mathbf{I}}\>,$
(17)

where $\mu_{\mathrm{B}}$ is the electron’s Bohr magneton, $\mathbf{g}$ is the
effective electronic Landè tensor, $\gamma_{\mathrm{N}}$ is the nuclear
gyromagnetic factor of 51V of and $\mathbf{A}$ is the hyperfine coupling
tensor. (1) crystallizes in a monoclinic unit-cell containing four molecular
units and eight tetraphenylphosphonium counter-ions. The entire unit-cell was
optimized with periodic-DFT (pDFT) and used for phonons calculations as
detailed in the Computational Methods. The tensors $\mathbf{g}$ and
$\mathbf{A}$ were computed with DFT after a structural optimization. Predicted
values of $(g_{xx},g_{yy},g_{zz})=(1.976,1.986,1.988)$ are in perfect
agreement with experimental onesAtzori  _et al._ (2016). Computed values of
$\mathbf{A}$ are found sensitive to the choice of scalar relativistic
corrections used in the DFT simulation, with ZORA underestimating the values
up to 23% and DKH overestimating them up to 38%. All the terms of the spin-
phonon coupling Hamiltonian of (1), including second order terms, are computed
by taking the numerical derivative of all the parameters in Eq. 17 that depend
on molecular coordinates and affect spin dynamics, namely $\mathbf{g}$ and
$\mathbf{A}$. Notably, here we perform the calculation of second-order spin-
phonon coupling coefficients fully ab initio and without using any machine-
learning interpolator as done in previous workLunghi and Sanvito (2020a). We
neglect the contributions of Eq. 6 to relaxation, which are known to lead to
Direct relaxationLunghi and Sanvito (2019a). The latter only overcome two-
phonon relaxation for $T<10$ KAtzori  _et al._ (2016); Lunghi and Sanvito
(2020a). In that regime experimental values of $\tau$ are affected by spin
diffusion, which would make the comparison with simulations not
straightforward. Eqs. 5 and 8 are solved for an initial condition of the spin
density matrix where a $\pi$ rotation of the electronic spin is applied to the
canonical equilibrium distribution. This initial state is chosen to mimic the
experimental condition of an inversion recovery electron paramagnetic
resonance (EPR) experiment used to measure T1. The time-dependence of the
magnetization of the electronic spin is computed to reach equilibrium
following a stretched-exponential decay $M_{z}(t)=(M_{z}(0)-M_{z}^{eq})$
exp$[-(t/\tau)^{\beta}]+M_{z}^{eq}$ with $\beta\sim 0.7-0.9$ for $A$-driven
relaxation and $\beta=1$ for $g$-driven relaxation. Fig. 2 reports the
comparison between the values of T1 measured with X- and Q-Band EPR on
magnetically diluted samples, together with the predicted values of $\tau$ as
a function of temperature by assuming the modulation of $\mathbf{A}$ (top
panel) and $\mathbf{g}$ (bottom panel) as the sole relaxation mechanisms.

Figure 2: Spin-phonon relaxation time for (1). Simulated values of $\tau$ for
(1) including the contribution of A-tensor (top panel) and g-tensor (bottom
panel) modulations to the Raman mechanism of relaxation. Results in the top
panel are the average between calculations carried out with the DKH and ZORA
scalar relativistic corrections and the sky-blue shade represents the two
limiting values. Experimental values correspond to the T1 times obtained from
inversion recovery EPR experiments at X- and Q-bands. Figure 3: Redfield
Transition Rates for (2). Transition rates among general elements of the
density matrix predicted by the Redfield theory employing the non-diagonal
secular approximation and diagonal secular approximation are reported in the
left and central panel, respectively. The right panel reports the same rates
for the non-diagonal secular approximation after rotating the molecular
geometry and applying a small external field. The transitions rates have been
computed at $T=40$ K. The first four elements correspond to the population
terms, and all the remaining corresponds to coherence terms.

Results obtained considering the modulation of the hyperfine interaction show
a good agreement with inversion recovery experiments, both in terms of
absolute values and external field dependence. Depending on the fine details
of the DFT calculations agreement up to a factor five and not exceeding one
order of magnitude is observed. Conversely, the modulation of the Zeeman
interaction is shown to overestimate spin lifetime by at least one order of
magnitude. More importantly, the relaxation rate due to the modulation of g is
found to depend on the magnitude of the external field as $B^{2}$, which does
not fit the generally observed experimental trend for Raman relaxation in this
class of complexes. The computational results for (1) are in qualitative
agreement with a previous study on the crystal of VO(acac)2Lunghi and Sanvito
(2020a). However, differently from the VO(acac)2 case, T1 data obtained on
diluted crystals for (1) are availableAtzori  _et al._ (2016), making it
possible to demonstrate that once cross-relaxation mediated by dipolar
interactions is excluded, the deviation between simulations and experiments is
below one order of magnitude and the $\tau$ vs $T$ profile is nicely
reproduced. Indeed, the residual deviation is consistent with the lack of
phonons outside the $\Gamma$-point and the inaccuracies of electronic
structure methods. In order to determine the most important phonons for
relaxation, we perform simulations by including only vibrations in the low-
energy window and starting from the first optical mode at $\Gamma$-point, here
computed at $\sim 12$ cm-1. We find that optical phonons up to $\sim 50$ cm-1
are the main responsible for spin relaxation even at room temperature. This
can be interpreted as an effect of thermal population, which is always larger
for lowest-energy phonons. As discussed previously, low-energy vibrations of 1
are an admixture of molecular rotations and delocalized intra-molecular
distortionsAlbino  _et al._ (2019).

Spin-phonon relaxation in $S>1/2$ systems.

Next we address the spin dynamics of (2) as an example of relaxation in a
mononuclear transition-metal complex with large zero-field splitting. This
molecule contains a high-spin Co2+ ion that can be described with an effective
$S=3/2$ Hamiltonian

$\hat{H}_{\mathrm{s}}=D\hat{S}_{z}+E(\hat{S}_{x}-\hat{S}_{y})\>,$ (18)

where the large axial anisotropy term $D\sim-115$ cm-1 removes the degeneracy
of the two Kramers doublets (KDs) $M_{s}=\pm 3/2$ and $M_{s}=\pm
1/2$Rechkemmer  _et al._ (2016); Rechkemmer (2016). This compound crystallizes
in an orthorhombic lattice with four molecules in the unit cell and eight
NHEt${}_{3}^{+}$ counter-ion moleculesRechkemmer  _et al._ (2016). All the
spin-phonon coupling coefficients and the phonons for this compound were
computed in a previous work employing CASSCF and DFT-based machine-learning
force fields (MLFF), respectivelyLunghi and Sanvito (2020b). In this previous
work, the time evolution of $\hat{\rho}^{s}$ was studied with the diagonal
terms of Eq. 6 and Eq. 16, therefore neglecting the role of coherence terms.
In order to investigate the relevance of these terms we compute all the terms
of $R2^{1-\mathrm{ph}}_{ab,cd}$. Fig. 3a shows the value of
$R2^{1-\mathrm{ph}}_{ab,cd}$ when only the terms $R2^{1-\mathrm{ph}}_{aa,bb}$
and $R2^{1-\mathrm{ph}}_{ab,ab}$ are retained. We will refer to this
approximation as diagonal secular approximation from now on. Instead, Fig. 3b
shows the value of $R2^{1-\mathrm{ph}}_{ab,cd}$ when all terms that fulfil the
condition $(\omega_{ac}+\omega_{db})=0$ are retained. We will refer to this
level of approximation as non-diagonal secular approximation from now on. As
it is evident from the comparison of Figs. 3a-b, the non-diagonal secular
approximation includes non-negligible terms that couples coherence and
population elements of $\rho^{\mathrm{s}}_{ab}$. Finally, Fig. 3c shows that
the diagonal secular approximation becomes fulfilled when Kramers degeneracy
is removed by applying a small external field of 0.01 Tesla along the
molecule’s easy axis. We note that the eigenvalues of the matrices reported in
Figs. 3a,c are virtually identical, and at the same time different from those
of the matrix reported in Fig. 3b. Moreover, we find that the spectrum of $R$
in the diagonal secular approximation is not rotationally invariant even in
zero-field, a clear indication of inconsistency.

Figure 4: Spin-phonon relaxation time for (2). Experimental relaxation times
measured at 1000 Oe with AC magnetometry are reported as black continuous
linesRechkemmer (2016). Simulated Orbach rates computed with the non-diagonal
(continuous line and triangles) and diagonal secular approximation (dashed
line and circles) are reported in red. Simulated Raman rates computed with the
diagonal secular approximation with external field $B_{z}=0.01$ Tesla
(continuous line and triangles) and $B_{z}=0.0$ Tesla (dashed line and
circles) are reported in green. Raman rates in non-zero external field are
computed by rotating the molecular geometry in the frame of the g-tensor’s
eigenvectors.

In order to determine the importance of including coherence-transfer terms
into the description of spin-phonon relaxation, we compute the value of spin
lifetime $\tau$ for (2) by taking the inverse of the first non-zero eigenvalue
of the matrix $R2^{1-\mathrm{ph}}_{ab,cd}$, as commonly done for Markov
processes. The latter is equivalent to studying the decay of $M_{z}(t)$ for a
molecule with the magnetization easy-axis parallel to $z$. Fig. 4 shows the
results for the two different levels of secular approximation, together with
experimental results obtained in the presence of a small external field. The
application of a small external field in experiments helps removing the effect
of dipolar relaxationRechkemmer (2016), not included in the simulations and
therefore makes the comparison between theory and experiments more reliable.
The relaxation time predicted with the diagonal secular approximation is the
same reported beforeLunghi and Sanvito (2020b) and despite the qualitative
agreement with experiments, it shows between one and two orders of magnitude
of underestimation of spin lifetime. We note that the prediction of Raman
contributions with Eq. 16 only slightly improves with respect to those
presented previouslyLunghi and Sanvito (2020b). Once the non-diagonal secular
approximation is introduced for one-phonon processes with the use of Eq. 6,
the agreement between experimental and simulated $\tau$ becomes virtually
exact in the high-$T$ relaxation regime, where the Orbach mechanism drives
relaxation. For temperature below 20 K, the two-phonon Raman relaxation
described by Eq. 16 takes over the Orbach one. Although a non-diagonal secular
approximation for this process is not available, we attempted to remove the
effect of coherence transfer by applying a small external field along the
molecular easy-axis of magnetization. As suggested by the study of
$R2^{1-\mathrm{ph}}_{ab,cd}$ matrix elements, the results for the diagonal and
non-diagonal secular approximations becomes identical under these conditions.
Applying this strategy to $R4_{aa,bb}^{2-\mathrm{ph}}$, results improve
significantly with only a negligible residual deviation left between
experiments and simulations. Finally, we note that simulations for (2) have
been found extremely robust with respect to the choice of basis set, scalar
relativistic effects, and dynamical correlation. Moreover, predictions
obtained with $\Gamma-$point phonons computed from MLFF or DFT agree very well
among them, thus validating the use of MLFFLunghi and Sanvito (2019b) to
integrate the phonons’ Brillouin zone and converge the values of $\tau$ below
10 K, where border-zone and acoustic phonons start contributing. As discussed
elsewhereLunghi and Sanvito (2020b) and in agreement with what is observed
here for (1), low energy phonons up to $\sim 50$ cm-1 are the main drive for
Raman spin relaxation in (2). Also in this case, it was found that these low
energy vibrations correspond to molecular rotations combined to small and
delocalized intra-molecular distortionsLunghi and Sanvito (2020b).

Spin-phonon relaxation in J=15/2 lanthanide systems.

In order to further test the ability of our models, we extended our study to
(3). This molecule exhibits a Dy3+ ion in almost perfectly axial crystal
field. This coordination stabilizes a $J=15/2$ ground state and imposes a
strong zero-field splitting that separates the ground- and highest excited KDs
($Mj=\pm 15/2$ and $Mj=\pm 1/2$, respectively) by $\sim 1500$ cm-1Goodwin  _et
al._ (2017). The spin states of (3) can be described with a generalized spin
HamiltonianUngur and Chibotaru (2017); Jung  _et al._ (2019)

$\hat{H}_{s}=\sum_{l=2,4,6}\sum_{m=-l}^{l}=B^{l}_{m}\hat{O}_{m}^{l}(\vec{\mathbf{J}})\>,$
(19)

where $\hat{O}_{m}^{l}(\vec{\mathbf{J}})$ are tesseral tensor operatorsTennant
_et al._ (2000). (3) crystallizes in a triclinic unit-cell with two molecular
units as well as two [B(C6F5)4]-Goodwin  _et al._ (2017) molecules as counter-
ions. Phonons and spin-phonon coupling coefficients were computed for (3)
following the strategy validated on (2) and that involves the computation of
$\Gamma-$point phonons with pDFT and linear spin-phonon coupling coefficients
with CASSCF. Results of spin relaxation for the one-phonon Orbach process in
(3) are reported in Fig. 5 together with experimental values. Once again, the
agreement between the simulated and experimental Orbach relaxation rates is
excellent, but only after the non-diagonal secular approximation is
introduced. The prediction of Orbach rates obtained using the diagonal secular
approximation is off by many orders of magnitude, depending on the value of
$T$. Such a large disagreement is mostly due to the fact that in absence of
the non-diagonal secular approximation, relaxation is predicted to be promoted
by the first excited KD, computed at $\sim$450 cm-1, instead of being driven
by absorption of phonons in resonance with the KD with energy $\sim 1300$
cm-1. We tested the effect of removing Kramers degeneracy by aligning the
Cartesian frame along the molecular easy-axis, as expressed by the
eigenvectors of ground-state KD’s g-tensor, and by applying an external field
of 0.01 Tesla along the same direction. In agreement with what observed for
(2), under these conditions, Orbach relaxation rates predicted with the
diagonal and the non-diagonal secular approximations are coincident among them
and with experimental relaxation rates. Exporting our computational strategy
to the case of Raman relaxation, we observe a massive improvement of results,
with a negligible factor of residual deviation between experiments and
simulations. Differently from (2), we note a slight dependence of results with
respect to basis set, the number of CASSCF’s solutions considered and the
inclusion of dynamical correlation, hinting to the possibility to further
improve results upon a systematic exploration of these technicalities.

Figure 5: Spin-phonon relaxation time for (3). Experimental relaxation times
measured in zero field with AC magnetometry and magnetization decay are
reported as a black continuous lineGoodwin  _et al._ (2017). Simulated Orbach
rates computed with the non-diagonal (continuous line and triangles) and
diagonal secular approximation (dashed line and circles) are reported in red.
Simulated Raman rates computed with the diagonal secular approximation with
external field $B_{z}=0.01$ Tesla (continuous line and triangles) and
$B_{z}=0.0$ Tesla (dashed line and circles) are reported in green. Raman rates
in the presence of non-zero external field are computed by rotating the
molecular geometry in the frame of the g-tensor’s eigenvectors of the ground-
state KD.

Turning to the analysis of phonon contributions to Raman spin relaxation in
(3), relaxation was found to receive contributions from several optical
phonons with energy up to 80 cm-1. A visual inspection of these modes shows
that similarly to (2) they are dominated by local molecular rotations
overlapped to small intra-molecular distortions. The latter mainly involve the
tBu groups and the position of Dy3+ ion within the two Cp- rings. We
determined the contribution of different excited KDs to Eq. 16 by including
only one of them at the time in the simulation of $\tau$ at 40 K. This
analysis shows that Raman relaxation is almost entirely determined by the
contribution of the first excited KD.

## Discussion

Although the theory of relaxation for $d$/$f$-block ions has been thoroughly
investigated for decades, the many parameters populating its equations have
prevented its use as a predictive tool. However, once relaxation theory is
combined with ab initio methods, it becomes possible to estimate spin
relaxation time without any input from experiments. Such a method stands out
as a key tool to benchmark our understanding of spin-phonon interaction and
relaxation in an unbiased way and to provide new strategies towards
optimization of relaxation times. Nonetheless, a computational method of such
complexity possesses many assumptions and approximations that require careful
validation.

In this work we have shown that a careful comparison with experimental data, a
consistent use of all fourth-order processes, and the inclusion of coherence
terms in the Redfield equations are key for obtaining quantitative
predictions. Most importantly we have shown how seemingly small details, such
as neglecting the coherence terms of the density matrix, can lead to many
orders of magnitude of deviation between predictions and experiments. As long
as a parametric approach to open-quantum-system dynamics is used, it is
impossible to detect such pitfalls, as errors are removed by rescaling
factors. These results thus provide a unique proof of concept that open
quantum system theory can be applied in a parameter-free fashion for complex
systems such as solid-state.

We envisage that further improvement of the agreement between experiments and
simulations could be achieved by carefully benchmarking both electronic
structure methods and the remaining theoretical approximations to the open-
quantum system dynamical equations. The former includes a careful study of how
vdW corrections to DFT affect phonons’ simulations and how much electronic
structure correlation and the choice of basis-set impact the prediction of
spin Hamiltonian and spin-phonon coupling coefficients. From a theoretical
point of view, the derivation of an equation able to predict the evolution of
the entire spin density matrix up to the fourth-order stands out as the next
fundamental challenge. Assessing the validity of the Markov approximation in
situations of phonon-bottleneckTesi  _et al._ (2016) and the inclusion of
spin-spin interactions into the formalism are two other interesting directions
for the field.

The unprecedented level of accuracy of simulations for such different
compounds validates the theoretical framework and the ab initio methods at the
same time, making it possible to confirm recent theoretical findings in terms
of spin relaxation mechanism in vdW crystals of molecular Kramers systemLunghi
and Sanvito (2019a); Chiesa  _et al._ (2020); Lunghi and Sanvito (2020a); Gu
and Wu (2020); Lunghi and Sanvito (2020b); Briganti  _et al._ (2021). Direct
and Orbach relaxation due to resonant phonons are correctly accounted for by
second-order contributions to density matrix dynamics and linear spin-
coupling. Two-phonon spin relaxation due to second-order contributions to
density matrix dynamics and quadratic spin-phonon coupling is the main source
of Raman relaxation for $S=1/2$, while fourth-order perturbation theory and
linear coupling is necessary to explain Raman relaxation in $S>1/2$ with large
uni-axial anisotropy. Interestingly, our results show that only THz vibrations
contribute to Raman relaxation, regardless of the energy of excited KDs.
Simulations thus support the use of the expression

$\tau^{-1}=\sum_{i}\frac{A^{1-\mathrm{ph}}_{i}}{(e^{\beta\omega_{i}}-1)}+\sum_{i}A^{2-\mathrm{ph}}_{i}\frac{e^{\beta\omega_{i}}}{(e^{\beta\omega_{i}}-1)^{2}}$
(20)

to fit one- and two-phonon contributions to experimental spin lifetime of vdW
crystals of magnetic moleculesLunghi and Sanvito (2020a). The first term of
Eq. 20 can either lead to Direct or Orbach relaxation depending on the energy
of resonant phonons. This approach has been used already on a few occasions,
providing additional insights with respect to a fitting done with simple
power-lawsSantanni  _et al._ (2020); de Camargo  _et al._ (2021); Kazmierczak
_et al._ (2021); Pfleger  _et al._ (2021).

In conclusion, we have studied spin-phonon relaxation in three molecular
compounds representing the most relevant classes of slow-relaxing Kramers
systems and shown that open-quantum system theory can be implemented in a
fully-ab initio fashion and used to accurately predict spin relaxation rates.
Now that ab initio spin dynamics is rapidly maturing into a quantitative
predictive tool it will be possible to deploy it to study other systems such
as poly-nuclear ion clustersFeltham and Brooker (2014), organic radicalsEaton
and Eaton (2002), nuclear spinsCorzilius  _et al._ (2014), paramagnetic
defects in semiconductorsJarmola  _et al._ (2012) and magnetic impurities in
solid-state hostsWolfowicz  _et al._ (2021) or adsorbed on surfacesDonati  _et
al._ (2016).

## Computational Methods

Electronic structure simulations. The unit-cell X-ray structure of (1)Atzori
_et al._ (2016) and (3)Goodwin  _et al._ (2017) were used as starting points
for a periodic DFT optimization with the software CP2KKühne  _et al._ (2020).
Density functional theory (DFT) with the PBE functionalPerdew  _et al._
(1996a), including Grimme’s D3 van der Waals corrections Grimme  _et al._
(2010), was used together with a double-zeta polarised (DZVP) MOLOPT basis
set. A plane-wave cutoff of 2500 Ry, and 1500 Ry was used for (1) and (3),
respectively. All the unit-cell lattice force constants were computed with a
two-point numerical differentiation with step of 0.01 Å. Geometry
optimization, phonons and generation of a MLFF for (2) were presented before.
The ORCA software Neese  _et al._ (2020) has been employed for the calculation
of all spin Hamiltonian terms. Simulations with ORCA were carried out on
molecular geometries optimized with pDFT. The tensors g and A for (1) were
computed with DFT employing the PBE0 functionalPerdew  _et al._ (1996b). The
spin Hamiltonian of (3) was computed with CASSCF employing a (9,7) active
space and by using all the solutions with multiplicity six, 128 solutions with
multiplicity four, and 130 solutions with multiplicity two. Spin-orbit
contributions were included through quasi-degenerate perturbation theory. ZORA
and DKH scalar relativistic corrections were used for (1) and (3),
respectively. The calculation of spin-phonon coupling coefficients for (1) was
also repeated with DKH scalar relativistic corrections. The RIJCOSX
approximation with GridX6 integration grid was used for both (1) and (3). The
basis sets DKH/ZORA-def2-TZVPP was used for V, S, C and H. SARC-DKH-TZVPP was
used for Dy instead. Spin-phonon coupling coefficients for (3) were computed
using a slightly cheaper setup: state-average CASSCF simulations only included
the roots with multiplicity six, and RIJCOSX was used with the GridX4 option.
The basis sets DKH-SARC-def2-QZVP, DKH-def2-TZVPP and DKH-def2-SVP were used
for Dy, C, and H, respectively. Spin Hamiltonian and spin-phonon coupling
coefficients for (2) were computed previously with CASSCFLunghi and Sanvito
(2020b).

Spin-phonon relaxation simulations. Phonons’ frequency and normal modes of
vibrations were computed by diagonalization of the dynamical matrix

$D_{ij}(\mathbf{q})=\sum_{l}\Phi^{0,l}_{ij}e^{i\mathbf{q}\cdot\mathbf{R}_{l}}\>,$
(21)

where $\mathbf{q}$ is a Brillouin zone vector and $\Phi^{0,l}_{ij}$ is the
mass-weighed lattice force constant between the degree of freedom $i$ in the
reference unit-cell $l=0$ and the degree of freedom $j$ in the unit-cell $l$.
First-order spin-phonon coupling coefficients are computed with the expression

$\Big{(}\frac{\partial\hat{H}_{\mathrm{s}}}{\partial
q_{\alpha\mathbf{q}}}\Big{)}=\sum_{i}^{3N}\sqrt{\frac{\hbar}{N_{q}\omega_{\alpha\mathbf{q}}m_{i}}}L^{\mathbf{q}}_{\alpha
i}\Big{(}\frac{\partial\hat{H}_{\mathrm{s}}}{\partial X_{i}}\Big{)}\>,$ (22)

where $L^{\mathbf{q}}_{\alpha i}$ and $\omega_{\alpha\mathbf{q}}^{2}$ are the
eigenvectors and eigenvalues of $D_{ij}(\mathbf{q})$, respectively, $N_{q}$ is
the total number of q-points used to integrate the Brillouin zone and the sum
over $i$ is extended to the $3N$ molecular degrees of freedom. A similar
expression is used for second-order spin-phonon couplingLunghi and Sanvito
(2020a). The derivatives of the spin Hamiltonian terms appearing in Eq. 22 are
computed by numerical differentiation. In the case of $(\partial
B^{l}_{m}/\partial X_{i})$, the parameters $B^{l}_{m}$ were computed for six
distortions between $\pm$ 0.05 Å for each molecular degree of freedom. Their
value was then interpolated with a third-order polynomial function, where the
linear coefficient corresponds to $(\partial B^{l}_{m}/\partial X_{i})$. For
the second-order differentiation of g and A, we used a four-point finite-
difference expression with step 0.01 Å. Dirac’s delta functions appearing in
the spin dynamics equations are smeared out with a Gaussian function with
$\sigma=10$ cm-1. A value of $i\sigma$ is also added to the denominators of
Eq. 16 to prevent possible divergence of calculated rates. Results are found
to be well converged with respect to $\sigma$ according to a strategy
previously testedLunghi and Sanvito (2020a, b). Low-temperature results for
Raman relaxation in (2) were also converged by sampling the vibrational
Brillouin zone thanks to supercell phonons calculationsLunghi and Sanvito
(2020b).

Acknowledgements
This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement No. [948493]). Computational resources were provided by the
Trinity College Research IT and the Irish Centre for High-End Computing
(ICHEC).

Conflict of interests
The authors declare that they have no competing interests.

## References

  * Wu  _et al._ (2010) MW Wu, JH Jiang, and MQ Weng, “Spin dynamics in semiconductors,” Physics Reports 493, 61–236 (2010).
  * Maehrlein  _et al._ (2018) Sebastian F Maehrlein, Ilie Radu, Pablo Maldonado, Alexander Paarmann, Michael Gensch, Alexandra M Kalashnikova, Roman V Pisarev, Martin Wolf, Peter M Oppeneer, Joseph Barker, _et al._ , “Dissecting spin-phonon equilibration in ferrimagnetic insulators by ultrafast lattice excitation,” Science advances 4, eaar5164 (2018).
  * Žutić  _et al._ (2004) Igor Žutić, Jaroslav Fabian, and S Das Sarma, “Spintronics: Fundamentals and applications,” Reviews of modern physics 76, 323 (2004).
  * Wolfowicz  _et al._ (2021) Gary Wolfowicz, F Joseph Heremans, Christopher P Anderson, Shun Kanai, Hosung Seo, Adam Gali, Giulia Galli, and David D Awschalom, “Quantum guidelines for solid-state spin defects,” Nature Reviews Materials 6, 906–925 (2021).
  * Atzori and Sessoli (2019) Matteo Atzori and Roberta Sessoli, “The second quantum revolution: role and challenges of molecular chemistry,” Journal of the American Chemical Society 141, 11339–11352 (2019).
  * Eaton and Eaton (2002) Sandra S Eaton and Gareth R Eaton, “Relaxation times of organic radicals and transition metal ions,” in  _Distance measurements in biological systems by EPR_ (Springer, 2002) pp. 29–154.
  * Breuer  _et al._ (2001) Heinz-Peter Breuer, Bernd Kappler, and Francesco Petruccione, “The time-convolutionless projection operator technique in the quantum theory of dissipation and decoherence,” Annals of Physics 291, 36–70 (2001).
  * Timm (2011) Carsten Timm, “Time-convolutionless master equation for quantum dots: Perturbative expansion to arbitrary order,” Physical Review B 83, 115416 (2011).
  * Koller  _et al._ (2010) Sonja Koller, M Grifoni, M Leijnse, and MR Wegewijs, “Density-operator approaches to transport through interacting quantum dots: Simplifications in fourth-order perturbation theory,” Physical Review B 82, 235307 (2010).
  * Waller (1932) I Waller, “Über die magnetisierung von paramagnetischen kristallen in wechselfeldern,” Zeitschrift für Physik 79, 370–388 (1932).
  * Van Vleck (1940) JH Van Vleck, “Paramagnetic relaxation times for titanium and chrome alum,” Physical Review 57, 426 (1940).
  * Redfield (1957) Alfred G Redfield, “On the theory of relaxation processes,” IBM Journal of Research and Development 1, 19–31 (1957).
  * Orbach (1961) R Orbach, “Spin-lattice relaxation in rare-earth salts,” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 264, 458–484 (1961).
  * Escalera-Moreno  _et al._ (2017) Luis Escalera-Moreno, Nicolas Suaud, Alejandro Gaita-Arino, and Eugenio Coronado, “Determining key local vibrations in the relaxation of molecular spin qubits and single-molecule magnets,” The journal of physical chemistry letters 8, 1695–1700 (2017).
  * Lunghi  _et al._ (2017a) Alessandro Lunghi, Federico Totti, Roberta Sessoli, and Stefano Sanvito, “The role of anharmonic phonons in under-barrier spin relaxation of single molecule magnets,” Nature communications 8, 14620 (2017a).
  * Goodwin  _et al._ (2017) Conrad AP Goodwin, Fabrizio Ortu, Daniel Reta, Nicholas F Chilton, and David P Mills, “Molecular magnetic hysteresis at 60 kelvin in dysprosocenium,” Nature 548, 439–442 (2017).
  * Lunghi and Sanvito (2019a) Alessandro Lunghi and Stefano Sanvito, “How do phonons relax molecular spins?” Science advances 5, eaax7163 (2019a).
  * Park  _et al._ (2020) Jinsoo Park, Jin-Jian Zhou, and Marco Bernardi, “Spin-phonon relaxation times in centrosymmetric materials from first principles,” Physical Review B 101, 045202 (2020).
  * Xu  _et al._ (2020) Junqing Xu, Adela Habib, Sushant Kumar, Feng Wu, Ravishankar Sundararaman, and Yuan Ping, “Spin-phonon relaxation from a universal ab initio density-matrix approach,” Nature communications 11, 2780 (2020).
  * Albino  _et al._ (2019) Andrea Albino, Stefano Benci, Lorenzo Tesi, Matteo Atzori, Renato Torre, Stefano Sanvito, Roberta Sessoli, and Alessandro Lunghi, “First-principles investigation of spin–phonon coupling in vanadium-based molecular spin quantum bits,” Inorganic chemistry 58, 10260–10268 (2019).
  * Mirzoyan and Hadt (2020) Ruben Mirzoyan and Ryan G Hadt, “The dynamic ligand field of a molecular qubit: decoherence through spin–phonon coupling,” Physical Chemistry Chemical Physics 22, 11249–11265 (2020).
  * Lunghi and Sanvito (2020a) Alessandro Lunghi and Stefano Sanvito, “The limit of spin lifetime in solid-state electronic spins,” The Journal of Physical Chemistry Letters 11, 6273–6278 (2020a).
  * Kazmierczak  _et al._ (2021) Nathanael P. Kazmierczak, Ruben Mirzoyan, and Ryan G. Hadt, “The impact of ligand field symmetry on molecular qubit coherence,” Journal of the American Chemical Society 143, 17305–17315 (2021).
  * Lunghi  _et al._ (2017b) Alessandro Lunghi, Federico Totti, Stefano Sanvito, and Roberta Sessoli, “Intra-molecular origin of the spin-phonon coupling in slow-relaxing molecular magnets,” Chemical science 8, 6051–6059 (2017b).
  * Moseley  _et al._ (2018) Duncan H Moseley, Shelby E Stavretis, Komalavalli Thirunavukkuarasu, Mykhaylo Ozerov, Yongqiang Cheng, Luke L Daemen, Jonathan Ludwig, Zhengguang Lu, Dmitry Smirnov, Craig M Brown, _et al._ , “Spin-phonon couplings in transition metal complexes with slow magnetic relaxation,” Nature communications 9, 2572 (2018).
  * Ullah  _et al._ (2019) Aman Ullah, Jesús Cerdá, José J Baldoví, Sergey A Varganov, Juan Aragó, and Alejandro Gaita-Ariño, “In silico molecular engineering of dysprosocenium-based complexes to decouple spin energy levels from molecular vibrations,” The journal of physical chemistry letters 10, 7678–7683 (2019).
  * Lunghi and Sanvito (2020b) Alessandro Lunghi and Stefano Sanvito, “Multiple spin–phonon relaxation pathways in a kramer single-ion magnet,” The Journal of Chemical Physics 153, 174113 (2020b).
  * Reta  _et al._ (2021) Daniel Reta, Jon GC Kragskow, and Nicholas F Chilton, “Ab initio prediction of high-temperature magnetic relaxation rates in single-molecule magnets,” Journal of the American Chemical Society 143, 5943–5950 (2021).
  * Blockmon  _et al._ (2021) Avery L. Blockmon, Aman Ullah, Kendall D. Hughey, Yan Duan, Kenneth R. O’Neal, Mykhaylo Ozerov, José J. Baldoví, Juan Aragó, Alejandro Gaita-Ariño, Eugenio Coronado, and Janice L. Musfeldt, “Spectroscopic analysis of vibronic relaxation pathways in molecular spin qubit [Ho(W5O18)2]9–: Sparse spectra are key,” Inorganic Chemistry 60, 14096–14104 (2021).
  * Briganti  _et al._ (2021) Matteo Briganti, Fabio Santanni, Lorenzo Tesi, Federico Totti, Roberta Sessoli, and Alessandro Lunghi, “A complete ab initio view of orbach and raman spin-lattice relaxation in a dysprosium coordination compound,” Journal of the American Chemical Society 143, 13633–13645 (2021).
  * Zabala-Lekuona  _et al._ (2021) Andoni Zabala-Lekuona, José Manuel Seco, and Enrique Colacio, “Single-molecule magnets: From Mn12-ac to dysprosium metallocenes, a travel in time,” Coordination Chemistry Reviews 441, 213984 (2021).
  * Moreno-Pineda and Wernsdorfer (2021) Eufemio Moreno-Pineda and Wolfgang Wernsdorfer, “Measuring molecular magnets for quantum technologies,” Nature Reviews Physics 3, 645–659 (2021).
  * Atanasov  _et al._ (2015) Mihail Atanasov, Daniel Aravena, Elizaveta Suturina, Eckhard Bill, Dimitrios Maganas, and Frank Neese, “First principles approach to the electronic structure, magnetic anisotropy and spin relaxation in mononuclear 3d-transition metal single molecule magnets,” Coordination Chemistry Reviews 289, 177–214 (2015).
  * Ungur and Chibotaru (2017) Liviu Ungur and Liviu F Chibotaru, “Ab initio crystal field for lanthanides,” Chem.-Eur. J 23, 3708–3718 (2017).
  * Neese  _et al._ (2019) Frank Neese, Mihail Atanasov, Giovanni Bistoni, Dimitrios Maganas, and Shengfa Ye, “Chemistry and quantum mechanics in 2019: give us insight and numbers,” Journal of the American Chemical Society 141, 2814–2824 (2019).
  * Atzori  _et al._ (2016) Matteo Atzori, Elena Morra, Lorenzo Tesi, Andrea Albino, Mario Chiesa, Lorenzo Sorace, and Roberta Sessoli, “Quantum coherence times enhancement in vanadium (IV)-based potential molecular qubits: the key role of the vanadyl moiety,” Journal of the American Chemical Society 138, 11234–11244 (2016).
  * Rechkemmer  _et al._ (2016) Yvonne Rechkemmer, Frauke D Breitgoff, Margarethe Van Der Meer, Mihail Atanasov, Michael Hakl, Milan Orlita, Petr Neugebauer, Frank Neese, Biprajit Sarkar, and Joris Van Slageren, “A four-coordinate cobalt (II) single-ion magnet with coercivity and a very high energy barrier,” Nature communications 7, 10467 (2016).
  * Rechkemmer (2016) Yvonne Rechkemmer,  _Spectroscopic investigations of the magnetic anisotropy of lanthanide- and cobalt-based molecular nanomagnets_ , Ph.D. thesis, University of Stuttgart (2016).
  * Tesi  _et al._ (2016) Lorenzo Tesi, Alessandro Lunghi, Matteo Atzori, Eva Lucaccini, Lorenzo Sorace, Federico Totti, and Roberta Sessoli, “Giant spin–phonon bottleneck effects in evaporable vanadyl-based molecules with long spin coherence,” Dalton Transactions 45, 16635–16643 (2016).
  * Eastham  _et al._ (2016) PR Eastham, P Kirton, HM Cammack, BW Lovett, and J Keeling, “Bath-induced coherence and the secular approximation,” Physical Review A 94, 012110 (2016).
  * Gösweiner  _et al._ (2020) Christian Gösweiner, Per-Olof Westlund, and Hermann Scharfetter, “Spin–spin relaxation of nuclear quadrupole resonance coherences and the important role of degenerate energy levels,” Molecular Physics 118, e1743888 (2020).
  * Garlatti  _et al._ (2020) E Garlatti, L Tesi, A Lunghi, M Atzori, DJ Voneshen, P Santini, S Sanvito, T Guidi, R Sessoli, and S Carretta, “Unveiling phonons in a molecular qubit with four-dimensional inelastic neutron scattering and density functional theory,” Nature communications 11, 1751 (2020).
  * Lunghi and Sanvito (2019b) Alessandro Lunghi and Stefano Sanvito, “A unified picture of the covalent bond within quantum-accurate force fields: From organic molecules to metallic complexes’ reactivity,” Science advances 5, eaaw2210 (2019b).
  * Jung  _et al._ (2019) Julie Jung, M Ashraful Islam, Vincent L Pecoraro, Talal Mallah, Claude Berthon, and Hélène Bolvin, “Derivation of lanthanide series crystal field parameters from first principles,” Chemistry–A European Journal 25, 15112–15122 (2019).
  * Tennant  _et al._ (2000) WC Tennant, CJ Walsby, RFC Claridge, and DG McGavin, “Rotation matrix elements and further decomposition functions of two-vector tesseral spherical tensor operators; their uses in electron paramagnetic resonance spectroscopy,” Journal of Physics: Condensed Matter 12, 9481 (2000).
  * Chiesa  _et al._ (2020) A Chiesa, F Cugini, R Hussain, E Macaluso, G Allodi, E Garlatti, M Giansiracusa, CAP Goodwin, F Ortu, D Reta, _et al._ , “Understanding magnetic relaxation in single-ion magnets with high blocking temperature,” Physical Review B 101, 174402 (2020).
  * Gu and Wu (2020) Lei Gu and Ruqian Wu, “Origins of slow magnetic relaxation in single-molecule magnets,” Physical Review Letters 125, 117203 (2020).
  * Santanni  _et al._ (2020) Fabio Santanni, Andrea Albino, Matteo Atzori, Davide Ranieri, Enrico Salvadori, Mario Chiesa, Alessandro Lunghi, Andrea Bencini, Lorenzo Sorace, Federico Totti, _et al._ , “Probing vibrational symmetry effects and nuclear spin economy principles in molecular spin qubits,” Inorganic chemistry 60, 140–151 (2020).
  * de Camargo  _et al._ (2021) Luana C de Camargo, Matteo Briganti, Francielli S Santana, Danilo Stinghen, Ronny R Ribeiro, Giovana G Nunes, Jaísa F Soares, Enrico Salvadori, Mario Chiesa, Stefano Benci, _et al._ , “Exploring the organometallic route to molecular spin qubits: The [CpTi (cot)] case,” Angewandte Chemie 133, 2620–2625 (2021).
  * Pfleger  _et al._ (2021) Rouven F Pfleger, Sören Schlittenhardt, Marcel P Merkel, Mario Ruben, Karin Fink, Christopher E Anson, Jesper Bendix, and Annie K Powell, “Terminal ligand and packing effects on slow relaxation in an isostructural set of [Dy (H2dapp) X2]+ single molecule magnets,” Chemistry (Weinheim an der Bergstrasse, Germany) 27, 15086 (2021).
  * Feltham and Brooker (2014) Humphrey LC Feltham and Sally Brooker, “Review of purely 4f and mixed-metal nd-4f single-molecule magnets containing only one lanthanide ion,” Coordination Chemistry Reviews 276, 1–33 (2014).
  * Corzilius  _et al._ (2014) Björn Corzilius, Vladimir K Michaelis, Susanne A Penzel, Enrico Ravera, Albert A Smith, Claudio Luchinat, and Robert G Griffin, “Dynamic nuclear polarization of 1H, 13C, and 59Co in a tris (ethylenediamine) cobalt(III) crystalline lattice doped with Cr(III),” Journal of the American Chemical Society 136, 11716–11727 (2014).
  * Jarmola  _et al._ (2012) A Jarmola, VM Acosta, K Jensen, S Chemerisov, and D Budker, “Temperature-and magnetic-field-dependent longitudinal spin relaxation in nitrogen-vacancy ensembles in diamond,” Physical review letters 108, 197601 (2012).
  * Donati  _et al._ (2016) Fabio Donati, S Rusponi, Sebastian Stepanow, C Wäckerlin, Aparajita Singha, Luca Persichetti, R Baltic, K Diller, Francois Patthey, E Fernandes, _et al._ , “Magnetic remanence in single atoms,” Science 352, 318–321 (2016).
  * Kühne  _et al._ (2020) Thomas D Kühne, Marcella Iannuzzi, Mauro Del Ben, Vladimir V Rybkin, Patrick Seewald, Frederick Stein, Teodoro Laino, Rustam Z Khaliullin, Ole Schütt, Florian Schiffmann, _et al._ , “Cp2k: An electronic structure and molecular dynamics software package-quickstep: Efficient and accurate electronic structure calculations,” The Journal of Chemical Physics 152, 194103 (2020).
  * Perdew  _et al._ (1996a) John P Perdew, Kieron Burke, and Matthias Ernzerhof, “Generalized gradient approximation made simple,” Physical review letters 77, 3865 (1996a).
  * Grimme  _et al._ (2010) Stefan Grimme, Jens Antony, Stephan Ehrlich, and Helge Krieg, “A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu,” The Journal of chemical physics 132, 154104 (2010).
  * Neese  _et al._ (2020) Frank Neese, Frank Wennmohs, Ute Becker, and Christoph Riplinger, “The ORCA quantum chemistry program package,” The Journal of Chemical Physics 152, 224108 (2020).
  * Perdew  _et al._ (1996b) John P Perdew, Matthias Ernzerhof, and Kieron Burke, “Rationale for mixing exact exchange with density functional approximations,” The Journal of chemical physics 105, 9982–9985 (1996b).