# Achieving Strong Magnon Blockade through Magnon Squeezing in a Cavity
Magnetomechanical System

M. Amazioug  LPTHE, Department of Physics, Faculty of sciences, Ibnou Zohr
University, Agadir, Morocco.  D. S. Dutykh  Department of Mathematics, Khalifa
University, Abu Dhabi 127788, United Arab Emirates.  B. Teklu  Department of
Mathematics, Khalifa University, Abu Dhabi 127788, United Arab Emirates.
Center for Cyber-Physical Systems (C2PS), Khalifa University, 127788 Abu
Dhabi, UAE.  M. Asjad<EMAIL_ADDRESS>Department of Mathematics, Khalifa
University, Abu Dhabi 127788, United Arab Emirates.

###### Abstract

We propose a scheme to achieve magnon (photon) blockade by using magnon
squeezing within a cavity magnomechanical system under weak pump driving.
Under ideal conditions, we observe a substantial magnon blockade effect, as
well as simultaneous photon blockade. Moreover, both numerical and analytical
results match perfectly, providing robust evidence of consistency. In addition
to calculating optimal parametric gain and detuning values, we can improve the
second-order correlation function. The proposed scheme will be a pioneering
approach towards magnon (photon) blockade in experimental cavity
magnomechanical systems.

##  I Introduction

Photon statistics are integral to quantum optics [1]. Initially unveiled by
Hanbury Brown [2, 3] and confirmed by Kimble et al. [4], photon antibunching
illustrates light’s intriguing properties. An enthralling phenomenon, photon
blockade (PB), arises within nonlinear media, where the presence of a single
photon obstructs the entry of another photon [5]. PB’s theoretical relevance
extends to diverse nonlinear optical systems, encompassing Kerr-type nonlinear
cavities [6] and optomechanical setups [7, 8, 9, 10, 11, 12, 13, 14].
Experimental manifestations of PB materialize in contexts like an optical
cavity coupled to a confined atom [15] and a photonic crystal cavity linked
with a quantum dot [16, 17, 18]. The ability to generate single-photon sources
is of paramount importance for quantum information technologies [19, 20]. The
photon blockade effect plays a crucial role in achieving this objective as it
produces strongly anti-bunched photons. This effect enables the generation of
sub-Poissonian light even when the system is driven by a classical light
field.

In recent years, the photon blockade effect has gained significant importance
in quantum information processing and quantum communication. Currently, there
are two well-known mechanisms for achieving a strong photon blockade effect:
unconventional photon blockade (UCPB) [21, 22, 23], which relies on
destructive quantum interference between two different quantum transition
paths from the ground state to a two-excited state, and conventional photon
blockade (CPB) [5, 24, 25] which relies on large nonlinearities to alter the
energy-level structure of the system. Experimental implementations of both
UCPB and CPB mechanisms have been reported in [15, 16, 26, 27, 28].

Cavity magno-mechanics has emerged as a promising research area, providing a
solid platform for studying ferrimagnetic crystals such as YIG spheres coupled
with microwave cavities [29, 30, 31]. This coupling is achieved through
magnetostrictive force, which combines the ferromagnet’s vibratory deformation
mode with the microwave cavity mode, and the magnetic dipole interaction. When
the mechanical mode frequency is significantly lower than the magnon
frequency, a dispersive interaction known as the magnetostrictive interaction
occurs, which can be compared to radiation pressure in larger ferromagnetic
systems [32, 33, 34]. Such interactions enable the coupling of mechanical and
magnetic degrees of freedom in nanoscale devices, facilitating the
manipulation and control of magnetic properties through mechanical means. By
exploiting this coupling, it is possible to control and manipulate the flow of
magnons within the system. Magnon blockade refers to the phenomenon where one
magnon inhibits additional magnons through a given channel.

In this paper, we present a scheme to investigate the magnon (photon) blockade
effect in an opto-magnomechanical system with magnon squeezing under the weak
driving limit, focusing on the second-order correlation function of an opto-
magnomechanical system. Analytically and numerically, we demonstrate that a
magnon (photon) blockade can be generated by controlling the nonlinear
parameter value and selecting appropriate system parameters to generate a
magnon (photon) blockade. Furthermore, we investigate the robustness of the
magnon blockade in the presence of pure dephasing in the system.

This paper is organized as follows: in Section II, we provide an overview of
the opto-magnomechanical system with magnon squeezing and present its
corresponding Hamiltonian. In Section III, we investigate the magnon (photon)
blockade using the second-order correlation function. We analyze the blockade
effect under optimal parameters using analytical and numerical methods.
Additionally, we explore system dynamics and discuss the obtained results.
Finally, we present our conclusions in Section III.1.

##  II The Model Hamiltonian

The microwave cavity and the YIG sphere, which include the magnon-photon and
magnon-phonon-induced interactions of the magnetic dipole and
magnetostriction, respectively, make up the magnomechanical system of the
cavity, as depicted in Figure 1. The YIG sphere’s geometry deforms due to the
magnetic change brought on by the magnon’s excitation. This creates vibratory
modes (phonons) in the sphere, and vice versa. In the YIG sphere position, the
external drive magnetic field is situated along the y axis, the microwave
cavity’s additional polarization magnetic field is positioned along the z
axis, and the cavity mode’s magnetic field is oriented along the x axis. The
magnetostrictive effect causes the photon and magnon to couple. One can change
the magneton-photon coupling force by using an extra magnetic field and moving
the YIG sphere. Due to the high spin density $(4.2\times 10^{21}cm^{-3})$ of
the magnon in the YIG sphere, strong coupling to the microwave cavity is
possible. The total Hamiltonian describing the system in the rotating frame
approximation is given by ($\hbar=1$)

$\displaystyle\mathcal{H}$ $\displaystyle=$
$\displaystyle\Delta_{c}c^{{\dagger}}c+\Delta_{m}m^{{\dagger}}m+\omega_{b}b^{{\dagger}}b-g_{mb}(b+b^{{\dagger}})m^{{\dagger}}m$
(1) $\displaystyle+$ $\displaystyle
g_{mc}(m^{{\dagger}}c+c^{{\dagger}}m)+{i}\lambda\big{(}(m^{\dagger})^{2}{e}^{{i}\theta}-m^{2}{e}^{-{i}\theta}\big{)}$
$\displaystyle+$ $\displaystyle\mathcal{E}(m^{{\dagger}}+m),$

Figure 1: A single-mode cavity contain inside a $YIG$ sphere with magnon
squeezing. In order to improve the magnomechanical coupling, the magnon mode
is directly driven by a microwave source (not shown), which causes a large
number of spins in a macroscopic ferrimagnet to move together in a coordinated
motion. Another force that drives the cavity is an electromagnetic field with
amplitude $\mathcal{E}$. Magnons and cavity photons are coupled through
magnetic dipole interaction, while phonons and magnons are connected through
magnetostrictive (radiation pressure-like) interaction.

where $\Delta_{c}=\omega_{c}-\omega_{l}$ and
$\Delta_{m}=\omega_{m}-\omega_{l}$ are the detunings of the microwave cavity
mode and magnon mode, respectively. When a unitary transformation
$\mathcal{U}=\exp[g_{mb}(b-b^{{\dagger}})/\omega_{b}]$ is applied to
$\mathcal{H}$ and in the case of $g_{mb}\ll\omega_{b}$, The unitary
transformation $\mathcal{U}$ decouples the two mechanical and the optical
modes in the special case of weak optomechanical coupling, i.e.
$g_{mb}/\omega_{b}\ll 1$ thus we can safely ignore the mechanical mode.
Therefore, the Hamiltonian can be further expressed as

$\displaystyle\mathcal{H}_{1}$ $\displaystyle=$
$\displaystyle\Delta_{c}c^{{\dagger}}c+\Delta_{m}m^{{\dagger}}m+g_{mc}(m^{{\dagger}}c+c^{{\dagger}}m)-\mu(m^{{\dagger}}m)^{2}$
(2) $\displaystyle+$
$\displaystyle{i}\lambda\big{(}(m^{\dagger})^{2}{e}^{{i}\theta}-m^{2}{e}^{-{i}\theta}\big{)}+\mathcal{E}(m^{{\dagger}}+m),$

where $\mu=g_{mb}^{2}/\omega_{b}$ is the Kerr-type nonlinear strength. the
Hamiltonian of non-Hermitian system is given by

$\mathcal{H}_{2}=\mathcal{H}_{1}-{i}\frac{\kappa_{c}}{2}c^{\dagger}c-{i}\frac{\kappa_{m}}{2}m^{\dagger}m.$
(3)

##  III Magnon (photon) blockade

We concentrate on the analytical solution of the non-Hermitian Schrödinger
equation as well as numerical solution using Master equation. In addition, we
talk about the weak coupling regime’s ability to realize a strong magnon
(photon) blockade ($g_{mb}\ll\omega_{b}$). The analytical expression of the
correlation function can be obtained by solving the Schrödinger equation in
the weak driving limit, i.e. ${i}\frac{\partial|\varphi(t)\rangle}{\partial
t}=\mathcal{H}_{2}|\varphi(t)\rangle$, where $|\varphi(t)\rangle$ denotes the
system’s current state.

The low-excitation subspace has a limit on the evolution space of up to 2,
i.e., $\mathcal{D}=\\{|q,r\rangle/q+r\leq 2\\}$. The state of the system with
the bare-state bases are : $|\varphi(t)\rangle=\sum_{q,r}^{q+r\leq
2}\mathcal{P}_{qr}(t)|q,r\rangle$, where $\mathcal{C}_{qr}(t)$ is the
amplitude of the bare state’s probability. The symbol $|q,r\rangle$ with $q$
denotes the presence of magnons in the cavity, while $r$ denotes the presence
of photons in the cavity. The Schrödinger equation generates the subsequent
set of linear differential equations for the probability amplitudes
$\mathcal{P}_{qr}(t)$

$\displaystyle{i}\frac{\partial\mathcal{P}_{00}}{\partial t}$ $\displaystyle=$
$\displaystyle\mathcal{E}\mathcal{P}_{10}-{i}\sqrt{2}\lambda\mathcal{P}_{20},$
$\displaystyle{i}\frac{\partial\mathcal{P}_{10}}{\partial t}$ $\displaystyle=$
$\displaystyle\mathcal{E}\mathcal{P}_{00}+(\Delta^{\prime}_{m}-\mu)\mathcal{P}_{10}+g_{mc}\mathcal{P}_{01}+\sqrt{2}\mathcal{E}\mathcal{P}_{20},$
$\displaystyle{i}\frac{\partial\mathcal{P}_{01}}{\partial t}$ $\displaystyle=$
$\displaystyle
g_{mc}\mathcal{P}_{10}+\Delta^{\prime}_{c}\mathcal{P}_{01}+\mathcal{E}\mathcal{P}_{11},$
$\displaystyle{i}\frac{\partial\mathcal{P}_{11}}{\partial t}$ $\displaystyle=$
$\displaystyle\sqrt{2}g_{mc}\mathcal{P}_{20}+\mathcal{E}\mathcal{P}_{01}+(\Delta^{\prime}_{c}+\Delta^{\prime}_{m}-\mu)\mathcal{P}_{11}$
$\displaystyle+$ $\displaystyle\sqrt{2}g_{mc}\mathcal{P}_{02},$
$\displaystyle{i}\frac{\partial\mathcal{P}_{02}}{\partial t}$ $\displaystyle=$
$\displaystyle\sqrt{2}g_{mc}\mathcal{P}_{11}+2\Delta^{\prime}_{c}\mathcal{P}_{02},$
$\displaystyle{i}\frac{\partial\mathcal{P}_{20}}{\partial t}$ $\displaystyle=$
$\displaystyle\sqrt{2}\mathcal{E}\mathcal{P}_{10}+2(\Delta^{\prime}_{m}-2\mu)\mathcal{P}_{20}+\sqrt{2}g_{mc}\mathcal{P}_{11}$
(4) $\displaystyle+$
$\displaystyle{i}\lambda\sqrt{2}{e}^{i\theta}\mathcal{P}_{00},$

where $\Delta^{\prime}_{c}=\Delta_{c}-{i}\kappa_{c}/2$ and
$\Delta^{\prime}_{m}=\Delta_{m}-{i}\kappa_{m}/2$. The above set of equations
(Eq. III) can be solved analytically in the steady-state by setting
$\partial\mathcal{P}_{qr}/\partial t$ to zero. Under the condition of weak
driving ($\mathcal{E}\ll\kappa$), the probability amplitudes satisfy the
relation $|\mathcal{P}_{00}|\simeq
1\gg|\mathcal{P}_{01}|,|\mathcal{P}_{10}|\gg|\mathcal{P}_{11}|,|\mathcal{P}_{02}|,|\mathcal{P}_{20}|$.
The corresponding steady-state values of the probability amplitudes, assuming
$\Delta_{c}=\Delta_{m}=\Delta$, $\kappa_{c}=\kappa_{m}=\kappa$, and
$\theta=0$, are given by $\mathcal{P}_{10}=\mathcal{E}\Delta^{\prime}/\chi$,
$\mathcal{P}_{01}=g_{mc}\mathcal{E}/\chi$,
$\mathcal{P}_{20}=-\mathcal{A}/\mathcal{B}$ and
$\mathcal{P}_{02}=-\mathcal{C}/\mathcal{B}$, where
$\chi=g^{2}_{mc}+\mu\Delta^{\prime}-\Delta^{\prime 2}$,

$\displaystyle\mathcal{A}$ $\displaystyle=$
$\displaystyle\mathcal{E}^{2}\Delta^{\prime
2}(2\Delta^{\prime}-\mu)-i\lambda\chi\left(\chi-\Delta^{\prime 2}\right)$
$\displaystyle\mathcal{B}$ $\displaystyle=$
$\displaystyle\sqrt{2}\chi\left[\Delta^{\prime}(\Delta^{\prime}-2\mu)(2\Delta^{\prime}-\mu)-2g^{2}_{mc}(\Delta^{\prime}-\mu)\right]$
$\displaystyle\mathcal{C}$ $\displaystyle=$ $\displaystyle
g^{2}_{mc}[2\mathcal{E}^{2}(\Delta^{\prime}-\mu)+i\lambda\chi].$ (5)

with $\Delta^{\prime}_{m}=\Delta^{\prime}_{c}=\Delta^{\prime}$. By examining
the magnon (photon) statistics, we can investigate the occurrence of magnon
(photon) blockade. The second-order correlation function in the steady state
provides valuable information in this regard. Specifically, we analyze
$g_{m}^{(2)}(0)$ and $g_{c}^{(2)}(0)$, which are respectively given by:

$g_{m}^{(2)}(0)=\frac{2|\mathcal{P}_{20}|^{2}}{|\mathcal{P}_{10}|^{4}}\quad;\quad
g_{c}^{(2)}(0)=\frac{2|\mathcal{P}_{02}|^{2}}{|\mathcal{P}_{01}|^{4}}.$ (6)

The condition $g_{j}^{(2)}(0)>1$ ($j=m,c$) corresponds to the magnon (photon)
bunching effect, indicating a higher probability of detecting multiple magnons
(photons) simultaneously. Conversely, when $g_{j}^{(2)}(0)<1$ ($j=m,c$), it
signifies magnon (photon) antibunching, where there is a reduced likelihood of
detecting multiple magnons (photons) together. Moreover, the magnon blockade
effect occurs in the cavity when $g_{m}^{(2)}(0)=0$, indicating the complete
suppression of magnon emission. Similarly, the strong photon blockade effect
can be observed in the cavity when $g_{c}^{(2)}(0)=0$, indicating the complete
suppression of photon emission. The Master equation approach, which enables us
to perform numerical calculations, can be expressed as follows:

$\frac{{\rm d}\rho}{{\rm
d}t}=-{i}[\mathcal{H}_{1},\rho]+\mathcal{L}_{c}(\rho)+\mathcal{L}_{m}(\rho),$
(7)

where $\mathcal{L}_{c}(\rho)=\frac{\kappa}{2}(2c\rho
c^{+}-c^{\dagger}c\rho-\rho c^{\dagger}c)$ and
$\mathcal{L}_{m}(\rho)=\frac{\kappa}{2}(2m\rho m^{+}-m^{\dagger}m\rho-\rho
m^{\dagger}m)$. $\mathcal{H}_{1}$ is the Hamiltonian given by Eq.(2). $\kappa$
is the corresponding decay rate of the mode.

Figure 2: Plot of the equal-time second-order correlation function
$g_{j}^{(2)}(0)$ $(j=m,c)$ versus the detuning $\Delta$ with
$\lambda=\lambda^{opt}$. The blue solid line represents the analytical
results, and the red dotted line represents the numerical results,
respectively. See the text for value of other parameters.

###  III.1 RESULTS AND DISCUSSIONS

The achieved magnon and photon blockade effects in this study correspond to
the Unconventional Cross-Phase-Modulation Blockade (UCPB) mechanism, which
arises from destructive quantum interference. The following parameter values
were used: $\kappa/2\pi=1\times 10^{6}$ Hz, $\omega_{b}=\kappa$,
$\mathcal{E}=0.01\kappa$, $g_{mb}=3\kappa$, and $g_{mc}=0.5\kappa$. The
optimal parameter pairs $\lambda$ and $\Delta$ can be obtained alongside the
other fixed parameters. The real solutions for the optimal parameter values
are as follows: (i) for magnon blockade, ${\Delta^{opt}\to
9.03\omega_{b},\lambda^{opt}\to 2\times 10^{-4}\omega_{b}}$, and (ii) for
photon blockade, ${\Delta^{opt}\to-0.03\omega_{b},\lambda^{opt}\to 4\times
10^{-4}\omega_{b}}$.

Figure 3: Plot of the analytical results of the equal-time second-order
correlation function $g_{j}^{(2)}(0)$ $(j=m,c)$ versus the detuning $\Delta$
for different values of the parameter $\lambda$. The dashed horizontal line
delimits the region under which represents the photon-magnon antibunching in
the cavity. See the text for value of other parameters.

In Fig. 2, we present the equal-time second-order correlation function
$g_{m}^{(2)}(0)$ and $g_{c}^{(2)}(0)$ as a function of the normalized detuning
$\Delta/\omega_{b}$ for optimal values obtained from both analytic and
numerical solutions. In Fig. 2(a), the second-order correlation function
exhibits a strong unconventional magnon blockade at the precise value of
$\Delta=9.03\omega_{b}$ for $\lambda=2\times 10^{-4}\omega_{b}$, as accurately
predicted by the analytically calculated optimal parameters. Similarly, in
Fig. 2(b), unconventional photon blockade is observed at
$\Delta=-0.03\omega_{b}$ for $\lambda=4\times 10^{-4}\omega_{b}$. This
behavior can be attributed to destructive quantum interference between
different excitation paths. Furthermore, in Fig. 2(a), another small dip is
observed with $g_{m}^{(2)}(0)>1$, indicating a bunching effect. The numerical
results are in excellent agreement with the analytical results, verifying
their consistency under weak driving conditions. It is noteworthy that the
optimal value leading to a complete photon-magnon blockade is achieved through
analytical calculations.

Figures 3(a)-3(b) display the second-order correlation function
$g_{j}^{(2)}(0)$ ($j=m,c$) as a function of the detuning $\Delta$ for
different values of $\lambda$, while keeping other parameters fixed. It is
observed that the values of $g_{m}^{(2)}(0)$ are significantly smaller for the
listed ${\Delta,\lambda}$ pairs than for others. Specifically, Figure 3(a)
demonstrates a strong unconventional magnon blockade effect at
$\Delta=9.03\omega_{b}$ with $\lambda=2\times 10^{-4}\omega_{b}$, while Figure
3(b) shows a strong unconventional photon blockade effect at
$\Delta=-0.03\omega_{b}$ with $\lambda=4\times 10^{-4}\omega_{b}$. These
results agree with the analytically calculated optimal parameters. We note
that when $\lambda=0$, implying the absence of magnon squeezing, the
unconventional magnon blockade (UMB) cannot be achieved, whereas the
conventional magnon blockade (CMB) is attainable at $\Delta\approx
4.5\omega_{b}$, as shown in Figure 3(a). Additionally, for instance, when
$\lambda=500$ Hz, a small dip appears corresponding to the conventional magnon
blockade, as depicted in Figure 3(a). On the other hand, in Figure 2(a), when
$\lambda=\lambda^{opt}$, another small dip exists with $g_{m}^{(2)}(0)>1$,
indicating a bunching effect. As anticipated in Figure 2(b), the different
curves lack a sharp dip, implying photon blockade cannot be achieved.

Figure 4: (a) Time evolution of the second-order correlation function,
$g_{m}^{(2)}(\tau)$ $(t=0)$, (b) Plot of $g_{m}^{(2)}(0)$ as a function of $t$
$(\tau=0)$. See the text for value of other parameters.

Figure 5: Plot of the equal-time second-order correlation function
$g_{j}^{(2)}(0)$ $(j=m,c)$ versus the detuning $\Delta$ for different pure-
dephasing rates $\gamma_{p}$. See the text for value of other parameters.

In Fig. 4(a), we depict the time evolution of the second-order correlation
function $g_{m}^{(2)}(\tau)$, which is calculated using the expression:

$g_{m}^{(2)}(\tau)=\frac{\langle
m^{\dagger}(t)m^{\dagger}(t+\tau)m(t+\tau)m(t)\rangle}{\langle
m^{\dagger}(t)m(t)\rangle^{2}}.$ (8)

This quantity is a significant measure of the photon statistical properties.
It is proportional to the joint probability of detecting one photon at $t=0$
and detecting the next photon at $t=\tau$. It can be observed that for delay
times $\tau>0$, $g_{m}^{(2)}(\tau)$ is greater than $g_{m}^{(2)}(0)$. This
indicates that the emitted photons exhibit blockade behavior and possess sub-
Poissonian statistics. Notably, as the delay time $\tau$ increases,
$g_{m}^{(2)}(\tau)$ also increases and approaches its maximum value (unity) at
$\tau=1.6\mu$s, which corresponds approximately to the lifetime of the photons
in the cavity. This signifies that with increasing time delay, photons become
more coherent and exhibit stronger correlation.

In Fig. 4(b), we present the time evolution of the second-order correlation
function, $g_{m}^{(2)}(0)$. It can be observed that $g_{m}^{(2)}(0)$ decreases
as time increases. Additionally, there is a significant oscillation in the
initial short time, which is a result of the classical driving field. However,
this oscillation gradually stabilizes, reaching a steady state. Throughout the
time evolution, the second-order correlation function remains below 1
($g_{m}^{(2)}(0)\ll 1$), ensuring the validity of a strong magnon blockade.

###  III.2 Pure dephasing effects

Pure dephasing can have a negative impact on the magnon (photon) blockade
characteristics and introduce unwanted decoherence in the system. Hence, we
investigate the influence of pure dephasing on cavity photon antibunching
characteristics. To model the effects of pure dephasing, an additional
Lindblad term of the form
$\mathcal{L}_{p}(\rho)=\frac{\gamma_{p}}{2}[2c^{{\dagger}}c\rho
c^{\dagger}c-(c^{{\dagger}}c)^{2}\rho-\rho(c^{{\dagger}}c)^{2}]$ is included
in the master equation, where $\gamma_{p}$ represents the pure dephasing rate
of the cavity mode.

In Fig. 5(a)-(b), we present the plots of the equal-time second-order
correlation function $g_{j}^{(2)}(0)$ $(j=m,c)$ against the detuning $\Delta$
for various values of the pure-dephasing rates $\gamma_{p}$. We observe that
even after including pure-dephasing induced losses in the system, photon
blockade and magnon blockade remain robust. When $\gamma_{p}=0$,
$g_{m}^{(2)}(0)$ exhibits strong sub-Poissonian quantum statistics at the
detuning $\Delta=\omega_{b}$. Furthermore, $g_{m}^{(2)}(0)$ increases with
increasing decay rates $\gamma_{p}$. In other words, for higher pure-dephasing
rates, $g_{m}^{(2)}(0)$ approaches classical Poissonian statistics similar to
a thermal source, as expected (Fig. 5(a)). Similar features can be observed
for $g_{c}^{(2)}(0)$ for this value of $\gamma_{p}$ (Fig. 5(b)).

##  IV Conclusion

In conclusion, we have proposed a scheme to achieve strong magnon (photon)
blockade in an opto-magnomechanical system under weak pump driving with magnon
squeezing. We investigated the magnon (photon) correlations in terms of the
second-order correlation function. The magnon (photon) blockade is realized
under optimal parameters, and we have found complete agreement between the
analytical and numerical results. Additionally, we have analyzed that the
magnon (photon) blockade effect can be achieved by adjusting the nonlinear
parameter and selecting appropriate parameters. We discussed both
unconventional and conventional magnon blockade. Furthermore, we have studied
the time-delay second-order correlation function in the presence of magnon
squeezing. Finally, we have shown that the magnon blockade remains robust even
in the presence of pure-dephasing induced losses in the system.

## References

  * Glauber [2006] R. J. Glauber, Nobel lecture: One hundred years of light quanta, Rev. Mod. Phys. 78, 1267 (2006).
  * Hanbury Brown  _et al._ [1952] R. Hanbury Brown, R. C. Jennison, and M. K. D. Gupta, Apparent Angular Sizes of Discrete Radio Sources: Observations at Jodrell Bank, Manchester, Nature (London) 170, 1061 (1952).
  * Brown and Twiss [1956] R. H. Brown and R. Q. Twiss, Correlation between Photons in two Coherent Beams of Light, Nature (London) 177, 27 (1956).
  * Kimble  _et al._ [1977] H. J. Kimble, M. Dagenais, and L. Mandel, Photon antibunching in resonance fluorescence, Phys. Rev. Lett. 39, 691 (1977).
  * Imamoḡlu  _et al._ [1997] A. Imamoḡlu, H. Schmidt, G. Woods, and M. Deutsch, Strongly interacting photons in a nonlinear cavity, Phys. Rev. Lett. 79, 1467 (1997).
  * Liao and Law [2010] J.-Q. Liao and C. K. Law, Correlated two-photon transport in a one-dimensional waveguide side-coupled to a nonlinear cavity, Phys. Rev. A 82, 053836 (2010).
  * Rabl [2011] P. Rabl, Photon blockade effect in optomechanical systems, Phys. Rev. Lett. 107, 063601 (2011).
  * Amazioug  _et al._ [2023a] M. Amazioug, S. Singh, B. Teklu, and M. Asjad, Coherent feedback control of quantum correlations in cavity magnomechanical system with magnon squeezing, arXiv e-prints , arXiv (2023a).
  * Sarma and Sarma [2017] B. Sarma and A. K. Sarma, Quantum-interference-assisted photon blockade in a cavity via parametric interactions, Phys. Rev. A 96, 053827 (2017).
  * Liao and Law [2013] J.-Q. Liao and C. K. Law, Correlated two-photon scattering in cavity optomechanics, Phys. Rev. A 87, 043809 (2013).
  * Liao and Nori [2013] J.-Q. Liao and F. Nori, Photon blockade in quadratically coupled optomechanical systems, Phys. Rev. A 88, 023853 (2013).
  * Wang  _et al._ [2015] H. Wang, X. Gu, Y.-x. Liu, A. Miranowicz, and F. Nori, Tunable photon blockade in a hybrid system consisting of an optomechanical device coupled to a two-level system, Phys. Rev. A 92, 033806 (2015).
  * Zhu  _et al._ [2018] G.-L. Zhu, X.-Y. Lü, L.-L. Wan, T.-S. Yin, Q. Bin, and Y. Wu, Controllable nonlinearity in a dual-coupling optomechanical system under a weak-coupling regime, Phys. Rev. A 97, 033830 (2018).
  * Zou  _et al._ [2019] F. Zou, L.-B. Fan, J.-F. Huang, and J.-Q. Liao, Enhancement of few-photon optomechanical effects with cross-kerr nonlinearity, Phys. Rev. A 99, 043837 (2019).
  * Birnbaum  _et al._ [2005] K. M. Birnbaum, A. Boca, R. Miller, A. D. Boozer, T. E. Northup, and H. J. Kimble, Photon blockade in an optical cavity with one trapped atom, Nature 436, 87 (2005).
  * Faraon  _et al._ [2008] A. Faraon, I. Fushman, D. R. Englund, N. G. Stoltz, P. M. Petroff, and J. Vucković, Coherent generation of non-classical light on a chip via photon-induced tunnelling and blockade, Nature Physics 4, 859 (2008).
  * Reinhard  _et al._ [2011] A. Reinhard, T. Volz, M. Winger, A. Badolato, K. J. Hennessy, E. L. Hu, and A. Imamoğlu, Strongly correlated photons on a chip, Nature Photonics 6, 93 (2011).
  * Müller  _et al._ [2015] K. Müller, A. Rundquist, K. A. Fischer, T. Sarmiento, K. G. Lagoudakis, Y. A. Kelaita, C. Sánchez Muñoz, E. del Valle, F. P. Laussy, and J. Vučković, Coherent generation of nonclassical light on chip via detuned photon blockade, Phys. Rev. Lett. 114, 233601 (2015).
  * Knill  _et al._ [2001] E. Knill, R. Laflamme, and G. J. Milburn, A scheme for efficient quantum computation with linear optics, Nature 409, 46 (2001).
  * Kok  _et al._ [2007] P. Kok, W. J. Munro, K. Nemoto, T. C. Ralph, J. P. Dowling, and G. J. Milburn, Linear optical quantum computing with photonic qubits, Rev. Mod. Phys. 79, 135 (2007).
  * Liew and Savona [2010] T. C. H. Liew and V. Savona, Single photons from coupled quantum modes, Phys. Rev. Lett. 104, 183601 (2010).
  * Bamba  _et al._ [2011] M. Bamba, A. Imamoğlu, I. Carusotto, and C. Ciuti, Origin of strong photon antibunching in weakly nonlinear photonic molecules, Phys. Rev. A 83, 021802 (2011).
  * Bajcsy  _et al._ [2013] M. Bajcsy, A. Majumdar, A. Rundquist, and J. Vučković, Photon blockade with a four-level quantum emitter coupled to a photonic-crystal nanocavity, New Journal of Physics 15, 025014 (2013).
  * Leoński and Tanaś [1994] W. Leoński and R. Tanaś, Possibility of producing the one-photon state in a kicked cavity with a nonlinear kerr medium, Phys. Rev. A 49, R20 (1994).
  * Tian and Carmichael [1992] L. Tian and H. J. Carmichael, Quantum trajectory simulations of two-state behavior in an optical cavity containing one atom, Phys. Rev. A 46, R6801 (1992).
  * Hoffman  _et al._ [2011] A. J. Hoffman, S. J. Srinivasan, S. Schmidt, L. Spietz, J. Aumentado, H. E. Türeci, and A. A. Houck, Dispersive photon blockade in a superconducting circuit, Phys. Rev. Lett. 107, 053602 (2011).
  * Snijders  _et al._ [2018] H. J. Snijders, J. A. Frey, J. Norman, H. Flayac, V. Savona, A. C. Gossard, J. E. Bowers, M. P. van Exter, D. Bouwmeester, and W. Löffler, Observation of the unconventional photon blockade, Phys. Rev. Lett. 121, 043601 (2018).
  * Vaneph  _et al._ [2018] C. Vaneph, A. Morvan, G. Aiello, M. Féchant, M. Aprili, J. Gabelli, and J. Estève, Observation of the unconventional photon blockade in the microwave domain, Phys. Rev. Lett. 121, 043602 (2018).
  * Lachance-Quirion  _et al._ [2019] D. Lachance-Quirion, Y. Tabuchi, A. Gloppe, K. Usami, and Y. Nakamura, Hybrid quantum systems based on magnonics, Applied Physics Express 12, 070101 (2019).
  * Amazioug  _et al._ [2023b] M. Amazioug, B. Teklu, and M. Asjad, Enhancement of magnon–photon–phonon entanglement in a cavity magnomechanics with coherent feedback loop, Scientific Reports 13, 3833 (2023b).
  * Yuan  _et al._ [2022] H. Yuan, Y. Cao, A. Kamra, R. A. Duine, and P. Yan, Quantum magnonics: When magnon spintronics meets quantum information science, Physics Reports 965, 1 (2022).
  * Zhang  _et al._ [2016] X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, Cavity magnomechanics, Science Advances 2, e1501286 (2016).
  * Asjad  _et al._ [2023] M. Asjad, J. Li, S.-Y. Zhu, and J. You, Magnon squeezing enhanced ground-state cooling in cavity magnomechanics, Fundamental Research 3, 3 (2023).
  * Fan  _et al._ [2022] Z.-Y. Fan, R.-C. Shen, Y.-P. Wang, J. Li, and J. Q. You, Optical sensing of magnons via the magnetoelastic displacement, Phys. Rev. A 105, 033507 (2022).


# Spin/Phonon Dynamics in Single Molecular Magnets: I. quantum embedding.

Nosheen Younas1,2,3, Yu Zhang2, Andrei Piryatinski2, Eric R. Bittner1,3 1
Department of Physics, University of Houston, Houston, Texas 77204, United
States
2Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico
87545, United States
3Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, New
Mexico 87545, United States

###### Abstract

Single molecular magnets (SMMs) and metal organic frameworks (MOFs) have
attracted significant interest because of their potential in quantum
information processing, scalable quantum computing, and extended lifetimes and
coherence times. The limiting factor in these systems is often the spin
dephasing caused by interactions and couplings with the vibrational motions of
the molecular framework. This work introduces a systematic
projection/embedding scheme to analyse spin-phonon dynamics in molecular
magnets. This scheme consolidates all spin/phonon couplings into a few
collective degrees of freedom. In this first of two papers, we establish the
projection scheme and show its equivalence to the complete set of modes within
second-order Redfield dynamics. Using parameters obtained from ab initio
methods for spin/phonon coupling via the Zeeman interaction, we apply this
approach to compute the electronic spin relaxation times for a single molecule
qubit VOPc(OH)8, which features a single unpaired electron located on the
central vanadium atom. However, our general embedding scheme can be applied to
any single-molecule magnet or qubit MOF with any coupling/interaction
Hamiltonian. This development offers a crucial tool for simulating spin
relaxation in complex environments with significantly reduced computational
complexity.

##  I Introduction

Recent efforts have explored various physical systems for qubit architectures
in the pursuit of practical quantum computing realisation, including
superconducting circuits [1], trapped ions [2], quantum dots [3], single
molecular magnets (SMMs) [4], and neutral atoms [5]. Recently, molecular
magnets and metal-organic frameworks (MOFs) have attracted attention for their
potential in quantum information processing, high-density magnetic data
storage, and as magnetic resonance contrast agents [6, 7, 8, 9, 10, 11]. These
materials allow for the tailoring of magnetic behaviour through the synthetic
versatility of molecular compounds. However, they face challenges because they
are susceptible to molecular vibrations, which can affect qubit performance.
The interaction between spins and molecular vibrations plays a crucial role in
molecular imaging, optoelectronics, and quantum technology, influencing
applications such as MRI, energy harvesting, and understanding qubit
decoherence [7]. This spin-phonon coupling is particularly significant in
SMMs, known for their magnetic bistability and memory retention at low
temperatures. MOF qubit engineering aims to slow the spin relaxation dynamics
to enhance magnetic memory and quantum coherence.

Figure 1: (left) Single Molecule qubit of VOPc(OH)8 (a VOPc derivative) with
carbon, nitrogen, hydrogen, and oxygen atoms shown in grey, blue, white, and
red, respectively. The spin-qubit electron is shown with a pink arrow. (right)
Energy level diagram for the d-orbitals for V+4. Oxygen with a short double
bond in vanadium IV causes $d_{xy}$ to be the lowest energy. Hence, due to
$d^{1}$ configuration of V+4, the spin-half qubit electron resides in $d_{xy}$
orbital. [12, 13]

Among various MOF architectures, vanadyl phthalocyanine (VOPc) is a prominent
candidate for room-temperature quantum computing [13, 14, 15, 16]. Fig. 1
shows a single VOPc molecule along with the d-orbital structure of the V^+4
centre hosting a single electron (in $d^{1}$ configuration) that is
responsible for the qubit nature of this molecule. First-principles
calculations, especially density functional theory (DFT), have become
essential in studying the electronic and magnetic properties of these
molecular magnets [17, 18, 19, 20, 21]. These methods facilitate the
investigation of spin-phonon dynamics without phenomenological parameters,
aiding in the design of materials with tailored properties. However,
challenges remain in accurately modelling these processes due to the high
computational cost and complexity of spin-phonon interactions, which involve
extensive degrees of freedom.

Highlighting these shortcomings underscores the value of exploring other
first-principles approaches in the theoretical and computational literature
that focus on phonon-assisted spin relaxation modelling. Lunghi et al.
recently used a first-principles approach to model the spin dynamics of
$\ce{VO(acac)2}$, focussing on first-order spin-phonon couplings through
phonon-assisted modulation of Zeeman and dipolar Hamiltonians.[22] This work
validates using a low-order Taylor expansion approach for Redfield dynamics.
However, even when modelling the complete crystal environment and phonons from
the entire Brillouin zone using $3\times 3\times 3$ supercells containing 1620
atoms at the DFT level, the estimates of spin-relaxation times remain orders
of magnitude different from experimental results [23]. Similarly, other first-
principles studies that incorporate Raman and Orbach mechanisms for phonon-
assisted spin relaxation in single-molecule magnets (and qubits) only capture
experimental relaxation values within a few orders of magnitude [24, 25, 22].
Despite quantitative disagreements, it is crucial to recognise that these
first-principles studies successfully capture the qualitative aspects of
phonon-assisted spin relaxation without relying on adjustable parameters. This
provides a powerful tool for understanding quantum spin degradation mechanisms
as well as engineering molecular architectures that protect these quantum
properties.

Beyond first-principles efforts, various phenomenological methods have been
employed to deduce the functional form of the $T_{1}$ variation as a function
of magnetic field and temperature by fitting curves to experimental data [26,
27]. Additionally, researchers have used ligand field methods to
systematically analyse the impact of each phonon on spin relaxation through
both direct and Raman processes. These analyses aim to uncover the role of
each phonon based on phonon symmetry. However, these approaches require
parameter adjustment to match the experimental values of $T_{1}$
quantitatively [12]. Another notable effort has been presented by Aruachan et
al. using a semi-empirical approach by constructing a parameterised Redfield
quantum master equation to describe the interaction of molecular spin qubits
with lattice phonons and electron/nuclear spin baths. Based on the Haken-
Strobl theory [28], this approach treats system-reservoir interactions as
stochastic fluctuations of the system Hamiltonian, with spin-lattice
interaction and spin-spin interaction modelled by a fluctuating gyromagnetic
tensor and a fluctuating local magnetic field, respectively. The master
equation is then derived up to the first order in fluctuations. This semi-
empirical procedure introduces fitting parameters into the bath spectral
densities of the Redfield tensors, which must be calibrated by fitting to
experimental measurements. [29] Although this methodology yields excellent
quantitative and qualitative agreement with the experimental data of $T_{1}$
and $T_{2}$ dynamics, it depends on prior experimental results for
calibration.

Our approach in this paper follows the steps of other first-principles studies
and extends them by isolating a handful of collective motions responsible for
the spin-phonon interaction. Our first-principles model presented in this work
is utilised in direct spin-phonon interaction mechanism. Although this
interaction is crucial at high fields and low temperatures ($B>5\ T$ and
$T<20\ K$), it does not prevail at low magnetic fields and intermediate
temperatures. Consequently, deviations arise when comparing our results with
the experimental data (Fig. 7). However, such deviations are expected in
first-principles calculations, which do not employ adjustable parameters. We
then introduce a quantum embedding approach using Singular Value Decomposition
(SVD) to simulate spin-qubit dynamics efficiently. This method is inspired by
mode projection schemes of Bittner et al. [30, 31, 32]. It reduces the number
of phonon modes required for simulating the spin dynamics in molecular
frameworks. Here we demonstrate the proficiency of this technique by
simulating the spin lifetime dynamics in VOPc(OH)8 as an example. This
approach also bridges the gap between SVD-based quantum embedding techniques
and mode projection, creating an isomorphism between the two methodologies.

The following section details the theoretical framework for spin-phonon
interaction in molecular qubits and the development of the SVD-based mode
projection scheme. We also present details of the implementation, simulation
results, and comparisons with previous studies. Finally, we summarises our
findings and discusses future directions. In our subsequent paper (Ref. [33])
we use the results presented herein to perform numerically exact quantum
dynamical calculations of the spin/phonon system and show how these projected
modes act as irreversible conduits of quantum information from the spin into
the thermal environment.

##  II Theory

###  II.1 Spin-Vibrational Hamiltonian

Figure 2: Schematic diagram of quantum embedding for spin-phonon interaction.
The SVD projects the original phonon modes into a primary set of phonons
directly coupled to central spins and a secondary phonon bath that has no
coupling to the central spin but is coupled to the primary phonon bath,
leading to the lifetime of the primary phonon modes.

Both SMMs and MOFs feature one or a few central spins surrounded by
environments of nuclear spins and molecular phonons (or vibrations), as
illustrated in Figure 1. Typically, in transition-metal systems, the Zeeman
term dominates. Disregarding the hyperfine and dipole-dipole terms, we express
the spin Hamiltonian as

$H_{s}=\beta{\vec{\mathbf{S}}}\cdot\bm{g}\cdot{\vec{\mathbf{B}}}$ (1)

where $\beta=\mu_{B}/\hbar$ represents the ratio of the Bohr magneton to the
reduced Planck constant,
${\vec{\mathbf{S}}}=(\hat{S}_{x},\hat{S}_{y},\hat{S}_{z})$ denotes the vector
spin, ${\vec{\mathbf{B}}}=(B_{x},B_{y},B_{z})^{T}$ represents the magnetic
field, and $\bm{g}$ is the $g$-tensor that characterises the coupling to the
external magnetic field. For VOPc(OH)8, we computed the g-tensor as

$\mathbf{g}=\begin{pmatrix}2.073&-1.64\times 10^{-2}&-9.25\times 10^{-3}\\\
-1.64\times 10^{-2}&2.037&3.50\times 10^{-3}\\\ -9.22\times 10^{-3}&3.58\times
10^{-3}&2.032\end{pmatrix}$ (2)

using density functional theory as discussed below. The fact that this is not
diagonal poses no real problem since it can be brought into diagonal form by
rotating the inertial frame of the molecule to the external (i.e. laboratory)
frame of the applied magnetic field giving
$(g_{x},g_{y},g_{z})=(2.081,2.031,2.030)$ in the rotated frame. Although these
values differ from the experimental values of
$(g_{x},g_{y},g_{z})=(1.966,1.989,1.989)$ reported in Ref. 12 for VoPc, we
made no effort to adjust the PBE exchange/correlation functional used in our
calculations to fit our computed values to the experimental values.

Although nuclear displacements do not appear directly in this expression,
$\bm{g}$ depends on the local molecular geometry via second-order perturbation
theory;

$\displaystyle
g_{\alpha\beta}=g_{o}+\sum_{e}\frac{\langle\psi_{g}|\hat{L}_{\alpha}|\psi_{e}\rangle\langle\psi_{e}|\hat{L}_{\beta}|\rangle}{E_{e}-E_{g}},$
(3)

where $g_{o}$ is the g value for a free electron and $\hat{L}_{\alpha}$ are
the orbital angular momentum operators coupling the ground electronic state
$|\psi_{g}\rangle$ to the unoccupied electronic states $|\psi_{e}\rangle$
taken at the equilibrium nuclear geometry. Of the terms entering this
expression, the energy gap is most sensitive to fluctuations in equilibrium
geometry. Further, we can use simple symmetry arguments to deduce which
orbital terms contribute to the tensor since the integrand of each of the
matrix elements must be in a totally symmetric irreducible representation,
i.e.

$\displaystyle\Gamma_{A_{1}}=\Gamma_{g}\otimes\Gamma_{op}\otimes\Gamma_{e}$
(4)

For a $C_{4v}$ molecule such as VOPc, the $d_{xy}$ orbital is in the $B_{2}$
irreducible representation and is coupled via $\hat{L}_{x}$ and $\hat{L}_{y}$
the orbital angular momentum operators to the $d_{xz}$ and $d_{yz}$ orbitals
($E$ irreducible representation) and similarly to the $d_{x^{2}-y^{2}}$
orbital ($B_{1}$ irreducible representation) via the $\hat{L}_{z}$ operator
($A_{2}$ irreducible representation) Recent work by Kazmierczak et al.
indicates that the couplings are the most sensitive to those vibrational
normal modes that are also in these irreducible representations. [12]

The phonon or vibrational degrees of freedom are described by the harmonic
oscillator Hamiltonian:

$\displaystyle H_{b}$
$\displaystyle=\sum_{k}^{n_{q}}\frac{1}{2}\left(\hat{p}^{2}_{k}+\omega_{k}^{2}\hat{x}_{k}^{2}\right)$
(5)

where $\hat{x}_{k}$ and $\hat{p}_{k}$ represent the mass-weighted normal
coordinates and their conjugate momenta, respectively, the summation runs over
$n_{q}$ phonon degrees of freedom. To capture the effect of the phonons on the
spin, we expand the spin Hamiltonian concerning the phonon displacements about
the ground-state geometry and truncate the expansion q5 first order. The
resulting Hamiltonian for the total system is of the form

$\displaystyle H$
$\displaystyle=\sum_{\alpha}^{n_{s}}h_{\alpha}\hat{S}_{\alpha}+\sum_{k}^{n_{q}}\frac{1}{2}(\hat{p}_{k}^{2}+\omega_{k}^{2}\hat{x}_{k}^{2})+\sum_{\alpha
k}^{n_{s},n_{q}}g_{\alpha k}\hat{S}_{\alpha}\hat{x}_{k},$ (6) $\displaystyle
H$ $\displaystyle=H_{s}+H_{b}+H_{sb}.$ (7)

Here, the first term sums over the complete set of operators describing the
$n_{s}$ spin degrees of freedom, where we have defined

$\displaystyle h_{\alpha}=\beta\sum_{j}^{3}\bm{g}_{\alpha j}\cdot B_{j}$ (8)

as the spins and external magnetic field coupling via spin operator
$\hat{S}_{\alpha}$. The last term defines the coupling between the spins and
the vibrational degrees of freedom, determined from the gradient of $H_{s}$
with respect to the nuclear displacements.

$\displaystyle g_{\alpha
k}=\beta\sum_{j}^{3}\left(\frac{\partial\bm{{g_{\alpha
j}}}}{\partial\vec{x}_{k}}\right)^{(0)}B_{j}.$ (9)

The system usually has a small number of discrete quantum levels compared to
the larger degrees of freedom (DOFs) in the bath, such that $n_{s}\ll n_{q}$.

The gradient term can be computed using Density Functional Theory (DFT)
calculations from software packages such as ORCA and CP2K. However, this
process presents significant challenges as it requires using a finite-
difference method, which is computationally demanding and must be performed
for each system individually. In this work, we used ORCA (v5.0) [34] to
compute the magnetic properties (g-tensor) of the SMMs in both equilibrium and
distorted configurations. We employ the def2-TZVP base set for metal and O
elements while using the def2-SVP basis set for other lighter elements,
complemented by an auxiliary def2-TZVP/C basis set for all elements. We
calculate the g-tensor using the DFT method with the PBE functional. [35]
Initially, we optimised the molecular system’s geometry using DFT and then
calculated the Hessian matrix to obtain the vibrational normal modes. We then
estimate the spin-phonon coupling coefficients by calculating the gradient of
$g$ with respect to the nuclear coordinates using the finite difference
method. Specifically, we approximate

$\displaystyle\frac{\partial g(\bm{R})}{\partial
X_{ix}}\approx\frac{g(\bm{R}+\Delta_{ix})-g(\bm{R}-\Delta_{ix})}{2\Delta_{ix}}$
(10)

where $X_{ix}$ denotes the $x$ Cartesian component of the $i$-th atom in the
optimized geometry and $\Delta_{ix}$ is the perturbation along that direction
, $\bm{R}$ represents the equilibrium grometry. To achieve this, we distort
the optimised geometry with small forward and backward displacements ($\pm
0.01$ Å) and compute the $g$-tensors of the distorted geometries.

For a molecule with $N$ atoms, there are $3N-6$ normal modes. However, not
every vibrational mode will contribute equally to the coupling. To understand
the effectiveness of each phonon displacement towards spin-phonon coupling, we
computed the $l_{2}$ norms of the resulting g-tensor gradient with respect to
dimensionless displacements $q_{k}$ for VOPc(OH)8 via

$||\partial g/\partial q_{k}||^{2}=\sum_{i,j}^{3}\sqrt{\left(\partial
g_{ij}/\partial q_{k}\right)^{2}}.$ (11)

Consequently, just as the $l_{2}$ norm of a vector gives its length, the
$l_{2}$ norm of the gradient tensor gives its overall amplitude and indicates
the participation of each normal mode $k$ in the spin-phonon coupling.
Although spin-phonon coupling linearly depends on the applied magnetic field,
the gradient still determines the participation of each individual phonon;
therefore, using the $g$-tensor gradient instead of the full spin-phonon
coupling highlights the contribution of each phonon in the spin-phonon
coupling.

An analysis of Figure 3 reveals that the few lowest frequency modes have the
highest norm of the $g$-tensor gradient, although the high-frequency modes
also play a role. Furthermore, the amplitude of the norm $l_{2}$ for all modes
is of the order of $10^{-3}$ (compared to the values of the g-tensors of
$\approx 2$). Therefore, spin-phonon coupling is relatively weak in VOPc(OH)8
MOFs, which explains their relatively long $T_{1}$ relaxation times and their
candidacy for molecular qubits.

Figure 3: $l_{2}$ norm of the g-tensor derivative ($\vec{\nabla}g$) with
respect to each phonon mode for VOPc(OH)8 as obtained from DFT calculations.
The normal mode frequencies ($cm^{-1}$) are on the x-axis, with $l_{2}$ norms
on the y-axis (unitless). The histograms show a distinct band structure, with
the first group of frequencies representing the bending modes of the The
second group, the phenylocyanine ring, belongs to the ring-breathing modes,
and the highest-frequency modes are the OH and CH stretch modes.

Defining the entire system/bath Hamiltonian is straightforward, but solving
its complete quantum dynamics is a significant challenge for any reasonably
sized system. This difficulty arises because quantum dynamics exhibits
superlinear scaling with the bath DOF which scales as $3N$, where $N$ is the
number of atoms involved. This challenge can be addressed by implementing a
series of systematic projections, as illustrated in Fig. 2. These projections
simplify the fully coupled (star) system into a series of independent
branches, where only a few bath degrees of freedom are directly coupled to the
system, while the rest act as a thermal reservoir. This approach effectively
transforms a high-dimensional and intractable quantum problem into a more
manageable one that can be solved using contemporary quantum dynamics
approaches.

###  II.2 Decomposition

The general model for a quantum system that interacts with an environment can
be represented as a tensor product of system operators ($\hat{S}$) and bath
($\hat{B}$) operators, with coupling constants $g$. The Hamiltonian of this
system, as shown in Eq. 6, includes the summation over all bath degrees of
freedom (DOFs) in the interaction term $H_{sb}$. As previously mentioned, Fig.
2 graphically presents this situation in its left diagram. The objective is to
reduce the number of DOFs that directly interact with the system through
projective transformations. These transformations lead to the creation of two
sets of new DOFs. One set, called _system-phonons_ , interacts directly with
the spin system, and the other, called _bath-phonons_ , interacts only with
the first set, as depicted in the right panel of Fig. 2. To achieve this,
Bittneret al. introduced a Mode Projection scheme, which determines the
optimal partitioning of the entire dynamical system into interacting
subcomponents [30, 36]. The approach has been used successfully to determine
electronic charge and energy transfer in small molecular systems with various
degrees of complexity[31, 32] using ab initio derived electron/phonon
couplings.

This work presents a single-value decomposition (SVD)–based mode projection
method for reducing the spin-phonon coupling’s dimension. Consider the
coupling matrix and its SVD, defined as

$\bm{g}=[g_{\alpha k}]_{n_{s}\times n_{q}}=\bm{U}_{n_{s}\times
n_{s}}\cdot\bm{\Sigma}_{n_{s}\times n_{q}}\cdot\bm{V}^{\dagger}_{n_{q}\times
n_{q}}$ (12)

where $\Sigma$ is a matrix of singular values with at most $\min(n_{s},n_{q})$
nonzero singular values. Taking column vectors of $V^{\dagger}$ corresponding
to non-zero singular values and creating projection operator

$\mathbb{P}=\sum_{r}|V_{r}^{\dagger}\rangle\langle V_{r}^{\dagger}|,$ (13)

projects the phonon modes into reduced dimensional subspace. Likewise its
complement

$\mathbb{Q}=\mathbb{I}-\mathbb{P},$ (14)

projects them into the complement subspace. Using these projection operators,
one could reduce the dimensionality of the spin-phonon coupling; however, to
convert this projection scheme into a canonical transformation, we proceed
with diagonalisation of the Hessian in the two projected sub-spaces. In the
original normal mode basis, the Hessian is a diagonal matrix defined by
$\Omega^{2}=[\omega_{ii}^{2}]$. Using the resolution of identity as a sum of
projection operators,

$\displaystyle\Omega^{2}$
$\displaystyle=\left(\mathbb{P}+\mathbb{Q}\right)\Omega^{2}\left(\mathbb{P}+\mathbb{Q}\right)$
(15)
$\displaystyle=\mathbb{P}\Omega^{2}\mathbb{P}+\mathbb{Q}\Omega^{2}\mathbb{Q}+\left(\mathbb{P}\Omega^{2}\mathbb{Q}+\mathbb{P}\Omega^{2}\mathbb{Q}\right)$
$\displaystyle=\Omega^{2}_{PP}+\Omega^{2}_{QQ}+\Omega^{2}_{PQ}+\Omega^{2}_{QP}$

Diagonalising $\Omega^{2}_{PP}$ and $\Omega^{2}_{QQ}$ provides the frequencies
$\omega_{k}^{2}$ for the new system and bath modes, respectively. Since the
number of DOF of phonons is much larger than the spin DOF, $\Omega^{2}_{PP}$
will have at most $n_{s}$ nonzero eigenvalues, while $\Omega^{2}_{QQ}$ will
have $n_{q}-n_{s}$ nonzero eigenvalues. The corresponding eigenvectors or
nonsingular eigenvalues for $\Omega^{2}_{PP}(\\{\tilde{x}^{P}\\})$ and
$\Omega^{2}_{QQ}(\\{\tilde{x}^{Q}\\})$ constitute the system and bath modes,
respectively. And $\Omega^{2}_{PQ}$ defines the couplings between the system
and the bath modes. Using projected modes, the phonon and interaction terms of
the total Hamiltonian Eq. 6 become

$\displaystyle H_{b}+H_{sb}$ $\displaystyle=\sum_{\alpha r}\tilde{g}_{\alpha
r}S_{\alpha}\tilde{x}_{r}+\sum_{r}^{n_{s}}\frac{1}{2}(\tilde{p}_{r}^{2}+\omega_{r}^{2}\tilde{x}_{r}^{2})+\sum_{j}^{n_{q}-n_{s}}\frac{1}{2}(\tilde{\pi}_{j}^{2}+\omega_{j}^{2}\tilde{y}_{j}^{2})+\sum_{r,j}^{n_{s},n_{q}-n_{s}}\gamma_{r,j}\tilde{x}_{r}\tilde{y}_{j}.$
(16)

The first term reflects the direct coupling between the $n_{s}$ spin degrees
of freedom and $n_{r}$ projected modes $\\{\tilde{x}_{r}\\}$ with momenta
$\\{\tilde{p}_{r}\\}$ (note, $n_{r}\leq n_{s}$). These are the primary modes
since they are directly coupled to the spin variables via projected forces
$\tilde{g}_{\alpha r}$. The second term represents the harmonic oscillator
Hamiltonian for these primary modes. The third term is a harmonic oscillator
representation of the residual modes $\\{\tilde{y}_{j}\\}$ with conjugate
momenta $\\{\tilde{\pi}_{j}\\}$. The residual modes are linearly coupled to
the primary modes via $\gamma_{r,j}$ in the fourth term.

Figure 4: Projected modes are linear combinations of original phonons. The
contribution of initial phonons in projected system modes 1, 2, and 3 is given
in the top, middle, and bottom panels, respectively.

##  III Results

Using this SVD mode projection approach and the computed spin-phonon coupling
matrix for VOPc(OH)8), we determined three primary modes at frequencies
170.05, 352.20, and 533.980 $cm^{-1}$, respectively. These projected modes are
linear combinations of the normal vibrational modes of the molecule. The
atomic displacements within the VOPc(OH)8 molecule associated with the primary
modes are shown in Fig. 4 (top). The primary mode displacements include the
bending and twisting of the phthalocyanine structure and the movement of the
central oxygen. The couplings of these modes with spin through spin operators
$\hat{S}_{x}$, $\hat{S}_{y}$, and $\hat{S}_{z}$ are given in the Supplemental
Information. Fig. 4 (bottom) shows the contribution of each original phonon to
projected modes 1, 2, and 3, at frequencies 170.05, 352.20, and 533.980
$cm^{-1}$, respectively. In particular, original modes with low frequencies
are dominant contributors, but even high-frequency modes make small
contributions to the projected modes. However, it should be pointed out that
replacing the terminal OH groups with methyl groups affects only the very high
frequencies assigned to the OH stretch. They are replaced by the CH stretch
modes. Because the spectral densities drop off quickly with phonon frequency,
replacing the OH group with the methyl group does not affect the spin
dynamics. This observation is supported by identical spin-phonon couplings and
identical phonon frequencies of the central molecular structure for the two
cases.

A useful measure of the decomposition efficiency can be obtained by
calculating the von Neumann entropy.

$\displaystyle S=-\sum_{n}r^{2}_{n}\ln r^{2}_{n}$ (17)

where $r_{n}=\lambda_{n}/||\lambda_{n}||$ are the components of the vectors
representing phonons. Using Eq. 17 and considering $r_{n}$ as elements of the
projected vectors, the von Neumann entropies for these vectors are $1.96$,
$3.22$, and $3.42$ for the first, second, and third projected modes,
respectively. Unlike the original normal modes, the projected modes exhibit an
increased entropy, highlighting the efficiency obtained by compressing the
information carried by the full environment into a small number of collective
degrees of freedom. In particular, the number of primary vectors that SVD
returns is proportional to the matrix’s rank, which is at most 3 for a single
spin-1/2 system. The original normal coordinates can be recovered by
diagonalizing the last three terms. As a result of this transformation, the
coupling between spin and phonons involves a sum over only the system modes;
since $n_{s}\ll n_{q}$ is generally satisfied in molecular magnets, the SVD-
based projection can significantly reduce the number of phonon modes in the
simulations.

Figure 5: Spectral density (Eq. 20) of the phonons specific to $\sigma_{x}$
(a), $\sigma_{y}$(b), and $\sigma_{z}$(c) (blue curve in the bottom, middle,
and top panels, respectively) at 300 K temperature with a magnetic field of 1
T along the molecular $z$-axis. The derivatives of the g-tensor elements are
shown in red. The spectral density for high-frequency modes is low even at
ambient temperature. Ultimately, this spectral density affects phonons’
participation in the spin system’s Redfield dynamics, penalising the higher
frequency components. Movies of the nuclear motions are given in the
Supplemental Information.

The couplings, both for the full set of modes and projected modes, contribute
to the bath spectral densities required by the Redfield master equation for
simulating the quantum dynamics of the spin system. Expressions for these
spectral densities are given in the Supplemental Information and take the
usual form of a sum over all bath modes. Fig. 5 presents the spectral density
for each of the three spin components without applying SVD and using full set
of original normal modes. Additionally, the g-tensor gradients projected along
each normal mode are also given. The spectral densities decay rapidly with
increasing phonon frequency, indicating that low-energy phonons predominantly
contribute to the bath spectral densities in the Redfield master equation.
Furthermore, the $x$ and $z$ components dominate the spin/phonon spectral
densities, whereas the components of the g-tensor gradient are almost
uniformly distributed throughout the vibrational spectrum.

Figure 6: Spectral density for mode projection at 300K and 1T (top), couplings
between system modes and bath modes (middle), and couplings between spin and
system modes (bottom).

After SVD mode projection, we compute the spectral density of the primary
modes using the couplings between the spin and primary modes and the phonon-
phonon couplings between the primary and residual modes. Fig. 6 shows the
spectral density along with the spin-phonon and phonon-phonon couplings. For
clarity, we show the spectral density and spin-phonon couplings through the
spin operators $\hat{S}_{x},\hat{S}_{y}$, and $\hat{S}_{z}$. The spectral
density exhibits a unique lineshape, indicating resonances at both the primary
and residual mode frequencies. This characteristic results from the driven-
dissipative harmonic oscillator nature of the primary modes, as detailed in
Eq. 20 and the Supplementary Information that accompanies this paper.

###  III.1 Spin dynamics with original modes

With the spin-phonon couplings at hand, we can investigate the spin relaxation
dynamics of vanadyl phthalocyanine (VOPc(OH)8) using the Bloch-Redfield
equations of motion for the reduced density matrix of the isolated spin,

$\displaystyle\partial_{t}\rho_{s}=-i[H_{s},\rho_{s}]+{\cal R}(\rho_{s})$ (18)

where $({\cal R}(\rho_{s})$ is the Redfield superoperator

$\displaystyle({\cal
R}(\rho_{s}))_{\alpha\beta}=\sum_{\gamma\delta}R_{\alpha\beta;\gamma\delta}[\rho_{s}]_{\gamma\delta}$
(19)

which depends upon the spectral density of the bath as implemented in the
QuTip (v.5) quantum dynamics package. [37, 38] This is given explicitly by

$S_{\alpha}(\omega)=\pi\hbar\sum_{k}^{n_{q}}\frac{g_{\alpha
k}^{2}}{\omega_{k}}\left[\left(n(\omega_{k})+1\right)\delta(\omega-\omega_{k})+n(\omega_{k})\delta(\omega+\omega_{k})\right]$
(20)

where $n(\omega)$ are the thermal phonon occupations given by Bose-Einstein
statistics. For computational purposes, we replace the delta
$\delta(\omega\pm\omega_{j})$ functions with the Gaussian functions i.e..
$\frac{1}{\sqrt{2\pi}\Delta\omega}e^{-\frac{\left(\omega\pm\omega_{j}\right)^{2}}{2\Delta\omega}}$
using $\Delta\omega=2\text{ meV}\approx 16\ cm^{-1}$ as per the FTIR spectra
of VOPc(OH) and VOPc(Me) compounds. The spectral densities and g-tensor
derivatives for the normal modes of single molecule VOPc(OH)8 are shown in
Fig. 5

Fig. 7 (a) shows the temperature dependence of the $T_{1}$ population
relaxation time and the $T_{2}$ pure dephasing time produced by our model. As
expected, $T_{1}$ shows a marked dependence on the system’s temperature,
whereas $T_{2}$ time is nearly constant over the same range due to its
dependence on nuclear hyperfine interaction that is ignored in the current
approximation. Similarly, Fig. 7 (b) compares $T_{1}$ using the classical
versus quantum versions of the spectral density functions. As expected, the
two diverge at low temperatures. These results provide a baseline for
comparison against any error introduced by the mode projection method. In
fact, over the entire temperature range, the mode-projection/reduction
approach gives identical results compared to the standard Redfield treatment
as shown in Fig. 7 (b) and discussed in the next section.

Finally, we compare our ab initio model with recent experimental data on the
VoPC system reported by Atzori et al. taken at 345 mT. In Fig. 7 (c), we
compare the ratio of $T_{2}/T_{1}$. Analytically, this ratio should equal 2 if
the spectral density takes a simple Lorenzian form. In general, the comparison
between the model and the experimental data gives The correct overall trend is
at low temperatures but differs by several orders of magnitude at higher
temperatures. This can be accounted for because the simulations only account
for a single molecule in a vacuum, thus neglecting the effects of low-energy
phonons in solid-state systems. Furthermore, the linear dependence on the
phonon displacement encoded in the spin-phonon interaction Hamiltonian does
not account for higher-order inelastic processes that become the dominant
mechanism at low magnetic fields and intermediate temperatures compared to
direct phonon coupling. [39]

Figure 7: (a) Temperature dependence of $T_{1}$ and $T_{2}$ (red and blue,
respectively) calculated from Redfield simulations using the full set of modes
at 1 Tesla Magnetic field. (b) $T_{1}$ as a function of temperature using a
complete set of modes with quantum and classical spectral density and
projected modes with classical spectral density. The green curve overlays the
orange curve. (c) Temperature dependence of the ratio between $T_{1}$ and
$T_{2}$ calculated from Redfield simulations (blue-dashed) and from
experimental observations (yellow-solid) reported by Atzori et al. at 345 mT.
[40]

###  III.2 Spin Dynamics with Projected Modes

Having defined the projected modes, we can perform numerically exact quantum
simulations within the tensor space defined by the central spin and the few
projected modes using the reduced Hamiltonian

$\displaystyle H_{red}=\sum_{\alpha}h_{\alpha}\hat{S}_{\alpha}+\sum_{\alpha
r}\tilde{g}_{\alpha,r}S_{\alpha}\tilde{x}_{r}+\sum_{r}^{n_{s}}\frac{1}{2}(\tilde{p}_{r}^{2}+\omega_{r}^{2}x_{r}^{2})$
(21)

and treating the $n_{q}$ residual modes as a thermal bath for the primary
($\tilde{x}_{r}$) degrees of freedom.

Using this, we evaluated the effectiveness of projected modes in capturing the
influence of the phonon bath on the spin system by simulating the spin
dynamics over the lifetime $T_{1}$. This simulation employs the spectral
density for the primary modes within the framework of the Redfield master
equation, as previously discussed. Fig. 7 shows the resulting temperature
dependence of $T_{1}$. Since we compute the spectral density of the projected
modes using a classical formalism, the $T_{1}$ results using projected modes
match those obtained using the full set of phonons classically. This complete
match indicates that the projected modes successfully capture the phonon bath
originally represented by 192 phonons. Additionally, in the Supplementary
Information, we compare the performance of projected modes against a selection
of original phonon modes using a trivial cutoff on spin-phonon coupling by
applying a threshold at $35\%$ of the maximum spin-phonon coupling. In doing
so, we can reproduce the relaxation $T_{1}$ as a function of temperature,
indicating that the SVD approach is highly effective in identifying an optimal
partitioning of the phonon environment.

##  IV Summary

This study presents an innovative quantum embedding approach for studying
spin-phonon dynamics in molecular magnets. The proposed method effectively
investigates the spin-lattice relaxation process with high precision. We
employed this approach to analyse the spin dynamics in single molecular
magnets. Through this methodology, we have gained a more profound
understanding of the collective intra-molecular motions governing spin
relaxation in molecular magnets, which is critical for optimising their
performance in various applications. Moreover, this research establishes a
foundation for rationalising new materials with tailored magnetic properties
suitable for quantum information processing and other technological
applications.

In addition to shedding light on the collective intramolecular motions that
dominate spin-phonon coupling, our embedding approach enables the quantum
treatment of phonons when their energies are comparable to spin splitting.
Such scenarios can occur in the high magnetic field limit and in solid-state
systems that feature low-frequency crystal phonons with energies similar to
the Zeeman splitting of the spin system. In our accompanying paper[33], we
explore the quantum entanglement between electronic spin and projected modes.
Here, the projected modes serve as conduits for channelling quantum
information from the spin to the remaining or “residual” degrees of freedom
corresponding to the modes spanning the $\mathbb{Q}$ subspace from Eq. 14.

In our future work, we intend to expand the quantum embedding method to
include spin-spin interactions in molecular magnets. This extension aims to
provide a deeper understanding of the interplay between spin-spin and spin-
phonon interactions, elucidating their roles in shaping the overall relaxation
dynamics of molecular qubits. By exploring spin-spin interactions, we
anticipate gaining enhanced insights into these complex systems, which will
contribute to developing advanced materials for quantum computing and other
magnetism applications. This proposed expansion is expected to significantly
broaden the scope and utility of the quantum embedding approach, thereby
augmenting its value as a research tool within the domain of molecular
magnets.

###### Acknowledgements.

The research presented in this article was supported by the LANL LDRD program
(number 20220047DR). LANL is operated by Triad National Security, LLC, for the
National Nuclear Security Administration of the U.S. Department of Energy
(contract no. 89233218CNA000001). The authors thank the LANL Institutional
Computing (IC) programme for access to HPC resources. The work at the
University of Houston was funded in part by the National Science Foundation
(CHE-2404788) and the Robert A. Welch Foundation (E-1337). This work was
performed, in part, at the Center for Integrated Nanotechnologies, an Office
of Science User Facility operated for the U.S. Department of Energy (DOE)
Office of Science by Los Alamos National Laboratory (Contract
89233218CNA000001) and Sandia National Laboratories (Contract DE-NA-0003525).

## Data Availability Statement

The data supporting this study’s findings are available from the corresponding
author upon reasonable request.

## Author Contribution Statement

The authors confirm their contribution to the paper: study conception and
design: ERB, YZ, and AP; data collection and simulations: NY; analysis and
interpretation of results: NY, YZ, AP, ERB; draft manuscript preparation: NY,
ERB. All authors reviewed the results and approved the final version of the
manuscript.

## References

  * [1] Krantz P, Kjaergaard M, Yan F, Orlando TP, Gustavsson S, Oliver WD. 2019 A quantum engineer’s guide to superconducting qubits. Applied Physics Reviews 6. (10.1063/1.5089550)
  * [2] Cirac JI, Zoller P. 2000 Quantum simulations with trapped ions. Nature 404, 579–586. (10.1038/35007038)
  * [3] Veldhorst M, Yang CH, Hwang JCC, Huang W, Dehollain JP, Muhonen JT, Simmons S, Laucht A, Hudson FE, Itoh KM, Morello A, Dzurak AS. 2015 A two-qubit logic gate in silicon. Nature 526, 410–414. (10.1038/nature15263)
  * [4] Bertaina S, Gambarelli S, Mitra T, Tsukerblat B, Müller A, Barbara B. 2008 Quantum oscillations in a molecular magnet. Nature 453, 203—206. (10.1038/nature06962)
  * [5] Lukin MD, Fleischhauer M, Cote R, Duan LM, Jaksch D, Cirac JI, Zoller P. 2001 Dipole Blockade and Quantum Information Processing in Mesoscopic Atomic Ensembles. Phys. Rev. Lett. 87, 037901. (10.1103/PhysRevLett.87.037901)
  * [6] Gaita-Ariño A, Luis F, Hill S, Coronado E. 2019 Molecular spins for quantum computation. Nat. Chem. 11, 301–309.
  * [7] Wasielewski MR, Forbes MDE, Frank NL, Kowalski K, Scholes GD, Yuen-Zhou J, Baldo MA, Freedman DE, Goldsmith RH, Goodson III T, Kirk ML, McCusker JK, Ogilvie JP, Shultz DA, Stoll S, Whaley KB. 2021 Exploiting chemistry and molecular systems for quantum information science. Nature Chemistry 13, 352–359. (10.1038/s41557-021-00675-3)
  * [8] Iqbal M, Tantivasadakarn N, Verresen R, Campbell SL, Dreiling JM, Figgatt C, Gaebler JP, Johansen J, Mills M, Moses SA, Pino JM, Ransford A, Rowe M, Siegfried P, Stutz RP, Foss-Feig M, Vishwanath A, Dreyer H. 2024 Non-Abelian topological order and anyons on a trapped-ion processor. Nature 626, 505–511. (10.1038/s41586-023-06934-4)
  * [9] Yu CJ, von Kugelgen S, Krzyaniak MD, Ji W, Dichtel WR, Wasielewski MR, Freedman DE. 2020 Spin and Phonon Design in Modular Arrays of Molecular Qubits. Chemistry of Materials 32, 10200–10206. (10.1021/acs.chemmater.0c03718)
  * [10] Yamabayashi T, Atzori M, Tesi L, Cosquer G, Santanni F, Boulon ME, Morra E, Benci S, Torre R, Chiesa M, Sorace L, Sessoli R, Yamashita M. 2018 Scaling Up Electronic Spin Qubits into a Three-Dimensional Metal-Organic Framework. J Am Chem Soc 140, 12090–12101. (10.1021/jacs.8b06733)
  * [11] Lavroff RH, Pennington DL, Hua AS, Li BY, Williams JA, Alexandrova AN. 2021 Recent Innovations in Solid-State and Molecular Qubits for Quantum Information Applications. The Journal of Physical Chemistry Letters 12, 10742–10745. PMID: 34758627 (10.1021/acs.jpclett.1c03269)
  * [12] Kazmierczak NP, Mirzoyan R, Hadt RG. 2021 The Impact of Ligand Field Symmetry on Molecular Qubit Coherence. J. Am. Chem. Soc. 143, 17305–17315. (10.1021/jacs.1c04605)
  * [13] Atzori M, Tesi L, Morra E, Chiesa M, Sorace L, Sessoli R. 2016 Room-Temperature Quantum Coherence and Rabi Oscillations in Vanadyl Phthalocyanine: Toward Multifunctional Molecular Spin Qubits. J Am Chem Soc 138, 2154–2157. (10.1021/jacs.5b13408)
  * [14] Bonizzoni C, Ghirri A, Atzori M, Sorace L, Sessoli R, Affronte M. 2017 Coherent coupling between Vanadyl Phthalocyanine spin ensemble and microwave photons: towards integration of molecular spin qubits into quantum circuits. Scientific Reports 7, 13096. (10.1038/s41598-017-13271-w)
  * [15] Malavolti L, Briganti M, Hänze M, Serrano G, Cimatti I, McMurtrie G, Otero E, Ohresser P, Totti F, Mannini M, Sessoli R, Loth S. 2018 Tunable Spin-Superconductor Coupling of Spin 1/2 Vanadyl Phthalocyanine Molecules. Nano Lett. 18, 7955–7961. (10.1021/acs.nanolett.8b03921)
  * [16] Nadeem S, Khan MN, Muhammad N, Ahmad S. 2019 Mathematical analysis of bio-convective micropolar nanofluid. Journal of Computational Design and Engineering 6, 233–242. (10.1016/j.jcde.2019.04.001)
  * [17] Neese F. 2006 Importance of direct spin spin coupling and spin-flip excitations for the zero-field splittings of transition metal complexes: a case study. J Am Chem Soc 128, 10213–10222. (10.1021/ja061798a)
  * [18] Timco GA, Carretta S, Troiani F, Tuna F, Pritchard RJ, Muryn CA, McInnes EJL, Ghirri A, Candini A, Santini P, Amoretti G, Affronte M, Winpenny REP. 2009 Engineering the coupling between molecular spin qubits by coordination chemistry. Nature Nanotechnology 4, 173–178. (10.1038/nnano.2008.404)
  * [19] Chilton NF, Anderson RP, Turner LD, Soncini A, Murray KS. 2013 PHI: A powerful new program for the analysis of anisotropic monomeric and exchange-coupled polynuclear d- and f-block complexes. Journal of Computational Chemistry 34, 1164–1175. (https://doi.org/10.1002/jcc.23234)
  * [20] Chilton NF, Goodwin CA, Mills DP, Winpenny RE. 2015 An electrostatic model for the determination of magnetic anisotropy in dysprosium complexes. Chem. Sci. 6, 3949–3955.
  * [21] Reta D, Kragskow JGC, Chilton NF. 2021 Ab Initio Prediction of High-Temperature Magnetic Relaxation Rates in Single-Molecule Magnets. Journal of the American Chemical Society 143, 5943–5950. (10.1021/jacs.1c01410)
  * [22] Lunghi A, Totti F, Sessoli R, Sanvito S. 2017 The role of anharmonic phonons in under-barrier spin relaxation of single-molecule magnets. Nat. Commun. 8, 14620.
  * [23] Lunghi A, Sanvito S. 2019 How do phonons relax molecular spins?. Science Advances 5, eaax7163. (10.1126/sciadv.aax7163)
  * [24] Lunghi A, Sanvito S. 2020a The Limit of Spin Lifetime in Solid-State Electronic Spins. The Journal of Physical Chemistry Letters 11, 6273–6278. PMID: 32667205 (10.1021/acs.jpclett.0c01681)
  * [25] Lunghi A, Sanvito S. 2020b Multiple spin–phonon relaxation pathways in a Kramer single-ion magnet. The Journal of Chemical Physics 153, 174113\. (10.1063/5.0017118)
  * [26] Atzori M, Morra E, Tesi L, Albino A, Chiesa M, Sorace L, Sessoli R. 2016 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. PMID: 27517709 (10.1021/jacs.6b05574)
  * [27] Santanni F, Albino A, Atzori M, Ranieri D, Salvadori E, Chiesa M, Lunghi A, Bencini A, Sorace L, Totti F, Sessoli R. 2021 Probing Vibrational Symmetry Effects and Nuclear Spin Economy Principles in Molecular Spin Qubits. Inorg Chem 60, 140–151. (10.1021/acs.inorgchem.0c02573)
  * [28] Čápek V. 1985 Generalized Haken-Strobl-Reineker model of excitation transfer. Zeitschrift für Physik B Condensed Matter 60, 101–105. (10.1007/BF01312648)
  * [29] Aruachan K, Colón YJ, Aravena D, Herrera F. 2023 Semi-empirical Haken–Strobl model for molecular spin qubits. New Journal of Physics 25, 093031. (10.1088/1367-2630/acf2bd)
  * [30] Pereverzev A, Bittner ER, Burghardt I. 2009 Energy and charge-transfer dynamics using projected modes. The Journal of Chemical Physics 131, 034104. (10.1063/1.3174447)
  * [31] Yang X, Keane T, Delor M, Meijer AJ, Weinstein J, Bittner ER. 2017 Identifying electron transfer coordinates in donor-bridge-acceptor systems using mode projection analysis. Nature Communications 8, 14554.
  * [32] Yang X, Pereverzev A, Bittner ER. 2018 Inelastic Charge Transfer Dynamics in Donor-Bridge-Acceptor Systems Using Optimal Modes. Advances in Chemical Physics, Volume 163 163, 167.
  * [33] Younas N, Zhang Y, Piryatinski A, Bittner ER. 2024 Systematic projections for quantum spin-phonon dynamics in a vanadium-based qubit: II. spin/phonon entanglement. Submitted, Proc. Roy. Soc. A.
  * [34] Neese F, Wennmohs F, Becker U, Riplinger C. 2020 The ORCA quantum chemistry program package. The Journal of Chemical Physics 152, 224108. (10.1063/5.0004608)
  * [35] Perdew JP, Burke K, Ernzerhof M. 1996 Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 77, 3865–3868. (10.1103/PhysRevLett.77.3865)
  * [36] Yang X, Bittner ER. 2014 Intramolecular Charge- and Energy-Transfer Rates with Reduced Modes: Comparison to Marcus Theory for Donor–Bridge–Acceptor Systems. The Journal of Physical Chemistry A 118, 5196–5203. (10.1021/jp503041y)
  * [37] Johansson J, Nation P, Nori F. 2012 QuTiP: An open-source Python framework for the dynamics of open quantum systems. Computer Physics Communications 183, 1760–1772. (https://doi.org/10.1016/j.cpc.2012.02.021)
  * [38] Johansson J, Nation P, Nori F. 2013 QuTiP 2: A Python framework for the dynamics of open quantum systems. Computer Physics Communications 184, 1234–1240. (https://doi.org/10.1016/j.cpc.2012.11.019)
  * [39] Lunghi A. 2022 Toward exact predictions of spin-phonon relaxation times: An ab initio implementation of open quantum systems theory. Science Advances 8, eabn7880. (10.1126/sciadv.abn7880)
  * [40] Atzori M, Tesi L, Morra E, Chiesa M, Torre R, Sorace L, Sessoli R. 2016 Room-temperature quantum coherence and Rabi oscillations in organometallic single-ion magnets: the role of anharmonic vibrations. Phys. Rev. Lett. 117, 227202.