Spin–phonon coupling and magnetic relaxation in single-molecule magnets

Electron–phonon coupling is important in many physical phenomena, e.g. photosynthesis, catalysis and quantum information processing, but its impacts are difficult to grasp on the microscopic level. One area attracting wide interest is that of single-molecule magnets, which is motivated by searching for the ultimate limit in the miniaturisation of binary data storage media. The utility of a molecule to store magnetic information is quantified by the timescale of its magnetic reversal processes, also known as magnetic relaxation, which is limited by spin–phonon coupling. Several recent accomplishments of synthetic organometallic chemistry have led to the observation of molecular magnetic memory effects at temperatures above that of liquid nitrogen. These discoveries have highlighted how far chemical design strategies for maximising magnetic anisotropy have come, but have also highlighted the need to characterise the complex interplay between phonons and molecular spin states. The crucial step is to make a link between magnetic relaxation and chemical motifs, and so be able to produce design criteria to extend molecular magnetic memory. The basic physics associated with spin–phonon coupling and magnetic relaxation was outlined in the early 20th century using perturbation theory, and has more recently been recast in the form of a general open quantum systems formalism and tackled with different levels of approximations. It is the purpose of this Tutorial Review to introduce the topics of phonons, molecular spin–phonon coupling, and magnetic relaxation, and to outline the relevant theories in connection with both the traditional perturbative texts and the more modern open quantum systems methods.


Lattice-dynamics calculations
Lattice-dynamics calculations were performed using pseudopotential plane-wave densityfunctional theory (DFT), as implemented in the Vienna Ab initio Simulation Package (VASP) code, 1 in conjunction with the supercell finite-displacement approach implemented in the Phonopy code. 2 Vibrational modes of CO 2 Calculations were performed on an isolated CO 2 molecule placed at the centre of a large cubic box, with an initial 15Å between periodic images, using the PBE exchange-correlation functional 3 with the DFT-D3 dispersion correction (i.e. PBE-D3). 4 The ion cores were modelled using projector augmented-wave (PAW) pseudopotentials 5,6 with the C and O 2s/2p electrons in the valence region. The valence Kohn-Sham wavefunctions were represented in a plane-wave basis with an 800 eV kinetic-energy cutoff and sampled at the Γ point (i.e. a single k-point). The atomic positions were optimised while keeping the cell shape and volume fixed, and the phonon frequencies and eigenvectors at the Brillouin-zone centre (q = Γ) were computed using the finite-displacement routines in VASP with a displacement step of 10 −2Å . To ensure accurate forces, the PAW projection was performed in reciprocal space and an enhanced charge-density grid with 8× as many points as the standard grids was used to represent the augmentation charges.
was fully optimised to tolerances of 10 −8 eV on the electronic total energy and 10 −2 eVÅ −1 on the forces. The second-order interatomic force constants (IFCs) were determined in a 4 × 4 × 4 supercell of the primitive cell (128 atoms) using the default displacement step of 10 −2Å , and the k-point sampling reduced accordingly. The PAW projection was performed in reciprocal space and an additional charge-density grid with 8× as many points as the standard grids was used to represent the augmentation charges. The atom-projected phonon density of states was obtained by interpolating the phonon frequenies and eigenvectors onto a uniform Γ-centered q-point mesh with 48×48×48 subdivisions. The phonon dispersion was obtained by interpolating the frequencies along strings of q-points passing through the q = L, Γ and X high-symmetry wavevectors in the F m3m Brillouin zone. Non-analytical corrections to the phonon dispersion to account for the splitting of the longitudinal and transverse modes close to q = Γ were included using the approach of Gonze et al., 9,10 with the required high-frequency dielectric constant ε ∞ and Born charge tensors Z * determined using the density-functional perturbation theory (DFPT) routines in VASP. 11 Phonon dispersion and density of states of crystalline NH 3 Calculations were performed on the cubic crystaline phase of NH 3 (spacegroup P 2 1 3) using PBE-D3. The ion cores were modelled using PAW pesudopotentials with the H 1s and N 2s/2p electrons in the valence region. The valence electronic structure was modelled using a plane-wave basis with an 800 eV cutoff and a Γ-centered 2×2×2 k-point grid. The initial structure was taken from the Materials Project database 12 (mp-29145) and has four NH 3 molecules (16 atoms) in the primitive cell. The atomic positions and unit-cell volume were optimised to tolerances of 10 −8 eV on the total energy and 10 −2 eVÅ −1 on the forces. The second-order IFCs were determined in a 4×4×4 supercell with 1,024 atoms using a 10 −2Å displacement step, and the k-point sampling reduced accordingly. The PAW projection was performed in reciprocal space and an enhanced charge-density grid with 8× as many points as the standard grids was used to represent the augmentation charges. The phonon density of states was evaluated on a 24 × 24 × 24 q-point grid, and the dispersion wa computed along a path including the high-symmetry X, Γ, M and R wavevectors in the P 2 1 3 Brillouin zone. In order to classify the phonon modes, the modes at q = Γ were visualised using the MolecularCrystalPhononAnimation code. 13

Spin-phonon coupling formalism
As noted in the main text, the phonon position operatorŝ q∈q * jV In the last line we have introduced a new set of phonon creation operatorŝ which are defined for q ∈ q * , and the Hermitean mode position operatorsX qj =b qj +b † qj .
The spin-phonon couplings ∂Ĥ S /∂X qj can be calculated in terms of the cartesian dis-placements as if q ∈ q + , and The Γ point and the q-points at the edge of the Brillouin zone need to be handled separately. In that case, Q qj is already Hermitean, and we defineX qj = 2ω qj h Q qj = a qj +â † qj , such that the spin-phonon coupling operator simply becomes if q / ∈ q * . Note that for the practical implementation of Equation 6 we adjust the arbitrary complex phase of the eigenvector W (qj) to obtain an overall real quantity for the displacement amplitudes W α κ (qj) exp (iq · r κl ). In analogy with the re-definitions above, second order coupling parameters ∂ 2Ĥ S ∂X qj ∂X q ′ j ′ can be expressed in the same basis of normal mode operatorsX qj .
Since the transformation defined in Equation (3) between the set of mode operatorsâ qj andb qj is unitary, the canonical commutation relations are preserved, e.g.
This allows us to easily calculate the bath correlation functions that we will use to write down the quantum master equations for the reduced density matrix of the spin system in the next section. In particular, the two-time bath correlation function at thermal equilibrium, which determines the one-phonon (Orbach) rates, is given by wheren qj = 1/(eh ω qj /k B T − 1) is the Bose-Einstein occupation number. It can be shown that the four-time equilibrium correlation functions entering the Raman-I and Raman-II rates can be decomposed in terms of two-time correlations functions as This fact only relies on the operatorsX qj being Hermitian linear combinations of creation/annihilation operators, and on the Gaussian property of thermal states.

Quantum Master Equations
In this section we apply the time-convolutionless (TCL) expansion 14 of the generator of the reduced dynamics to obtain 2 nd and 4 th order quantum master equations describing the dynamics of a spin system coupled to a very large vibrational bath in thermal equilibrium.
We consider both linear and quadratic couplings to a harmonic environment, in order to obtain expressions for the Orbach, Raman-I and Raman-II magnetic relaxation rates. The

system-bath coupling Hamiltonian is thuŝ
are the first and second derivatives of the spin HamiltonianĤ S with respect to modes i and j at the equilibrium geometry. For ease of notation, we seth = 1 and we drop the q-point dependence of the modes, letting a single mode index i run over both phonon bands and q-points within the first Brillouin zone. We also assume, without loss of generality, that every single term in the sums in Equation (9) is Hermitian, based on the results shown in the previous section.
Up to 4 th order in the system-bath coupling, the spin reduced density matrix ρ evolves according to the quantum master equation The 4 th order term consists of two contributions K 4 = K I 4 + K II 4 , describing the effect of the linear spin-phonon coupling to 4 th order (K I 4 ) and the 2 nd order contribution of the quadratic spin-phonon coupling (K II 4 ).
The matrix element of the generator K ab,cd describes the influence of ⟨ψ c | ρ |ψ d ⟩ on the time evolution of ⟨ψ a | ρ |ψ b ⟩. If we are only interested in relaxation rates for spin populations, we can focus just on the terms K f f,ii , which represents the transition rate between two eigenstates of the system Hamiltonian H S at energies E i and E f . In general, population dynamics is influenced by the presence of coherences, i.e. K f f,ab ̸ = 0 for a ̸ = b. However, one often invokes the secular approximation 14 (or rotating wave approximation), which leads to a decoupling of populations from coherences, provided the system has a non-degenerate spectrum (in this context, "non-degenerate" implies ω if ≫ 2πτ −1 , which is easily satisfied for magnetic molecules in small magnetic fields of a few Oe). If the spectrum is degenerate, it can be shown that the only coherence-to-population terms K f f,ab surviving the secular approximation are the ones involving coherences between degenerate states, i.e. E a = E b .
The secular approximation amounts to neglecting oscillatory terms in the interaction picture representation of the generator K, based on the idea that energy gaps between system eigenstates correspond to frequencies that are much faster than the relaxation dynamics, and therefore average to zero on the long timescale of magnetic relaxation.

Orbach rates
To lowest order in the spin-phonon coupling, the rate of change of the spin reduced density matrix is given by the 2 nd order TCL expansion (TCL2) in the linear system-bath coupling, which we call K 2 ρ. In the interaction picture, whereρ(t) = e iĤ S t ρ(t)e −iĤ S t , this is given by (h.c. means Hermitian conjugate) which is widely known as the Redfield Equation.
The function c i (t − t 1 ) represents the two-time bath correlation function calculated on a thermal equilibrium state, i.e.
where ω i is the mode frequency andn i = 1/(e ω i /k B T −1) is the occupation number at thermal equilibrium.
We introduce the jump operators between system energy levelŝ which allow us to decompose the time evolution of system operators aŝ This decomposition allows us to separate all the time-dependent factors in Equation (11) and perform the time integral. Transforming back to the Schrödinger picture, we obtain where the rate R i (ω 1 ) is given by Note that we have introduced two approximations. First, we have taken the upper limit of integration to go to infinity, i.e. t → ∞. This corresponds to only focusing on the long-time limit of the Markovian dynamics. Second, we are deliberately neglecting the imaginary part of the rates. This is because we are only interested in the dissipative relaxation dynamics induced by the real part, as opposed to the oscillatory dynamics described by the imaginary part. The limit t → ∞ gives rise to selection rules for transitions between spin states depending on the phonon energies, based on the relation Substituting Equation (12) into the expression for R i (ω 1 ), and calculating the contribution to ⟨ψ f | K 2 ρ |ψ f ⟩ stemming from ⟨ψ i | ρ |ψ i ⟩, we obtain the Orbach relaxation rates where we introduced the notation V ab i = ⟨ψ a |V i |ψ b ⟩ and ω ab = E a − E b . Extension to modes with finite linewidth is straightforward, and simply amounts to integrating the rate expression over a continuous distribution of mode energies described by a lineshape function (see main text).

Raman-II rates
The Raman-II relaxation mechanism is captured by the TCL2 expansion in the quadratic spin-phonon coupling in Equation (9). In this case, we need to account for the non-vanishing expectation value of the quadratic spin-phonon coupling at equilibrium, i.e. ⟨X iXj ⟩ eq = δ ij c i (0) ̸ = 0. We do so by adding and subtracting this expectation value to the bath operators and incorporating it into a redefined system HamiltonianĤ S + 1 2 ijŴ ij ⟨X iXj ⟩ eq . We rewrite the residual quadratic spin-phonon coupling as where the index µ runs over all mode pairs (i, j) (including i = j) and we have introduced V µ =Ŵ ij /2 andX µ =X iXj −⟨X iXj ⟩ eq . These definitions allow us to exploit the analogy with the TCL2 expansion for the linear coupling carried out in the previous section. Replacinĝ V i →V µ in Equation (15) andX i →X µ in Equation (12), yields the Raman-II rates.
The two-time correlation function for the operatorsX µ can be readily obtained from Equation (8), yielding Following the same steps outlined in the previous section for Orbach rates, we obtain the Raman-II rates
We can now calculate the Raman-I rate for a transition between initial and final states ψ i and ψ f by taking the expectation value ⟨ψ f | K 4 ρ |ψ f ⟩ and collecting all terms proportional to ⟨ψ i | ρ |ψ i ⟩. Upon expanding all the nested commutators in Equation (23), the TCL4 generator takes the form of a sum of terms where four different spin-phonon coupling operators act on the spin density matrix, from the left, right, or any possible combination, i.e. terms of the formVVVV ρ,VVV ρV ,VV ρVV ,V ρVVV , and ρVVVV , and with all possible permutations of mode indices i and j. Grouping together all terms with the same dependence on the operatorsV i , we obtain ij (ω af , ω ba , ω ib ) + R (2) ji (ω ab , ω bi , ω if ) + R ji (ω f i , ω ba , ω ib ) + R ij (ω af , ω ba , ω ib ) − R (2) ji (ω bi , ω ab , ω if ) − R ji (ω f i , ω ib , ω ba ) − R (4) ji (ω bi , ω if , ω ab ) (37) r (4) ij (ω ai , ω bf , ω ib ) − R (4) ij (ω f b , ω ia , ω bi ) (38) ij (ω bi , ω ia , ω af ) + R (4) ij (ω bi , ω ia , ω f b ). (39) The first three lines in Equation (34) correspond to terms of the formVVV ρV andV ρVVV , while lines four and five correspond toVV ρVV , while there is no contribution from terms of the form ρVVVV orVVVV ρ. All five terms in Equation (34) can contribute to population dynamics via the effective single-phonon transitions arising from Equations (28) and (29), and thus can be seen as higher order corrections to the Orbach process. However, it can be shown that only the last two lines contribute to genuine two-phonon Raman transitions.
Focussing only on these terms, we can finally write the Raman-I rate as