Introduction

The two main pillars of the current research in molecular magnetism are the possibilities of exploiting molecules as classical bits for high-density magnetic memories1,2,3,4,5 or as quantum bits (qubits) for quantum information processing6,7,8,9,10,11,12. The use of molecular nanomagnets (MNMs) as classical bits relies on their slow relaxation of the magnetisation, which is undermined by spin-phonon interactions13. Molecular vibrations may also play an important role in determining the magnitude and temperature dependence of coherence times in molecular qubits14, which need to be very long to implement quantum algorithms.

Despite their importance, very limited experimental investigations on phonons in MNMs have been performed so far. Indeed, only integrated information on their low-energy spectrum or \(\Gamma\)-point energies have been determined in some compounds through specific heat measurements or Raman and THz spectroscopy techniques13,15,16,17. As also pointed out in a recent critical perspective on the topic18, the key to construct a reliable model of phonon-induced relaxation dynamics in MNMs is to have access to phonon dispersions. In fact, several mechanisms take part into relaxation dynamics of MNMs (direct, Orbach, Raman and quantum tunnelling processes5,19,20,21,22,23), and the interplay between them depends on the phonon spectrum. Since experimental results on phonon dispersions were not available so far, many relaxation theories on MNMs have relied on the simple Debye model24,25,26. Even if they have been successful in understanding nuclear magnetic resonance (NMR) or dynamical susceptibility signals in some systems26,27,28,29, the fundamental limitations of the Debye model prevent a quantitative description of relaxation dynamics in the new generations of MNMs5,14,16,20,30. Moreover, phonon eigenvectors are necessary for the evaluation of spin-phonon coupling coefficients, which requires a full description of molecular vibrations, including atomic displacements due to accessible phonon modes. In particular, the possible presence of low-lying optical phonons producing sizeable intra-molecular distortions can be crucial for magnetic relaxation31. Furthermore, anti-crossings (ACs) between these low-lying optical branches and acoustic phonons, which can be experimentally detected only by measuring phonon dispersions, are also important for magnetic relaxation. While single molecules could be the long-term goal for applications, the majority of experimental studies on the relaxation dynamics of MNMs are performed in crystals or polycrystalline samples. In addition, the use of magnetically diluted single crystals32 can bring advantages in quantum simulation protocols, where the ensemble measurements immediately yield expectation values of the observables. Hence, the combination of an experimental technique directly addressing phonon dispersions and eigenvectors of a MNM crystal with state-of-the-art ab initio calculations is the starting point to investigate the physics behind relaxation mechanisms and benchmark theoretical models.

Here we exploit the four-dimensional inelastic neutron scattering approach (4D-INS)21,33 and the LET spectrometer34 to directly measure for the first time phonon dispersions in a molecular qubit prototype. Indeed, INS is a very powerful technique to study phonons, because it enables a direct access to phonon eigenvalues and eigenvectors. Moreover, the recent advent of spectrometers combining the time-of-flight technique with position-sensitive detectors makes measuring the four-dimensional scattering function \(S({\bf{Q}},\omega )\) in large portions of the reciprocal space possible8,21,35,36. This experimental method is largely applied to study phonon dispersions in inorganic systems with extended structures37,38,39,40,41, but it has never been applied to complex molecular crystals. This work focuses on VO-acetylacetonate ([VO(acac)\({}_{2}\)]), a prototypical complex embedding the vanadyl (VO) unit, archetype of a new generation of molecular qubits with long coherence times up to high temperatures14,16,30. The results of this challenging 4D-INS experiment are compared to density functional theory (DFT) calculations. Simulations reproduce important qualitative features of the 4D-INS data, such as the presence of low-lying optical modes and ACs between acoustic and optical phonons, providing insights on their effect on spin-phonon couplings.

Results

The [VO(acac)\({}_{2}\)] molecular qubit

VO-based molecules contain a penta-coordinated \({\mathrm{V}}^{\mathrm{{IV}}}\) metal centre lying above the basal plane of a distorted square pyramidal coordination geometry. The double bond of the vanadyl VO\({}^{2+}\) ion is much shorter than single V–O bonds, leading to a strong axial distortion of the ligand field acting on the metal centre. (It should be highlighted that, despite the VO2+ unit is commonly described by a double bond, a triple bond has also been proposed42,43 and supported by theoretical calculations44. The resulting splitting of the V\({}^{\mathrm{{IV}}}\)\(d\)-orbitals yields a \({d}_{xy}\) orbital lying the lowest in energy, singly occupied and well separated from the other orbitals. This also quenches the orbital contribution to the ground state, making these \(S=1/2\) systems ideal candidates as molecular spin qubits. The large splitting between the ground state doublet and the first excited levels also ensures that no magnetic excitations are present in the energy range up to more than 1500 meV44. Recent studies on this family of compounds suggest that the rigidity of the VO moiety is also at the basis of their long quantum coherence times up to room temperature14,30,45.

[VO(acac)\({}_{2}\)] is the simplest molecule among all vanadyl complexes investigated so far (see Fig. 1), thus it is an ideal model system to investigate phonons in MNMs from both computational and experimental point of views. [VO(acac)\({}_{2}\)] crystallises in the triclinic centrosymmetric space group, with the unit cell containing two molecules, symmetric with respect to the inversion centre. The absence of magnetic excitations and the characteristics of its low-energy phonon spectrum, as predicted by DFT calculations, make it possible to measure both acoustic and optical modes with neutron incident energies which ensure a sufficiently high resolution. Moreover, it is also possible to grow sufficiently large single crystals of [VO(acac)\({}_{2}\)] of several \({\mathrm{{mm}}}\) in size, required to perform INS experiments, with a high percentage of deuteration (ca. 98%).

Fig. 1: The [VO(acac)\({}_{2}\)] molecular qubit.
figure 1

Molecular structure of VO-acetylacetonate (green, V; red, O; grey, C; white, H). The double bond typical of the VO\({}^{2+}\) moiety is the apical one.

Unveiling phonons with 4D-INS

We exploited the cold-neutron LET spectrometer at ISIS34 to measure the four-dimensional scattering function \(S({\bf{Q}},\omega )\) of [VO(acac)\({}_{2}\)] (more details about \(S({\bf{Q}},\omega )\) in the Methods sections), since the remarkable features of this instrument fit the requirements of this challenging experiment. Indeed, high-energy resolutions are required to resolve the phonon spectrum and, at the same time, it is important to explore wide portions of the reciprocal space. Moreover, we also need to measure at relatively small \(Q\) even at high energies, because in this condition it is easier to focus on the coherent scattering signal. LET offers high-energy resolutions combined with position-sensitive detectors covering a wide solid angle, thus a wide \({\bf{Q}}\) range. In addition, its high neutron flux and low background enable one to obtain high-quality data also from weakly scattering samples. Large crystals and a high percentage of deuteration are indeed required for these INS experiments, but these two requirements are not easily fulfilled in MNMs. The investigated sample was composed of 40 [VO(acac)\({}_{2}\)] co-aligned deuterated single crystals, in order to obtain a final mass of about 1 g (more details in the Methods sections). Measurements were performed at the temperature \(T\) = 5 K to minimise the phonon line broadening, with two different incident neutron energies \({E}_{i}\) = 7.3 and 13 meV (see Methods). The [VO(acac)\({}_{2}\)] 4D scattering function \(S({\bf{Q}},\omega )\) was then reconstructed using the software HORACE46, which allowed us to slice the 4D datasets into 1D curves to visualise excitations intensities as a function of energy around desired \({\bf{Q}}\) values and into 2D slices along specified trajectories in the \(({\bf{Q}},E)\)-space to visualise phonon dispersions along different directions.

An important feature of the INS technique is the possibility to experimentally distinguish between phonon modes with longitudinal or transverse polarisation, characterised respectively by lattice vibrations parallel or perpendicular to the direction of propagation. This is done by properly selecting the relative direction between the neutron scattering vector \({\bf{Q}}\) and the phonon polarisation vector \({{\boldsymbol{\sigma }}}_{j}^{d}({\bf{q}})\), thanks to the dot-product in Eq. (4) (see Methods). For instance, longitudinal modes are extracted by inspecting data obtained with \({\bf{Q}}\) vectors parallel to the desired \({{\boldsymbol{\sigma }}}_{j}^{d}({\bf{q}})\): e.g., longitudinal modes along \(\Gamma\)Z are obtained by selecting \({\bf{Q}}\) vectors of the type (0, 0, \({Q}_{z}\)). Conversely, \({\bf{Q}}\)-trajectories with sizeable components perpendicular to \({{\boldsymbol{\sigma }}}_{j}^{d}({\bf{q}})\) have been chosen to probe transverse phonon modes. Given the very low symmetry of [VO(acac)\({}_{2}\)], we do not expect phonon branches to be purely transverse or longitudinal (e.g., \({{\boldsymbol{\sigma }}}_{j}^{d}({\bf{q}})\) may not contain displacements only perpendicular or parallel to \({\bf{Q}}\))47. Anyway, as far as acoustic branches are concerned, it is often convenient to call them “longitudinal” and “transverse”, although they are not pure.

Phononic excitations intensities of [VO(acac)\({}_{2}\)] are reported as a function of energy in Fig. 2 for some representative \({\bf{Q}}\) values along \(\Gamma\)X, \(\Gamma\)Y and \(\Gamma\)Z symmetry directions, for both \({E}_{i}\) = 7.3 meV (Fig. 2a–c) and 13 meV (Fig. 2d–g) incident energies. The former show one or two excitation peaks, due to longitudinal and transverse acoustic modes, while \({E}_{i}\) = 13 meV data show phononic excitations up to 11 meV from both acoustic and optical modes.

Fig. 2: Phononic excitations intensities vs energy.
figure 2

Black scatters: Representative examples of [VO(acac)\({}_{2}\)] phononic excitations as a function of energy, obtained from LET data with \({E}_{i}\) = 7.3 meV (ac) and \({E}_{i}\) = 13 meV (dg) at \(T\) = 5 K. The 1D cuts have been obtained by integrating the 4D datasets over the three components of the neutron scattering vector \({\bf{Q}}=[\eta ,\zeta ,\xi ]\) (expressed in terms of reciprocal lattice vector units) around the desired values, with error bars representing the SE. For data obtained with \({E}_{i}\) = 7.3 meV the intensity is normalised for the maximum in each panel, while for data obtained with \({E}_{i}\) = 13 meV the intensity is normalised in order to saturate acoustic modes and give prominence to optical ones. Red lines: intensity curves as a function of energy calculated with Eq. (3) (see Methods) and with phonon energies and polarisation vectors obtained from DFT.

1D cuts give a detailed insight onto excitations, whereas 2D slices of the 4D datasets provide an immediate visualisation of phonon dispersions, also demonstrating the capabilities of the 4D-INS technique. For instance, the data obtained with \({E}_{i}\) = 7.3 meV in Fig. 3 provide a picture of [VO(acac)\({}_{2}\)] acoustic phonons dispersions up to 6 meV. Representative examples of the measured dispersion curves along \(\Gamma\)Z, \(\Gamma\)X and \(\Gamma\)Y symmetry directions are reported in Fig. 3. Along \(\Gamma\)Z, the longitudinal acoustic mode is the only intense one for data along the \([0,0,\xi ]\) longitudinal direction (Fig. 3a). Indeed, it corresponds to the steepest among the three modes along \(\Gamma\)–Z (see below). High-intensity modes in Fig. 3c and d are instead due to the two transverse acoustic branches along \(\Gamma\)X and \(\Gamma\)Y, respectively, which reach the Brillouin zone boundary at only about 4 meV. In addition, low-lying optical modes are visible around 5 meV along the transverse \(\Gamma\)Z direction reported in Fig. 3b. In order to visualise [VO(acac)\({}_{2}\)] low-energy optical modes, which could play an important role in magnetic relaxation, we also report phonon dispersions obtained with \({E}_{i}\) = 13 meV in Fig. 4. Figure 4a confirms the longitudinal acoustic mode along \(\Gamma\)Z shown in Fig. 3a and, in addition, is now possible to observe optical modes around 8 and 10 meV. Furthermore, optical modes up to 11 meV are clearly visible also in transverse configuration (\({\bf{Q}}\nparallel {\bf{q}}\)) along \(\Gamma\)Z, \(\Gamma\)X and \(\Gamma\)-Y directions (Fig. 4b–d, respectively).

Fig. 3: Phonon dispersions up to 6 meV.
figure 3

In ad we show representative examples of [VO(acac)\({}_{2}\)] phonon dispersions measured on LET with \({E}_{i}\) = 7.3 meV and \(T\) = 5 K. This incident energy allowed us to resolve acoustic modes up to 6 meV. Dispersions curves have been obtained integrating over two out of the three components of the neutron scattering vector \({\bf{Q}}=[\eta ,\zeta ,\xi ]\) (expressed in terms of reciprocal lattice vector units) around different Bragg peaks. a Shows longitudinal acoustic modes along the symmetry direction \(\Gamma\)Z around [0, 0, −2]. bd Report transverse acoustic modes along three principal symmetry directions: \(\Gamma\)Z around [−1, 0, −3], \(\Gamma\)X around [0, 0, −2] and \(\Gamma\)Y around [0, 0, −2] and [0, 1, −2]. The colour bar reports the intensity normalised for the maximum in each panel. Simulated dispersion curves along the same directions are reported in eh and were calculated with Eq. (3) (see Methods), with phonon energies and polarisation vectors obtained from DFT.

Fig. 4: Phonon dispersions up to 11 meV.
figure 4

In ad we show representative examples of [VO(acac)\({}_{2}\)] phonon dispersions measured on LET with \({E}_{i}\) = 13 meV and \(T\) = 5 K. Dispersions curves have been obtained integrating over two out of the three components of the neutron scattering vector \({\bf{Q}}=[\eta ,\zeta ,\xi ]\) (expressed in terms of reciprocal lattice vector units) around different Bragg peaks. a Shows longitudinal acoustic modes along the symmetry direction \(\Gamma\)Z around [0, 0, −2], as in Fig. 3a. With this incident energy we have been able to resolve not only the acoustic but also optical phonon modes up to 11 meV. bd Display transverse modes along three principal symmetry directions: \(\Gamma\)Z around [0, 2, −5] and [0, 2, −4], \(\Gamma\)X around [−1, 0, −3] and \(\Gamma\)Y around [0, 2, −5]. The colour bar reports the intensity normalised for the maximum in each panel. Simulated dispersion curves along the same directions are reported in eh and were calculated with Eq. (3) (see Methods), with phonon energies and polarisation vectors obtained from DFT.

Phonon ACs involve branches belonging to the same symmetry representation, which therefore cannot cross each other. In particular, INS data show that low-lying optical branches display ACs with acoustic ones in [VO(acac)\({}_{2}\)] (see Fig. 5). A strong mixing between the modes is expected close to ACs, which should significantly affect atomic displacements and thus crystal-field modulations involved in magnetic relaxation (see Discussion). Moreover, ACs cause a reduction of phonon lifetimes48,49,50, again possibly affecting the relaxation dynamics23. We have therefore performed an analysis of acoustic-optical phonon ACs, considering transverse \(\Gamma\)Z and \(\Gamma\)X directions data. Results of this analysis are reported in Fig. 5. The position and intensity of each phononic excitation as a function of energy has been extracted from LET data over a significant portion of the Brillouin zone. Figure 5 highlights the typical features of phonon ACs: acoustic modes bend and lose linearity well before reaching the zone boundary and they exchange intensity with the optical branch, which tends to drift apart after the AC. ACs are also evident along the \(\Gamma\)Y direction from INS data in Fig. 4d, but are more complex to unravel. Indeed, more than one low-energy AC occur around the same \(q\) along \(\Gamma\)Y, as confirmed by DFT calculations (see the following section).

Fig. 5: Phonon anti-crossings.
figure 5

Panels a and c report 1D cuts as a function of energy of the LET data, obtained integrating around different values of Q for a \(\Gamma\)Z around [−1, 0, −3] with \({E}_{i}\) = 7.3 meV (Fig. 3b) and c \(\Gamma\)X around [−1, 0, −3] with \({E}_{i}\)  = 13 meV (Fig. 4c). Panels b and d report the position of the peaks obtained from the 1D cuts of the data for both acoustic (black scatters) and optical modes (red and blue scatters). The intensity of each peak is translated into scatter dimension and the dashed green lines interpolate the linear behaviour of the acoustic modes. The data show the typical features of phonon ACs: bending of the acoustic mode, intensity exchange between acoustic and optical branches, with the latter drifting apart after the AC.

DFT simulations of 4D-INS data

Preliminary force field calculations derived from DFT results had been used in ref. 30 for a first theoretical evaluation of the phonon dispersion curves of [VO(acac)\({}_{2}\)]. More accurate results have now been obtained by employing DFT with the finite-displacements method in a supercell. It is worth noting that calculations of phonons for a generic Brillouin zone point by DFT is a daunting task, as the number of atoms that needs to be included is usually very large and approaches the computational limits of this method. Moreover, the accurate description of molecular crystals is particularly challenging due to the known shortcomings of DFT in describing van der Waals (vdW) interactions51,52,53. The latter are expected to be particularly weak in [VO(acac)\({}_{2}\)], which is highly volatile as are most \(\beta\)-diketonate neutral complexes.

A periodic 3 × 3 × 3 supercell of [VO(acac)\({}_{2}\)] comprising 1620 atoms was used to determine all the translational symmetry-inequivalent reticular force constants. This calculation gives access to vibrational properties for the entire Brillouin zone and is here used to evaluate phonon frequencies and eigenvectors along the three directions of interest. The results reported in Fig. 6 show acoustic modes up to 5 meV and, most importantly, optical branches lying at very low energies and displaying ACs with the acoustic ones, in agreement with the experimental results. Indeed, Fig. 6 also shows the energies of some representative phononic excitations extracted from 1D cuts as in Fig. 2, superimposed to DFT calculations. This comparison shows that calculated phonon energies reproduce important qualitative features of the experimental findings, such as the presence of low-lying optical modes and ACs with acoustic modes. However, a rescaling of 13\(\%\) has been uniformly applied to all phonon modes, in order to lower the calculated phonon energies \({\omega }_{j}({\bf{q}})\) and better reproduce INS data. The observed discrepancies are consistent with other results in literature for molecular crystals54,55,56 and are mainly due to the limitations of DFT methods in describing vdW interactions in these systems51,52,53, especially at low energies. However, the comparison of phonon frequencies alone is not enough to guarantee the accuracy of DFT simulations. Indeed, it is crucial to verify the normal modes composition, since it enters in the definition of spin-phonon coupling coefficients and deeply affects the calculated spin dynamics. Our experiment probes also this composition, because the INS cross-section is determined by phonon eigenstates.

Fig. 6: DFT calculations of phonon dispersions compared to experimental excitation energies.
figure 6

DFT calculations of the low-energy phonon dispersion curves of [VO(acac)\({}_{2}\)] up to 15 meV along the path X(0.5, 0, 0)–\(\Gamma\)(0, 0, 0)–Y(0, 0.5, 0)–\(\Gamma\)(0, 0, 0)–Z(0, 0, 0.5) (black lines), where the coordinates are expressed in the reciprocal lattice basis units. Red scatters are the energies of experimental excitations extracted by fitting 1D cuts of the data as in Fig. 2 (error bars as SD are of the order of scatters dimension). Not all the modes are visible in the considered configurations. Blue circles highlight ACs between acoustic and low-lying optical branches along \(\Gamma\)X and \(\Gamma\)Z. Red and green dashed lines indicate the q vectors at which we have investigated the atomic displacements far from and at the ACs, respectively.

To investigate this aspect, we have simulated the coherent scattering function \(S({\bf{Q}},\omega )\) of Eq. (3) (see Methods), starting from phonon energies \({\omega }_{j}({\bf{q}})\) and normal modes polarisation vectors \({{\boldsymbol{\sigma }}}_{j}^{d}({\bf{q}})\) calculated by DFT. The simulations of 4D-INS spectra are reported in Fig. 2 (red lines) for intensity vs energy 1D plots and in Figs. 3e–h and 4e–h for 2D colour-maps. Our calculations show an overall qualitative agreement with the observed excitations and dispersions. In particular, it is evident from Figs. 3 and 4 that, while phonon dispersions are symmetric with respect to the Brillouin zone centre, neutron intensity patterns are not, as experimentally observed. However, quantitative discrepancies are evident especially from Fig. 2.

Given that DFT results are able to reproduce important qualitative features of experiments, we can further analyse them in order to extract the atomic motions associated with specific phonons. Of particular interest is the nature of acoustic and low-lying optical phonons near the AC points. For instance, Fig. 7 shows atomic displacements associated to a transverse acoustic mode before and at the AC with a low-lying optical branch along the \(\Gamma\)X direction. It is evident that, while before the AC the displacements have a translational character (Fig. 7a), at the AC molecular vibrations are profoundly changed and the rigid translation of the molecule typical of acoustic modes is now rather small. Indeed, Fig. 7b shows that the predominant motions induced by the acoustic phonon at the AC are a rigid translation strongly mixed with bending and torsions of acetylacetonate ligands and with methyl groups rotations, causing the VO\({}^{2+}\) moiety to move with respect to the basal plane. Thus, both distances and angles between the V centre and single-bonded oxygen ions are modified by this mixed acoustic mode, while, as expected, the double bond of the VO\({}^{2+}\) moiety remains rigid. We have also verified that these distortions are the same induced by the optical phonon, thus demonstrating a strong admixing of the two modes at the AC. Similar conclusions can be drawn for the AC involving the transverse acoustic mode and the low-lying optical branch along the \(\Gamma\)Z direction.

Fig. 7: DFT atomic displacements.
figure 7

a The DFT equilibrium molecular structure (blue ball-and-stick) is superimposed with the distorted one according to the nature of the second transverse acoustic mode along \(\Gamma\)X at q = (0.145, 0.0, 0.0) (red ball-and-stick, red dashed line in Fig. 6). b The DFT equilibrium molecular structure (blue ball-and-stick) is overlapped with the distorted one according to the nature of the second acoustic mode along \(\Gamma\)X at q = (0.345, 0.0, 0.0) (green ball-and-stick, green dashed line in Fig. 6). Arrows sketch the small rigid translation of the molecule which is now strongly mixed with the bending of the acetylacetonate ligands and their torsion around the dashed axis.

Discussion

Having direct access to phonon dispersions, and in particular for the first time to ACs between acoustic and optical modes, it is now interesting to focus on their effect on the relaxation dynamics. Indeed, the here-calculated phonons energies and eigenvectors enable the evaluation of the relaxation dynamics in this prototype molecular qubit. A first ab initio simulation of direct relaxation processes has just been performed57. In that contribution the authors showed that in high external magnetic field the dominant mechanism of spin relaxation comes from the modulation of the \({\bf{g}}\) tensor of the Spin Hamiltonian Zeeman term by low-energy acoustic modes. Starting from those findings, we here investigate the effect of ACs on this mechanism by calculating for the first time the q-resolved spin-phonon coupling coefficients. In particular, we report the squared norm of each coefficient for each phonon mode \(j\) as a function of the phonon wave vector \({\bf{q}}\):

$${V}_{\,\text{SPh}\,}^{2}({\omega }_{j},{\bf{q}})=\sum _{\alpha ,\beta = x,y,x}{\left(\frac{\partial {g}_{\alpha ,\beta }}{\partial {Q}_{j,{\bf{q}}}}\right)}^{2},$$
(1)

where \({g}_{\alpha ,\beta }\) are the components of the g tensor and \({Q}_{j,{\bf{q}}}\) are those of the normal modes57. Results for the AC along \(\Gamma\)X are reported in Fig. 8 (see also Supplementary Fig. 2 for results along \(\Gamma\)Z) and demonstrate that, in general, optical modes are characterised by stronger spin-phonon couplings coefficients with respect to the acoustic ones. However, at the ACs the stronger spin-phonon coupling of the almost flat optical modes is partially transferred to the dispersive acoustic modes, which display a maximum of \({V}_{\,{\text{SPh}}\,}^{2}\). This is evident for both the AC involving the transverse acoustic mode analysed in Figs. 5 and 7 and for the AC involving the longitudinal one at higher energies. Figure 8 shows that the mixing of the modes is not negligible even well before the AC at lower \(| {\bf{q}}|\) values and lower energies, where the acoustic mode is still strongly dispersive. Indeed, the enhancement of spin-phonon coupling for the nominally acoustic mode starts at energies significantly smaller than the AC one. It is well known that low-energy dispersive modes are crucial for magnetic relaxation. Hence, the small optical component present in the nominally acoustic mode at low energies can significantly influence magnetic relaxation, by modulating the g tensor with larger probability and in a wider energy range than without the mixing of modes.

Fig. 8: DFT spin-phonon couplings.
figure 8

Colours map spin-phonon couplings squared norm \({V}_{\,\text{SPh}\,}^{2}({\omega }_{j},{\bf{q}})\) calculated as in Eq. (1) for the low-energy phonon modes along \(\Gamma\)X, as a function of \(| {\bf{q}}|\). At the ACs the strong coupling of the optical modes is transferred to the acoustic ones.

These results about the effect of ACs on the spin-phonon couplings further confirm the importance of a complete description of phonons to understand relaxation dynamics in MNMs18. Recently, ab initio calculations have been used to predict the contribution of individual vibrational modes to the spin-phonon coupling23,31,58. The agreement with experimental results was limited because only vibrations at the \(\Gamma\)-point were included31 or due to the gas-phase approximation5. For comparison, gas-phase calculations of the magneto-vibrational couplings for the present molecule (see Methods) have been performed and are reported in Supplementary Note 1. At last, it was shown in thermoelectric materials with low-energy optical branches48,49,50 that the strong mixing of the two modes involved in the ACs decreases acoustic phonons group velocity and also phonon lifetimes. Indeed, the presence of low-lying optical modes is also expected to enhance phonon–phonon scattering processes. This decrease of phonon lifetimes could further affect magnetic relaxation23.

The scenario evidenced by this study is expected to be common to other MNMs. Indeed, the presence of soft coordination bonds and vdW interactions strongly reduces the energy of rotational and intra-molecular vibrations, leading to a non-negligible admixing between acoustic and optical branches. This effect is already visible in crystals of small-size molecules such as arenes54,55. Furthermore, being able to correlate the molecular structure of these systems with their experimentally tested vibrational and phonon spectra, will also enable the development of new strategies for the design of new and optimised molecular qubits and bits. Our results point out that long coherence times or slow magnetic relaxation are undermined by the presence of acoustic-optical phonons ACs. Thus, future synthetic strategies should, for instance, reduce the presence of soft coordination bonds with the aim to shift optical branches to higher energies and thus remove some of the ACs. Few ACs could be also eliminated by increasing the crystal symmetry. However, this would only affect crossings along symmetry directions, while ACs would still be present along generic \({\bf{q}}\) directions.

At last, the present study represents an important test for vdW-corrected DFT in describing phonon dispersions. We showed that while qualitative features of experimental results are reproduced, a full quantitative reproduction of both eigenvalues and eigenvectors will require further improvements of DFT to treat vdW interactions. The DFT community involved in this field commonly use molecular crystals’ cohesive energy59,60 as benchmark test, whose direct comparison with experimental estimation is known to represent a challenge, due to the presence of multiple effects contributing to the measured values (e.g., finite-temperature and zero-point-energy effects61). Conversely, here we have a direct comparison between the practically zero-temperature experimental data and the calculated phonon energies, dispersions and polarisation vectors. This comparison for molecular crystals of such a complexity is almost unique and it could represent also a new standard in the benchmarking of vdW correction schemes.

Methods

Synthesis and crystal structure of the samples

[VO(acac)\({}_{2}\)] is a prototypical molecular qubit containing a single magnetic ion: We use the term molecular nanomagnets (MNMs) also for this new generation of molecular bits and qubits.

[VO(acac-d7)2] was obtained by reaction of a basified (NaOD 99\(\%\) D) deuterated water solution (D\({}_{2}\)O 99.8\(\%\) D) of acetylacetone-d\({}_{8}\) (98\(\%\) D) and VOSO\({}_{4}\). The precipitate of [VO(acac-d7)2] was filtered and recrystallised in inert atmosphere with deuterated acetone-d\({}_{6}\) (99.95\(\%\) D). Single crystals were obtained by slow evaporation of the solvent (2 weeks). The effective percentage of deuteration of [VO(acac-d7)2] was estimated by electron spray ionisation mass spectrometry, and it was found to be 98\(\pm 1 \%\). All crystals were indexed using a single-crystal X-ray diffractometer. The investigated sample included about 40 single crystals for a total mass of about 1 g, co-aligned one-by-one using the ALF crystal alignment facility at ISIS. Since [VO(acac)\({}_{2}\)] crystallises in the triclinic centrosymmetric space group, this procedure represented an additional challenge for the experiment, but the final mosaicity of the sample was enough to obtain sufficiently focused Bragg peaks (see also Supplementary Fig. 1).

4D-INS experiment

The 4D-INS experiments on [VO(acac)\({}_{2}\)] were performed on the LET cold neutrons time-of-flight spectrometer at the ISIS Neutron and Muon source34. LET is a direct geometry spectrometer with \(\pi\)-steradians position-sensitive detectors bank divided in about 10\({}^{5}\) pixels, covering 180° of azimuthal angle and ±30° out-of-plane.

Forty [VO(acac)\({}_{2}\)] single crystals were placed on five different aluminium plates using a fluorinated oil, then stacked on a custom-made sample holder. Having aligned the crystals with the reciprocal vector \({{\bf{a}}}^{* }\) perpendicular to the scattering plane, we were able to explore large portions of the [VO(acac)\({}_{2}\)] reciprocal space in the \(\Gamma\)Z and \(\Gamma\)Y directions and a smaller portion in the \(\Gamma\)X one (see Supplementary Fig. 1). Measurements were performed by rotating the crystals (in steps of 1°) about the vertical axis.

Two incident neutron energies of 7.3 and 13 meV were selected in repetition rate multiplication mode, with an energy resolution (full-width at half-maximum) of 0.23 and 0.55 meV at the elastic line, respectively. All measurements were performed at \(T\) = 5 K.

DFT calculations

All the structural optimisation and Hessian calculations were performed with the CP2K software62 at the level of DFT with the Perdew–Burke–Ernzerhof functional63 including Grimme’s D3 vdW corrections64,65. Indeed, the requirement of vdW-corrected functionals for the modelling of molecular crystals is by now confirmed and settled by several results in literature51,52,53,55. A double-zeta polarised (DZVP) MOLOPT basis set and a 600 Ry of plane-wave cutoff were used for all the atomic species. All the translational symmetry independent force constants were computed by finite difference approach with a 0.01 Å atomic displacements.

The optimisation of the isolated molecule has been done with a plane-wave cutoff of 850 Ry and by turning off the periodic boundary conditions. An integration step of 0.01 Å was used for the calculation of the isolated molecules force constants.

The coefficients \(\left(\partial {g}_{\alpha \beta }/\partial {Q}_{j{\bf{q}}}\right)\) were calculated as

$$\left(\frac{\partial {g}_{\alpha \beta }}{\partial {Q}_{j{\bf{q}}}}\right)=\sum _{l}^{{N}_{\mathrm{{cells}}}}\sum _{is}^{N,3}\sqrt{\frac{\hslash }{{N}_{q}{\omega }_{j{\bf{q}}}{m}_{i}}}{\mathrm{e}}^{i{\bf{q}}\cdot {{\bf{R}}}_{l}}{L}_{is}^{j{\bf{q}}}\left(\frac{\partial {H}_{{\rm{s}}}}{\partial {X}_{is}^{l}}\right)\ ,$$
(2)

where \({X}_{is}^{l}\) is the \(s\) Cartesian coordinate of the ith atom of \(N\) with mass \({m}_{i}\), inside the unit-cell replica at position \({{\bf{R}}}_{l}\) and \({N}_{q}\) is the number of \(q\)-points used. Finally, \({\omega }_{j{\bf{q}}}\) and \({L}_{is}^{j{\bf{q}}}\) are the jth eigenvalue and eigenvector of the Hessian matrix at the q point. The Cartesian derivatives \(\left(\partial {g}_{\alpha \beta }/\partial {X}_{is}\right)\) were obtained from ref. 57.

Data analysis and simulations

The 4D-INS technique is named after the possibility to measure the four-dimensional scattering function \(S({\bf{Q}},\omega )\), i.e. as a function of the transferred energy and of the three components of the transferred momentum \({\bf{Q}}\). The one-phonon coherent scattering function \(S({\bf{Q}},\omega )\) is defined as66

$$S({\bf{Q}},\omega ){\,}\propto {\,} \sum _{j,{\bf{G}},{\bf{q}}}{\left|{F}_{j,{\bf{q}}}({\bf{Q}})\right|}^{2} \frac{1}{2{\omega }_{j}({\bf{q}})}\left[{n}_{j}({\bf{q}})\delta (\omega +{\omega }_{j}({\bf{q}}))\delta ({\bf{Q}}+{\bf{q}}-{\bf{G}})\right.\\ +\left.\ ({n}_{j}({\bf{q}})+1)\delta (\omega -{\omega }_{j}({\bf{q}}))\delta ({\bf{Q}}-{\bf{q}}-{\bf{G}})\right],$$
(3)

where \({\omega }_{j}({\bf{q}})\) are the phonons’ frequencies. The summation in Eq. (3) is extended to all phononic branches \(j\). The momentum conservation law of the scattering events involves the neutron scattering vector \({\bf{Q}}\), the reciprocal lattice vector \({\bf{G}}\) and the phonon wave vector \({\bf{q}}\), while the energy conservation depends on the neutron energy transfer \(\omega\). The term \({n}_{j}({\bf{q}})={(\exp (\beta {\omega }_{j}({\bf{q}}))-1)}^{-1}\) is the phonon Bose factor with \(\beta ={({k}_{B}T)}^{-1}\).

The structure factor \({F}_{j,{\bf{q}}}({\bf{Q}})\) takes into account interference effects between the different atoms \(d\) in the unit cell:

$${F}_{j,{\bf{q}}}({\bf{Q}})=\sum _{d}\frac{{b}_{d}}{\sqrt{{m}_{d}}}{e}^{-{W}_{d}({\bf{Q}})}{e}^{i{\bf{Q}}\cdot {\bf{R}}_{d}}({\bf{Q}}\cdot {{\boldsymbol{\sigma }}}_{j}^{d}({\bf{q}})),$$
(4)

where \({e}^{-{W}_{d}({\bf{Q}})}\) is the Debye-Waller factor, \({b}_{d}\) is the coherent scattering length, \({m}_{d}\) the mass and \({\bf{R}}_{d}\) the position vector in the real space of each atom \(d\) and \({{\boldsymbol{\sigma }}}_{j}^{d}({\bf{q}})\) are the phonons’ polarisation vectors.

Time-of-flight data reduction were performed with Mantid analysis suite67 and measurements for different rotation angles were combined using HORACE46, yielding four-dimensional \(S({\bf{Q}},\omega )\) datasets for the selected incident energies. HORACE also allows to slice these 4D datasets into 2D slices along specified trajectories in the \(({\bf{Q}},E)\)-space, in order to visualise phonon dispersion along different directions, and into 1D curves, to visualise excitations intensities as a function of energy around desired \({\bf{Q}}\) values. In order to reduce the noise-to-signal ratio, data along each symmetry direction were integrated over the other two within a \(\pm {\!}\)0.1 range (in reciprocal lattice vector units) centred around the Bragg peak (for longitudinal \(\Gamma\)Z and transverse \(\Gamma\)X and \(\Gamma\)Y modes with \({E}_{i}\) = 7.3 meV it is necessary to integrate within a \(\pm {\!}\)0.2 range). Empty sample holder data were subtracted, corrections for multiple-scattering flat backgrounds and a 7 bins smoothing were also applied. Low-\(Q\) regions of the reciprocal space were selected for our in-depth analysis, in order to better focus on coherent scattering. 1D cuts as a function of energy were obtained by integrating 2D slices over the surviving direction in a \(\pm {\!}\)0.1 interval centred around the desired \({\bf{Q}}\) value and by applying a seven bins smoothing, with error bars representing the SE. They were then fitted with Gaussian functions in order to extract excitation energies in Fig. 6.

Data simulations were performed by calculating the scattering cross-section in Eq. (3) with DFT calculated phonon energies \({\omega }_{j}({\bf{q}})\) and normal modes polarisation vectors \({{\boldsymbol{\sigma }}}_{j}^{d}({\bf{q}})\). To reproduce INS data, a small rescaling of about 13\(\%\) was uniformly applied to all the calculated phonon energies along all the symmetry directions. We then assumed a Gaussian line-shape with \(\sigma\) = 0.6 meV. For the simulation of 1D cuts in Fig. 2 we also added the contribution of the elastic signal.