How do phonons relax molecular spins?

See allHide authors and affiliations

Science Advances  27 Sep 2019:
Vol. 5, no. 9, eaax7163
DOI: 10.1126/sciadv.aax7163


The coupling between electronic spins and lattice vibrations is fundamental for driving relaxation in magnetic materials. The debate over the nature of spin-phonon coupling dates back to the 1940s, but the role of spin-spin, spin-orbit, and hyperfine interactions has never been fully established. Here, we present a comprehensive study of the spin dynamics of a crystal of Vanadyl-based molecular qubits by means of first-order perturbation theory and first-principles calculations. We quantitatively determine the role of the Zeeman, hyperfine, and electronic spin dipolar interactions in the direct mechanism of spin relaxation. We show that, in a high magnetic field regime, the modulation of the Zeeman Hamiltonian by the intramolecular components of the acoustic phonons dominates the relaxation mechanism. In low fields, hyperfine coupling takes over, with the role of spin-spin dipolar interaction remaining the less important for the spin relaxation.


Spin 1/2 systems represent the fundamental prototype of magnetic materials, and the understanding of their properties is pivotal for the rationalization of more complex magnetic compounds. Their study has deep roots in the early days of quantum mechanics, and despite their basic nature, they still represent a very rich quantum playground with nontrivial and elusive dynamical properties. The spin lifetime in insulating materials is essentially limited by the interaction between spins and lattice vibrations, namely, the spin-phonon coupling. At finite temperature, spins can absorb/emit one or multiple phonons from/into the lattice and relax to an equilibrium state. The detailed understanding of this process at the first-principles level represents a long-standing question in physics and chemistry and goes well beyond the fundamental theory aspect.

Spin-lattice relaxation is key in several fields. For instance, the efficiency of magnetic resonance imaging contrast agents (1) and the spin lifetime of single-molecule magnets (2) are determined by the magnetization decay rate of paramagnetic elements. Turning to the main focus of this work, in both molecular and solid-state qubits, spin-lattice relaxation sets the upper limit for the coherence time (36), and the engineering of this interaction is a primary challenge in the spin-based quantum computing field. The synthetic versatility of molecular compounds offers an intriguing route to the tailoring of spin-phonon coupling but needs to be supported by a rational understanding of the physical principles governing spin relaxation.

The debate over the role of phonons in the relaxation of electronic spins can be traced back to the 1940s when it was discussed for the first time in the context of transition-metal chemistry. Relaxation pathways occurring through the modulation of spin dipolar interactions or the modulation of the d-electrons’ crystal field were first pointed out by Waller (7) and Van Vleck (8, 9), respectively. More recently, the role of hyperfine interaction has also been discussed (10). Van Vleck’s mechanism remains to date the most commonly accepted explanation for the microscopic origin of spin-lattice relaxation of electronic spins. Several questions, however, remain unanswered. In particular, early models, being phenomenological, could not entirely address the nature of the coupling between phonons and spins and simply ascribe acoustic vibrations as responsible for this interaction. Molecular compounds are inherently complex, and a rational description of spin relaxation in terms of molecular motions is still to be developed.

First-principles calculations represent the perfect ground to provide unbiased answers to these fundamental questions. The possibility to predict the spin relaxation times in magnetic materials without introducing any phenomenological parameter is also of fundamental importance from a materials-design perspective. Quantum chemistry and related disciplines have been already proven invaluable tools for providing insights into the physics of new chemical systems (11), and they have been used as a screening method to predict new materials with tailored properties (12). These computational strategies represent a rich opportunity for the field of quantum informatics. The development of a first-principles framework able to predict the spin relaxation time is the first step in that direction. In this work, we offer a theoretical and computational explanation of the nature of the atomistic processes that lead to the relaxation of a molecular electronic S = 1/2 system.

Our approach builds on previous contributions from ourselves (13, 14) and others (1518) and substantially extends the state of the art in the field. We include in the formulation the contribution of phonons from the entire Brillouin zone and the hyperfine and dipolar spin-spin interactions, both at a static and dynamical level. The description of all the fundamental spin relaxation channels occurring at the first order of perturbation theory is therefore here complete. Our method is purely first principles and provides a rationalization of spin-lattice relaxation without any previous knowledge of the system’s properties other than its crystal structure. The S = 1/2 system investigated is the crystal of the Vanadyl complex VO(acac)2, with acac=acetylacetonate (19). Vanadyl-based molecular qubits represent archetypal systems for room temperature quantum-computing applications (20).

We demonstrate that Van Vleck’s mechanism is the dominant first-order relaxation channel for S = 1/2 molecular spins in a high magnetic field, while hyperfine contributions become dominant in a low-field regime. Our calculations show that acoustic phonons are not rigid molecular translations, shedding light on the origin of spin-phonon coupling in solid-state spin 1/2 compounds and thus solving an 80-year-old controversy. Last, we discuss the implications of our results for the field of molecular qubits and magnets and discuss a possible way forward for the field. In particular, we suggest how a tighter agreement with experiments can be achieved by developing the theory beyond the first order.


First-principles spin dynamics

The quantum dynamics of even a single spin is, in principle, entangled with the dynamics of all the other spin and lattice degrees of freedom that it is interacting with. More explicitly, it is driven by the total Hamiltonian Ĥ=Ĥs+Ĥph+Ĥsph, where the three terms are the spins and phonons Hamiltonian and the spin-phonon coupling, respectively. It is convenient to think of the problem as composed of two parts: the simulation of the sole spin degrees of freedom as an isolated system and the interaction of this ensemble with a thermal bath, namely, with the phonons.

The spin Hamiltonian of Eq. 1 describes the fundamental interactions taking place within an ensemble of Ns spinsĤs=iNsβiBg(i)S(i)+12ijNsS(i)D(ij)S(j)(1)where the i-th spin, S(i), interacts with an external magnetic field B through the gyromagnetic tensor, βig(i). Equation 1 can account for both electronic and nuclear spins on the same footing. Thus, the spin-spin interaction tensor, D, may coincide with the point-dipole interaction, Ddip, or the hyperfine tensor, A, depending on whether the interaction is among electronic and/or nuclear spins. The state of the spin system can be described in term of the spin density matrix, ϱ̂s, whose dynamics is regulated by the Liouville equationdϱ̂sdt=i[Ĥs,ϱ̂s](2)

Once the eigenstates and eigenvalues, {Ea}, of the spin Hamiltonian are known, it is possible to integrate Eq. 2 to obtain an expression for the evolution of ρabs in the so-called configuration interactionρabs(t)=eiωab(tt0)ϱ̂abs(t0)(3)where ωab = (EaEb)/ℏ. Despite its simple form, Eq. 3 hides a high level of complexity, which originates from the high dimension of the Hilbert space it acts on. Several numerical schemes exist to reduce the computational costs associated with this problem (21, 22), and these will be the subject of future investigations. Here, we have decided to retain the full Hilbert space description and to restrict the study to a small number of interacting spins.

The basic first-principles theory for phonon-driven spin relaxation has been derived before (13, 14). Here, we extend it by including acoustic phonons and the phonon reciprocal space dispersions. For a periodic crystal, defined by a set of reciprocal lattice vectors, q, and by N atoms in the unit cell, the lattice’s dynamics can be described in terms of periodic displacement waves (phonons), Q̂αq, with frequency ωαq and obeying to a harmonic HamiltonianĤph=αqωαq(n̂αq+12)(4)

Here, n̂αq is the phonon’s number operator. The spin-phonon interaction, responsible for the energy exchange between the spins and the lattice, is generally assumed to operate in the weak coupling regime and at the first-order level of approximation. Thus, the spin-phonon coupling Hamiltonian is written asĤsph=αq(ĤsQαq)Q̂αq(5)

Building on the weak coupling approximation, a compact expression for the spin quantum equation of motion can be derived. As explained in detail in the first section of the Supplementary Materials, the Born-Markov approximation can be used to derive the nonsecular Redfield equations. Besides the weak coupling regime, this level of approximation assumes the phonon dynamics to develop at a much shorter time scale than that of the spin relaxation, as it is usually the case. A possible shortcoming of this approximation is generally encountered when a spin-phonon bottleneck is present. This effect mainly arises from the small thermal conductivity of phonons, which are unable to dissipate the excess of energy absorbed from the relaxing spin degrees of freedom (19). However, if no spin-phonon bottleneck is at play, the Markov approximation can be assumed to be valid, and the phonon dynamics are always considered at the thermal equilibrium. From these assumptions follow the nonsecular Redfield equations that can be used to study the spin dynamics under the effect of a phonon bath (23). In this framework, the spin density matrix ρabs(t) evolves according to the equationdρabs(t)dt=cdei(ωac+ωdb)tRab,cdρcds(t)(6)

The transition rates between the elements of the density matrix are of the formRab,cdπ22αqVacαqVdbαq(G(ωdb,ωαq)+G(ωca,ωαq))(7)where Vabαq=aĤsQαqb is the matrix element of the spin-phonon Hamiltonian in the spin Hamiltonian eigenfunctions basis and Gij, ωαq) is the Fourier transform of the single phonon correlation function at finite temperature. The entire expression for Rab, cd appearing in Eq. 7 is provided in the Supplementary Materials. For harmonic lattices, Gij, ωαq) is defined asG(ωij,ωαq)=n¯αqδ(ωαqωij)+(n¯αq+1)δ(ωαq+ωij)(8)where n¯αq is the Bose-Einstein population at a temperature T. When lattice anharmonic interactions are important, the finite lifetime of phonons 2πΔαq−1 must be included in the model, and Gij, ωαq) assumes the usual Lorentzian line shapeG(ωij,ωαq)=1π[ΔαqΔαq2+(ωijωαq)2n¯αq+ΔαqΔαq2+(ωij+ωαq)2(n¯αq+1)](9)

Once Eq. 6 is solved, the magnetization dynamics for the i-th spin can be computed by the canonical expression M(i)=Tr{ρ̂s(t)S(i)}.

Spin phonon coupling in a molecular 1/2 electronic spin

To make the physics of Eq. 7 more transparent, we hereafter provided a study of the spin-phonon coupling terms, Vab, for the molecular qubit VO(acac)2. This V4+ complex crystallizes with a triclinic primitive cell containing two inversion symmetry related molecular units (24), as reported in Fig. 1. Each molecule bears a single electronic S = 1/2 spin in addition to the I = 7/2 nuclear spin of 51V. Let us consider one electronic spin Si interacting with the rest of the crystal electronic spins Sj and the local vanadium nuclear spin Ii. Given the definition in Eq. 5, the spin-phonon coupling Hamiltonian will contain three distinguished contributions: a first intramolecular spin-phonon coupling due to the modulation of the Landé tensor, g; a second intramolecular coupling coming from the modulation of the hyperfine interaction, A; and an interelectronic spin interaction originating from the modulation of the dipolar terms, Ddip, namelyĤs(i)Qαq=βiBg(i)QαqS(i)+S(i)A(ii)QαqI(i)++jNsS(i)Ddip(ij)QαqS(j)(10)

Fig. 1 VO(acac)2 structure and spin phonon coupling distributions.

(A) The geometrical structure of the two VO(acac)2 molecular units inside the crystal’s unit cell. Vanadium atoms are represented in pink, oxygen in red, carbon in green, and hydrogen in white. (B) The spin-phonon coupling distribution relative to the Zeeman energy as function of the phonons’ frequency. (C) The spin-phonon coupling distribution relative to the dipolar spin-spin energy as function of the phonons’ frequency. (D) The spin-phonon coupling distribution relative to the hyperfine energy as function of the phonons’ frequency.

A 3 × 3 × 3 super cell containing 1620 atoms is optimized at the density functional theory (DFT) level, and it is used to compute the lattice vibrational properties. The molecular optimized structure has been further used for the calculation of all the spin-phonon coupling coefficients appearing in Eq. 10. The details about the protocol used for the phonons and spin-phonon calculations can be found in the Materials and Methods. All these interactions break the single-spin time-reversal symmetry and are potentially active in intra–Kramer doublet spin relaxation, but the physics beyond the three processes is quite different. To make these differences more evident, we have calculated the spin-phonon coupling squared norms Vsph2 (defined in the Materials and Methods) associated to each phonon mode and their corresponding distributions. The results, reported in Fig. 1 (B to D), show notable differences both in qualitative and quantitative terms.

The anisotropy of the spin Hamiltonian is the fingerprint of the dependence of spin degrees of freedom on atomic positions, and the same interactions contributing to magnetic anisotropy are also contributing to the derivatives of Eq. 10. The tensors g and A have anisotropic components coming from spin-orbit coupling, and other interactions which are localized on the Vanadium center. These short-ranged interactions are only influenced by localized intramolecular vibrations and local rotations affecting the first coordination sphere. In this case, high-energy phonons are still operative. Conversely, dipolar interactions act between different molecules and are both nonlocal and long wavelength in nature. These are expected to be prominently modulated by molecular translations. Their interaction thus vanishes at high frequency.

VO(acac)2 crystal spin dynamics

The spin-phonon coupling coefficients discussed in the previous section have been used together with Eq. 6 to simulate the longitudinal relaxation time τ of three different systems: one single electronic spin, one electronic spin interacting with the V nuclear spin, and two interacting electronic spins. Tests including more than two coupled electronic spins are reported in the Supplementary Materials and show a very small dependence of the longitudinal relaxation time, τ, on spins farther than the first-neighbor ones. In all simulations, spins are interacting with all the phonons of the periodic crystal, calculated by integrating the Brillouin zone with homogeneous grids up to 64 × 64 × 64 k-points. Note that this system exhibits spin-phonon bottleneck effects, and therefore all the comparisons to experimental results have been done with respect to the finest ground crystals studied by Tesi et al. (19). In this condition, the spin-phonon bottleneck is mostly quenched.

Figure 2 shows the spin relaxation time as function of the external magnetic field B at 20 K for a single electronic spin coupled to the V nuclear spin and a purely harmonic phonons bath. At this temperature, spin-phonon bottleneck effects are already vanished (19). The combined effect of the Zeeman and hyperfine modulation, described by the black line and dots, shows a field dependence approaching B4 in the high-field limit (B > 5 T), where τ starts to converge to the experimental value (19). The green curve in Fig. 2, relative to the single electronic spin contribution, shows that in this regime the leading mechanism is the modulation of the Zeeman energy. In the low-field regime, the relaxation time as a function of field approaches a plateau where it is about three orders of magnitude longer than the experimental value. The red curve in Fig. 2 shows that the rate-determining mechanism in this regime is the modulation of the hyperfine Hamiltonian. Numerical tests concerning the effects of the numerical noise in the calculation of the frequencies and the spin-phonon coupling parameters are reported in the Supplementary Materials and prove the robustness of our results. In particular, we note that the slight overestimation of τ in the high-field regime (B ∼ 10 T) might come from the overall overestimation of the calculated vibrational frequencies with respect to the experimental ones. It is also important to remark that the energies at play in the low-field regime (B < 1 T) are extremely small if compared with the phonons’ frequencies and the simulations. Thus, a more sophisticated Brillouin zone integration scheme might be needed to obtain a more robust estimate of the relaxation times in this field regime. Nonetheless, the calculated relaxation time for low fields is orders of magnitude longer than the experimental one, suggesting the presence of higher-order relaxation mechanisms at play in this regime. Figure 3 displays the simulated and experimental temperature dependence of the spin relaxation time in a field of 5 T.

Fig. 2 Spin relaxation time as function of the B field for one electronic spin coupled to one nuclear spin.

The relaxation time, τ, in milliseconds and T = 20 K as a function of the external field in Tesla is reported for the simulations of one electronic spin coupled to one nuclear spin and relaxing due to the phonon modulation of the Zeeman and hyperfine energies (black line and dots). The contribution coming from the sole hyperfine energy modulation is represented by a red line and dots, while the sole Zeeman contribution is reported by the green line and dots. The experimental relaxation time as extracted from AC magnetometry (19) is also reported (blue dots and line).

Fig. 3 Spin relaxation time as function of the temperature T.

The relaxation time, τ, in milliseconds at B = 5 T as a function of the temperature is reported for the simulations of one electronic spin coupled to one nuclear spin and relaxing due to the modulation of the Zeeman and hyperfine energies by harmonic phonons (black line and dots), phonons with 1 cm−1 of line width (red line and dots), and phonons with a linearly T-dependent phonon lifetime (green line and dots). A coefficient γ = 0.1 cm−1/K is chosen. The experimental relaxation time as extracted from AC magnetometry (19) is also reported (blue dots and line).

Assuming a perfectly harmonic model for the vibrational bath, we simulate a T−1 behavior. Our result is in agreement with what is expected from a first-order approximation to the spin-phonon coupling and a harmonic lattice dynamics but deviates from the experimental results, which instead show τ ∝ Tn (with n > 1) (19). Deviations from the T−1 power law can be considered as fingerprints of higher-order processes taking place. Among the several multiphonon processes that could influence the spin dynamics, here we investigate the role of a finite and T-dependent phonon lifetime. In molecular magnets, the phonons’ lifetime is mainly limited by phonon-phonon scattering events, and it can be estimated within perturbation theory. However, the computational requirements for this kind of calculation are extremely high, and the first-principles simulation of this quantity is only possible for very simple systems. While we postpone a detailed study of these processes to future work, here we provide a phenomenological estimation of their effects. First, we show the effects of a finite phonon lifetime by using Eq. 9 instead of Eq. 8 in evaluating the spin relaxation rate. By choosing a Lorentzian phonon linewidth of 1 cm−1, the relaxation rate significantly drops. This is not unexpected, given the fact that a Lorentzian profile has longer tails than the Gaussian one used to simulate a δ-like function (see Materials and Methods). This means that the Lorentzian is expected to converge to a Dirac’s ∆ function for a smaller line width. The phonon line width can also be T dependent. We simulate this effect by assuming the linear relation ℏΔαq = γT between the line width and the temperature. This model is justified by perturbation theory arguments (25). Figure 3 shows that this model has the effect of increasing the τ versus T power law of one unit and, depending on the size of the coefficient γ, speeding up the overall spin relaxation rate. Although the introduction of a T-dependent phonon lifetime offers a better agreement with experimental results, only a rigorous estimation of these quantities by first principles will make a final assessment of the role of phonons’ lifetime and of other multiphonon processes possible. However, note that the harmonic lattice case yields an upper limit to the spin relaxation time. Our simulations also show a residual T-independent process in the T → 0 limit. In this regime, due to the absence of populated phonon states, the only possible relaxation pathway is provided by the T-independent spontaneous phonon emission from a spin excited state. This phenomenon has been recently observed in Nitrogen-Vacancy (N-V) centers (17).

Figure 4 presents the spin relaxation time as a function of the external magnetic field at 20 K for two coupled electronic spins. At high fields (B > 1 T), the simulated relaxation dynamics, including both inter- and intraspin direct relaxation mechanisms (black curve and dots in Fig. 3), shows no significant difference from that of two isolated spins relaxing through the modulation of the Zeeman interaction (green curve and dots in Fig. 3). The relaxation time due to the sole dipolar contribution (red curve and dots in Fig. 3) becomes predominant only at low fields, and it is found to be two orders of magnitude slower than that associated to the hyperfine coupling.

Fig. 4 Spin relaxation time as function of external field for two coupled electronic spins.

The relaxation time, τ, in milliseconds as a function of the external field in Tesla is reported for the simulations of two electronic spins relaxing due to the phonon modulation of the Zeeman and dipolar energies (black line and dots). Green line and dots represent the relaxation time of two isolated spins, where only the Zeeman energy is modulated by phonons. The contribution coming from the sole dipolar energy modulation is represented by a red line and dots. The experimental relaxation time as extracted from AC magnetometry (19) is also reported (blue dots and line).

To understand the nature of the interaction between S = 1/2 spins and the lattice, it is now necessary to look at the nature of the phonons involved in the process. It has been shown that, in the absence of spin-spin interactions, spin relaxation can only occur through the modulation of the spin Hamiltonian by local rotations and intramolecular distortions (14). However, for a spin 1/2 in reasonable external fields, the Zeeman spin splitting (up to a few cm−1) is much smaller than the first Γ-point optical mode, here approximately 50 cm−1. Energy conservation (see Eq. 8) leaves only acoustic phonons as candidates for the spin-phonon interaction, suggesting that no energy exchange between lattice and spin is possible under these conditions. The solution of this conundrum is provided by the analysis of the phonons by means of their decomposition into local molecular translation, local molecular rotations, and intramolecular distortions (14). The results of this analysis, carried out over the phonons in the entire Brillouin zone, are summarized in Fig. 5, where the total and decomposed phonon density of states are reported.

Fig. 5 Phonon density of states.

The total phonon density of states (DOSs) as a function of the frequency is reported in black. The total phonon density of states was also decomposed in a pure translational contribution (red line), a rotational contribution (blue line), and an intramolecular contribution (green line), all relative to a single molecule inside the unit cell. The inset shows the details of the density of states in the low-energy part of the spectrum. The Brillouin zone was integrated with a uniform mesh of 643 points. The reported density was smeared with a Gaussian function with breadth of 1.0 cm−1. a.u., arbitrary units.

The total phonon density of states shows the typical ~ω2 dependence at low frequency, where acoustic phonons dominate the vibrational spectra. However, the decomposition of these modes shows that at low frequency, the nature of the acoustic modes is far from being that of a pure molecular translation. A significant rotational and intramolecular contribution is present at frequencies corresponding to the energy levels’ Zeeman splitting considered in this work. These contributions provide an efficient relaxation pathway even for a single spin isolated from other magnetic centers.


The results of the first-principles calculations presented here are in agreement with Van Vleck’s interpretation of the spin relaxation in high fields, where a direct relaxation mechanism due to the Landé tensor modulation is the relevant relaxation pathway. Here, we have demonstrated the mechanism underlying the energy transfer between phonons and spins. The presence of an intramolecular contribution inside the low-energy acoustic phonons is of fundamental importance to open a relaxation channel that otherwise would be completely inactive due to the translational invariance of the spin-orbit coupling interaction.

The presence of intramolecular components in acoustic modes can be understood as a mixing between rigid reticular translations and soft molecular modes, and it can be used as a rationale to engineer solid-state qubits. More rigid molecular modes are expected to diminish sensibly the contamination of low-energy modes, therefore extending the spin lifetime. Such a synthetic strategy has been recently attempted (26, 27) on the basis of similar observations for high-spin molecular magnets (14).

We expect our findings to be of general importance for molecular systems. The relaxation mechanism due to the intramolecular components of the acoustic phonons depends on the presence of low-lying optical vibrations and on the locality of the A and g tensors. These features are common to any molecular crystal. The magnetic field intensity at which the different interactions overcome each other and the overall relaxation rate clearly depend on the specific features of the system under investigation. Nevertheless, here we show that it is possible to build a quantitative theory of the direct relaxation in electronic spins on the sole basis of first-principles calculations. The adaptation of the Redfield equations to the first-principles treatment of the spin-phonon interaction represents a completely general framework and can be readily applied to any solid-state magnetic system, including molecules in frozen solutions, surface-adsorbed molecules, and nonmolecular qubits, as defects in Si and color centers in diamond. Moreover, the same formalism works for many spins systems and giant spins (S > 1/2).

Regarding the possibility to quantitatively predict the spin relaxation time completely from first principles, further work is needed. The direct relaxation mechanism, dominant in the high-field regime, is well captured by the present approach. Given the large number of quantities that one needs to compute and the complete absence of any adjustable parameter, an agreement with experiments within approximately one order of magnitude has to be considered as an outstanding result. An even better agreement with experimental data in this regime can probably be achieved through the use of more accurate first principles methods. In particular, we have shown the importance of an accurate determination of the low-energy phonon frequencies. Details concerning the spin dynamics at low fields and the nature of the polynomial dependence of the relaxation rate at high T (19) remain partially elusive and suggest that additional higher-order processes may take place. Capturing these probably represents the most important challenge in this nascent field and that the inclusion of second-order spin-phonon coupling and phonon-phonon scattering processes will bring the present approach on a fully quantitative ground.

To conclude, we have presented a general and fully first-principles method to study spin-lattice dynamics in magnetic materials. We have predicted the correct spin relaxation time and field dependence for a solid-state molecular qubit in high external fields. In particular, we have given a full microscopic rationale of the spin relaxation in solid-state molecular qubits by ranking the three fundamental interactions at play among electronic and nuclear spins at the first order of perturbation theory. The method presented here can be readily extended to include higher-order processes such as Raman relaxation mechanism and phonon-phonon interactions. This will be the subject of future work.


Lattice dynamics

All the structural optimization and Hessian calculations were performed with the CP2K software (28) at the level of DFT with the Perdew-Burke-Ernzerhof (PBE) functional including Grimme’s D3 van der Waals corrections (29, 30). A double zeta-polarized MOLOPT basis set and a 600 Ry of plane-wave cutoff were used for all the atomic species. The comparison between simulated and experimental lattice parameters is available in the Supplementary Materials. All the translational symmetry–independent force constants have been computed by a finite difference approach with a 0.01-Å step. Φij(lm) is the force constant that couples the i-th atomic degrees of freedom in the lattice cell at the position Rl and the j-th atomic degrees of freedom in the lattice cell of position Rm. The dynamical matrix, K(q), at the q point, is built asKij(q)=lΦij(l0)eiqRl(11)

The eigenvalues of K(q) are ω2(q), while the eigenvectors L(q) defines the normal modes of vibration. The calculated phonons are in good agreement with those previously performed at the sole Γ-point (19).

Spin-phonon coupling coefficients

The ORCA software (31) was used for the computation of the g and A tensors for both equilibrium and distorted geometries. We used the basis sets def2-TZVP for V and O, def2-SVP for C and H, and a def2-TZVP/C auxiliary basis set for all the elements. For the calculation of the A tensors, the entire basis set was decontracted. The calculations of the g tensors were carried out at the CASSCF + NEVPT2 level of theory with a (1,5) active space and spin-orbit contributions included through quasi-degenerate perturbation theory. The calculations of the A tensors were performed at the DFT level with the PBE functional (29).

The spin phonon coefficients relative to the g and A tensors were calculated as numerical derivatives. Ten Cartesian displacements ranging from ±0.01 Å were used to estimate ∂g/∂Xis and ∂A/∂Xis, where Xis refers to the s Cartesian component of the i-th atom in the DFT-optimized unit cell. The g versus Xis and A versus Xis profiles were fitted with a fourth-order polynomial expression and set to zero if the fitting error on the linear term exceeded 7%. The spin-phonon coupling coefficients relative to the point dipole-dipole interaction were obtained by analytical differentiation. The Cartesian derivatives Ĥs/Xis were then projected onto the normal modes by means of the expression(ĤsQαq)=lNcellsisN,3NqωαqmieiqRlLisαq(ĤsXisl)(12)where Xisl is the s Cartesian coordinate of the i-th atom of N with mass mi, inside the unit-cell replica at position Rl, and Nq is the number of q points used. All the data regarding spin Hamiltonian parameters and their differentiation are reported in the Supplementary Materials. The spin-phonon coupling distributions were calculated starting from the spin-phonon coupling coefficients squared norm, defined asVsph2(ωα)=1Nqqij3(gijQαq)2(13)where gij are the matrix element of the g tensor. Vsph2(ωα) is analogously defined for the A and DDip tensors. The Dirac’s ∆ function appearing in Eq. 8 was evaluated as a Gaussian function in the limit for infinite q points and vanishing Gaussian breadth. A grid of 643 q points and a Gaussian breadth of 1 cm−1 were estimated to accurately reproduce this limit for all the temperature and field values investigated. Some results about these convergence tests are reported in the Supplementary Materials.


Supplementary material for this article is available at

Note S1. Derivation of the nonsecular Redfield equations

Note S2. Lattice parameters and spin Hamiltonian

Note S3. K points and phonon’s line width convergence

Note S4. Number of spins convergence

Note S5. Effect of numerical noise

Table S1. Lattice parameters for the VO(acac)2 molecular crystal.

Table S2. Spin Hamiltonian parameters for VO(acac)2.

Fig. S1. Spin relaxation time as function of external fields for different values of k points.

Fig. S2. Spin relaxation time as function of the Gaussian breadth for different external fields.

Fig. S3. Spin relaxation time as function of the number of coupled spins.

Fig. S4. Effect of numerical noise on spin relaxation time.

Spin-phonon coupling coefficients and ab initio calculation' input files

Reference (32)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We acknowledge L. Tesi and R. Sessoli for providing the original experimental data and the stimulating scientific discussions. We also acknowledge the MOLSPIN COST action CA15128. Funding: This work was sponsored by the Science Foundation Ireland (grant 14/IA/2624). Computational resources were provided by the Trinity Centre for High-Performance Computing (TCHPC) and the Irish Centre for High-End Computing (ICHEC). Author contributions: A.L. and S.S. conceived the project and discussed the results. A.L. developed the theory, the computational tools and ran the simulations. A.L. wrote the manuscript with inputs from S.S. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article