# The role of electronic excited states in the spin-lattice relaxation of
spin-1/2 molecules

Lorenzo A. Mariano1 Vu Ha Anh Nguyen1 Jonatan B. Petersen2 Magnus Bjørnsson3
Jesper Bendix3 Gareth R. Eaton4 Sandra S. Eaton4<EMAIL_ADDRESS>Alessandro Lunghi1<EMAIL_ADDRESS>1School of Physics, AMBER and CRANN
Institute, Trinity College, Dublin 2, Ireland  2Department of Chemistry, The
University of Manchester, Manchester M13 9PL, United Kingdom  3Department of
Chemistry, University of Copenhagen, DK-2100 Copenhagen, Denmark  4Department
of Chemistry and Biochemistry, University of Denver, Denver, CO 80210, USA

###### Abstract

Magnetic resonance is a prime method for the study of chemical and biological
structures and their dynamical processes. The interpretation of these
experiments relies on considering the spin of electrons as the sole relevant
degree of freedom. By applying ab inito open quantum systems theory to the
full electronic wavefunction, here we show that contrary to this widespread
framework the thermalization of the unpaired electron spin of two Cr(V)
coordination compounds is driven by virtual transitions to excited states with
energy higher than 20,000 cm-1 instead of solely involving low-energy spin
interactions such as Zeeman and hyperfine ones. Moreover, we found that a
window of low-energy THz phonons contributes to thermalization, rather than a
small number of high-energy vibrations. This work provides a drastic
reinterpretation of relaxation in spin-1/2 systems and its chemical control
strategies, and ultimately exemplifies the urgency of further advancing an ab
initio approach to relaxometry.

## Introduction

The ubiquitous presence of spin in chemical compounds, either because of the
presence of spin-active nuclei or unpaired electrons, makes magnetic resonance
a central experimental technique for the study of chemical and biological
processes. In particular, the time it takes for a spin to reach its thermal
equilibrium state, namely the spin relaxation time, provides unique insights
into the interactions between the spin and its chemical environment. For
instance, the study of spin relaxation time provides insights into the
dynamics of the molecules they are embedded in, which can be exploited to
determine the aggregation state of bio-moleculeseaton2000relaxation , their
folding dynamicsparigi2021unveiling , the presence of analytes in
solutionsli2019all , and enhance in-vivo imaging
techniquesalzola2022comprehensive . In recent years, relaxometry experiments
have also become a central tool of quantum sensingdegen2017quantum , where the
relaxation of highly coherent spin states is used to probe magnetic fields and
magnetic states with unprecedented sensitivitymzyk2022relaxometry .

$\mathbf{A}$$\mathbf{C}$$\mathbf{B}$ Figure 1: (A) Molecular structures.
CrN(trop)2 (top) and CrN(pyrdtc)2 (bottom). Color code: cyan for Cr, blue for
N, grey for C, red for O, orange for S, white for H. (B) Absorption spectra.
Simulated absorption peaks for CrN(pyrdtc)2 (left) and CrN(trop)2 (right).
CASSCF molecular orbitals for each state are also reported using an isovalue
of $0.02$ $e$/Å3. Positive (negative) values of the wavefunction are plotted
in purple (yellow). (C) EPR spectra. Simulated (solid line) and measured
(dashed line) cw-EPR spectum for CrN(pyrdtc)2 in green, CrN(trop)2 in blue,
and CrN(trop-d4)2 in orange.

However, the information that can be extracted from these experiments is only
as good as the theoretical framework used to interpret experimental evidence.
The canonical quantum theory of spin relaxation is based on a perturbative
treatment of how the thermal motion of atoms inside a crystalline lattice,
i.e. the phonons, affects the crystal field of unpaired electrons located in
the 3$d$ or 4$f$ shell of an ionvan1940paramagnetic . Following Orbach’s work
in the ’60sorbach1961spin , this complex quantum mechanical framework is
generally reformulated in terms of a phenomenological parametric description
of the sole ground-state multiplet, namely the $2S+1$ possible orientations
associated with a spin of magnitude $S$. This effective ground-state approach
to spin-lattice relaxation theory still represents the cornerstone of the
interpretation of experiments.

Despite its success, such a phenomenological theory of spin-lattice relaxation
can only partly assist the interpretation of relaxometry experiments due to
the fact that it is not a predictive tool and does not make it possible to
unequivocally determine the leading relaxation mechanism, nor the interaction
at the origin of relaxation itself. Only recently, these striking limitations
have been overcome by combining spin-lattice relaxation theory with ab initio
methodslunghi2023spin , where all the coefficients appearing in the equations
describing the dynamics of spin are fully determined from first principles.
Such an ab initio theory of spin relaxation has already provided a successful
prediction of both one- and two-phonon spin relaxation rates for a diverse
range of systems, going from coordination compounds of 3$d$ and 4$f$ ions with
long relaxation timeslunghi2022toward ; mondal2022unraveling ; reta2021ab ;
nabi2023accurate , known as single-molecule magnets, to solid-state quantum
sensors based on color centers, such as defects in diamond and hexagonal
boron-nitridemondal2023spin ; cambria2023temperature . A common denominator
for all these compounds is a spin moment larger than 1/2 and energy levels
defined by zero-field splitting.

Despite the momentous success of ab initio spin dynamics, the description of
spin relaxation in two-level spin-1/2 compoundslunghi2019phonons ;
lunghi2020limit , arguably the most fundamental type of system, has proved to
be the toughest challenge, and fundamental questions about its nature remain
unanswered. In particular, none of the fundamental interactions involving a
typical spin-1/2 system seems to be able to fully explain its dynamics at
temperatures above $\sim 20$ K (Raman relaxation), with the Zeeman and the
hyperfine interactions failing to explain i) the time-scale of relaxation and
the absence of a correlation with the external magnetic field
intensitylunghi2020limit ; lunghi2022toward , and ii) the angular dependence
of relaxation rateskazmierczak2022illuminating  and their correlation with
spin-orbit coupling strengthchakarawet2021effect , respectively. In the
attempt to reconcile some of these aspects of theory and experiments, the
involvement of electronic excited states has been
invokedkazmierczak2022illuminating , but a full ab initio assessment of their
role to Raman relaxation and a comparison with other proposed relaxation
mechanisms are still missing. Overall, such a state of affairs casts serious
doubts on our understanding of this fundamental physical process and must be
urgently addressed.

Here, we show that this conundrum can be solved by lifting the fundamental
assumption of Orbach’s canonical approach to spin-lattice relaxation theory,
namely by going beyond considering the sole spin-1/2 molecular ground state
and by considering virtual transitions to high-energy electronic excited. We
show that this novel theoretical framework can reproduce the relaxation times
of two prototypical spin-1/2 Cr(V) coordination compounds as measured by
electron paramagnetic resonance (EPR) inversion recovery experiments,
advancing the establishment of a complete ab initio theory of spin relaxation.
Moreover, we review recent results on other spin-1/2 coordination compounds
and critically discuss the need to include all possible mechanisms in a
comprehensive approach to spin-lattice relaxation theory.

## Results

The compounds nitrido bis(pyrrolidine dithiocarbamate)chromium(V)
(CrN(pyrdtc)2) and nitrido bis(tropolone)chromium(V) (CrN(trop)2) are chosen
for our investigation of the role of the various contributions to Raman
relaxation in spin-1/2. In this work we adopt a very inclusive definition of
Raman relaxation as any process that simultaneously involves two
phononslunghi2023spin . Among the many molecules with spin-1/2, Cr(V) systems
are particularly advantageous for our study because samples with natural
nuclear isotopic abundance are comprised of 90.5% 52Cr with nuclear spin
moment $I=0$ and 9.5% 53Cr with $I=3/2$. This permits the evaluation of the
impact of nuclear hyperfine interaction on relaxation in the same sample by
selectively measuring one or the other nuclear isotope. CrN(trop)2 and its
deuterated counterpart also allow the study of the effect of a rigid auxiliary
coordination sphere and nuclear spin-free ligands on spin relaxation time,
which have been flagged as possible sources of
decoherencezadrozny2015millisecond ; lunghi2020limit .

Synthesis, electronic structure and magnetism. The synthesis follows the
general protocol developed for nitrido-chromium(V) complexes, which involves
atom-transfer (formally transfer of N-) from N,N’-ethylene
bis(salicylideniminato)nitrido manganese(V), MnN(salen), to semi-labile
CrCl3(THF)3 in poorly coordinating solvents. Good yields over 70% are obtained
consistently for synthesis carried out in dry acetonitrile. The synthesis of
CrN(pyrdtc)2 was first reported in ref birk2003atom , whilst CrN(trop)2 is new
among tropolonate complexes. The solid-state structural parameters determined
from single crystal X-ray diffraction follow similar trends among nitrido
complexes. The molecular structures are depicted in Fig. 1A. The complex with
deuterated tropolonate (CrN(trop-d4)2), is also synthesized following the same
procedure. Full deuteration of all positions except the most distant
5-position in the tropolone ring is achieved as previously reported for other
derivativessakakibara2021direct . Finally, the compound ReN(pyrdtc)2,
isomorphous to the chromium complex, is synthesizedDT9720001079  to serve as a
diamagnetic solid-state host for CrN(pyrdtc)2. Dilution of paramagnetic
species into a diamagnetic host removes dipolar spin-spin interactions that
interfere with the measurement of spin-lattice relaxation.

The electronic excitations for CrN(pyrdtc)2 and CrN(trop)2 are computed with
multireference methods and reported in Fig. 1B together with their description
in terms of molecular orbitals, revealing a single unpaired electron located
in the $d_{xy}$ orbital of their Cr(V) ion. Excited states corresponding to
intra-band $d$ transition are present at high energy, starting from $\sim$
30,000 cm-1 for CrN(pyrdtc)2 and $\sim$ 20,000 cm-1 for CrN(trop)2, revealing
a well-isolated spin-1/2 Kramers doublet ground state. Simulations are
consistent with spectroscopic studiesbirk2003atom ;
bendix1998nitridocyanometalates , apart from an overall overestimation of the
transitions’ energies.

$\mathbf{A}$$\mathbf{B}$$\mathbf{C}$$\sim$$\sim$Energy$\Gamma_{I}$$|a\rangle$$|b\rangle$$\bar{n}_{\alpha}$,
$\bar{n}_{\beta}$$\bar{n}_{\alpha}$ \+ 1 $\bar{n}_{\beta}$ \+ 1
$\Gamma_{II}$$|a\rangle$$|b\rangle$$|c\rangle$$\bar{n}_{\alpha}$,
$\bar{n}_{\beta}$$\bar{n}_{\alpha}$ \+ 1 $\bar{n}_{\beta}$ \+ 1  Figure 2: (A)
Spin-lattice relaxation time. The relaxation times for CrN(pyrdtc)2 and
CrN(trop)2 are reported in the left and right panels, respectively. Red lines
and triangles correspond to experimental results obtained with inversion
recovery at the $g_{\perp}$ field and X-band frequency. Blue lines and dots
report simulation results obtained with $\Gamma_{II}$. Green lines and squares
report simulation results obtained with $\Gamma_{I}$. (B) $\mathbf{T_{1}}$
angular dependency. The relaxation time computed at 100 K as a function of the
direction of the applied magnetic field is reported. Triangles and circles
correspond to CrN(pyrdtc)2 and CrN(trop)2, respectively. Continuous red lines
show the fitting of the computed points using $f(\theta)=sin^{2}(\theta)$. (C)
Schematic representation of spin-phonon transitions. The left panel describes
a spin transition due to a resonant two-phonon process, where one phonon is
emitted and one is simultaneously absorbed. The right panel describes a spin
transition among the same, where a two-phonon process occurs thanks to virtual
transitions (red dashed lines) to an excited electronic state $|c\rangle$.

Thanks to the large energy separation between ground and excited states, these
compounds represent the prototypical systems where their magnetic properties
can be very accurately described by an effective spin-1/2 Hamiltonian, and all
excited states are assumed to play negligible roles up to room temperature.
The X-band continuous wave (CW) spectra of CrN(pyrdtc)2, CrN(trop)2 and
CrN(trop-d4)2, reported in Fig. 1C, were interpreted with a canonical spin
Hamiltonian

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

where the electron spin operator $\vec{\mathbf{S}}$ is coupled to the external
magnetic field $\vec{\mathbf{B}}$ and the 53Cr’s nuclear spin
$\vec{\mathbf{I}}$ through the Bohr’s magneton, $\mu_{\mathrm{B}}$, the Lande’
tensor, $\mathbf{g}$, and the hyperfine coupling tensor, $\mathbf{A}$,
respectively. The dominant features are from the $I=0$ 52Cr isotope which is
about 90.5% of natural abundance. Weaker, but well-defined lines are observed
for the $I=3/2$ 53Cr isotope. The eigenvalues of $\mathbf{g}$ and $\mathbf{A}$
obtained by simulating the EPR spectra are summarized in Table 1 and exhibit
characteristic axial anisotropy. Such symmetry axis roughly coincides with the
CrN chemical bond and will be referred to as the $g_{\parallel}$ direction.
The $g_{\perp}$ values for CrN(pyrdtc)2 and CrN(trop)2 are similar, but the
$g_{\parallel}$ value for CrN(trop)2 is significantly lower than for
CrN(pyrdtc)2. The comparison with values for other Cr(V) nitrido complexes
that have been reported previously in the literaturekonda1994electron ;
du1995electron  shows substantial similarities. The eigenvalues of
$\mathbf{g}$ and $\mathbf{A}$ for CrN(pyrdtc)2 in 3:7 CH2Cl2:toluene and doped
into the Re analog are the same within experimental uncertainty, demonstrating
that the two different environments do not significantly change the molecular
structure.

|  | X-Band | DFT | NEVPT2
---|---|---|---|---
CrN(pyrdtc)2 | $g_{iso}$ | 1.9875 | 1.984 | 1.983
$g_{\parallel}$ | 1.9673 | 1.969 | 1.955
$g_{\perp}$ | 1.9976 | 1.992 | 1.996
$A_{iso}$ | 74.20 | 79.31 | –
$A_{\parallel}$ | – | 125.71 | –
$A_{\perp}$ | 49.84 | 57.33 | –
CrN(trop)2 | $g_{iso}$ | 1.9678 | 1.976 | 1.971
$g_{\parallel}$ | 1.9415 | 1.953 | 1.925
$g_{\perp}$ | 1.9958 | 1.986 | 1.994
$A_{iso}$ | 74.76 | 82.61 | –
$A_{\parallel}$ | – | 128.66 | –
$A_{\perp}$ | 49.84 | 62.48 | –
Table 1: Spin Hamiltonian parameters. The table reports the eigenvalues of A
and g from Eq. 1 obtained by fitting the CW-EPR spectra and simulated with DFT
and NEVPT2. Hyperfine couplings are reported in MHz. Experimental values for
CrN(pyrdtc)2 and CrN(trop)2 were obtained at 60-70 K for the compound diluted
in the Re analog and in glassy toluene:CH2Cl2, respectively.

Table 1 also reports the eigenvalues of the tensors appearing in Eq. 1
calculated by two different electronic structure theory methods. As observed
in other studies, Density Functional Theory (DFT) methods are well suited to
recover electronic dynamical correlation in systems with a single unpaired $d$
electrons and achieve a good agreement for both $\mathbf{g}$ and A
tensorssingh2018challenges ; garlatti2023critical . NEVPT2 is found to
overestimate the $g_{\perp}$ shifts but still correctly capture the overall
experimental trends. The calculation of $\mathbf{A}$ at this level of theory
has not been attempted, but recent work suggests this is
possiblebirnoschi2024relativistic . In conclusion, both methods adequately
describe the electronic structure of these compounds.

Spin-lattice relaxation measurement. Spin-lattice relaxation for CrN(pyrdtc)2
and CrN(trop)2 are measured by 3-pulse inversion recovery as detailed in the
Materials and Methods section. In these experiments, the direction of the
electron spin is inverted with respect to the direction of the magnetic field,
and the time it takes to relax back to the thermal equilibrium, namely
$T_{1}$, is monitored. Emphasis is placed on data obtained at temperatures
above about 20 K, where the Raman process is proposed to
dominateeaton2000relaxation ; eaton2018relaxation . Results obtained at the
X-band frequency (see Fig. 2A) reveal a sharp decrease of $T_{1}$ as a
function of temperature, in agreement with the vast literature on spin-1/2
transition metal complexeskonda1994electron ; bader2014room ;
zadrozny2015millisecond ; atzori2016room . Faster relaxation for CrN(trop)2
than for CrN(pyrdtc)2 is consistent with the general trend that sees compounds
with larger $g$-shift (deviation from the free electron value of 2.0023) relax
faster due to a larger effective spin-orbit couplingsato2007impact . The
differences in $T_{1}$ between the two compounds can also be partly attributed
to the different host matrix, Re-doped solid and glassy solution,
respectively. The latter environment likely presents softer lattice
vibrations, which are known to speed up spin relaxation sato2007impact .
Values of $T_{1}$ for CrN(trop)2 in 3:7 CH2Cl2:toluene are indistinguishable
from CrN(trop-d4)2 in 3:7 CD2Cl2:toluene-d8 which shows that deuteration plays
a negligible role.

Experiments were repeated for different values of the external static magnetic
field, isolating the contribution of molecules with different orientations,
different field intensities, and different nuclear isotopes. The comparison
between experiments conducted at X-band and Q-band conclusively shows that the
overall intensity of the static magnetic field up to $\sim 1.2$ T does not
substantially affect $T_{1}$ for the same orientation between field and
molecule. Although it is not feasible to obtain data that define the full
orientation dependence of $T_{1}$ for randomly oriented molecules in the
glassy or powdered samples, it is possible to probe molecules oriented
parallel or perpendicular to the external field. For CrN(pyrdtc)2 in the Re
analog the ratio of $T_{1}$ for the parallel and perpendicular directions is
1.4 at X-band and 1.8 at Q-band, therefore pointing to longer values of
$T_{1}$ observed when the field is along $g_{\parallel}$. In glassy
toluene:CH2Cl2 the ratio is 1.2. For CrN(trop)2 the ratio of $T_{1}$
(parallel/perpendicular) is 1.6 at X-band and 1.5 at Q-band. Finally, the
comparison between $T_{1}$ for the two Cr’s isotopes does not reveal any
substantial difference, supporting the negligible role of hyperfine
interactions in the spin relaxation process for these complexes in rigid
matrices.

The experimental $T_{1}$ vs $T$ curves are interpreted with conventional
empirical models. According to the latter, above about 10 or 20 K the dominant
contributions for most spin-1/2 species are Raman processes due to a Debye-
like distribution of phonons plus a local-mode phonon with a discrete energy,
usually above the Debye cutoff temperatureeaton2000relaxation ;
eaton2018relaxation . The temperature dependence of $T_{1}$ for CrN(pyrdtc)2
in ReN(pyrdtc)2 and for CrN(trop)2 in 3:7 CH2Cl2 could be fit equally well
with two models: i) two local modes with energies of 55 and 200 cm-1 for
CrN(pyrdtc)2 and 42 and 145 cm-1 for CrN(trop)2, or ii) replacement of the
lower energy local modes with a Raman process due to a Debye distribution of
phonons with Debye temperatures of 70 and 55 cm-1 for CrN(pyrdtc)2 and
CrN(trop)2, respectively. Both models indicate that a substantial range of
phonon energies are required to match the data. The ability to fit the data
with two rather different models of the phonon distribution is a caution that
a match with the temperature dependence of the experimental data does not
constitute proof of the validity of the proposed phonon energy distribution
responsible for the relaxation. A comparison of the temperature dependence of
$T_{1}$ for CrN(pyrdtc)2 and CrN(trop)2 with other Cr(V) shows substantial
similarities in the temperature dependence of $T_{1}$ across compounds and
lattices suggest that the calculations and interpretations in this paper have
important implications for other Cr(V) complexes.

Spin-lattice relaxation theory and simulations. According to the canonical
picture of spin-lattice relaxation the interactions appearing in the spin
Hamiltonian of Eq. 1 are modulated by the activity of phonons,
$q_{\alpha}$orbach1961spin . Historically, this physical model was based on
the ideal picture of extended ionic crystal lattices described by a Debye
distribution. However, ab initio simulationslunghi2019phonons ;
lunghi2020multiple  and experimentsmoseley2018spin ; garlatti2020unveiling ;
chiesa2020understanding ; garlatti2023critical  have recently made it possible
to describe the real distribution of phonons in molecular crystals, showing
large deviations from this ideal picture. In particular, very low-energy
optical phonons are present at energies well below the common Debye
temperature of ionic compounds. The latter are characterized by a strong
admixing of intra-molecular and rotational components which can efficiently
modulate molecular spin-orbit interactionslunghi2019phonons , ultimately
making it possible for energy to flow between spin and lattice. This physical
picture of phonons in molecular crystals also helps to understand why the
relaxation rate of a molecule is not drastically affected by different
environments, e.g. crystals and frozen solutions, where optical phonons retain
much of the molecular normal modes of vibration and rotational motions are not
very environment-specific. This coupling between spin and phonons is formally
represented by the operators

$\hat{V}^{\alpha}=\left(\frac{\partial\hat{H}_{\mathrm{s}}}{\partial
q_{\alpha}}\right)\quad\text{and}\quad\hat{V}^{\alpha\beta}=\left(\frac{\partial^{2}\hat{H}_{\mathrm{s}}}{\partial
q_{\alpha}\partial q_{\beta}}\right)\>,$ (2)

where one or two phonons are modulating the spin’s properties, respectively.

In previous worklunghi2020limit ; lunghi2022toward ; garlatti2023critical ,
the second-order terms of Eq. 2 have been used to derive the Raman relaxation
rate as a function of temperature through the expression

$\Gamma_{I}=\frac{\pi}{2\hbar^{2}}\sum_{\alpha\geq\beta}|V^{\alpha\beta}_{ab}|^{2}\>G(\omega_{ab},\omega_{\alpha},\omega_{\beta})\>,$
(3)

where

$G(\omega_{ba},\omega_{\alpha},\omega_{\beta})=\delta(\omega_{ba}-\omega_{\alpha}+\omega_{\beta})\bar{n}_{\alpha}(\bar{n}_{\beta}+1)\>,$
(4)

and $V^{\alpha\beta}_{ab}$ is the matrix elements of the quadratic coupling
operator in Eq. 2 among the eigenstates of $\hat{H}_{\mathrm{s}}$. The latter
states approximately correspond to the parallel and antiparallel orientation
of the spin along the magnetic field direction. $E_{a}$ is the energy of said
spin states, $\hbar\omega_{ab}=E_{a}-E_{b}$, $\hbar\omega_{\alpha}$ is the
energy of the $\alpha$-phonon, and
$\bar{n}=[\mathrm{exp}(\hbar\omega/\mathrm{k_{B}}T)-1]^{-1}$ is the Bose-
Einstein thermal phonon population. For a two-level system, the relaxation
time is connected to Eq. 3 as $T_{1}=(2\Gamma_{I})^{-1}$. Fig. 2C provides a
graphical representation of this process. A transition between two spin
states, e.g. a spin flip, is induced by two simultaneous phonon processes,
where one is absorbed by the spin and the other emitted. To fulfil energy
conservation, introduced by the presence of a Dirac delta in Eq. 4, the energy
difference between the two phonons must match the spin energy splitting.

Eq. 3 is used to simulate $T_{1}$ as a function of temperature, field
intensity, and field orientation to compare with the experimental data from
inversion recovery experiments. When the separated contributions of A and g to
Eq. 2 are considered, the modulation of the former is predicted to be the
leading relaxation mechanism at all fields investigated in this study ($\sim$
0.33 and 1.2 T) and to lead to values of $T_{1}$ that are independent of the
magnitude of the external magnetic field. The modulation of g instead leads to
$T_{1}\propto B^{-2}$, which is inconsistent with the observed similarity of
$T_{1}$ values at X-band ($\sim$ 0.33 T) and Q-band (1.2 T). When studying the
dependence of $T_{1}$ on the magnetic field orientation, we observe that both
the modulation of A and g lead to orientation anisotropy. All these results
are in agreement with former reports for several Vanadyl compounds at both X-
and Q-band (vide infra Fig. 4)lunghi2020limit ; garlatti2023critical .
However, differently from previous studies, a relaxation mechanism due to the
modulation of A in these compounds is inconsistent with the similarity in
$T_{1}$ for the two nuclear isotopes of Cr ($I=0$ and $I=3/2$). Moreover, a
large deviation between predictions of $T_{1}$ and experiments is observed for
the Cr(V) compounds, with the predictions based on Eq. 3 overestimating the
$T_{1}$ by two orders of magnitude (see Fig. 2A). In conclusion, neither terms
of Eq. 1 is fully consistent with experimental data.

We shall now consider an alternative mechanism of Raman relaxation. When the
linear coupling of Eq. 2 is combined with fourth-order time-dependent
perturbation theorylunghi2022toward ; lunghi2023spin , Raman relaxation is
described by a rate population transfer among spin states that reads

$\displaystyle\Gamma_{II}=\frac{\pi}{2\hbar^{2}}\sum_{\alpha\geq\beta}\Big{|}T^{\alpha\beta,+}_{ab}+T^{\beta\alpha,-}_{ab}\Big{|}^{2}G(\omega_{ab},\omega_{\alpha},\omega_{\beta})\>,$
(5)

where

$T^{\alpha\beta,\pm}_{ab}=\sum_{c}\frac{V^{\alpha}_{ac}V^{\beta}_{cb}}{E_{c}-E_{b}\pm\hbar\omega_{\beta}}\>.$
(6)

All the terms in Eqs. 5 and 6 carry the same meaning as in Eq. 3, with the
additional presence of a sum over the excited states, $c$. The latter commonly
take the name of virtual states to stress that although they correspond to
real molecular electronic states, the process described by Eq. 5 does not
involve their population. This relaxation mechanism is graphically described
in Fig. 2C, where the same spin and phonon transition are concerned, but the
process is now mediated by virtual transitions.

Eq. 5 has been successfully used to simulate spin relaxation in molecules with
spin larger than 1/2lunghi2022toward ; mondal2022unraveling . However, for
two-level systems such as the Cr(V) compounds, the use of Eq. 5 is prevented
by the absence of any excited state $c$. It is then apparent that Eq. 5 can
only contribute to the relaxation rate of a spin-1/2 system if electronic
excited states beyond the spin Hamiltonian formalism are included. We
therefore apply Eq. 5 by reinterpreting the matrix elements of the first-order
derivatives of Eq. 2 as

$\hat{V}^{\alpha}_{ab}=\langle\varphi_{a}|\left(\frac{\partial\hat{H}_{\mathrm{el}}}{\partial
q_{\alpha}}\right)|\varphi_{b}\rangle\>,$ (7)

where the state vectors $\\{|\varphi_{a}\rangle\\}$ represent the eigenstates
of the full electronic Hamiltonian $\hat{H}_{\mathrm{el}}$, inclusive of spin-
orbit coupling. These include the spin-1/2 ground-state doublet and all
electronic excited states represented in Fig. 1B. These new matrix elements of
Eq. 7 are computed with a multireference ab initio method recently developed
by some of the authorsmariano2023manual .

$\mathbf{A}$$\mathbf{B}$ Figure 3: (A) High-energy phonon energy cutoff
dependence of relaxation time. The value of $T_{1}$ at 100 K is reported for
CrN(pyrdtc)2 (left) and CrN(trop)2 (right). $T_{1}$ is computed with
$\Gamma_{I}$ (green squares and line) and $\Gamma_{II}$ (blue squares and
line) by excluding phonons with $\omega>\Omega_{c}$. (B) Low-energy phonon
cutoff dependence of relaxation time. Phonon density of states are reported
for CrN(pyrdtc)2 (left) and CrN(trop)2 (right). The curves are shaded
according to the contribution to the computed relaxation time $T_{1}$ at 100
K, with exclusion of phonons with $\omega<\omega_{c}$, where the latter
corresponds to the value reported on the x-axis. Coloring of the phonon
distributions above the center line and reflected below the center line
correspond to the relaxation times obtained with $\Gamma_{II}$ and
$\Gamma_{I}$, respectively, and normalised to the value obtained including all
the phonons in the simulation.

As depicted in Fig. 2A, simulated values of $T_{1}$ using Eq. 5 and the
redefinition of Eq. 7 are in excellent agreement with experimental values for
both compounds, with the only notable slight difference being the slope of the
curves. A slightly larger deviation between simulations and experiments is
observed for CrN(trop)2, which we attribute to comparing simulations of
crystal phonons and measurements in frozen solution. Fig. 2B also reports the
value of $T_{1}$ as a function of the angular dependence of the external
magnetic field. We found excellent agreement between simulations and
experimental results for compound CrN(pyrdtc)2. In the case of CrN(trop)2, we
observe a reverse dependence, with longer $T_{1}$ for the perpendicular
orientation. However, we note that the orientation dependency of $T_{1}$ in
CrN(trop)2 is very small in both simulations and experiments, making a direct
comparison less straightforward. Simulations also correctly reproduce the
$\sin^{2}(\theta)$ orientational dependency of $T_{1}$ recently observed in
the literature for other spin-1/2 moleculeskazmierczak2022illuminating ;
kazmierczak2023t . Finally, this mechanism of Raman relaxation is found to be
independent of the magnitude of the external magnetic field and the nuclear
isotope, therefore fitting all experimental evidence.

Finally, we address the important question of identifying the relevant phonons
for spin relaxation. To determine the most relevant vibrational energy window,
the simulation of $T_{1}$ at 100 K is performed including all phonons up to a
cutoff $\Omega_{c}$. The full results of Fig. 2A are recovered for
$\Omega_{c}\rightarrow\infty$. Fig. 3A reports the results of this analysis
for both compounds clearly showing that predicted values of $T_{1}$ are
converged by $\Omega_{c}\sim 70$ cm-1 for both mechanisms considered
($\Gamma_{I}$ and $\Gamma_{II}$), meaning that phonons of higher energy do not
contribute to the observed values of relaxation time. This result is in
agreement with previous considerations based on calculations, where the
optical phonons with the lowest energy were identified as the dominant
contribution to relaxationlunghi2020limit ; lunghi2022toward ;
garlatti2023critical . This result is understood by noticing two key facts: i)
the thermal population of phonons, $\bar{n}$, decreases rapidly as a function
of their energy for a given $T$, making high-energy phonons less effective,
and ii) that a large density of states is observed in the low energy part of
the vibrational spectrum of the lattice (reported in Fig. 3B). We now turn to
the question of how much $T_{1}$ would increase if we removed the lowest
energy phonons, e.g. by synthetic design. To answer this question we perform
the simulation of $T_{1}$ at 100 K including all phonons down to a lower
cutoff $\omega_{c}$. When $\omega_{c}$ is set to a lower value than the first
available phonon, the full results of Fig. 2A are recovered. By shifting
$\omega_{c}$ to higher values, we exclude the low energy phonons from the
computation and observe how $T_{1}$ is affected. Fig. 3B reports the results
of this analysis for both compounds. Through the color scale, we can observe
that the value of $T_{1}$ predicted through $\Gamma_{I}$ increases drastically
as $\omega_{c}$ is increased. When we consider the contribution of phonons
through $\Gamma_{II}$ we observe a slightly different scenario, with a slower
increase of $T_{1}$ as $\omega_{c}$ is increased, effectively pointing to the
fact that even though the lowest energy phonons are still the ones that
determine $T_{1}$ through a winner-takes-it-all effect, a large window of
vibrations provides only slightly less effective alternative relaxation
pathways and that phonons up to $\sim$ 120 cm-1 would need to be removed to
extend $T_{1}$ by one order of magnitude 100 K. Interestingly, the thermal
population of phonons enters the expression of $\Gamma_{I}$ and $\Gamma_{II}$
in identical ways and the explanation for such a difference between the two
mechanisms must be found in the different pre-factors of the function $G$.
Indeed, the expression of $\Gamma_{II}$ shows that phonons at high energy are
favoured by minimizing the denominators of the functions
$T^{\alpha\beta,\pm}_{ab}$. Such condition is however weaker than the one
imposed by the phonon thermal population in virtue of the fact that the former
is a quadratic expression and the latter is exponential. As a result, phonons
with energy higher than the lowest ones offer more efficient alternative
relaxation pathways in $\Gamma_{II}$ than in $\Gamma_{I}$.

## Discussion

Despite their conceptual simplicity, the description of spin relaxation in
spin-1/2 molecules has proven to be the hardest nut to crack among a variety
of magnetic systemslunghi2023spin . Spin lattice relaxation times for a
variety of organic radicals and transition metal complexes have been observed
to decrease as the $g$-shift values increasesato2007impact ;
eaton2000relaxation ; eaton2018relaxation . Early modeling of $T_{1}$ and
based on a spin Hamiltonian ansatz elevated this experimental correlation to
causation and attempted to quantitatively predict spin relaxation in spin-1/2
in terms of the modulation of the $g$-tensoralbino2019first ;
mirzoyan2020dynamic . Both previouslunghi2019phonons ; lunghi2020limit ;
lunghi2022toward  and present computational results have however shown that
the modulation of the $g$-tensor cannot account for spin relaxation in
spin-1/2 molecules, due to its spurious field dependenceeaton2001frequency ;
atzori2016quantum  and leading to overall too long values of $T_{1}$. Without
abandoning the conceptual framework of the spin Hamiltonian, the hyperfine
interaction stood out as the only interaction able to explain spin relaxation
in spin-1/2 systems, and despite the independence of $T_{1}$ on Cr’s nuclear
spin is not consistent with this mechanism, the modulation of A has managed to
adequately explain all experimental evidence in the case of Vanadyl compounds
(see Fig. 4)lunghi2020limit ; lunghi2022toward ; garlatti2023critical

As we have shown here, the solution to this puzzle can be found by overcoming
the limits of Steven’s effective spin Hamiltonian formalism to interpret the
dynamical properties of spin systemselliott1952theory . The pervasive use of
this theoretical framework can be traced back to Orbach’s seminal paper on the
use of Stevenson effective operators to interpret spin-lattice relaxation in
rare-earth saltsorbach1961spin  and has become the fundamental language of
magnetic resonance ever since. The results obtained in this work for
CrN(pyrdtc)2 and CrN(trop)2 demonstrate that molecular electronic excited
states instead contribute to the Raman virtual transitions described by
$\Gamma_{II}$, dominating two-phonon spin relaxation at a temperature above
$\sim$ 20 K. Surprisingly, this is true despite the high energy of these
electronic excited states, a fact that has probably inhibited considering this
mechanism before. The theory proposed in this work naturally explains the
experimental correlation between the $g$-shift and relaxation time. Indeed,
spin-orbit coupling is at the same time at the origin of the $g$-shift and a
fundamental ingredient of $\hat{H}_{el}$ entering Eq. 7. An increased value of
such interaction automatically leads to a stronger coupling between lattice
and spin and to deviations to the free electron $g$-value. It is also
important to remark that $\Gamma_{I}$ and $\Gamma_{II}$ mechanisms are
operative at the same time. Depending on the relative strength of the crystal
field and hyperfine coupling, either mechanism could become the leading one
for a given molecule. To stress the importance of accounting for all the
relaxation mechanisms in the interpretation of relaxometry experiments, we
show in Fig. 4 the simulation of Raman relaxation time for the compound
VO(dmit)2atzori2016quantum  using both $\Gamma_{I}$ and $\Gamma_{II}$. Some of
the present authors previously simulated the spin-lattice relaxation time for
this molecule based on the $\Gamma_{I}$ mechanism and found excellent
agreement with experimentslunghi2022toward . The inclusion of $\Gamma_{II}$
contributions shows that the latter also accounts for similar relaxation
rates, making it hard to attribute spin relaxation to a single mechanism,
especially once numerical uncertainties are considered.

Figure 4: VO(dmit)2 Spin-lattice relaxation time. Red lines and triangles
correspond to experimental results obtained with inversion recovery at X-band
frequency. Blue lines and dots report simulation results obtained with
$\Gamma_{II}$. Green lines and squares report simulation results obtained with
$\Gamma_{I}$.

Our study also shines new light on the identification of the phonons
responsible for relaxation and highlights the limitations of using the common
local-mode model for the interpretation of $T_{1}$ v.s. $T$ curves. The
fitting of relaxation data often leads to the conclusion that modes with
energy above 100 cm-1 are necessary to explain spin-lattice relaxation. This
interpretation is in disagreement with ab initio simulations of $T_{1}$, which
instead points to low energy modes of tens of cm-1 as the sole relevant. Our
study confirms the latter interpretation, with the caveat that a larger window
of low energy vibrations than previously noted should be removed to sensibly
extend $T_{1}$. Since both ab initio methods and the phenomenological local-
mode model descend from the same underlying assumptionslunghi2023spin , e.g.
weak coupling, Markovian dynamics, harmonic phonons, etc., the disagreement
between the two suggests that conclusions based on the latter must be treated
with considerable care and that the proposed higher energies for a local mode
currently lacks an explanation. We also note that recent attempts to
complement local-mode models with computed vibrational spectra have led to the
reinforcement that high-energy modes are responsible for relaxation in
spin-1/2 moleculeskazmierczak2021impact ; kazmierczak2024determining .
However, such attempts include adjustable coefficients and are invariably
applied to gas-phase normal modes, which arguably differ from the real crystal
vibrational density of states and provide a poor description of low-energy
vibrationslunghi2019phonons . Indeed, a very recent experimental study has
shown that the modulation of the low-energy vibrational spectrum of vanadyl-
based phthalocyanines deeply affects their $T_{1}$lohaus2024quantifying , as
predicted by previous simulations accounting for the full crystal vibrational
spectrum garlatti2023critical . At the same time, we acknowledge that although
overall trends and relaxation rates are well reproduced, the calculated slope
of $T_{1}$ v.s. $T$ deviates from the experimental one. Whilst discrepancies
at low $T$ can be easily interpreted as the lack of acoustic phonons in the
simulations, the high-$T$ behaviour cannot be easily reconciled with the
current theoretical framework. The latter inevitably leads to the high-$T$
limit of $T_{1}\propto T^{-2}$, thus underestimating the experimental $T_{1}$
in the temperature range in which empirical fittings invoke a high-energy
local mode. We perform numerical tests where we artificially increase the
coupling to spin of high-energy phonons in the range individuated by the
local-mode analysis and found that at least a factor four of increase in the
coupling would be necessary to have a visible effect of the predicted $T_{1}$.
Although a systematic investigation of the accuracy of electronic structure
methods has not been performed yet, successful applications of our method to
other compoundslunghi2022toward ; mondal2022unraveling  make it improbable
that numerical inaccuracies could lead to such a drastically different
qualitative picture of relaxation. On the other hand, the accuracy of our
methods in the prediction of low-energy vibrations in crystals has been
thoroughly benchmarked against experimentsgarlatti2020unveiling ;
garlatti2023critical . As such, we envision that the sophistication of the
underlying open quantum systems methods is a priority for further improving
the agreement between theory and experiments and conclusively assessing the
role of high energy phonons in the spin relaxation of spin-1/2 systems.
Interestingly, as this manuscript was being prepared, Shushkov reported a DFT-
based computational method to predict the relaxation rate of Cu(II) spin-1/2
moleculesShushkov2024novel . Differently from the method presented here,
Sushkov described an additional contribution to $\Gamma_{I}$ due to non-
adiabatic effects. The latter contribution was found to dominate $T_{1}$ over
the entire temperature range with respect to the modulation of g and to be
driven by high-energy gas-phase molecular modes. Although this study also
potentially suffers from a poor representation of low energy modes, the
proposed relaxation mechanism stands out as an additional competing effect to
those accounted for by $\Gamma_{II}$. To the best of our understanding, the
same non-adiabatic effects would be fully included in $\Gamma_{I}$ by
extending Eq. 7 to the second order, making the possibility to compare all
these methods within the same ab initio framework within
reachmariano2023manual .

Thanks to these new results on the interpretation of Raman relaxation, we are
now in the position to critically review the list of chemical strategies at
our disposal to control $T_{1}$. Indeed, achieving a large separation among
ground and excited states as well as reducing the density of states in the
low-energy part of the vibrational spectrum point to longer values of $T_{1}$.
The different relaxation rates of CrN(pyrdtc)2 and CrN(trop)2 can indeed be
interpreted following these considerations, with both lower excited states and
lower energy phonons for the tropolonate. The overall remarkable properties of
these square pyramidal complexes can also now be easily interpreted through
the presence of the strongly $\sigma$ and $\pi$ donor nitrido group, which
induces a very large splitting of $d$ orbitals . In light of these
considerations, one would be tempted to employ heavier elements like 4$d$ or
5$d$ transition metals, which boast much larger crystal field splitting
lever1968inorganic . This however would also lead to an increased value of
effective spin-orbit coupling, which is well documented as being detrimental
for $T_{1}$chakarawet2021effect . Interestingly, despite having provided a
strong reinterpretation of Raman relaxation mechanism in spin-1/2 and the
justification for its chemical control guidelines, the latter does not
drastically shift from those previously inferred from considering the
modulation of the $g$-tensor as the leading relaxation
mechanismalbino2019first ; mirzoyan2020dynamic . Incidentally, the $g$-shift
and the residual orbital angular momenta of a molecule, often claimed to
determine $T_{1}$albino2019first ; ariciu2019engineering ; mirzoyan2020dynamic
, also depend on the energy of excited states and have thus suggested large
crystal field splitting as a design rule towards long-lived spin excitations.
Similarly, for what concerns the phonons, both Raman mechanisms receive the
largest contribution from low-energy vibrations, making previous advice still
validgarlatti2023critical , albeit with reduced expected effectiveness, as
suggested by Fig. 3.

In conclusion, our theoretical and experimental study has demonstrated that
the interpretation of spin relaxation in spin-1/2 compounds requires going
beyond considering only the ground state of molecules and including high-
energy electronic virtual excitations. The inclusion of these contributions
not only made it possible to reproduce experiments to an unprecedented level
of accuracy but also provided novel conceptual tools to study and chemically
engineer this important class of compounds. We envision that further
development of ab initio open quantum system theory in synergy with novel
experiments will soon make it possible to achieve a long-sought full
reconciliation between theory and experiments.

## Materials and Methods

Syntheses and samples preparation. All solvents were HPLC grade from VWR
chemicals, except for dichloromethane (DCM) which was from Fischer Chemical.
Mn(N)Salen, and mer-Cr(Cl)3(THF)3 were synthesized as previously
reportedbirk2003atom . Tropolone (98%) was obtained commercially from Combi
Blocks, and recrystallized from isopropanol before use. CF3COD (Sigma
Aldrich), platinum on carbon, 10% (Sigma Aldrich), and ammonium pyrrolidine
dithiocarbamate (Acros) were used as received. CrN(pyrdtc)2 and ReNCl2(PPh3)2
were synthesized as described in the literaturebirk2003atom ; DT9720001079 .

Synthesis of tropolone-d4: The deuteration of tropolone followed literature
protocols for aromatic compounds. Tropolone (1 g; 8.19 mmol) and 10% Pt/C (0.2
g) were added to deuterium oxide (15 g; 75 mmol)sakakibara2021direct . The
reaction was kept acidic by addition of two drops of d-trifluoroacetic acid to
subdue rearrangement of the tropolonate under the forcing conditions. The
reaction was done in a Teflon lined steel autoclave at 145 C∘ for three days.
Aqueous workup and recrystallization from isopropanol yielded 0.390 g (37.8%)
of Tropolone-d4. The deuteration degree was inferred to be 80% from the HRMS
of the resulting Cr(N)(trop-d4)2 in agreement with the previously reported
difference in reactivity of the 5-position compared to the remaining ring
positionssakakibara2021direct .

Synthesis of Cr(N)(trop)2 and Cr(N)(trop-d4)2: The procedure follows the
general procedure described by one of the authorsbirk2003atom . mer-
Cr(Cl)3(THF)3 (0.153 g; 0.409 mmol) and Mn(N)(salen) (0.275 g; 0.821 mmol)
were mixed in acetonitrile (10 mL) under a nitrogen atmosphere. The reaction
was stirred at RT for 20 min. and the precipitated Mn(Cl)(salen) was filtered
off. To the filtrate, tropolone (0.105 g; 0.859 mmol) in acetonitrile (3 mL)
was added dropwise, while stirring under under N2 upon formation of a brownish
orange precipitate. After completion of the precipitate, the product was
filtered off and washed with acetonitrile (2 times 30 mL) before being dried
in a stream of N2. The yield was 89.5 mg (71.1%) of a brownish orange
microcrystalline product. The polycrystalline sample was void of Cr(trop)3
based on the comparison with the PXRD (Cu-K$\alpha$), of the latter. Crystals
suitable for SC-XRD were obtained by recrystallization from acetonitrile/DCM
(1:2). The deuterated version Cr(N)(trop-d4)2 was synthesized by the same
protocol yielding 70.5% of isolated product with similar peak positions in the
PXRD as the non-deuterated sample. HRMS, m/z=317.0604 (100%), 316.0540 (92%)
corresponding to $>$ 80% deuteration.

Synthesis of ReN(pyrdtc)2, Nitridobis (pyrrolidinedithiocarbamato)rhenium(V):
ReNCl2(PPh3)2 (0.60 g, 0.75 mmol ) and NH4(pyrdtc) (0.50 g, 3.0 mmol) were
stirred in acetone (60 ml) overnight at 40∘C under N2. The brownish solution
containing a yellow precipitate was evaporated to dryness in a stream of N2.
The resulting solid was extracted with DCM/MeOH (9:1, 100 ml), filtered, and
the extract was reduced in volume to ca. 10 ml yielding a bright yellow
crystalline product, which was washed with ether, and dried in vacuum (0.24 g,
48%). Crystals suitable for single crystal X-ray diffraction were obtained by
slow evaporation of a DCM solution.

X-ray structure determination: Single crystals were obtained as described
under the synthesis. They were coated with mineral oil, mounted on kapton
loops, and transferred to the nitrogen cold stream of the diffractometer.
Diffraction data were recorded at 100(2) K on a Bruker D8 VENTURE
diffractometer equipped with a Mo K$\alpha$ high-brilliance I$\mu$S radiation
source ($\lambda$ = 0.71073 Å), a multilayer X-ray mirror and a PHOTON 100
CMOS detector, and an Oxford Cryosystems low-temperature device. The
instrument was controlled with the APEX3 software package using SAINT. Final
cell constants were obtained from least squares fits of several thousand
strong reflections. Intensity data were corrected for absorption using
intensities of redundant reflections with the program SADABS. The structures
were solved in Olex2 using SHELXT and refined using SHELXL.

Electron Paramagnetic Resonance Spectroscopy. Solutions of 0.6 mM
Cr(N)(pyrdtc)2 and 0.5 mM CrN(trop)2 were prepared by dissolution in CH2Cl2
followed by dilution with toluene to achieve a 3:7 CH2Cl2:toluene ratio. A
solution of 0.5 mM CrN(trop-d4)2 was prepared analogously in 3:7
CD2Cl2:toluene-d8. Approximately 200 $\mu$L of solution or the doped $\sim$
0.05% Cr(N)(pyrdtc)2 in ReN(pyrdtc)2 solid was placed in 4 mm OD quartz tubes.
Solutions were deoxygenated by multiple freeze-pump-thaw cycles. All EPR
samples were backfilled with $\sim$ 200 mtorr of helium prior to flame sealing
the tube. The low-pressure He provides thermal contact between the sample and
the walls of the tube to facilitate thermal equilibration.

EPR spectroscopy and electron spin relaxation time measurements were performed
on a Bruker E580 spectrometer with an Oxford ESR935 cryostat. X-band data were
obtained with an ER4118X-MS5 split ring resonator and Q-band data were
obtained with an ER5107 resonator. For pulse measurements the X-band resonator
was over-coupled to a ringdown 1/e time of $\sim$ 2.5 – 3 ns which corresponds
to $Q\sim 150$. The Q-band resonator was overcoupled to permit signal
acquisition with less than 200 ns deadtime. Temperature was controlled with a
ColdEdge/Bruker Stinger closed-cycle helium cooling system. Temperature was
monitored with a calibrated Lakeshore Cernox sensor and an Oxford Mercury iTC
readout. Continuous wave (CW) spectra were measured at 60 or 70 K using low
microwave power and 0.5 G modulation amplitude at 10 kHz to minimize passage
effects. The spectra were simulated using the Bruker software Anisospinfit.
$T_{1}$ was measured by 3-pulse inversion recovery with sequence
$\pi-t-\pi/2-\tau-\pi-\tau-echo$, with 2-step phase cycling, and the value of
$t$ was incremented to record the recovery curve. The constant $\tau$ was set
to 360 or 380 ns. The shot repetition time (SRT) was about 10 times the value
of $T_{1}$ calculated with a stretched exponential. For the relatively high
spin concentration in the 0.05% CrN(pyrdtc)2 in the Re analog the time
constants were dependent on lengths of the $\pi/2$ pulse lengths between 16
and 160 ns. To minimize the impact of spectral diffusion on $T_{1}$ the
inversion recovery data for 0.05% CrN(pyrdtc)2 were acquired with $\pi/2$
pulse lengths of 16 ns. For 0.6 mM Cr(N)(pyrdtc)2 in 3:7 CH2Cl2:toluene the
values of $T_{1}$ were independent of pulse length, but for consistency with
the doped solid the same pulse lengths were used. For 0.5 mM CrN(trop)2 and
CrN(trop-d4)2 in 3:7 CH2Cl2:toluene or 3:7 CD2Cl2:toluene-d8 $\pi/2$ pulse
lengths of 40 ns were used. Inversion recovery curves were not single
exponentials, indicating the presence of distributions of relaxation times.
Inversion recovery data were fit with a stretched exponential

$Y(t)=e^{-(t/T_{1})^{\beta}}\>,$ (8)

and values of $\beta$ ranged from about 0.65 at 10 K to about 0.95 at 70 K.
All fitted values are reported in the Supporting Information Tables S2 and S3.

Electronic structure simulations. Starting from X-ray structures of
CrN(pyrdtc)2 and CrN(trop)2, cell and geometry optimization and simulations of
$\Gamma$-point phonons have been performed with periodic density functional
theory (pDFT) using the software CP2Kkuhne2020cp2k . Cell optimization was
performed employing a very tight force convergence criteria of 10-7 a.u. and
SCF convergence criteria of 10-10 a.u. for the energy. A plane wave cutoff of
1000 Ry, DZVP-MOLOPT Gaussian basis sets, and Goedecker-Tetter Hutter
pseudopotentials were employed for all atoms. The Perdew-Burke-Ernzerhof (PBE)
functionalperdew1996generalized  and DFT-D3 dispersion
correctionsgrimme2010consistent  were used.

Electronic structure and magnetic properties were computed on the pDFT-
optimized structures using the software ORCAneese2020orca . The decontracted
basis set DKH-def2-TZVPP, along with the Douglas-Kroll-Hess (DKH) scalar
relativistic correction to the electronic Hamiltonian, was employed. Picture-
change effects were included in the calculations. The computation of
$\mathbf{g}$ and $\mathbf{A}$ tensors at the DFT level of theory was performed
using a PBE0 exchange-correlation functionaladamo1999toward  with the RIJCOSX
approximation for Coulomb and Hartree-Fock exchange. Multireference
calculations were performed using Complete Active Space Self Consistent Field
(CASSCF) in combination with N-Electron Valence Perturbation Theory (NEVPT2).
The active space used to build the CASSCF wavefunction is (1,5), i.e., one
electron in five 3$d$-orbitals. All possible states with multiplicity two were
included in the state-average procedure. Mean-field spin-orbit coupling
operator, along with Quasi Degenerate Perturbation Theory (QDPT), was employed
to account for the mixing of spin-free states.

Spin-phonon relaxation simulations. First-order vibronic coupling matrix
elements of Eq. 7 have been evaluated within the framework described
beforemariano2023manual . Briefly, starting from the ab initio CASSCF
wavefunctions $|\varphi_{a}\rangle$ and energies $E_{a}$, the matrix elements
of the $\nabla\hat{H}_{el}$ operator are expressed using the Hellmann-Feynman
theorem as

$\langle\varphi_{a}|\nabla\hat{H}_{el}|\varphi_{b}\rangle=\nabla
E_{b}\delta_{ab}+(E_{b}-E_{a})\langle\varphi_{a}|\nabla\varphi_{b}\rangle\>.$
(9)

To evaluate the non-adiabatic coupling (NAC) terms
$\langle\varphi_{a}|\nabla\varphi_{b}\rangle$, numerical differentiation of
wavefunction overlap has been performed using an in-house Python code
interfaced with the program WFOVERLAPplasser2016efficient . The step for the
numerical differentiation has been set to 0.001 Å along the Cartesian
coordinates of the system. A systematic convergence study of the $T_{1}$
relaxation time as a function of the differentiation step was performed. Once
the vibronic coupling matrix elements have been computed in the Cartesian
coordinates $\vec{R}$, they are transformed into the normal mode reference
framework by using

$\left(\frac{\partial\hat{H}_{\mathrm{el}}}{\partial
Q_{{\alpha}}}\right)=\sum_{a}^{3N}\sqrt{\frac{\hbar}{\omega_{\alpha}m_{a}}}L_{\alpha
a}\left(\frac{\partial\hat{H}_{\mathrm{el}}}{\partial R_{a}}\right)\>,$ (10)

where $m_{i}$ is the $i$-th atomic mass and $L_{\alpha i}$ the Hessian matrix
eigenvectors. The latter is computed with a two-step numerical differentiation
of forces. The vibronic coupling vectors from CASSCF, together with NEVPT2
energies in Eq. 6, are finally employed to simulate spin relaxation using Eq.
5 with the software MolForgelunghi2022toward . The use of CASSCF vibronic
coupling, in conjunction with NEVPT2 vertical excitation energies, is
necessary here, as orbital relaxation at the NEVPT2 level is not supported in
ORCA. The Dirac delta functions appearing in the spin dynamics equations are
broadened with a Gaussian smearing of 30 cm-1 after a convergence study on the
computed $T_{1}$.

The evaluation of the second-order terms of Eq. 2 requires computing the
second-order derivatives of the tensors $\mathbf{A}$ and $\mathbf{g}$, namely
$(\partial^{2}\mathbf{A}/\partial q_{\alpha}\partial q_{\beta})$ and
$(\partial^{2}\mathbf{g}/\partial q_{\alpha}\partial q_{\beta})$. In this
study, we employed the machine-learning-based strategy to calculate the
second-order derivatives of the spin Hamiltonian in Eq. 1lunghi2020limit ;
garlatti2023critical . These coefficients are subsequently employed to
simulate the time-dependency of the $z$-component of the simulated
magnetization $\vec{\mathbf{M}}$. Finally, $T_{1}$ is extracted by fitting
$M_{z}$ with a double exponential function.

Calculations of relaxation rates and magnetization time evolution were
obtained with the development branch T4 of the software
MolForgelunghi2022toward , available at github.com/LunghiGroup/MolForge.

Acknowledgements and Funding
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).

Authors Contributions
A.L., S.S.E., and G.R.E. conceived the project. L.M. and A.L. developed the
theoretical and computational models. A.L., L.M., and V.H.A.N carried out the
calculations and their analysis. G.R.E. and S.S.E. performed the EPR
measurements. J.B., J.B.P. and M.B. synthesized the compounds. All authors
contributed to the discussion of results, and to writing the manuscript.

Conflict of interests The authors declare no competing interests.

## References

  * (1) Eaton, S. S. & Eaton, G. R.  Relaxation times of organic radicals and transition metal ions.  In Berliner, L. J., Eaton, S. S. & Eaton, G. R. (eds.) _Distance measurements in biological systems by EPR_ , 29–154 (Springer US, 2000).
  * (2) Parigi, G., Ravera, E., Fragai, M. & Luchinat, C.  Unveiling protein dynamics in solution with field-cycling nmr relaxometry.  _Progress in Nuclear Magnetic Resonance Spectroscopy_ 124, 85–98 (2021).
  * (3) Li, C., Chen, M., Lyzwa, D. & Cappellaro, P.  All-optical quantum sensing of rotational brownian motion of magnetic molecules.  _Nano Letters_ 19, 7342–7348 (2019).
  * (4) Alzola-Aldamizetxebarria, S., Fernández-Méndez, L., Padro, D., Ruíz-Cabello, J. & Ramos-Cabrer, P.  A comprehensive introduction to magnetic resonance imaging relaxometry and contrast agents.  _ACS Omega_ 7, 36905–36917 (2022).
  * (5) Degen, C. L., Reinhard, F. & Cappellaro, P.  Quantum sensing.  _Reviews of Modern Physics_ 89, 035002 (2017).
  * (6) Mzyk, A., Sigaeva, A. & Schirhagl, R.  Relaxometry with nitrogen vacancy (NV) centers in diamond.  _Accounts of Chemical Research_ 55, 3572–3580 (2022).
  * (7) Van Vleck, J.  Paramagnetic relaxation times for titanium and chrome alum.  _Physical Review_ 57, 426 (1940).
  * (8) 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).
  * (9) Lunghi, A.  Spin-phonon relaxation in magnetic molecules: Theory, predictions and insights.  In Rajaraman, G. (ed.) _Computational Modelling of Molecular Nanomagnets_ , 219–289 (Springer, 2023).
  * (10) Lunghi, A.  Toward exact predictions of spin-phonon relaxation times: An ab initio implementation of open quantum systems theory.  _Science Advances_ 8, eabn7880 (2022).
  * (11) Mondal, S. & Lunghi, A.  Unraveling the contributions to spin–lattice relaxation in kramers single-molecule magnets.  _Journal of the American Chemical Society_ 144, 22965–22975 (2022).
  * (12) 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).
  * (13) Nabi, R. _et al._ Accurate and efficient spin–phonon coupling and spin dynamics calculations for molecular solids.  _Journal of the American Chemical Society_ 145, 24558–24567 (2023).
  * (14) Mondal, S. & Lunghi, A.  Spin-phonon decoherence in solid-state paramagnetic defects from first principles.  _npj Computational Materials_ 9, 120 (2023).
  * (15) Cambria, M. _et al._ Temperature-dependent spin-lattice relaxation of the nitrogen-vacancy spin triplet in diamond.  _Physical Review Letters_ 130, 256903 (2023).
  * (16) Lunghi, A. & Sanvito, S.  How do phonons relax molecular spins?  _Science Advances_ 5, eaax7163 (2019).
  * (17) Lunghi, A. & Sanvito, S.  The limit of spin lifetime in solid-state electronic spins.  _The Journal of Physical Chemistry Letters_ 11, 6273–6278 (2020).
  * (18) Kazmierczak, N. P. & Hadt, R. G.  Illuminating ligand field contributions to molecular qubit spin relaxation via T1 anisotropy.  _Journal of the American Chemical Society_ 144, 20804–20814 (2022).
  * (19) 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).
  * (20) 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).
  * (21) Birk, T. & Bendix, J.  Atom transfer as a preparative tool in coordination chemistry. synthesis and characterization of Cr (V) nitrido complexes of bidentate ligands.  _Inorganic Chemistry_ 42, 7608–7615 (2003).
  * (22) Sakakibara, Y., Osada, K., Uraki, Y., Ubukata, M. & Shigetomi, K.  Direct deuteration of hinokitiol and its mechanistic study.  _Bioscience, Biotechnology, and Biochemistry_ 85, 215–222 (2021).
  * (23) Fletcher, S. R. & Skapski, A. C.  Crystal and molecular structure of the square-pyramidal complex nitridobis(NN-diethyldithiocarbamato)rhenium(V).  _Dalton Transactions_ 10, 1079–1082 (1972).
  * (24) Bendix, J. _et al._ Nitridocyanometalates of CrV, MnV, and MnVI.  _Inorganic Chemistry_ 37, 1767–1775 (1998).
  * (25) Konda, R., Du, J. L., Eaton, S. & Eaton, G.  Electron spin relaxation rates for nitridochromium (V) tetratolylporphyrin and nitridochromium (V) octaethylporphyrin in frozen solution.  _Applied Magnetic Resonance_ 7, 185–193 (1994).
  * (26) Du, J.-L., Eaton, G. & Eaton, S.  Electron-spin-lattice relaxation in natural abundance and isotopically enriched oxo-chromium (V) bis (2-hydroxy-2-ethylbutyrate).  _Journal of Magnetic Resonance, Series A_ 115, 236–240 (1995).
  * (27) 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).
  * (28) Garlatti, E. _et al._ The critical role of ultra-low-energy vibrations in the relaxation dynamics of molecular qubits.  _Nature Communications_ 14, 1653 (2023).
  * (29) Birnoschi, L., Oakley, M. S., McInnes, E. J. & Chilton, N. F.  Relativistic quantum chemical investigation of actinide covalency measured by electron paramagnetic resonance spectroscopy.  _Journal of the American Chemical Society_ 146, 14660–14671 (2024).
  * (30) Eaton, S. S. & Eaton, G. R.  Relaxation mechanisms.  In Goldfarb, D. & Stoll, S. (eds.) _EPR Spectroscopy: Fundamentals and Methods_ , 175–192 (John Wiley & Sons, 2018).
  * (31) Bader, K. _et al._ Room temperature quantum coherence in a potential molecular qubit.  _Nature Communications_ 5, 5304 (2014).
  * (32) 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).
  * (33) 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).
  * (34) Lunghi, A. & Sanvito, S.  Multiple spin–phonon relaxation pathways in a kramer single-ion magnet.  _The Journal of Chemical Physics_ 153, 174113 (2020).
  * (35) Moseley, D. H. _et al._ Spin–phonon couplings in transition metal complexes with slow magnetic relaxation.  _Nature Communications_ 9, 2572 (2018).
  * (36) 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).
  * (37) Chiesa, A. _et al._ Understanding magnetic relaxation in single-ion magnets with high blocking temperature.  _Physical Review B_ 101, 174402 (2020).
  * (38) Mariano, L. A., Mondal, S. & Lunghi, A.  Spin-vibronic dynamics in open-shell systems beyond the spin hamiltonian formalism.  _Journal of Chemical Theory and Computation_ 20, 323–332 (2024).
  * (39) Kazmierczak, N. P., Luedecke, K. M., Gallmeier, E. T. & Hadt, R. G.  T1 anisotropy elucidates spin relaxation mechanisms in an S= 1 Cr (IV) optically addressable molecular qubit.  _The Journal of Physical Chemistry Letters_ 14, 7658–7664 (2023).
  * (40) Albino, A. _et al._ First-principles investigation of spin–phonon coupling in vanadium-based molecular spin quantum bits.  _Inorganic Chemistry_ 58, 10260–10268 (2019).
  * (41) 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).
  * (42) Eaton, S., Harbridge, J., Rinard, G., Eaton, G. & Weber, R.  Frequency dependence of electron spin relaxation for three s= 1/2 species doped into diamagnetic solid hosts.  _Applied Magnetic Resonance_ 20, 151–157 (2001).
  * (43) 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).
  * (44) Elliott, R. J. & Stevens, K.  The theory of the magnetic properties of rare earth salts: cerium ethyl sulphate.  _Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences_ 215, 437–453 (1952).
  * (45) 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).
  * (46) Kazmierczak, N. P., Lopez, N. E., Luedecke, K. M. & Hadt, R. G.  Determining the key vibrations for spin relaxation in ruffled Cu (II) porphyrins via resonance raman spectroscopy.  _Chemical Science_ 15, 2380–2390 (2024).
  * (47) Lohaus, S. H. _et al._ Quantifying the effects of peripheral substituents on the spin-lattice relaxation of a vanadyl molecular quantum bit.  _Journal of Porphyrins and Phthalocyanines_ (2024).  URL https://doi.org/10.1142/S1088424624500329.
  * (48) Shushkov, P.  A novel non-adiabatic spin relaxation mechanism in molecular qubits.  _The Journal of Chemical Physics_ 160, 164105 (2024).
  * (49) Lever, A. B. P.  _Inorganic electronic spectroscopy_ (Elsevier Publishing, 1968).
  * (50) 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).
  * (51) 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).
  * (52) Perdew, J. P., Burke, K. & Ernzerhof, M.  Generalized gradient approximation made simple.  _Physical Review Letters_ 77, 3865 (1996).
  * (53) 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).
  * (54) Neese, F., Wennmohs, F., Becker, U. & Riplinger, C.  The orca quantum chemistry program package.  _The Journal of Chemical Physics_ 152, 224108 (2020).
  * (55) Adamo, C. & Barone, V.  Toward reliable density functional methods without adjustable parameters: The PBE0 model.  _The Journal of chemical physics_ 110, 6158–6170 (1999).
  * (56) Plasser, F. _et al._ Efficient and flexible computation of many-electron wave function overlaps.  _Journal of Chemical Theory and Computation_ 12, 1207–1219 (2016).