Temperature Dependence of Spin–Phonon Coupling in [VO(acac)2]: A Computational and Spectroscopic Study

Molecular electronic spins are good candidates as qubits since they are characterized by a large tunability of their electronic and magnetic properties through a rational chemical design. Coordination compounds of light transition metals are promising systems for spin-based quantum information technologies, thanks to their long spin coherence times up to room temperature. Our work aims at presenting an in-depth study on how the spin–phonon coupling in vanadyl-acetylacetonate, [VO(acac)2], can change as a function of temperature using terahertz time-domain spectroscopy and density functional theory (DFT) calculations. Powder THz spectra were recorded between 10 and 300 K. The temperature dependence of vibrational frequencies was then accounted for in the periodic DFT calculations using unit-cell parameters measured at two different temperatures and the optimized ones, as usually reported in the literature. In this way, it was possible to calculate the observed THz anharmonic frequency shift with high accuracy. The overall differences in the spin–phonon coupling magnitudes as a function of temperature were also highlighted showing that the computed trends have to be ascribed to the anisotropic variation of cell parameters.


■ INTRODUCTION
The quantum bit or qubit 1 is the basic element of quantum information theory (QIT); 2 it differs from its classical equivalent, the bit, since it can exploit the quantum superposition of the two states 0 and 1. Several systems are under investigation as potential qubits, e.g., electronic defects in diamond, 3 photons, 4 ion traps, 5 superconducting circuits, 6 and spins. 7−9 Molecular electronic spin qubits are good candidates since they are characterized by a large tunability of their electronic and magnetic properties through a rational chemical design. Among them, coordination compounds of light transition metals, such as Ti(III), V(IV), and Cu(II), 10−12 have demonstrated to be promising systems at least from the point of view of the superposition state's lifetime. This quantity is determined by the spin−spin relaxation time, T 2 , but it is usually quantified through the phase memory time, T m , that is the measurable lower limit of T 2 . 13 Another fundamental parameter in the evaluation of the qubit performance is the spin−lattice relaxation time, T 1 . 14 If it is too short, it limits T m , 11 preventing the implementation of more complex algorithms. In solid-state qubits, T 1 is closely connected with lattice vibrations, i.e., phonons. Vibrations perturb spin degrees of freedom through the modulation of orbital contributions and the presence of spin−orbit coupling. Therefore, the knowledge and characterization of phonons in solid-state qubit candidates are of crucial importance. A macroscopic crystal can present a huge number of atomic displacement degrees of freedom and, consequently, of vibrational modes. The use of space-group symmetry and the introduction of the reciprocal lattice allows a fundamental rationalization of this huge number of states. At the center of the Brillouin zone (BZ), the so-called Γ-point, i.e., for wavevector k = 0, the degrees of freedom are reduced to 3N, N being the number of atoms in the crystal unit cell so that the vibrational modes are 3N − 3 (3 are the rigid displacements of the unit cell, corresponding to zero frequency). Only the k = 0 vibrations can be measured with optical techniques at equilibrium; all of the remaining vibrational modes are classified, throughout the entire BZ, in 3N dispersion bands and are optically inaccessible. A peculiar feature of most molecular crystals is that intramolecular forces are much stronger than intermolecular ones so that the intramolecular vibrations in the crystal are almost unchanged from those of the isolated molecule. If the crystal unit cell contains n molecules, there are then 6n − 3 intermolecular phonons at k = 0, corresponding to rigid translations and rotations of the molecules. These "external" phonons have frequencies much lower (typically, less than 200 cm −1 ) than intramolecular vibrations. However, in some cases, where molecular vibrations of unusually low frequency are present, the separation of "internal" and "external" phonons is impossible and lowfrequency phonon modes have mixed character. 15 In the harmonic approximation, phonon modes are orthogonal to each other and their frequencies depend on temperature indirectly, as the thermal expansion of the lattice induces an elastic strain in the solid, which results in a shift (generally negative) of the phonon frequencies. The thermal expansion of crystals is a consequence of the anharmonicity of the crystal potential; however, it has been shown 16,17 that the effect of thermal strain can be dealt with in a quasiharmonic approximation, adopting an effective crystal Hamiltonian that maintains the quadratic dependence on the atomic displacements. On the other hand, the anharmonic terms of the potential energy expansion enter directly into the description of the phonon dynamics as they mix different phonons and make phonon relaxation possible, thus determining the homogeneous linewidths of vibrational spectral lines. Furthermore, the same anharmonic terms cause additional temperature dependence of the frequencies through the temperature-dependent occupation number of phonons.
From a computational point of view, the harmonic treatment of molecular vibrations cannot explain all of the fine details of the vibrational spectra and the physical properties of these crystals that depend on the temperature. 18 Density-functional perturbation theory or molecular dynamics simulations (from several ps to ns) 19 are out-of-reach. The direct method of calculation of phonons 20 here employed can be considered as a cheaper alternative. Unfortunately, even with this approach, the extension of lattice dynamics calculations beyond the second-order of expansion of the potential energy surface (PES) is severely limited by the enormous number of terms appearing in the expression of a large molecule; therefore, it is mandatory to assume some approximations when aiming to calculate higher orders of interaction.
In such a framework, an extensive characterization of the vibrational properties as a function of temperature becomes feasible and it is of paramount importance for the derivation of key chemical insights useful in the engineering of efficient molecular spin qubits. Indeed, as shown in other works, 10,21−27 the study of phonons provides insights into the spin relaxation pathways pointing toward the rationalization of new synthetic strategies to achieve more performing systems.
The aim of this work is an in-depth study of the spin− phonon (hereafter, Sp-Ph for brevity) coupling in vanadylacetylacetonate, [VO(acac) 2 ], both from a spectroscopic and theoretical point of view. Thanks to its simplicity (only 60 atoms in the unit cell), [VO(acac) 2 ] has got the physique du role as a perfect training ground for extensive DFT simulations to investigate the rich vibrational properties of a molecular crystal. 26 [VO(acac) 2 ], already studied by inelastic neutron scattering, 28 is here investigated by means of terahertz radiation in crystalline powder by using of an in-house built setup performing THz time-domain spectroscopy (THz-TDS). 29−31 In our computational approach, the temperature dependence of the vibrational frequencies is accounted for using the unit-cell parameters measured at different temperatures. The results of this study confirm the paramount importance of lowenergy phonons, where intramolecular and lattice characters are admixed. 15,25,26,32 ■ MATERIALS AND METHODS Crystal Structure. The [VO(acac) 2 ] metal center comprises a vanadyl ion VO 2+ coordinated by two chelating bidentate β-diketonate ligands ( Figure 1). The coordination geometry around the penta-coordinated V IV metal center is a distorted square pyramid, as already reported in previous studies. 33 This molecule crystalizes in a P1̅ space group (no. 2), for which the only point symmetry element is inversion. The asymmetric unit contains one molecule and the primitive unit cell contains two molecules, for a total of 60 atoms (see Table 1).
The crystallographic structure has been already solved using high-resolution data in 1995 34 at 294 K. Here, we present new data collected at 300 K, 1 300K , and 100 K, 1 100K , both on the same sample (Table S1). A good agreement between 300 and 294 K datasets is found, while the shrinking of cell parameters at 100 K is highlighted in Table 1.
Experimental Apparatus and Data Analysis. THz spectra were measured by time-domain spectroscopy using a table-top experimental setup working in a transmission configuration. It allows one to achieve a very high signal-tonoise ratio, higher than what is usually obtained in the standard far-infrared investigation (FTIR). Furthermore, the nontrivial data analysis allows one to decouple the signals arising from the multiple reflections inside the sample, i.e., the Fabry−Peŕot effect. The THz-TDS spectrometer is a home-built system based on optical laser pulses (T-light 780 nm fiber laser, MenloSystems) and low-temperature GaAs photoconductive antennas. The emitted THz radiation (pump) is collimated and focused on the sample by a couple of off-axis parabolic mirrors. The signal transmitted through the sample is again collimated and focused on the detector, a second photoconductive antenna, biased by a second optical pulse (probe). It results in a photocurrent that is measured by a lock-in amplifier and an acquisition board. The amplitude of the THz electric field is measured by varying the time delay between the pump and probe pulses so as to reconstruct its time profile. The homemade software in Matlab code was implemented to acquire the processed signal together with the reading of the delay line encoder and retraces the final time-dependent THz field. The whole THz setup is placed in a nitrogen-purged chamber to eliminate the strong absorption of water vapor. To improve the data quality, the effects of external perturbations during the acquisition (e.g., temperature fluctuations) measurements must be taken into account. For this purpose, measurements with and without a sample (reference) are cyclically repeated. We also implemented a new method of data analysis, described in detail in ref 29, that enables the extraction of the index of refraction, n, absorption coefficient, α, and thickness, d, of the material. This method is based on an iterative fitting process of the transmission parameters and it takes into account the Fabry−Peŕot effect. More details on the description of the experimental setup and the data analysis procedure are reported elsewhere. 29 Sample Preparation. THz measurements were performed on two samples: (1 a ) a pellet of [VO(acac) 2 ] and high-density polyethylene (HDPE) (1:1 mass ratio) prepared by milling vanadyl-acetylacetonate with HDPE powder (Merck CAS 42.901-5), the powder was pressed by means of a manual hydraulic press (about 0.8 GPa); (1 b ) a pure HDPE pellet on which we deposited a few drops of [VO(acac) 2 ] dichloromethane solution of several microns by the drop-casting method.
Computational Methods. Vibrational Properties. For the density functional theory-based calculations, we used the [VO(acac) 2 ] crystallographic unit-cell coordinates. The temperature effects on the crystal structure were tested exploiting the X-ray data collected at 100 and 300 K. The availability of these data allowed to optimize the molecular structure keeping the experimental cell parameters fixed, 35 in such a way the temperature effects can be indirectly accounted. The simulation of cell parameters at finite temperature can be accomplished only via dynamical DFT methods (ab initio molecular dynamics). Vibrational properties at 0 K were also computed by fully relaxing the electronic degrees of freedom, lattice constants, and atomic positions; the resulting crystal structure is hereafter referred to as 1 0K,cell−opt . A tight DFT optimization is mandatory to avoid imaginary eigenvalues in the Hessian matrix. The different atoms in the unit cell are displaced by small amounts, and the forces on all other atoms are recorded to calculate the second derivative of the potential energy surface (PES) through numerical differentiation. Any is the equilibrium position, generates forces on all other atoms of the determined cell according to the relationship where n and m are the primitive unit-cell indices, and ν and μ are the atomic indices. In the present study, n = 1 because only one primitive cell was included in the computation. This relates the forces generated to the force constant matrix Φ(n, ν, m, μ) and atomic displacement U(m, μ). We imposed a set of constraint equations on the force constant matrix, known as the acoustic sum rule, to ensure that the crystal energy is invariant under the global translation of the whole crystal (ω acoustic ). The dynamical matrix is defined as 16 Here, M μ and M ν are masses of atoms, and k is the wave vector. The solution of the eigenproblem (eq 3) gives eigenvalues ω 2 (k, j) and eigenvectors e(k, j) of 3N orthogonal normal modes of vibrations The calculation is restricted to one primitive cell, giving access only to 3N − 3 zone-center phonons (at the Γ-point). 37 The Γ-point approximation comes with no loss of information. The full inclusion of the Brillouin zone integration, important to quantitatively estimate the optical phonon bands, is not part of the main focus of the present work. 25 A more than fair approximation in the analysis of phonons' decay mechanisms is achieved limiting the calculation to the DFT Γ-point optical modes, along with the acoustic phonons in the whole Brillouin zone. The latter were already computed by some of us at the force field (FF) level, 33 at the DFT level with the finite displacement method at Γ-point, 25 and in the Brillouin zone. 28 Periodic ab initio calculations of energy and forces were performed using the Quickstep module 38 present in the CP2K atomistic simulation package. 39 GTH pseudopotentials 40−42 for the core electrons and valence pseudo-wavefunctions are expanded in Gaussian-type orbitals, and the density is represented in a plane wave auxiliary basis set. In particular, for V, O, C, and H elements, the DZVP-MOLOPT-SR-GTH basis set was used. 38  47,48 was employed to carry out the structure optimization. The convergence of the plane wave basis cut-off was reached for 10 000 Ry with a convergence threshold of 1.0 × 10 −9 au for the SCF energy and 5.0 × 10 −8 au/Å for the forces. A cutoff of 10 000 Ry is cumbersome and the best compromise between accuracy and cpu-time consumption was chosen to be 5000 Ry (see Figures S1 and S2). Such a choice was taken also on the basis of the convergence of vibrational energies in the low-energy range (0−100 cm −1 ). Indeed, the error between the frequency values computed for 10 000 and 5000 Ry is just a few decimal digits of the wavenumber. The high level of accuracy chosen precluded to calculate the phonons' dispersions at the DFT level for all of the three structures investigated.
Magnetic Properties and the Spin−Phonon Coupling. The method to investigate the Sp-Ph coupling has been developed in previous works 15 and is based on a perturbative formalism. Calculation of the Sp-Ph Hamiltonian parameters was carried out with the ORCA package. 49 The level of theory chosen is DFT and the PBE0 hybrid functional 50 was used throughout the computation. A balanced polarized triple-ζ basis set was adopted (def2-TZVP) for V and S atoms, while a polarized single-ζ basis set (def2-SVP) was used for light elements such as C and H. The geometrical equivalence between the two molecules present in the crystal cell after p-DFT optimization with CP2K was ascertained. An RMSD of 0.12 was computed on V, O, and C atoms (0.94 on V, O, C, and H atoms) indicating that two molecules can be considered equivalent even after a symmetry-unconstrained optimization. For such a reason, only one molecule was used for the Sp-Ph calculations. The unperturbed Hamiltonian is evaluated with respect to external perturbation, sampling six points along each coordinate (see Figure S3). The g vs q α plot is then fitted with a polynomial whose first-order coefficient directly yields the The 0 for a S = 1/2 spin system, neglecting the hyperfine term that plays a key role only at the low magnetic field, 11 contains only the Zeeman term , where the main parameter depending on ionic coordinates is the secondorder correction to the Landéfactor g, μ b is the Bohr magneton, and B is the magnetic field.
The Sp-Ph coupling Hamiltonian connects molecular vibrations to the spin via the spin−orbit interaction, 51,52 and, in the weak coupling assumption, it is linear in the ionic The dynamics of the whole system (spin and phonons) can be studied by the Liouville−von Neumann equation for the density operator, within the approximations described by Redfield 53 in the framework of the open quantum system theory. 54 The spin−phonon coupling coefficients constitute the main ingredients to describe the interplay of the principal system (spin degrees of freedom) and the ancillary system (surrounding environment of the phonons), as shown in previous works. 26 ■ RESULTS AND DISCUSSIONS THz Spectra of Powder Samples.
(1 a ) and (1 b ) have been prepared according to two different procedures to check that the applied pressure, essential to make the pellet, does not induce any phase transition or changes in the crystal structure. The THz spectra of (1 a ) and (1 b ) are shown in Figure S4.
The two spectra are identical, except for the intensity that is related to the different concentrations of vanadyl-acetylacetonate in the samples. No evident effects of the applied pressure occur; the same result is observed by PXRD investigation on (1 b ) compared with the vanadyl-acetylacetonate powder and pure HDPE, as shown in Figure S4. Therefore, hereafter, the discussion will be focused on (1 a ). Figure 2 shows the THz spectra measured in the temperature range 10−300 K. At the lowest temperature, five peaks (p1−5) are detected. As the sample temperature is increased, all of them become broader and most of them shift to lower energy. It is worth noting that the lowest phonon detected falls at 43.9 cm −1 (p1), which is exactly the value predicted by the Brons−van Vleck model applied to the spin relaxation rate in our previous work. 21 The extraction of peaks' parameters, i.e., central frequency and linewidth (Γ lw ), has been performed by fitting each of them with a pseudo-Voigt function. Moreover, the baseline correction has been taken into account. Figure 3 summarizes the temperature dependence of the five bands. The error analysis on the five peaks' position led to an average standard deviation of less than 1%. The weak peak at ca. 71 cm −1 (p2) is not detectable above 120 K, as it becomes too broad and the strong peaks at higher frequency get closer and broader. The general trend is a decrease of frequency when the temperature increases and it is ascribed primarily to the effect of the thermal expansion of the lattice, which results in the softening of the intermolecular forces. 16 The peak at ca. 77 cm −1 (p3) is the only one that deviates from the general trend: its temperature dependence is very weak and not monotonous. The nature of the computed mode which should correspond to p3 (vide infra) is characterized by an in-phase rigid vibration of the external acac ligands of the two molecules in the cell with the ones belonging to the molecules in the neighbor cell replica (see Figure S7). Therefore, the variation of cell parameters as a function of temperature is expected to have a very limited impact on this mode. The Journal of Physical Chemistry C pubs.acs.org/JPCC Article Figure 3 shows how the linewidth of each peak, in terms of pseudo-Voigt, varies by increasing the temperature from 10 to 300 K. A deconvolution process was carried out to separate the Lorentzian contribution (Γ L ) and the Gaussian one (Γ G ). The first contribution takes into account the homogeneous broadening, connected to the phonon−phonon and phonon−lattice interactions, through the anharmonic terms of the crystal potential expansion, whereas the second one takes into account the inhomogeneous broadening, which is predominantly attributed to the variation of frequency with the phonons wave vector (k) direction of the phonons active in IR. Indeed, in common spectroscopic techniques, e.g., IR and Raman, the conservation of the moment is fulfilled in a small surrounding of the Γ-point (k ≃ 0), making the phonons' frequency dependent on the k direction.
The results collected in Figure 4 show that, for the temperature range 10−100 K, the homogeneous broadening of peaks p1, p3, and p4 increases almost linearly with temperature, while the inhomogeneous width is practically constant. The relative standard deviation is of the order of 10− 20% since Γ G and Γ L are particularly sensitive to the irregularity of the peak tails. In addition, the loss of a welldefined shape of p2 above 70 K and the high noise that affects p5 do not allow to deduce a clear temperature trend for these peaks.
Calculated Vibrational Spectra. In the Introduction section, we pointed out that the availability of X-ray data collected at cryogenic and room temperature (100 and 300 K) allowed us to fix the experimental cell parameters in the DFT energy minimization, prompting the estimation of the errors introduced in the vibrational analysis when the temperature effects on the structure are neglected, i.e., when cell parameters are free parameters in the energy minimization procedure. We denote 1 0K,cell−opt , 1 100K,opt , and 1 300K,opt as the DFT optimized structures; in 1 0K,cell−opt , both cell parameters and atomic coordinates are optimized while, in 1 100K,opt and 1 300K,opt , only atomic coordinates were left to relax keeping the cell parameters fixed to experimental values. Lowering the temperature from 300 to 0 K, a contraction occurred along the b and c axes, while a nonmonotonous expansion along the a-axis was observed. The overall variation of the cell parameters led to a unit-cell volume contraction. 1 0K,cell−opt was found to be more stable than 1 100K,opt of 2.80 kcal/mol, while 1 100K,opt was more stable than 1 300K,opt of 0.92 kcal/mol, as expected for the most enthalpic favored structure. The full relaxation of 1 0K,cell−opt led to optimized cell parameters that show a maximum deviation with respect to 1 100K of 4.8 and 2.8% for lengths and angles of lattice constants, respectively. The volume underwent a shrinking of ca. 5%. Even more, it is worth highlighting that a, β, and γ parameters computed for 1 0K,cell−opt do not completely follow the expected trend indicated by the experimental parameters obtained for 1 100K and 1 300K unit cells. These results have also a general relevance, as they show that neglecting the temperature effects on cells, i.e., using static DFT cell parameters optimized at 0 K, a poorer agreement in reproducing the vibrational spectra at finite temperatures can be attained. See, for instance, the computed p1 C , which is shifted to almost 20 cm −1 from 0 to 300 K calculated phonons.
The computed IR spectra for the three cells are reported in Figure 5. The frequencies of the lowest 20 Raman-active (9) and IR-active (11) modes are reported in Table 2 with the trends for the first IR modes gathered in Figure 6. The full spectral data are reported in Table S3.
The computed redshifts of peaks in the range 30−130 cm −1 are in an overall good agreement with those observed for the THz spectra measured at 10, 100, and 300 K. A one-to-one correspondence with the experimental peaks is valid for p1−4;  The Journal of Physical Chemistry C pubs.acs.org/JPCC Article the comparison is impossible for p5, as this band results, according to our calculation, as the superposition of four peaks (p5 C −8 C ). The calculated relative intensities are also in good agreement with the experimental ones with the only exception for p2 and p3, for which an inversion of their intensities can be claimed. This interpretation is more plausible than an energy swap of the two modes because p3 C shows a larger shift as a function of the temperature than p2 C as, indeed, experimentally observed. The comparison of the calculated linewidths at different temperatures is impracticable, as in our approach the effects of temperature are taken into account only using different cell parameters. In Figure 7, we report the decomposition procedure outlined by Neto et al. 55 and developed by Lunghi et al. 15 to compute the amplitude of local translation, local rotation, and internal vibrations of a single [VO(acac) 2 ] molecule inside its crystal. The nature of the lowest eight peaks shows that the reticular contributions decrease as the energy of modes increases, as expected. Such a contribution is still detectable up to the 20th mode with a 20−30% of external character. The modes can be sketched as the following: (1) in p1 C , two [VO(acac) 2 ] molecules rigidly tilt one with respect to the other ( Figure S6). (2) p2 C −3 C can be considered quasidegenerate modes where the two acac ligands of a [VO(acac) 2 ] molecule tilt out-of-phase (p2 C ) or in-phase (p3 C ) with respect to one acac belonging to a neighboring [VO(acac) 2 ] molecule: in p2 C , the neighboring molecule is inside the unit cell; in p3 C , in an adjacent one (Figures S6 and  S7). Phonon Decay Mechanisms from Linewidth Analysis. The analysis of the temperature dependence of the homogeneous linewidth allows the rationalization of the decay mechanism of the phonons and the estimation of the third-order anharmonic coefficients.
The temperature-dependent anharmonic coupling of phonons leads to relaxation via multiphonon processes. 56 In the simplest case (predominant at low temperature), line broadening results from two-photon decay, involving downand up-conversion mechanisms. In the first case, the optical phonon decays into two phonons of lower frequencies (downconversion); in the second case, the optical phonon and a phonon of the thermal bath add up to produce a more energetic phonon (up-conversion). 56 Energy and momentum conservation are required for all processes. In the limit kT > hν (where ν is the frequency of the optical phonon), a linear trend of the Lorentzian linewidth indicates that two-phonon decay processes are the dominant relaxation mechanism.  α is the normal mode index, ν α is the frequency, int. is the IR intensity, and Sp-Ph is the magnitude of the spin−phonon coupling coefficient.
The Journal of Physical Chemistry C pubs.acs.org/JPCC Article The decay pathway identified for p1, p3, and p4 corresponds to the down-and up-conversion processes. The simplest downconversion process fulfilling the conservation of energy and momenta requires that the energy of generated phonons is half of the starting one, 2hν = hν 0 . If we assume, for instance, that the optical phonon ν 0 decays into two phonons of the same frequency ν, and the phonons ν 1 and ν 2 sum up to create a ν 0 optical phonon, we can write the following equation 56,57 δ ν ν δ ν ν ν where Γ 0 is the residual linewidth at zero temperature, n = [exp(hν/kT) − 1] −1 is the occupation number of the phonon involved in the decay at temperature T, and B is the third-order anharmonic coefficient; the Kronecker δ functions ensure energy conservation. In principle, the anharmonic lattice dynamics theory provides the recipes for the calculation of phonon linewidth, which requires knowledge of the phonon frequencies in the entire Brillouin zone. This goal goes far beyond the limits of this work; we adopt here a very simplified approach that   The Journal of Physical Chemistry C pubs.acs.org/JPCC Article involves drastic approximations. We verified that the upconversion mechanism gives a small contribution for p1 C , while larger contributions are found for p3 C and p4 C . In addition, we assume that the down-conversion mechanism in which two phonons of the same frequencies are created can be taken as representative of all the down-conversion processes. Within these limits, for p1 C up-conversion, might be allowed considering the interaction with the Raman-active mode 7 (57 cm −1 ) to create a Raman-active phonon at 100 cm −1 (mode 13). p3 C interacts with p1 C and creates a phonon p6 C at 116 cm −1 . The matching is fulfilled considering that the tail of p6 C extends up to 125 cm −1 . Phonon p4 C can combine with Raman mode 5 (48 cm −1 ) to generate a phonon at 146 cm −1 . Concerning the down-conversion processes, the two lowest optical phonons p1 C and p3 C can decay only into acoustic phonons, which, as previously reported, 28 extend from 0 to ca. 40 cm −1 at the Brillouin zone boundary. The relaxation of p4 C involves two phonons of 48 cm −1 , a frequency value at the lowest limit of the optical phonon branches. This schematic description of the relaxation processes is summarized in Figure  8.
The fitting to the experimental linewidth in the temperature range 10−100 K gives for p1 C a third-order anharmonic coefficient B of ca. 2 cm −1 , with a negligible contribution from the up-conversion term. Conversely, the B value for p3 C and p4 C is ca. 2 cm −1 (2.5 and 4 cm −1 , respectively, considering only the down-conversion process).
Spin−Phonon Coupling Analysis. A first qualitative fingerprint of interaction between magnetism and vibrations in each molecule requires the calculation of the phonon spectrum and the amplitude of Sp-Ph coupling for all possible modes,  Figure 9 for all of the three considered cells. Unfortunately, the lack of the inversion center in the molecule did not allow to fully exploit symmetry considerations in the rationalization of the normal modes, as evidenced in previous works. 32 As stated in the previous section, some cell parameters computed for 1 0K,cell−opt do not completely follow the expected trend indicated by the experimental parameters obtained for 1 100K and 1 300K unit cells. The different cell parameters led to an alteration of the compositions of the normal modes and, therefore, to their Sp-Ph couplings. A closer similarity of 1 100K vs 1 300K Sp-Ph coefficients is indeed found, in contrast to 1 100K vs 1 0K,cell−opt and 1 300K vs 1 0K,cell−opt . . For the first three modes, by increasing the temperature, a decrease of the Sp-Ph magnitude was observed. An opposite behavior was observed for the fourth, instead. For several other modes, a decrease or increase in their Sp-Ph coupling was observed though, even if the variation in temperature was not so pronounced. These differences can be ascribed to the modified perturbation of the first coordination sphere caused by the variation of the primitive cell parameters observed for the three considered temperatures (see Table 1). The nonmonotonous (a, β, γ) and monotonous (b, c, α) trends of the cell parameters can differently alter the composition of the normal modes and modify their Sp-Ph coupling efficiency.
Focusing our attention on the IR modes in the experimental THz range (30−130 cm −1 ), the differences between the computed Sp-Ph dispersions of 1 0K,cell−opt and 1 100K /1 300K are more evident.
The computed Sp-Ph coefficients are reported in Figure 9, and only 1 100K,opt and 1 300K,opt will be discussed for homogeneity.
The Sp-Ph coupling can decrease its magnitude or remain unaffected with the increase of temperature; its variation is not constant for each peak going from a slight variation of the Sp-Ph coupling for p1 C (rigid tilting of molecules) to a significant one for p3 C and p5 C (first coordination sphere distortion).
In a recent paper, some of us have shown for the [Dy(acac) 3 (H 2 O) 2 ] 58 complex that a large Sp-Ph coupling is found for those normal modes that can induce large   Figure  S11), which involve the largest charges variation, show the largest Sp-Ph coupling (see Figure 10).
The indications from the computed results are twofold. From a strictly computational point of view, our analysis shows that the use of static approaches for the description of the phonon distribution should be performed when the temperature effects can be indirectly accounted for the use of crystallographic cells collected at different temperatures. We have shown that in such a framework, a quantitative agreement with the THz experiment was obtained, while the ordinary approach of optimizing the geometries along with or without the crystallographic cell parameters can lead only to qualitative results for the vibrational spectrum and for the Sp-Ph coupling. From a structural point of view, it is possible to explicitly show how the Sp-Ph coupling associated with a single normal mode can vary as a function of temperature as a direct consequence of the variation with the temperature of the crystal cell parameters. It is also evident that the Sp-Ph relaxation mechanism can give rise to a more complex behavior than the polynomial or exponential functions usually considered in the literature 59,60 for the T 1 temperature dependence.

■ CONCLUSIONS
The harmonic approximation of molecular vibrations is highly valuable for interpretation of the main features of the vibrational spectra but it cannot explain all their fine structures, which depend on the temperature. At finite temperatures, anharmonic motions of the molecules must be considered since they allow for the possibility of interactions among the harmonic phonons, resulting in slight changes in the energy and composition of characteristic vibrations. The most straightforward parameters to be considered as modified by anharmonicity are the cell parameters. The crystal structures obtained from X-ray diffraction collected at cryogenic and room temperatures allowed us to calculate the observed THz anharmonic frequency shift with high accuracy.
The overall differences in the spin−coupling magnitudes as a function of temperature were highlighted. It turned out that Sp-Ph coupling can decrease or increase its magnitude with the increase of temperature. The reason for such trends has to be ascribed to variations of primitive vectors and unit-cell angles observed for the two considered temperatures. Indeed, the cell volume tunes the normal modes modifying the perturbation to the first coordination sphere of the vanadium ion. The overall results presented in this paper have also a general relevance allowing the quantification of the errors introduced by neglecting the temperature effects on the cell parameters. The cell parameters obtained by full optimization at 0 K do not follow the trends outlined by 100 and 300 K experimental structures, suggesting the possibility that a constrained optimization employing the experimental cells could yield more accurate results.
Moreover, the presented results can have a general relevance regarding the calculation of the T 1 values, too. Indeed, even if significant results have been recently achieved in its computation at the ab initio level, 25,26 the T 1 temperature dependence was only due to the phonon thermal distribution, keeping the Sp-Ph coupling as temperature-independent. However, our calculations suggest that even though the computed Sp-Ph coupling variations are small, the overall calculated variation becomes mandatory for a more rigorous reproduction of T 1 along a wide temperature range.
Therefore, the results of this study represent further pieces of the puzzle of the role played by vibrations in determining the relaxation properties of potential molecular qubits at the quantitative level. In particular, the deeper knowledge of the Sp-Ph coupling might allow one to exploit phonons for the accurate control, e.g., initialization and readout of the qubits. This can be performed by using THz pulses by means of timeresolved ultrafast spectroscopy.   Table S5, while the phonons' composition is depicted in Figure S10.