in liquid solutionlunghi2020insights . A schematic representation of the
former method, based on Ridge regression and bispectrum components as
molecular fingerprints is reported in Fig. 5. In a nutshell, a supercell
approach to the calculation of phonons is used in combination with DFT. The
optimized molecular structure is then distorted in the range of 1000 times by
applying small random perturbations of the atomic positions. The spin
Hamiltonian coefficients of each distorted structure are determined with
electronic structure methods and form the basis for the generation of a ML
model. A part of this set is usually not used for training the ML model but
just for testing purposes. The ML model used in Ref. lunghi2020limit
predicted the spin Hamiltonian coefficients by means of the knowledge of each
atom’s coordination environment, as graphically sketched in Fig. 5B. After a
successful training of the ML, it is possible to very quickly scan the spin
Hamiltonian coefficients along any pair of Cartesian
directionslunghi2020surfing . Fig 5D reports just an example out the many 2D
plots needed to computed all the derivatives needed in Eq. 178.

### Spin-Phonon Coupling Translational Invariance

Similarly to what noted for the calculation of force constants, also spin-
phonon coupling coefficients need to fulfill certain sum rules that arise from
the symmetries of the problem. For instance, let us assume that the
$\mathbf{D}$ tensor of a molecule depends linearly with respect to atomic
positions

$D_{\alpha\beta}(X_{is})=D_{\alpha\beta}+\sum_{is}\left(\frac{\partial
D_{\alpha\beta}}{\partial X_{is}}\right)(X_{is}-X_{is}^{0}).$ (179)

If the entire system is rigidly translated in space by an amount $a$ Å, then
$(X_{is}-X_{is}^{0})=a$, leading to

$D_{\alpha\beta}(X_{is})=D_{\alpha\beta}+a\sum_{is}\left(\frac{\partial
D_{\alpha\beta}}{\partial X_{is}}\right).$ (180)

Finally, since the value of $\mathbf{D}$ is independent on the particular
position of the molecule in space,
$D_{\alpha\beta}(X_{is}^{0}+a)=D_{\alpha\beta}(X_{is}^{0})$. Therefore we
obtain the condition

$\sum_{i}\left(\frac{\partial D_{\alpha\beta}}{\partial
X_{is}}\right)=0\>,\quad\text{for each }s\>.$ (181)

Similarly to the method used to regularize the force constants in phonons
calculations, translational invariance conditions can be enforced on the
calculated derivatives of the spin Hamiltonian tensors by rescaling the values
by a mean deviation of this condition. For instance, for the elements
$\alpha\beta$ of the tensor D, the correction for the $s$ Cartesian direction
reads

$\Big{(}\frac{\partial D_{\alpha\beta}}{\partial
X_{is}}\Big{)}=\Big{(}\frac{\partial D_{\alpha\beta}}{\partial
X_{is}}\Big{)}-Dev_{s}\>,\quad
Dev_{s}=\frac{1}{N_{at}}\sum_{i}^{N_{at}}\Big{(}\frac{\partial
D_{\alpha\beta}}{\partial X_{is}}\Big{)}\>,$ (182)

where $N_{at}$ is the number of atoms in the molecule. An equivalent
expression can be used for the components of the second-order derivatives that
shares the same Cartesian degrees of freedomlunghi2020limit .

### Spin-Phonon Relaxation Time

The calculation of the Redfield matrices is now the last step required to
compute the spin-phonon relaxation time.

The calculation of the elements of $R_{ab,cd}$ is relatively straightforward
once phonon frequencies and spin-phonon coupling coefficients are know, except
for two important technical details: the integration of the Brillouion zone
and dealing with the presence of a Dirac deltalunghi2019phonons ;
lunghi2020multiple . These two aspects are strongly connected and need to be
discussed simultaneously. Let us start by noting that the Dirac delta in not
strictly speaking a function, but a distribution in the functional analysis
sense, and that it only makes sense when appearing under the integral sign.
This does not lead to any issue as all the expressions of $R_{ab,cd}$ contain
a sum over $\mathbf{q}$ points, which for a solid-state periodic system extend
to an infinite summation, namely an integral. We thus have

$R_{ab,cd}\propto\sum_{\mathbf{q}}\delta(\omega-\omega_{\mathbf{q}})\rightarrow\int_{\mathrm{BZ}}\delta(\omega-\omega_{\mathbf{q}})d\mathbf{q}\>.$
(183)

In fact, differently from the sparse spectra of an isolated cluster of atoms,
the vibrational density of states $P(\omega)$ of a solid-state system is a
continuous function of $\omega$ in virtue of an infinite number of degrees of
system. Let us now discuss how to numerically implement an object such as the
one in Eq. 183. The integration with respect to the q-points can be simply
dealt with by discretizing the first Brillouin zone in a finite mesh of
uniformly sampled q-points. The value of $R_{ab,cd}$ is then computed for each
q-point, which are then summed up. By increasing the finesse of this mesh one
eventually reaches convergence to the desired accuracy and the Redfield matrix
will be numerically indistinguishable from the one obtained by integration.
However, as already discussed, the Dirac delta loses sense when the sign of
integral is removed and we need a numerical way to address this issue. The
common way to deal with the Dirac delta is to replace it with a real function
that mimics its properties. A common choice is to use a Gaussian function with
width $\sigma$ in the limit of $\sigma\rightarrow 0$. One can then substitute
the Dirac delta functions with a Gaussian

$\delta(x)\sim\frac{1}{\sigma\sqrt{\pi}}e^{-x^{2}/\sigma^{2}}\>,$ (184)

and recompute the Redfield matrix for decreasing values of $\sigma$ until
convergence is reached. To summarize one needs to evaluate the Redfield matrix
in the limit of infinite q-points and vanishing Gaussian smearing

$R_{ab,cd}\propto\int_{\mathrm{BZ}}\delta(\omega-\omega_{\mathbf{q}})d\mathbf{q}=\lim_{\sigma\rightarrow
0}\lim_{q\rightarrow\infty}\sum_{\mathbf{q}}\frac{1}{\sigma\sqrt{\pi}}e^{-(\omega-\omega_{\mathbf{q}})^{2}/\sigma^{2}}\>,$
(185)

One final important care in evaluating the limits of Eq. 185 is needed: the
limit with respect to $\mathbf{q}$ must always be evaluated first than the one
with respect to $\sigma$. This is required because in a simulation we are
always dealing with a finite number of phonon states, regardless the number of
$N_{q}$, and if $\sigma$ is reduced arbitrarily the elements of $R$ will
invariably converge to zero. This is because exact degeneracy between two
discrete spectra (the spin and the phonons ones) is never achieved and without
any smearing, the Dirac delta will always invariably result in a null value.
However, we know that the discrete nature of the phonon spectrum is only due
to the use of a finite mesh of q-points and it is not intrinsic of a solid-
state system. For this reason, $\sigma$ can never be chosen to be much smaller
than the finesse of the vibrational spectrum.

To summarize, the correct numerical procedure to evaluate Eq. 185 is the
following: i) $R$ is evaluated for a large smearing, e.g. $\sigma=10$ cm-1,
and one q-point, namely the $\Gamma$ point; ii) Keeping $\sigma$ fixed, $R$ is
recalculated for increasing values of $N_{q}$ until the desired level of
convergence is reached; iii) $\sigma$ is reduced and $N_{q}$ is increased
further until new convergence is reached; iv) step iii) is repeated until
convergence with respect to both $\sigma$ and $N_{q}$ is obtained.

If this procedure is employed, the elements of the Redfield matrix converge to
the right value without leading to any divergence of the relaxation
timelunghi2019phonons ; lunghi2020multiple . Until now we have considered a
Gaussian smearing to reproduce the Dirac delta function but that is not the
only option. Another common choice is the use of a Lorentzian profile, with
line-width $\Delta$. Both Gaussian and Lorentzian functions converge to the
Dirac delta for their width tending to zero, however, the rate of convergence
of the two functions is quite different. The tails of the Gaussian function
indeed drop to zero exponentially, while those of a Lorentzian function they
do so quadratically. Therefore, in order for a Gaussian and Lorentzian
smearing to give results in agreement, the latter must be much smaller than
the former. Finally, it is important to remark that in this discussion we have
assumed that smearing is only used to represent the Dirac delta, which
correspond to a perfectly harmonic phonons bath. In section II we have also
contemplated the use of Lorentzian smering, but for a different reason. A
Lorentzian line-width is in fact the natural line-shape of the Fourier
transform of the phonon correlation function in the presence of phonon-phonon
dissipationlunghi2017role . The use of Lorentzian smearing in these two
circumstances should not be confused. In the former case $\Delta$ must tend to
zero in order to converge to a Dirac delta, while in the latter case $\Delta$
is set by the phonon bath properties and will in general have a finite value
in the order of cm-1 or fractions.

Once the Redfield matrix is computed, all is left to do is to solve the
equation of motion for the density matrix

$\frac{d\rho_{ab}(t)}{dt}=\sum_{cd}R_{ab,cd}\rho_{cd}(t)\>.$ (186)

In order to treat a mathematical object with four indexes we introduce the
Liuville space, $\mathcal{L}$ which is defined as the tensor product of the
original Hilbert space, $\mathcal{L}=\mathcal{H}\otimes\mathcal{H}$. The space
$\mathcal{L}$ is therefore spanned by vectors of size $size(\mathcal{H})^{2}$
that lists all the possible pairs of vectors of the Hilbert space. What is
expressed as a matrix in the Hilbert space, then becomes a vector in the
Liuville space. For instance, the density matrix $\rho_{ab}$ itself would read
$\rho_{i}$, where $i=ab$. Accordingly, the Redfield matrix $R_{ab,cd}$ becomes
$R_{ij}$, where $i=ab$ and $j=cd$. An operator in matrix form in the Liuville
space is called a super-operator and indicated with a double hat, e.g.
$\hat{\hat{R}}$, to highlight that it acts on normal operators to give another
operator. Finally, the equation of motion of the density matrix in the
interaction picture reads

$\frac{d\rho_{i}(t)}{dt}=\sum_{j}R_{ij}\rho_{j}(t)\>,\quad\text{or
equivalently,}\quad\frac{d\hat{\rho}}{dt}=\hat{\hat{R}}\hat{\rho}(t)$ (187)

This differential equation can now be solved with common numerical strategies.
Once the Redfield matrix has been computed and diagonalized, one can construct
the propagator as

$L_{ij}(t)=\sum_{k}V_{ik}e^{i\lambda_{k}t}V^{-1}_{kj}\>,$ (188)

where $V_{ij}$ and $\lambda_{k}$ are the eigenvectors and eigenvalues of
$R_{ij}$, respectively, and $V^{-1}_{ij}$ are the elements of the inverse
matrix of the eigenvectors of $R_{ij}$. The evolution of the density matrix in
the interaction picture therefore is

$\rho_{i}(t)=\sum_{j}L_{ij}(t)\rho_{j}(t=0)\>.$ (189)

Importantly, given the properties of $R$, $\lambda_{k}$ should always have one
null eigenvalue and all the others of negative sign. Moreover, the eigenvector
associated with the null eigenvalue will correspond to the equilibrium
distribution of spin states’ population. These properties of the propagator
enforce that the time evolution of $\hat{\rho}$ in fact corresponds to a
relaxation process, where out-of-diagonal elements goes to zero an the
diagonal terms tend to their equilibrium value for $t\rightarrow\infty$.

Once the evolution of $\hat{\rho}$ has been determined, it is possible to
compute the time evolution of the magnetization with the expression
$\vec{\mathbf{M}}(t)=\mathrm{Tr}\left(\vec{\mathbf{S}}\hat{\rho}(t)\right)$,
where the trace is generally performed in the basis of the spin Hamiltonian
eigenvectors. The profile of $\vec{\mathbf{M}}(t)$ can be studied as commonly
done in experiments and in general different initial conditions for
$\hat{\rho}(t=0)$ and different orientations of the molecule in the external
field will lead to different demagnetization profiles. A common choice is to
orient the molecule with the easy/hard axis along the external field and
initialize $\hat{\rho}$ in such a way that the molecule is fully magnetized
along the same direction. In this scenario, for a single-spin system, an
exponential decay of $M_{z}(t)$ is generally observed. The latter can be
fitted with the common expression

$M_{z}(t)=\left[M_{z}(t=0)-M_{z}^{eq}\right]e^{-t/T_{1}}+M_{z}^{eq}$ (190)

where $M_{z}^{eq}$ is the equilibrium value of the magnetization and $T_{1}$
(or $\tau$) is the relaxation time. Under this conditions, $1/T_{1}$ coincide
with the value of the second-smallest eigenvalue of $R_{ij}$ (the first non-
zero one).

##  V The origin of spin-phonon coupling in magnetic molecules

Now that we have established a robust computational protocol to implement the
theory of spin-phonon relaxation, we are well positioned to discuss
applications of this method for different molecular systems. In particular we
are interested in understanding how spin-phonon coupling is influenced by
different molecular motions and by the chemical nature of the metal ions,
ligands, and their bond.

### Normal modes composition

Let us start our journey form a series of simple $S=1/2$ mono-nuclear
coordination compounds based on a V4+ metal center. Vanadium in this oxidation
state has a $[$Ar$]d^{1}$ electronic configuration and for the common square
pyramidal or distorted octahedral coordination geometries the unpaired
electron occupies a non-degenerate $d$ orbital. This leads to a fully quenched
angular momentum and a well defined $S=1/2$ ground state. The first electronic
excited state in V+4 compounds is often found at energies above 5000
cm-1albino2019first , thus much higher energies than any vibrational frequency
of molecular crystals. This electronic configuration makes it possible to use
Vanadium complexes to study spin-phonon relaxation in a prototypical two-level
systems.

The spin Hamiltonian for an isolated V4+ molecule reads

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

where the Lande’ tensor $\mathbf{g}$ and the hyperfine coupling tensor
$\mathbf{A}$ determine the spin spectrum of the compound. Additionally, when
the molecule is embedded in a crystal, the dynamics of the single spin is
influenced by the dipolar interactions with all the surrounding molecules.
This interaction is captured by a spin Hamiltonian of the form

$\hat{H}_{s}=\sum_{ij}\vec{\mathbf{S}}_{i}\cdot\mathbf{D}_{ij}^{Dip}\cdot\vec{\mathbf{S}}_{j}\>,$
(192)

where $\mathbf{D}^{Dip}$ depends on the molecule-molecule distance as $r^{-3}$
and the magnitude and orientation of the molecular magnetic
momentsbencini2012epr .

Figure 6: Molecular structure and spin-phonon coupling in VO(acac)2. The unit
cell of VO(acac)2 is reported in panel a. Vanadium atoms are in pink, Carbon
atoms in green, Oxygen atoms in red and Hydrogen atoms in white. The the
distribution of spin-phonon coupling intensity as function of phonon
frequencies for the g-tensor, the dipolar D-tensor and the the A-tensor are
reported in panels b, c, and d, respectively. Reprinted from ref.
lunghi2019phonons . $\copyright$ The Authors, some rights reserved; exclusive
licensee AAAS. Distributed under a CC BY-NC 4.0 License.

Besides shaping the static spin spectrum, the tensors $\mathbf{g}$,
$\mathbf{A}$ and $\mathbf{D}^{Dip}$ are also the interactions that get
modulated by molecular motion and that are at the origin of spin-phonon
coupling. In Refs. lunghi2019phonons ; lunghi2020limit ; garlatti2020unveiling
; albino2021temperature , the compound VO(acac)2 was studied following the
methods detailed in the previous section, shedding light on the nature of
molecular vibrations, and spin-phonon coupling in $S=1/2$ molecular crystals.
For this purpose a $3\times 3\times 3$ supercell of the VO(acac)2 crystal was
optimized and used to sample the crystal phonons across the entire Brillouin
zonelunghi2019phonons . The unit cell of the crystal is reported in Fig. 6A.
Moreover, linear spin-phonon coupling coefficients were computed with Eq. 176
for all the interactions in Eqs. 191 and 192. Fig 6 reports the norm of spin-
phonon coupling resolved as function of the phonon’s angular frequency and
integrated over all the Brillouin zone. This analysis reveals that the
modulation of dipolar interactions reaches a maximum around 50 cm-1 (Fig. 6c),
coincidentally very close to the energy of the first optical mode at the
$\Gamma$-point. On the other hand, the contribution of hyperfine and Lande’
tensors to spin-phonon coupling (Fig. 6b,c) maintain a finite value for all
energies. This observation can be explained by considering the nature of the
tensors $\mathbf{g}$, $\mathbf{A}$ and $\mathbf{D}^{Dip}$. The latter only
depends on the inter-molecular distance and is therefore most effectively
modulated by rigid molecular translations. This sort of molecular
displacements are expected to be associated with acoustic modes, which density
of states increases with $\omega$ but up to a finite cut-off value that
depends on the inter-molecular force constants and molecular mass. Fig 6 is
compatible with such an interpretation with the exception of the presence of a
sharp cutoff frequency. On the other hand, hyperfine and Lande’ tensors
represent interactions mediated by the electronic structure of the metal ion
and it is reasonable to expect them to be influenced by those optical modes
that modulate the metal’s coordination shell.

A full decomposition of the phonons density of states in terms of molecular
motions was also performed in Ref. lunghi2019phonons  and reported here in
Fig. 7. The total vibrational density of states is decomposed in terms of
rigid molecular translations, rotations and intra-molecular distortions,
revealing several important features of molecular crystals’ vibrational
properties. Low-energy phonons are dominated by rigid translations or
rotations up to frequencies of 100 cm-1, after which optical modes are fully
characterized by intra-molecular displacements. Not surprisingly, the dipolar
contribution to spin-phonon coupling closely mimics the translation
contribution to phonons which peaks at 50 cm-1. Although intra-molecular
contributions become dominant only at high-frequencies, it is important to
note that they are present at virtually any frequency. As shown by the inset
of Fig. 7, a tiny amount of intra-molecular contribution to the phonons is
always present due to the finite rigidity of molecular structures. Since rigid
translations cannot lead to the modulation of local molecular properties such
as $\mathbf{g}$ and $\mathbf{A}$, this intra-molecular contribution to
acoustic phonons must be at the origin of spin-phonon coupling at frequencies
much smaller than optical modes.

Figure 7: Mode composition and spin-phonon coupling in VO(acac)2. The left
panel reports the total density of states of VO(acac)2 (purple) and its
decomposition into translational molecular contributions (green), molecular
rotations (blue) and intra-molecular distortions (yellow). The right panel
reports the spin-phonon coupling intensity due to the modulation of the
g-tensor as function of phonon frequency and Brillouin zone vector
$\mathbf{q}$ along the path $X-\Gamma$. Left panel is reprinted from ref.
lunghi2019phonons . $\copyright$ The Authors, some rights reserved; exclusive
licensee AAAS. Distributed under a CC BY-NC 4.0 License. Right panel is
reprinted from garlatti2020unveiling . Copyright 2020, The Author(s).

The right panel of Fig. 7 provides an analysis of spin-phonon coupling in
VO(acac)2 from a slightly different perspective. In this case, the spin-phonon
coupling due the modulation of the $g$-tensor was projected on the vibrational
modes along a path in the reciprocal spacegarlatti2020unveiling . Fig. 7 makes
it possible to appreciate the three acoustic modes with vanishing frequencies
for $q\rightarrow\Gamma$. The same phonons also has a vanishing spin-phonon
coupling, as acoustic modes at $q=\Gamma$ correspond to a collective
translation of the entire crystal, which cannot lead to any coupling with the
spin due to the isotropic nature of space. The spin-phonon coupling norm of
acoustic modes reported in Fig. 7 increases linearly moving away from the
$\Gamma$-point up to an energy of $\sim 2$ meV (not shown explicitly), and
abruptly changes for higher energies when the admixing with optical phonons
becomes important. This is also highlighted by the presence of avoided
crossing points in the phonon dispersions. Considering that the norm used in
Fig. 7 is the square of spin-phonon coupling, the latter is then found to vary
as $\sqrt{q}$ for small values of $q$, or equivalently as $\sqrt{\omega}$ in
virtue of a pseudo-linear dispersion relation. Although this behaviour is in
agreement with the predictions of the Debye model, its origin is qualitatively
different. The Debye model assumes that the magnetic properties of an atoms
are directly influenced by the modulation of its distance with the atoms in
the adjacent unit lattice cell as prescribed by acoustic phonons. As already
discussed, the rigid translation of molecules in the crystal cannot lead to
any sizeable spin-phonon coupling due to the fact that the number of atoms per
unit cell is very large and that the magnetic properties of the metal ions are
very well screened from inter-molecular interactions. However, since the
molecules are not infinitely rigid, when molecules translate in space they
also slightly deform, leading to a finite value of spin-phonon coupling. From
the results of Ref. garlatti2020unveiling  we can deduce that the amount of
this intra-molecular contribution to acoustic modes depends linearly with $q$
(for small $q$), thus leading to results in qualitative agreement with the
Debye model at frequencies much smaller that the optical modes.

The study of linear spin-phonon coupling provides important insights on the
origin of this interaction in terms of molecular displacements and the same
principles can expected to be applicable also to higher order
couplingslunghi2020limit .

A similar analysis to the one just presented has also been performed for
transition metal coordination compounds with S$>$1/2. Ref. lunghi2017intra
reported a study on mono-nuclear compound based on a Fe2+ ion with a nominal
$S=1$ ground state, namely [(tpaPh)Fe]−
(H3tpaPh=tris((5-phenyl-1H-pyrrol-2-yl)methyl)amine). The molecular structure
is reported in Fig. 8.

Figure 8: Spin-phonon coupling in [(tpaPh)Fe]-. The top-left panel reports the
relative spin-phonon coupling intensity due to the modulation of the
$\mathbf{D}$ by the displacement of the Cartesian coordinates of each atom.
The molecular structure is reported for reference. The bottom-left panel
reports the vibrational density of states due to the motion of the first
coordination shell. The top-right panel reports the composition of the normal
modes as function of frequency and the bottom-right panel the spin-phonon
coupling intensity decomposed across translational, rotational and intra-
molecular contributions. This figure is adapted from ref lunghi2017intra .

The vibrations of this system were only studied at the level of the unit-cell,
therefore only giving access at the $\Gamma$-point vibrations. Despite this
limitation, the study showed a consistent picture of spin-phonon coupling in
molecular compounds to the one just presented. For instance, a comparison
between the left and right panels of Fig. 8, shows that local distortions of
the molecular units are responsible for the largest spin-phonon coupling and
that the spin-phonon coupling due to the motion of atoms outside the first
coordination shell rapidly decreases with the distance from the metal atom.
Indeed, the intra-molecular distortions that were found to be strongly coupled
to spin are those involving the first coordination shell.

The right panel of Fig. 8 also reports the analysis of the molecular motion
associated to each phonon and shows that optical modes are completely
characterized by intra-molecular distortions for energies above 100-200 cm-1,
where the maximum of spin-phonon coupling is achieved.

Finally, we present the case of the high-spin Co2+ complex [CoL2]2−
(H2L=1,2-bis(methanesulfonamido)benzene)rechkemmer2016four . [CoL2]2− exhibit
a tetrahedral coordination and a large uni-axial magnetic anisotropy with
$D\sim 115$ cm-1. [CoL2]2− represents the only case other than VO(acac)2 where
phonon dispersions and their role in spin-phonon coupling have been studied
with ab initio methods. Fig. 9 reports the results of Ref. lunghi2020multiple
, where the spin-phonon coupling intensity was reported as function of
frequency and $\mathbf{q}$-vector along an arbitrary direction.

Figure 9: Spin-phonon coupling in CoL2. Panel A reports the spin-phonon
coupling intensity for the high-energy window where vibrations are resonant
with the excited KD energy. Panel B reports the low energy window proper of
acoustic and the first optical modes. Panels C-F depict the molecular
vibrations with energy indicated by the white arrows in panels A and B.
Reprinted from ref. lunghi2020multiple , with the permission of AIP
Publishing.

Also in this case, optical modes are found to be significantly more coupled
than acoustic phonons, except at border-zone, where they have similar energies
to the optical ones. The left panel of Fig. 9 reports the atomic displacements
associated with some of the phonons and confirm that local distortions of the
coordination shell strongly affect the magnetic degrees of freedom. Here we
note a typical scenario, where low energy optical modes are an admixture of
rigid rotations and delocalized intra-molecular translations, while high-
energy optical modes are strongly localized on specific intra-molecular
degrees of freedom. When the latter involve the first-coordination shell of
the metal, they posses a large coupling with the spin.

Finally, it is worth mentioning that the experimental estimation of spin-
phonon coupling is also possible. Far infra-red or Raman spectroscopies
recorded in the presence of external magnetic fields have been used to extract
the spin-phonon coupling coefficients for the modes in resonance with spin
transitionsmoseley2018spin ; stavretis2019spectroscopic ; moseley2020inter ;
blockmon2021spectroscopic ; marbey2021 . Once combined with ab initio
simulations, these experimental techniques provide one of the most direct
insights on spin-phonon coupling.

### The coordination bond and ligand field contributions

The studies discussed so far highlight important features of spin-phonon
coupling in crystals of magnetic molecules but do not address the critical
point of how the chemistry of a compound influences this interaction and the
nature of the lattice vibrations.

Let us start once again from $S=1/2$ systems. In the work of Albino et
al.albino2019first , four complexes were studied: two penta-coordinated
vanadyl (VO2+) and two hexa-coordinated V4+ molecules, where the coordination
is obtained by catecholate and dithiolate ligands, namely [VO(cat)2]2-,
[V(cat)3]2-, [VO(dmit)2]2-, and [V(dmit)3]2-, where cat=catecholate and
dmit=1,3-dithiole-2-thione-4,5-dithiolate. The analysis of spin-phonon
coupling strength distribution as function of the frequency is reported in
Fig. 10 for all four compounds. The notable result is that the molecules with
cat ligands present higher frequencies of vibrations with respect to the
corresponding dmit ligands, but at the same time overall larger spin-phonon
coupling coefficients. The first observation can be traced back to the lower
mass and harder nature of the Oxygen donor atom of Cat with respect to the
Sulphur of dmit. To understand the effect ligand’s chemical nature on spin-
phonon coupling coefficients it was instead proposed a correlation with the
static $g$-shifts, which depend on the position of the electronic excited
states, in turns correlated to ligand’s donor power and covalency of the
coordination bond. These observations also advances the hypothesis that the
$g$-tensor could be used as a proxy for spin-phonon coupling, pointing to
isotropic ones associated with low spin-phonon coupling and vice versa.

Figure 10: Spin-phonon coupling in V4+ complexes. The spin-phonon coupling
density due to the g-tensor anisotropy is computed as function of vibrational
modes’ frequency for [VO(cat)2]2- (1), [V(cat)3]2- (2), [VO(dmit)2]2- (3), and
[V(dmit)3]2- (4). Reprinted with permission from ref. albino2019first .
Copyright 2019 American Chemical Society.

These results for the spin-phonon coupling associated to the modulation of the
$g$-tensor have also been formalized in the context of ligand field
theorymirzoyan2020dynamic  and multi-reference perturbation
theorylunghi2019ligand . Following the formalism of dynamical ligand field
proposed by Mirzoyan et al.mirzoyan2021deconvolving , we express the
$g$-tensor in terms of second-order perturbation theory over $d$ orbitals

$g_{ij}=g_{e}-\sum_{v}\frac{\zeta}{S\Delta_{0v}}\langle\psi_{0}|\hat{l}_{i}|\psi_{v}\rangle\langle\psi_{v}|\hat{l}_{j}|\psi_{0}\rangle\>,$
(193)

where $\hat{l}_{i}$ is a component of the electronic angular momentum of the
metal ion, $\Delta_{0v}$ is the energy separation between the ground state and
the $v$-excited state characterized by the wave functions $\psi_{0}$ and
$\psi_{v}$, respectively. Finally, $\zeta$ is the free-ion spin-orbit coupling
constant and $S$ is the value of the molecular ground-state spin. Considering
once again the case of [V(cat)3]2-, where the V4+ ion experience a distorted-
octahedral crystal field, the ground-state and the first excited-state wave
functions $\psi_{0}$ and $\psi_{1}$ can be written as

$\displaystyle|\psi_{0}\rangle=\alpha|d_{z^{2}}\rangle+\alpha^{\prime}|\phi_{0}\rangle$
(194)
$\displaystyle|\psi_{1}\rangle=\beta|d_{yz}\rangle+\beta^{\prime}|\phi_{1}\rangle$
(195)

where $\alpha$ and $\beta$ describe the contribution of V4+’s $d$-like
orbitals that participate to $g_{xx}$, while $\alpha^{\prime}$ and
$\beta^{\prime}$ represent the contribution of the ligands’ orbitals. Assuming
that only the portion of the orbitals localized on the V4+ centre can
contribute to the integrals of Eq. 193 and substituting Eq. 195 into Eq. 193,
one obtains that

$\delta
g_{xx}=g_{xx}-g_{e}=-6\zeta\frac{\alpha^{2}\beta^{2}}{\Delta_{01}}=-\zeta_{\mathrm{eff}}\Delta_{01}^{-1}\>,$
(196)

It is then possible to derive an analytical expression for spin-phonon
coupling by taking the derivative of Eq. 196 with respect to the atomic
positions $X_{is}$

$\displaystyle\Big{(}\frac{\partial g_{xx}}{\partial X_{is}}\Big{)}_{0}$
$\displaystyle=6\zeta\frac{\alpha^{2}\beta^{2}}{\Delta_{01}^{2}}\Big{(}\frac{\partial\Delta_{01}}{\partial
X_{is}}\Big{)}_{0}-6\zeta\Delta_{01}^{-1}\Big{(}\frac{\partial(\alpha^{2}\beta^{2})}{\partial
X_{is}}\Big{)}_{0}$ (197) $\displaystyle=-\frac{\delta
g_{xx}}{\Delta_{01}}\Big{(}\frac{\partial\Delta_{01}}{\partial
X_{is}}\Big{)}_{0}-6\zeta\Delta_{01}^{-1}\Big{(}\frac{\partial(\alpha^{2}\beta^{2})}{\partial
X_{is}}\Big{)}_{0}\>.$ (198)

Eq. 198 shows that there are two possible contributions to spin-phonon
coupling due to the modulation of the $g$-tensor. The first one is due to the
modulation of the first excited state’s energy and the second one is due to
the modulation of the metal-ligand hybridization due to the coordination bond.
Both DFT and CASSCF results confirm that the modulation of covalency is not
the dominant contribution to Eq. 198mirzoyan2020dynamic ; lunghi2019ligand ,
thus confirming the important role of $\Delta_{01}$. These observations also
lead to the conclusion that spin-phonon coupling in $S=1/2$ systems is
proportional to the static $g$-tensor, which can be used as a simple rule-of-
thumb to rank potential slow relaxing molecules. This hypothesis has also been
recently tested on a series of iso-structural $S=1/2$ complexes with V4+,
Nb4+, and Ta4+, where it was possible to correlate relaxation time to the
experimental $g$-shiftchakarawet2021effect .

These studies on the contributions to spin-phonon coupling also make it
possible to rationalize a body of recent literature that proposes a variety of
apparently orthogonal strategies to improve spin-lattice relaxation time in
magnetic molecules. For instance, in light of Eq. 198, several experimental
works reporting respectively $g$-tensor isotropyariciu2019engineering , large
molecular covalencyfataftah2019metal ; briganti2019covalency  and molecular
rigiditybader2014room ; atzori2016room  as primary ingredients for long spin-
lattice relaxation times in molecular qubits, can be reconciled by recognizing
that they all point to different contributions to spin-phonon coupling and
that they all must be considered at the same time in analyzing spin-lattice
relaxation times. This analysis holds for $S=1/2$ in high-field, were the
$g$-tensor is likely to drive relaxation. The study of $S>1/2$ transition
metals or lanthanides needs to follow a different strategy. However, the
analysis of $g$-tensor contributions to spin-phonon coupling brings forward an
important general point: strong spin-phonon coupling arises from the
modulation of those molecular degrees of freedom that more intensely
contribute to the magnitude of the static spin Hamiltonian parameters.
Interestingly, this analysis can also be extended to the discussion of
molecular symmetrysantanni2020probing ; Kazmierczak2021 .

Similarly to the case of $g$-tensor, magnetic anisotropy $\mathbf{D}$ in
transition metals can be described within a ligand field formalism and second-
order perturbation theory. Taking analytical derivatives of $\mathbf{D}$ one
obtains an expression equivalent to Eq. 196, pointing once again to the energy
of excited electronic states as the key variable influencing spin-phonon
couplinghigdon2020spin . In the case of transition metals with very large
anisotropy, electronic excited states are generally so close to the ground
state that the perturbation theory used to describe $D$ in the ligand field
approach breaks down. This electronic configuration leads to a pseudo-
degenerate ground state which maximise the effect of spin-orbit coupling
beyond perturbative regime, similarly to what happens in
lanthanidesatanasov2011detailed . The molecular vibrations that contribute to
the removal of such electronic degeneracy will then be those with large spin-
phonon coupling. In the case of [(tpaPh)Fe]− it was in fact observed that the
vibrations that maximally couple to spin are those Jahn-Teller-like
distortions that remove the C3v symmetry axis of the molecule by bending the
in-plane $\widehat{NFeN}$ angles, thus quenching the residual angular
momentumlunghi2017intra .

The spin properties of single-ion Lanthanide complexes are largely dominated
by electrostatic contributions, as discussed by Rinehart et
alrinehart2011exploiting . Depending on the element of the Ln series, the
electronic density associated to different $M_{J}$ states can either be
prolate, oblate or spherical. Crystal field interactions would then stabilize
certain shapes of electronic density depending on the spatial orientation of
the ligands. For instance, the states with $M_{J}=\pm 15/2$ and $M_{J}=\pm
1/2$ of Dy3+ ions are associate with an oblate and prolate electronic density,
respectively. A very axial distribution of ligands around the Dy ion would
then stabilize in energy the $M_{J}=\pm 15/2$ states with respect to
$M_{J}=\pm 1/2$, and vice versa for equatorial crystal fieldsungur2011magnetic
; chilton2015first ; chilton2015design ; ungur2016strategies . As a
consequence, it is natural to expect that vibrations able to break the axial
symmetry of the Dy ion’s crystal field will be able to strongly couple to its
magnetic moment. It has also been show that a second coupling mechanism can
arisebriganti2021 . This is due to the fact that the effective charge of the
first coordination sphere, determining the crystal field felt by the ion, is
also modulated during molecular motion. This adds a second channel of spin-
phonon coupling and in Ref. briganti2021  it was shown that even vibrations
that little affect the shape of the first coordination shell are able to
strongly couple to spin if the second-coordination shell motion is able to
modulate the ligands’ local charges. Fig. 11 provides a schematic
representation of these two different spin-phonon coupling mechanisms for the
molecule Dy(acac)3(H2O)2jiang2010mononuclear .

Figure 11: Spin-phonon coupling in Dy(acac). Panel A schematically shows two
different types of vibrations similarly coupled to spin. Panel B reports the
variation of local charge density on the Oxygen atoms during the molecular
motion. Reprinted with permission from ref.briganti2021 . Copyright 2021
American Chemical Society.

Vib1, schematically represented in the left panel of Fig. 11, is a molecular
vibration where the first coordination shell of Dy is little distorted, but at
the same time the charges of the ligands show large fluctuations due to the
vibrations involving the second coordination shell. Vice versa, Vib2
represents the typical case where local charges remain largely constant during
vibration but the first coordination shell gets significantly distorted. Both
types of vibrations were found to be strongly coupled to spin.

##  VI The mechanism of spin-phonon relaxation in magnetic molecules

In this last section of the chapter we will finally address the mechanism of
spin-phonon relaxation and the prediction of its rate.

### Direct and Raman Relaxation in $S=1/2$ systems

Figure 12: Relaxation mechanism in VO(acac)2. The spin levels are reported in
the section A of the scheme. The part B of the scheme instead report the
vibrational levels involved in direct and Raman relaxation. Reprinted with
permission from ref. lunghi2020limit . Copyright 2020 American Chemical
Society.

The molecule VO(acac)2 will serve once again as benchmark system to present
the main features of spin-phonon relaxation in $S=1/2$ systems. In the case of
a two-level system, one- and two-phonon transition rates are determined by
second-order perturbation theory, and first- and second-order coupling
strength, respectively. In Refs lunghi2019phonons  and lunghi2020limit  the
numerical strategies of section IV were used to compute the Redfield matrices
of Eqs. 58 and 85 and determine the time evolution of the magnetization
moment. The main results of that studies are reported in Figs. 12 and 13

Let us begin our analysis from direct spin relaxation. This relaxation
mechanism involves transitions between spin states due to resonant one-phonon
absorption or emission, as depicted in Fig. 12. Ab initio spin dynamics
simulations were carried out as function of magnetic field intensity and
showed that two possible relaxation regimes are possible: the first one is
active at low field, where the modulation of hyperfine coupling drives
relaxation, and the second one is active at high field, where the modulation
of the $g$-tensor is instead the main source of relaxation. For VO(acac)2 the
transition field between the two regimes was estimate to be around $B^{*}\sim$
1 T. Importantly, the two regimes have a different field dependence, with the
former going as $\tau\sim B^{-2}$ and the latter as $\tau\sim B^{-4}$.
Regarding the temperature dependence of relaxation, it was observed that
$\tau\sim T^{-1}$ at high $T$, and $\tau\sim T^{0}$ for $T\rightarrow 0$.

Figure 13: Spin-phonon relaxation in VO(acac)2. Panel A reports the simulated
Raman relaxation time as function of field and temperature. Panel B and C
reports both Raman (black) and direct (red) simulated relaxation time together
with experimental values from AC measurements (blue) and inversion recovery on
diluted sample of similar molecules (black star). Panel B refers to 20 K and
Panel C to 5 T. Panel D reports the molecular displacements associated to the
first two optical modes at the $\Gamma$-point. Reprinted with permission from
ref. lunghi2020limit . Copyright 2020 American Chemical Society.

Interestingly, ab initio spin dynamics simulations reveal a picture of direct
relaxation in agreement with the conclusions of section III. As anticipated in
section V, this is because at very low vibrational frequencies typical of
Zeeman/hyperfine splittings, the acoustic modes of a molecular crystal
qualitatively behave as prescribed by the Debye model. Therefore, for Zeeman
splittings much lower than the energy of the optical phonons, the qualitative
$\tau$ vs $B/T$ behaviour expected from classical models based on the Debye
model holds. It is important to note, however, that the energies of
VO(acac)2’s optical modes are on the high-end side of this class of compounds,
with several examples known where the first optical vibrations fall as low as
10 cm-1albino2019first , as opposite to 40-50 cm-1 in
VO(acac)2lunghi2019phonons . In these compounds severe deviations from the
conclusions of the Debye model should be expected for fields as low as a few
Tesla.

Let us now turn to the Raman relaxation, which it was found to be driven by
simultaneous absorption and emission of two phonons, as represented
schematically in Fig. 12. Also in the case of Raman relaxation, it was showed
that two field regimes are potentially active, with hyperfine driving
relaxation as $\tau\sim B^{0}$ and $g$-tensor as $\tau\sim
B^{-2}$lunghi2020limit . Once again, this is in agreement with the conclusions
reached in section III. The temperature dependence of Raman relaxation was
found to behave as $T\sim T^{-2}$ for $T>20$ K and then rapidly increase for
lower $T$. In section III we derived a similar behaviour of $\tau$ vs $T$
coming from the acoustic phonons as accounted in the Debye model. However, as
we have noted in the previous section, optical phonons are much more strongly
coupled with spin, and since they fall at surprisingly low energies, it is
natural to look at them as source of spin-phonon relaxation instead of
acoustic ones. The results for VO(acac)2 were then interpreted by noticing
that the Fourier transform of the two-phonon correlation function
$G^{\mathrm{2-ph}}$, reported in Fig. 5C for the process of simultaneous
absorption and emission of two phonons at $T=20$ K, shows a sharp maximum at
some frequency close to the first $\Gamma$-point optical frequency. Since only
the phonons in a small energy windows seem to contribute to relaxation, we can
treat them as an effective localized mode. According to the structure of Eq.
85, a pair of quasi-degenerate phonons contributes to $\tau$ as

$\tau^{-1}\sim\frac{e^{\beta\omega}}{(e^{\beta\omega}-1)^{2}}$ (199)

which gives the correct simulated behaviour in the limit of k${}_{B}T>\omega$.
A relaxation process that follows a single contribution such as Eq. 199 is
often associated with the local-mode relaxation mechanismeaton2002relaxation .
The latter was originally proposed to origin in the presence of local low-
energy optical modes in a magnetically doped crystal due to the mismatch in
vibrational frequency between the magnetic impurities and the diamagnetic
hostklemens1962localized . However, in the case of magnetic molecules there is
no reason to invoke such a mechanism as plenty of low-lying optical modes are
always present, even in defect-free crystalsgarlatti2020unveiling .

The $T^{-2}$ limit closely matches the experimental observation of relaxation
following a low-exponent power law. However, deviations from $T^{-2}$ are
often observed in favour of power-law dependencies $T^{-n}$ with
$n>2$atzori2016room . Three possible explanations are likely to apply to this
scenario. When several strongly-coupled low-energy optical modes are present,
multiple contributions in the form or Eq. 199 overlap, resulting in a more
complex $T$ dependency. On the other hand, when relaxation time is fast and
rapidly exits the window of measurable values $\tau$, it is hard to
demonstrate that the high-$T$ limit has been fully established. In this
scenario, the profile of $\tau\>\text{vs}\>T$ might still be affected by the
pseudo-exponential regime and therefore results might erroneously be fitted
with a power-law with exponent larger than $n=2$. Finally, at this stage it is
impossible to exclude the contribution of spin-phonon coupling strength beyond
the quadratic order and that should be considered in future works. Despite the
limitations of Eq. 199, we suggest that it should be preferred to power laws
when fitting experimental data. The latter, unless supported by a physical
picture as the one presented in section III, do not provide insights on the
spin relaxation mechanism. On the other hand, Eq. 199 contains the frequency
of the phonons contributing to relaxation and therefore is able to provide a
clearer interpretation of relaxation experiments. Such approach has been
successfully applied to a series of spin 1/2 systemssantanni2020probing ;
de2021exploring ; Kazmierczak2021 ; pfleger2021terminal .

It is important to remark that the agreement between experiments and
simulations become acceptable only in the regimes of high-$T$/high-$B$. This
is due to the fact that in the opposite regime, dipolar spin-spin coupling,
not included in the simulations, becomes the dominant source of relaxation.
Indeed, when diluted samples are measured with EPR inversion recovery, a much
better agreement is observedlunghi2021towards .

Ab initio spin dynamics simulations also make it possible to determine the
phonons responsible for relaxation. In the case of direct relaxation, only
phonons in resonance with the transition are able to contribute. As depicted
in Figs. 7, at fields lower than 10 T only acoustic phonons are present, and
their intra-molecular component is the main drive for relaxation. In the case
of Raman relaxation instead, all the phonons in the spectrum are potentially
able to contribute. However, only low-energy optical phonons significantly
contribute to spin dynamics as they are at the same time thermally populated
and strongly-coupled. Indeed, the population of high-energy phonons diminish
exponentially with their frequency, while the acoustic phonons, despite being
thermally populated, have a low-density of states. This can also be
appreciated from Fig. 7, where the density of states rapidly drops below 40-50
cm-1. Fig. 13D shows typical molecular distortions associated with the first
few optical modes at the $\Gamma$-point.

### Orbach and Raman Relaxation in $S>1/2$ systems with uni-axial anisotropy

Our second case-study is a high-spin Co2+ complex ($S=3/2$) with large uni-
axial magnetic anisotropyrechkemmer2016four , namely with $D=-115$ cm-1 and
vanishing rhombicity $|E/D|$. In the absence of an external field, the direct
transitions within Kramers’ doublets $M_{S}=\pm n/2$ are prohibited by time-
reversal symmetry and one-phonon relaxation must occur through the Orbach
mechanism. Interestingly, the same selection-rule applies to two-phonon
transitions due to second-order perturbation theory and quadratic spin-phonon
coupling, as the matrix elements of Eqs. 58 and 85 have the same dependency on
spin operators and $\hat{H}_{s}$ eigenfunctions. The source of two-phonon
relaxation is therefore to be looked for in Eq. 119, which uses fourth-order
perturbation theory.

Figure 14: Orbach and Raman Relaxation in a Co2+ complex. Left Panel: dashed
lines corresponds to the measured relaxation time in zero external field. The
red line and triangles is the computed Orbach relaxation and the blue line and
triangles is the computed Raman relaxation. Right Panel: (A) The spin states
of Co(pdms)2 are reported as function of their energy and the nominal $S_{z}$
value. The red and blue arrows represent possible relaxation pathways. (B) An
electronic excitation from the ground state to an excited state (red arrow in
panel A) can be accompanied by the absorption of a phonon with energy in
resonance with the spin transition (Orbach process). (C) The direct transition
$S_{z}$=3/2 $\rightarrow$ $S_{z}$=-3/2 (blue arrow in panel A) can be induced
by the simultaneous absorption and emission of two phonons with the same
energy (Raman process). Reprinted from ref. lunghi2020multiple , with the
permission of AIP Publishing.

The main results of Ref. lunghi2020multiple  are presented in Fig. 14. At high
temperature the Orbach mechanisms becomes the most favourable relaxation
pathway and it is dominated by one-phonon absorption of phonons resonant with
the $M_{S}=\pm 3/2\rightarrow M_{S}=\pm 1/2$ transition, as depicted in the
right panel of Fig. 14. At low temperature, the Raman mechanism instead drives
relaxation and show a pseudo-exponential trend with respect to $T$. From the
analysis of the transition rates computed with Eq. 119, it was possible to
interpret the Raman relaxation as a direct spin flip $M_{S}=\pm 3/2\rightarrow
M_{S}=\mp 3/2$ due to the simultaneous absorption and emission of two
degenerate phonons. This process is represented in the right panel of Fig. 14.
Overall, the temperature behaviour of Orbach and Raman spin-relaxation were
found to follow the expression

$\tau^{-1}=\tau_{0}^{-1}e^{-\beta
U_{\mathrm{eff}}}+\tilde{V}_{\mathrm{2-sph}}\frac{e^{\beta\tilde{\omega}}}{(e^{\beta\tilde{\omega}}-1)^{2}}\>,$
(200)

In Eq. 200, $U_{\mathrm{eff}}$ is the energy of the excited Kramers doublet
and $\tilde{\omega}$ is a value of energy compatible with the lowest-lying
optical modes, computed around $\sim 20$ cm-1. Similarly to the case of
$S=1/2$, Raman relaxation receives contributions from all the vibrational
spectrum. Among all the phonons, those contemporaneously populated and coupled
to spin will drive relaxation. However, in the case of $S=3/2$, the matrix
element of Eq. 119 will also favour phonons in resonance with the spin excited
states. In highly anisotropic compounds such as this, the excited Kramers
doublet lies at a much higher energy than the first optical modes, suggesting
that phonons well above the first optical transitions should be involved in
Raman relaxation. However, it should be noted that the resonant condition with
the excited states is only imposed at the power of two (see Eq. 119), while
the population of phonons decreases exponentially with their energy, as
dictated by the Bose-Einstein population. At low temperature the latter
condition overcomes the former, leading to the conclusion that low-energy
modes are the main drive for Raman relaxation in both $S=1/2$ and $S>1/2$. As
discussed for $S=1/2$, acoustic modes do not generally contribute
significantly to Raman relaxation in molecules as they have a very low density
of states. However, at very low temperature ($T<5$ K), when even optical modes
are strongly un-populated, then acoustic phonons might become relevant.
Similarly to what discussed for relaxation in $S=1/2$, the use of Eq. 200
should be preferred to power laws, as commonly done.

As it can be appreciated from Fig. 200, the accuracy reached in Ref.
lunghi2020multiple  is not fully satisfactory and the error between
experiments and predictions is around one and two orders of magnitude. As
discussed in Ref. lunghi2021towards , this is almost entirely due to the
misuse of the secular approximation, which does not totally decouple coherence
and population terms of the density matrix for Kramers systems in zero
external field. Once the secular approximation is correctly carried out, the
agreement between experiments and theory becomes virtually
exactlunghi2021towards .

Figure 15: Orbach and Raman Relaxation in Dy3+ complexes. Experimental
relaxation time for Dy(acac) (black dots) together to the simulated Orbach
(red line and symbols) and Raman (green line and symbols) ones. Adapted with
permission from ref. briganti2021 . Copyright 2021 American Chemical Society.

Dy3+ complexes with axial symmetry have received a large attention from the
molecular magnetism community in virtue of their slow spin
relaxationluzon2012lanthanides . Several molecules have then been studied by
means of ab initio methods to elucidate the role of spin-phonon coupling and
predict spin relaxation rates. Ref. goodwin2017molecular  represents the first
attempt at computing the relaxation times in Ln complexes. The system studied
was [DyCp${}^{\textrm{ttt}}_{2}$]+
(Cp${}^{\textrm{ttt}}$=[C5H${}_{5}^{\;\;\textrm{t}}$Bu3-1,2,4]), DyCp in
short, i.e. the first single-molecule magnet with a blocking temperature
approaching Nitrogen boiling point. DyCp exhibits a Dy3+ ion with a strongly
anisotropic ground state $J=15/2$ due to an almost perfectly axial
coordination geometry. As mentioned in section V, this situation is ideal to
create very large energy gaps between the KDs with $M_{J}=\pm 15/2$ and
$M_{J}=\pm 1/2$, which is estimated to be around 1500 cm-1 in DyCp. The
simulations presented in in Ref. goodwin2017molecular  were able to capture
the experimental results in the high-$T$ regime ($T>60$ K) with only one order
of magnitude of error. Simulations shown that the modulation of the effective
crystal field operator was responsible for transitions involving up to the 6th
excited Kramers doublet, supporting the use of the first term of Eq. 200 to
fit the Orbach relaxation regime and determining an effective anisotropy
barrier $U_{\mathrm{eff}}\sim 1225$ cm-1. The prediction of Orbach rates with
ab initio spin dynamics has been further validated through the study of a
series of Dy3+ molecules with long relaxation time, revealing interesting
correlations between the molecular motion of high-energy modes in resonance
with spin Kramers doublet and spin relaxationevans2019monophospholyl ;
yu2020enhancing ; reta2021ab . In particular, these studies showed that
relaxation may occurs though complex relaxation pathways, often involving
several of the 16 states of the Dy3+ ground-state multiplet. Moreover, it was
shown that, although the performance of slow-relaxing single-molecule magnets
is primarily due to the axiality and strength of the crystal field splitting,
among molecules showing similar values of $U_{\mathrm{eff}}$, the energy
alignment and the careful design of the vibrational density of states can play
a fundamental role in bringing down relaxation rates even furtherreta2021ab .
These series of studies were performed including only gas-phase vibrations as
well as considering only one-phonon processes. As a consequence they were not
able to predict spin relaxation at low temperature, where a second pseudo-
exponential relaxation regime appears, similarly to what observed in Co2+ (see
Fig. 14).

The compound Dy[(acac)3(H2O)2], Dyacac in short, was also studied with ab
initio spin dynamics, this time including all the phonons of the crystal’s
unit-cell and both one- and two-phonon transitionsbriganti2021 . Similarly to
what discussed for Co2+, relaxation occurs through Raman mechanism at low
temperature while Orbach is the dominating mechanism at high temperature. Once
again the relaxation time depends on temperature as reported in Eq. 200 and
depicted in Fig. 15. Differently from DyCp, the Orbach relaxation in Dyacac
was observed to be mediated by the first two excited Kramers doublets, due to
the low symmetry of the coordination of the Dy ion and the consequently large
admixing between different $M_{J}$. Regarding the Raman mechanism it was
instead observed that only the first few optical modes at the $\Gamma$-point
contribute to relaxation in virtue of a larger thermal population at low
temperature.

Simulations for DyCp were repeated using a similar approach to the one
employed for Dyacaclunghi2021towards , confirming that a two-phonon Raman
mechanism is able to explain spin relaxation at low-$T$. A similar
interpretation to the dynamics in Dyacac also applies to DyCp. Despite the
much higher energy of the spin transitions in DyCp respect to Dyacac ( the
first spin excited state falls at $\sim$ 450 cm-1 in DyCp and $\sim$ 100 cm-1
in Dyacac), simulations confirm that low-energy vibrations up to a few THz are
the only contribution to Raman spin-phonon relaxation. On the other hand, only
high-energy optical phonons in resonance with spin transitions are responsible
for Orbach relaxation. Importantly, the results on DyCp’s spin dynamics in
zero-field showed a remarkable dependence on the inclusion of coherence terms
in the spin reduced density matrix, with deviations of up to 9 orders of
magnitude from the experimental results when neglected. Vice versa, as
observed for Co2+ system discussed above, once the secular approximation is
correctly implemented, simulations become in perfect agreement with
experiments over the entire range of temperaturelunghi2021towards , therefore
providing a conclusive proof that spin-phonon relaxation in magnetic molecules
of Kramers ions presenting easy-axis magnetic anisotropy is fully understood.

## Outlook

In this chapter we have illustrated how the theory of open quantum systems can
be adapted to the problem of spin-phonon interaction and provide a
quantitative ground for the study of spin dynamics. Despite the handful of
systems studied so far, ab initio simulations have already provided deep
insights into the nature of spin-phonon coupling and spin relaxation in
crystals of magnetic molecules. Here we will attempt to summarize the most
important points.

Most of our understanding of spin relaxation in solid-state comes from the
application of the Debye model and perturbation theory. However, differently
from the solid-state doped materials initially studied, molecular crystals
posses many low-lying optical modes. Moreover, as one of the main goals
pursued in the field of molecular magnetism is the enhancement of magnetic
anisotropy, the energy splitting among spin states have been increasing by
orders of magnitudes along the yearszabala2021single ; duan2021data . This
situation made it so that the assumptions of section III have lost their
general validity. Ab initio simulations showed that for molecular materials,
acoustic modes are not particularly relevant and that optical ones are instead
responsible for relaxation. Moreover, it was found that optical modes appear
at surprisingly low energies, so that even at very low temperature they remain
sensibly populated. Thanks to ab initio spin dynamics, it has now been
possible to demonstrate the nature of spin relaxation mechanism in Kramers
systems. Direct and Orbach mechanisms are due to the absorption and emission
of resonant phononsgoodwin2017molecular ; lunghi2019phonons . The origin of
Raman relaxation is attributed to different mechanisms for $S=1/2$ and $S>1/2$
systems. In the former case, quadratic spin-phonon coupling is responsible for
two-phonon transitionslunghi2020limit , while in the latter case, fourth-order
perturbation theory is necessary to explain spin-relaxation in zero-external
fieldchiesa2020understanding ; gu2020origins ; lunghi2020multiple ;
briganti2021 . In both cases, Raman is mediated by low-energy optical phonons
exhibiting local molecular rotations admixed to delocalized intra-molecular
distortions.

These findings have major consequences for molecular magnetism practices. In
addition to well established approaches to increase the axiality of the
crystal fieldrinehart2011exploiting  and exchange coupling multiple
ionsGould2022ultrahard , ab initio spin dynamics results point to the design
of rigid molecular structures as key to reduce spin-phonon relaxation
rateslunghi2017intra . Indeed, by increasing molecular rigidity, optical modes
are shifted up in energy with a two-fold effect: i) they become less admixed
with acoustic modeslunghi2019phonons , so reducing the rate of direct
relaxation, and ii) they become less populated so that Raman relaxation rate
would also slow downlunghi2020limit ; lunghi2020multiple ; briganti2021 . A
second approach suggested by ab initio spin dynamics involves the careful
disalignement of spin and vibrational spectra in order to reduce the effect of
resonant Orbach processgoodwin2017molecular ; ullah2019silico ;
yu2020enhancing .

Finally, although ab initio spin dynamics is not yet a fully quantitative
computational approach, we have shown that it is rapidly maturing into one.
Since its first proposals in 2017escalera2017determining ; lunghi2017role ;
lunghi2017intra ; goodwin2017molecular , a level of accuracy down to one order
of magnitude has now been proved for different systems across different
relaxation regimeslunghi2021towards . Although further benchmarking is needed,
this level of accuracy already hints at the possibility to rank the spin
lifetime of different compounds in advance of experimental characterization.
This is an important milestone for the field of molecular magnetism, as it
opens the gates to a long-sought fully in silico screening of new compounds.

Let us now attempt to provide an overview of future developments as well as a
brief discussion of the potential impact of ab initio spin dynamics for other
adjacent fields.

Ab initio theory of spin relaxation now includes all the terms up to two-
phonon processes. This is essentially the same level of theory employed in the
past to derive the canonical picture of spin-phonon relaxation and we can
therefore claim that the theoretical formulation of such theories in terms of
ab initio models is available. However, this is far from claiming that no
further work will be necessary to obtain a conclusive and universal picture of
spin relaxation in solid-state. On the contrary, this theoretical advancements
form the stepping stone for an in-depth analysis of relaxation phenomenology,
for instance going beyond the classical Born-Markov approximation and the
inclusion of other decoherence channels such as magnetic dipolar interactions.
These two avenues in fact represent some of the most urgent development of the
proposed theoretical framework and have the potential to account for phonon-
bottlenecktesi2016giant ; ullah2021insights  and magnetization quantum
tunnelling effectsding2018field , respectively. These contributions to
relaxation are indeed ubiquitous in experiments and must find their place into
a quantitative theory of spin dynamics. Moreover, it is important to remark
that so far ab initio spin dynamics has been applied to only a handful of
systems and there is large scope for investigating the large and varied
phenomenology of molecular magnetism. Examples includes the study of exchange-
coupled poly-nuclear magnetic moleculeszabala2021single , surface-adsorbed
systemsmannini2009magnetic  and molecules in glassy solutionssato2007impact .

In order to reach a fully quantitative ground, ab initio spin dynamics now
requires a careful testing of the electronic structure methods underlying it.
It is well known that the inclusion of electronic correlation in the
simulations is key to achieve accurate prediction of spin Hamiltonian
parametersatanasov2015first ; ungur2017ab , and the same effect can be
expected for spin-phonon coupling coefficients. On the other hand, the
simulation of phonon modes requires significant advancements. The introduction
of sophisticated vdW corrections to DFTtkatchenko2012accurate  and the
inclusion of anharmonic lattice effectsruggiero2020invited  stand out as two
important avenues for future development. The development of efficient
simulation strategies is another necessary step towards the widespread use of
ab initio spin dynamics. On the one hand, the ab initio methods underlying
spin relaxation theory are particularly heavy in computational terms. Magnetic
molecules indeed often crystallize in large unit cells and the simulations of
phonons beyond the $\Gamma$-point is a daunting task. Similarly, the
calculation of spin-phonon coupling coefficients requires the quantum
chemistry simulations for many similarly distorted molecular geometries,
rapidly leading to immense computational overheads. In section V we have
hinted to machine learning as a game-changer in this area. Machine learning,
learning the underlying features of a distribution from examples drawn out it,
makes it possible to interpolate complex functions such as the potential
energy surface of a molecular crystallunghi2019unified ; lunghi2020multiple
or the relation between spin Hamiltonian and molecular
structurelunghi2020surfing . After the training stage, ML models can be used
to make predictions at very little computational cost, therefore offering a
significant speed-up of these simulations. Further adaptation of machine
learning schemes to this tasks is needed, but encouraging proof-of-concept
applications have recently been presentedlunghi2020surfing ; lunghi2020limit ;
lunghi2020multiple ; lunghi2020insights ; nandy2021computational ;
zaverkin2021thermally . Another interesting strategy involve the use of
parametrized Hamiltonians to compute the spin-phonon coupling coefficients.
When an accurate parametrization is possible, extensive speed-ups can be
achievedescalera2020design .

Although this chapter has been focusing on crystals of magnetic molecules,
many more systems are within the reach of the methods we have illustrated.
Indeed, the only strong assumptions of ab initio spin dynamics are that the
system must be well described with a spin Hamiltonian formalism and that it
must exist an accurate ab initio method to predict lattice and spin
properties. The computational scheme we have discussed should readily apply to
any kind of system exhibiting localized magnetism. Indeed, solid-state defects
or impurities in solid-state semiconductors are commonly described with the
same theoretical tools used in the field of molecular magnetism, and DFT-based
schemes for the simulation of their properties already existsivady2018first .
Many of these systems are currently under scrutiny for quantum sensing
technologies and there is a large interest in their relaxation
propertieswolfowicz2021quantum . Another type of system that we believe is
easily within the reach of the presented methods is nuclear spin. Indeed,
electronic structure methods are routinely used to predict nuclear spin
Hamiltonian’s parametershelgaker1999ab , and when such nuclear spins are part
of a molecular system, there is virtually no difference from the treatment of
electron spin dynamics.

In conclusion, we have provided a comprehensive overview of the state-of-the-
art in ab initio spin dynamics as well as a practical guide to its application
to crystals of magnetic molecules. We have shown that this novel computational
method is rapidly maturing into a quantitative tool able to provide unique
insights into spin relaxation, thus holding great promise for the field of
molecular magnetism and beyond. As such, we hope that this work will serve as
a starting point for the readers interested in adopting this strategy as well
as a review of the strives so far for expert readers.

Acknowledgements

We acknowledge funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No. [948493])

## References

  * (1) Bencini, A. & Gatteschi, D.  _EPR of exchange coupled systems_ (Courier Corporation, 2012).
  * (2) Berliner, L. J., Eaton, S. S. & Eaton, G. R.  _Distance measurements in biological systems by EPR_ , vol. 19 (Springer Science & Business Media, 2006).
  * (3) Zhang, S., Winter, P., Wu, K. & Sherry, A. D.  A novel Europium (III)-based MRI contrast agent.  _Journal of the American Chemical Society_ 123, 1517–1518 (2001).
  * (4) Sessoli, R., Gatteschi, D., Caneschi, A. & Novak, M.  Magnetic bistability in a metal-ion cluster.  _Nature_ 365, 141–143 (1993).
  * (5) Zabala-Lekuona, A., Seco, J. M. & Colacio, E.  Single-molecule magnets: From mn12-ac to dysprosium metallocenes, a travel in time.  _Coordination Chemistry Reviews_ 441, 213984 (2021).
  * (6) Žutić, I., Fabian, J. & Sarma, S. D.  Spintronics: Fundamentals and applications.  _Reviews of Modern Physics_ 76, 323 (2004).
  * (7) Bogani, L. & Wernsdorfer, W.  Molecular spintronics using single-molecule magnets.  _Nature Materials_ 7, 179–186 (2008).
  * (8) Sanvito, S.  Molecular spintronics.  _Chemical Society Reviews_ 40, 3336–3355 (2011).
  * (9) Aromí, G., Aguila, D., Gamez, P., Luis, F. & Roubeau, O.  Design of magnetic coordination complexes for quantum computing.  _Chemical Society Reviews_ 41, 537–546 (2012).
  * (10) Gaita-Ariño, A., Luis, F., Hill, S. & Coronado, E.  Molecular spins for quantum computation.  _Nature Chemistry_ 11, 301–309 (2019).
  * (11) Atzori, M. & Sessoli, R.  The second quantum revolution: role and challenges of molecular chemistry.  _Journal of the American Chemical Society_ 141, 11339–11352 (2019).
  * (12) Wasielewski, M. R. _et al._ Exploiting chemistry and molecular systems for quantum information science.  _Nature Reviews Chemistry_ 4, 490–504 (2020).
  * (13) Heinrich, A. J. _et al._ Quantum-coherent nanoscience.  _Nature Nanotechnology_ 16, 1318–1329 (2021).
  * (14) Mannini, M. _et al._ Magnetic memory of a single-molecule quantum magnet wired to a gold surface.  _Nature Materials_ 8, 194–197 (2009).
  * (15) Urdampilleta, M., Klyatskaya, S., Cleuziou, J.-P., Ruben, M. & Wernsdorfer, W.  Supramolecular spin valves.  _Nature Materials_ 10, 502–506 (2011).
  * (16) Bader, K. _et al._ Room temperature quantum coherence in a potential molecular qubit.  _Nature Communications_ 5, 5304 (2014).
  * (17) Zadrozny, J. M., Niklas, J., Poluektov, O. G. & Freedman, D. E.  Millisecond coherence time in a tunable molecular electronic spin qubit.  _ACS Central Science_ 1, 488–492 (2015).
  * (18) Atzori, M. _et al._ 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).
  * (19) Ariciu, A.-M. _et al._ Engineering electronic structure to prolong relaxation times in molecular qubits by minimising orbital angular momentum.  _Nature Communications_ 10, 3330 (2019).
  * (20) Atzori, M. _et al._ A two-qubit molecular architecture for electron-mediated nuclear quantum simulation.  _Chemical Science_ 9, 6183–6192 (2018).
  * (21) Godfrin, C. _et al._ Operating quantum states in single magnetic molecules: implementation of Grover’s quantum algorithm.  _Physical Review Letters_ 119, 187702 (2017).
  * (22) Bayliss, S. _et al._ Optically addressable molecular spins for quantum information processing.  _Science_ 370, 1309–1312 (2020).
  * (23) Moreno-Pineda, E., Godfrin, C., Balestro, F., Wernsdorfer, W. & Ruben, M.  Molecular spin qudits for quantum algorithms.  _Chemical Society Reviews_ 47, 501–513 (2018).
  * (24) Carretta, S., Zueco, D., Chiesa, A., Gómez-León, Á. & Luis, F.  A perspective on scaling up quantum computation with molecular spins.  _Applied Physics Letters_ 118, 240501 (2021).
  * (25) Biard, H. _et al._ Increasing the Hilbert space dimension using a single coupled molecular spin.  _Nature Communications_ 12, 4443 (2021).
  * (26) Gimeno, I. _et al._ Enhanced molecular spin-photon coupling at superconducting nanoconstrictions.  _ACS Nano_ 14, 8707–8715 (2020).
  * (27) Urtizberea, A. _et al._ Vanadyl spin qubit 2D arrays and their integration on superconducting resonators.  _Materials Horizons_ 7, 885–897 (2020).
  * (28) Atzori, M. _et al._ Radiofrequency to microwave coherent manipulation of an organometallic electronic spin qubit coupled to a nuclear qudit.  _Inorganic Chemistry_ 60, 11273–11286 (2021).
  * (29) Deutsch, J. M.  Eigenstate thermalization hypothesis.  _Reports on Progress in Physics_ 81, 082001 (2018).
  * (30) Nandkishore, R. & Huse, D. A.  Many-body localization and thermalization in quantum statistical mechanics.  _Annu. Rev. Condens. Matter Phys._ 6, 15–38 (2015).
  * (31) Breuer, H.-P., Petruccione, F. _et al._ _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * (32) Breuer, H.-P., Kappler, B. & Petruccione, F.  The time-convolutionless projection operator technique in the quantum theory of dissipation and decoherence.  _Annals of Physics_ 291, 36–70 (2001).
  * (33) Timm, C.  Time-convolutionless master equation for quantum dots: Perturbative expansion to arbitrary order.  _Physical Review B_ 83, 115416 (2011).
  * (34) Waller, I.  Über die magnetisierung von paramagnetischen kristallen in wechselfeldern.  _Zeitschrift für Physik_ 79, 370–388 (1932).
  * (35) Van Vleck, J.  Paramagnetic relaxation times for titanium and chrome alum.  _Physical Review_ 57, 426 (1940).
  * (36) Orbach, R.  Spin-lattice relaxation in rare-earth salts.  _Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences_ 264, 458–484 (1961).
  * (37) Klemens, P.  Localized modes and spin-lattice interactions.  _Physical Review_ 125, 1795 (1962).
  * (38) Singh, A. & Shrivastava, K.  Optical-acoustic two-phonon relaxation in spin systems.  _Physica Status Solidi (b)_ 95, 273–277 (1979).
  * (39) Gill, J.  The establishment of thermal equilibrium in paramagnetic crystals.  _Reports on Progress in Physics_ 38, 91 (1975).
  * (40) Stevens, K.  The theory of paramagnetic relaxation.  _Reports on Progress in Physics_ 30, 189 (1967).
  * (41) Shrivastava, K.  Theory of spin–lattice relaxation.  _Physica Status Solidi (b)_ 117, 437–458 (1983).
  * (42) Eaton, S. S. & Eaton, G. R.  Relaxation times of organic radicals and transition metal ions.  In _Distance measurements in biological systems by EPR_ , 29–154 (Springer, 2002).
  * (43) Sato, H. _et al._ Impact of molecular size on electron spin relaxation rates of nitroxyl radicals in glassy solvents between 100 and 300 k.  _Molecular Physics_ 105, 2137–2151 (2007).
  * (44) Ding, Y.-S. _et al._ Field-and temperature-dependent quantum tunnelling of the magnetisation in a large barrier single-molecule magnet.  _Nature Communications_ 9, 3134 (2018).
  * (45) Moseley, D. H. _et al._ Spin–phonon couplings in transition metal complexes with slow magnetic relaxation.  _Nature Communications_ 9, 2572 (2018).
  * (46) Chiesa, A. _et al._ Understanding magnetic relaxation in single-ion magnets with high blocking temperature.  _Physical Review B_ 101, 174402 (2020).
  * (47) Neese, F., Atanasov, M., Bistoni, G., Maganas, D. & Ye, S.  Chemistry and quantum mechanics in 2019: give us insight and numbers.  _Journal of the American Chemical Society_ 141, 2814–2824 (2019).
  * (48) Atanasov, M. _et al._ 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).
  * (49) Ungur, L. & Chibotaru, L. F.  Ab initio crystal field for lanthanides.  _Chemistry–A European Journal_ 23, 3708–3718 (2017).
  * (50) Maurice, R. _et al._ Universal theoretical approach to extract anisotropic spin hamiltonians.  _Journal of Chemical Theory and Computation_ 5, 2977–2984 (2009).
  * (51) Neese, F.  Spin-hamiltonian parameters from first principle calculations: theory and application.  In _High Resolution EPR_ , 175–229 (Springer, 2009).
  * (52) Garlatti, E. _et al._ Unveiling phonons in a molecular qubit with four-dimensional inelastic neutron scattering and density functional theory.  _Nature Communications_ 11, 1751 (2020).
  * (53) Tennant, W., Walsby, C., Claridge, R. & McGavin, D.  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).
  * (54) Qian, K. _et al._ Does the thermal evolution of molecular structures critically affect the magnetic anisotropy?  _Chemical Science_ 6, 4587–4593 (2015).
  * (55) Lunghi, A. & Sanvito, S.  How do phonons relax molecular spins?  _Science advances_ 5, eaax7163 (2019).
  * (56) Lunghi, A. & Sanvito, S.  The limit of spin lifetime in solid-state electronic spins.  _The Journal of Physical Chemistry Letters_ 11, 6273–6278 (2020).
  * (57) Lunghi, A. & Sanvito, S.  Multiple spin–phonon relaxation pathways in a kramer single-ion magnet.  _The Journal of Chemical Physics_ 153, 174113 (2020).
  * (58) Lunghi, A.  Towards exact predictions of spin-phonon relaxation times: an ab initio implementation of open quantum systems theory.  _arXiv preprint arXiv:2112.09236_ (2021).
  * (59) Kühne, T. D. _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).
  * (60) Neese, F., Wennmohs, F., Becker, U. & Riplinger, C.  The orca quantum chemistry program package.  _The Journal of Chemical Physics_ 152, 224108 (2020).
  * (61) Kresse, G., Furthmüller, J. & Hafner, J.  Ab initio force constant approach to phonon dispersion relations of diamond and graphite.  _EPL (Europhysics Letters)_ 32, 729 (1995).
  * (62) Wang, Y. _et al._ A mixed-space approach to first-principles calculations of phonon frequencies for polar materials.  _Journal of Physics: Condensed Matter_ 22, 202201 (2010).
  * (63) Neese, F.  Prediction of molecular properties and molecular spectroscopy with density functional theory: From fundamental theory to exchange-coupling.  _Coordination Chemistry Reviews_ 253, 526–563 (2009).
  * (64) Cramer, C. J. & Truhlar, D. G.  Density functional theory for transition metals and transition metal chemistry.  _Physical Chemistry Chemical Physics_ 11, 10757–10816 (2009).
  * (65) Perdew, J. P., Burke, K. & Ernzerhof, M.  Generalized gradient approximation made simple.  _Physical Review Letters_ 77, 3865 (1996).
  * (66) Perdew, J. P., Ernzerhof, M. & Burke, K.  Rationale for mixing exact exchange with density functional approximations.  _The Journal of Chemical Physics_ 105, 9982–9985 (1996).
  * (67) Tkatchenko, A. & Scheffler, M.  Accurate molecular van der waals interactions from ground-state electron density and free-atom reference data.  _Physical Review Letters_ 102, 073005 (2009).
  * (68) Grimme, S., Antony, J., Ehrlich, S. & Krieg, H.  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).
  * (69) Chan, G. K.-L. & Sharma, S.  The density matrix renormalization group in quantum chemistry.  _Annual Review of Physical Chemistry_ 62, 465–481 (2011).
  * (70) Ruggiero, M. T.  Invited review: Modern methods for accurately simulating the terahertz spectra of solids.  _Journal of Infrared, Millimeter, and Terahertz Waves_ 41, 491–528 (2020).
  * (71) Lunghi, A., Totti, F., Sanvito, S. & Sessoli, R.  Intra-molecular origin of the spin-phonon coupling in slow-relaxing molecular magnets.  _Chemical Science_ 8, 6051–6059 (2017).
  * (72) Albino, A. _et al._ Temperature dependence of spin–phonon coupling in [vo (acac) 2]: A computational and spectroscopic study.  _The Journal of Physical Chemistry C_ 125, 22100–22110 (2021).
  * (73) Ullah, A., Baldoví, J. J., Gaita-Ariño, A. & Coronado, E.  Insights on the coupling between vibronically active molecular vibrations and lattice phonons in molecular nanomagnets.  _Dalton Transactions_ 50, 11071–11076 (2021).
  * (74) Jung, J. _et al._ Derivation of lanthanide series crystal field parameters from first principles.  _Chemistry–A European Journal_ 25, 15112–15122 (2019).
  * (75) Chibotaru, L. F. & Ungur, L.  Ab initio calculation of anisotropic magnetic properties of complexes. I. unique definition of pseudospin hamiltonians and their derivation.  _The Journal of Chemical Physics_ 137, 064112 (2012).
  * (76) Singh, S. K., Atanasov, M. & Neese, F.  Challenges in multireference perturbation theory for the calculations of the g-tensor of first-row transition-metal complexes.  _Journal of Chemical Theory and Computation_ 14, 4662–4677 (2018).
  * (77) Briganti, M. _et al._ Covalency and magnetic anisotropy in lanthanide single molecule magnets: the DyDOTA archetype.  _Chemical Science_ 10, 7233–7245 (2019).
  * (78) Escalera-Moreno, L., Suaud, N., Gaita-Arino, A. & Coronado, E.  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).
  * (79) Lunghi, A., Totti, F., Sessoli, R. & Sanvito, S.  The role of anharmonic phonons in under-barrier spin relaxation of single molecule magnets.  _Nature Communications_ 8, 14620 (2017).
  * (80) Goodwin, C. A., Ortu, F., Reta, D., Chilton, N. F. & Mills, D. P.  Molecular magnetic hysteresis at 60 kelvin in dysprosocenium.  _Nature_ 548, 439–442 (2017).
  * (81) Harman, W. H. _et al._ Slow magnetic relaxation in a family of trigonal pyramidal Iron (II) pyrrolide complexes.  _Journal of the American Chemical Society_ 132, 18115–18126 (2010).
  * (82) Lunghi, A. & Sanvito, S.  A unified picture of the covalent bond within quantum-accurate force fields: From organic molecules to metallic complexes’ reactivity.  _Science Advances_ 5, eaaw2210 (2019).
  * (83) Lunghi, A. & Sanvito, S.  Surfing multiple conformation-property landscapes via machine learning: Designing single-ion magnetic anisotropy.  _The Journal of Physical Chemistry C_ 124, 5802–5806 (2020).
  * (84) Bartók, A. P., Kondor, R. & Csányi, G.  On representing chemical environments.  _Physical Review B_ 87, 184115 (2013).
  * (85) Lunghi, A.  Insights into the spin-lattice dynamics of organic radicals beyond molecular tumbling: A combined molecular dynamics and machine-learning approach.  _Applied Magnetic Resonance_ 51, 1343–1356 (2020).
  * (86) Albino, A. _et al._ First-principles investigation of spin–phonon coupling in Vanadium-based molecular spin quantum bits.  _Inorganic Chemistry_ 58, 10260–10268 (2019).
  * (87) Rechkemmer, Y. _et al._ A four-coordinate cobalt (ii) single-ion magnet with coercivity and a very high energy barrier.  _Nature Communications_ 7, 10467 (2016).
  * (88) Stavretis, S. E. _et al._ Spectroscopic studies of the magnetic excitation and spin-phonon couplings in a single-molecule magnet.  _Chemistry–A European Journal_ 25, 15846–15857 (2019).
  * (89) Moseley, D. H. _et al._ Inter-kramers transitions and spin–phonon couplings in a lanthanide-based single-molecule magnet.  _Inorganic Chemistry_ 59, 5218–5230 (2020).
  * (90) Spectroscopic analysis of vibronic relaxation pathways in molecular spin qubit [Ho(W5O18)2]9-: Sparse Spectra Are Key, author=Blockmon, Avery L and Ullah, Aman and Hughey, Kendall D and Duan, Yan and O’Neal, Kenneth R and Ozerov, Mykhaylo and Baldoví, José J and Aragó, Juan and Gaita-Ariño, Alejandro and Coronado, Eugenio and others, journal=Inorganic Chemistry, volume=60, number=18, pages=14096–14104, year=2021, publisher=ACS Publications .
  * (91) Marbey, J. _et al._ Vibronic coupling in a molecular 4f qubit.  _10.26434/chemrxiv.13725904.v1_ (2021).
  * (92) Mirzoyan, R. & Hadt, R. G.  The dynamic ligand field of a molecular qubit: decoherence through spin–phonon coupling.  _Physical Chemistry Chemical Physics_ 22, 11249–11265 (2020).
  * (93) Lunghi, A.  Ligand-field contributions to spin-phonon coupling in a family of vanadium molecular qubits from multi-reference electronic structure theory.  _arXiv preprint arXiv:1912.04545_ (2019).
  * (94) Mirzoyan, R., Kazmierczak, N. P. & Hadt, R. G.  Deconvolving contributions to decoherence in molecular electron spin qubits: A dynamic ligand field approach.  _Chemistry: a European Journal_ 27, 9482–9494 (2021).
  * (95) Chakarawet, K. _et al._ Effect of spin–orbit coupling on phonon-mediated magnetic relaxation in a series of zero-valent Vanadium, Niobium, and Tantalum isocyanide complexes.  _Inorganic Chemistry_ 60, 18553–18560 (2021).
  * (96) Fataftah, M. S. _et al._ Metal–ligand covalency enables room temperature molecular qubit candidates.  _Chemical Science_ 10, 6707–6714 (2019).
  * (97) Atzori, M. _et al._ Room-temperature quantum coherence and rabi oscillations in vanadyl phthalocyanine: toward multifunctional molecular spin qubits.  _Journal of the American Chemical Society_ 138, 2154–2157 (2016).
  * (98) Santanni, F. _et al._ Probing vibrational symmetry effects and nuclear spin economy principles in molecular spin qubits.  _Inorganic Chemistry_ 60, 140–151 (2020).
  * (99) Kazmierczak, N. P., Mirzoyan, R. & Hadt, R. G.  The impact of ligand field symmetry on molecular qubit coherence.  _Journal of the American Chemical Society_ 143, 17305–17315 (2021).
  * (100) Higdon, N. J., Barth, A. T., Kozlowski, P. T. & Hadt, R. G.  Spin–phonon coupling and dynamic zero-field splitting contributions to spin conversion processes in Iron (II) complexes.  _The Journal of Chemical Physics_ 152, 204306 (2020).
  * (101) Atanasov, M., Ganyushin, D., Pantazis, D. A., Sivalingam, K. & Neese, F.  Detailed ab initio first-principles study of the magnetic anisotropy in a family of trigonal pyramidal Iron (II) pyrrolide complexes.  _Inorganic Chemistry_ 50, 7460–7477 (2011).
  * (102) Rinehart, J. D. & Long, J. R.  Exploiting single-ion anisotropy in the design of f-element single-molecule magnets.  _Chemical Science_ 2, 2078–2085 (2011).
  * (103) Ungur, L. & Chibotaru, L. F.  Magnetic anisotropy in the excited states of low symmetry lanthanide complexes.  _Physical Chemistry Chemical Physics_ 13, 20086–20090 (2011).
  * (104) Chilton, N. F., Goodwin, C. A., Mills, D. P. & Winpenny, R. E.  The first near-linear bis (amide) f-block complex: a blueprint for a high temperature single molecule magnet.  _Chemical Communications_ 51, 101–103 (2015).
  * (105) Chilton, N. F.  Design criteria for high-temperature single-molecule magnets.  _Inorganic Chemistry_ 54, 2097–2099 (2015).
  * (106) Ungur, L. & Chibotaru, L. F.  Strategies toward high-temperature lanthanide-based single-molecule magnets.  _Inorganic Chemistry_ 55, 10043–10056 (2016).
  * (107) Briganti, M. _et al._ 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).
  * (108) Jiang, S.-D., Wang, B.-W., Su, G., Wang, Z.-M. & Gao, S.  A mononuclear dysprosium complex featuring single-molecule-magnet behavior.  _Angewandte Chemie_ 122, 7610–7613 (2010).
  * (109) de Camargo, L. C. _et al._ Exploring the organometallic route to molecular spin qubits: The [CpTi (cot)] case.  _Angewandte Chemie (International ed. in English)_ 60, 2588–2593 (2021).
  * (110) Pfleger, R. F. _et al._ Terminal ligand and packing effects on slow relaxation in an isostructural set of [Dy(H2dapp)X2]+ single molecule magnets.  _Chemistry–A European Journal_ 27, 15086–15095 (2021).
  * (111) Luzon, J. & Sessoli, R.  Lanthanides in molecular magnetism: so fascinating, so challenging.  _Dalton Transactions_ 41, 13556–13567 (2012).
  * (112) Evans, P., Reta, D., Whitehead, G. F., Chilton, N. F. & Mills, D. P.  Bis-monophospholyl dysprosium cation showing magnetic hysteresis at 48 K.  _Journal of the American Chemical Society_ 141, 19935–19940 (2019).
  * (113) Yu, K.-X. _et al._ Enhancing magnetic hysteresis in single-molecule magnets by ligand functionalization.  _Chem_ 6, 1777–1793 (2020).
  * (114) Reta, D., Kragskow, J. G. & Chilton, N. F.  Ab initio prediction of high-temperature magnetic relaxation rates in single-molecule magnets.  _Journal of the American Chemical Society_ 143, 5943–5950 (2021).
  * (115) Duan, Y. _et al._ Data mining, dashboard and statistical analysis: a powerful framework for the chemical design of molecular nanomagnets.  _arXiv preprint arXiv:2103.03199_ (2021).
  * (116) Gu, L. & Wu, R.  Origins of slow magnetic relaxation in single-molecule magnets.  _Physical Review Letters_ 125, 117203 (2020).
  * (117) Gould, C. A. _et al._ Ultrahard magnetism from mixed-valence dilanthanide complexes with metal-metal bonding.  _Science_ 375, 198–202 (2022).
  * (118) Ullah, A. _et al._ 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).
  * (119) Tesi, L. _et al._ Giant spin–phonon bottleneck effects in evaporable vanadyl-based molecules with long spin coherence.  _Dalton Transactions_ 45, 16635–16643 (2016).
  * (120) Tkatchenko, A., DiStasio Jr, R. A., Car, R. & Scheffler, M.  Accurate and efficient method for many-body van der Waals interactions.  _Physical Review Letters_ 108, 236402 (2012).
  * (121) Nandy, A. _et al._ Computational discovery of transition-metal complexes: from high-throughput screening to machine learning.  _Chemical Reviews_ 121, 9927–10000 (2021).
  * (122) Zaverkin, V., Netz, J., Zills, F., Köhn, A. & Kästner, J.  Thermally averaged magnetic anisotropy tensors via machine learning based on Gaussian moments.  _Journal of Chemical Theory and Computation_ 18, 1–12 (2021).
  * (123) Escalera-Moreno, L., Baldoví, J. J., Gaita-Arino, A. & Coronado, E.  Design of high-temperature f-block molecular nanomagnets through the control of vibration-induced spin relaxation.  _Chemical Science_ 11, 1593–1598 (2020).
  * (124) Ivády, V., Abrikosov, I. A. & Gali, A.  First principles calculation of spin-related quantities for point defect qubit research.  _npj Computational Materials_ 4, 76 (2018).
  * (125) Wolfowicz, G. _et al._ Quantum guidelines for solid-state spin defects.  _Nature Reviews Materials_ 6, 906–925 (2021).
  * (126) Helgaker, T., Jaszunski, M. & Ruud, K.  Ab initio methods for the calculation of NMR shielding and indirect spin-spin coupling constants.  _Chemical Reviews_ 99, 293–352 (1999).