Unveiling phonons in a molecular qubit with four-dimensional inelastic neutron scattering and density functional theory

Phonons are the main source of relaxation in molecular nanomagnets, and different mechanisms have been proposed in order to explain the wealth of experimental findings. However, very limited experimental investigations on phonons in these systems have been performed so far, yielding no information about their dispersions. Here we exploit state-of-the-art single-crystal inelastic neutron scattering to directly measure for the first time phonon dispersions in a prototypical molecular qubit. Both acoustic and optical branches are detected in crystals of [VO(acac)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}_{2}$$\end{document}2] along different directions in the reciprocal space. Using energies and polarisation vectors calculated with state-of-the-art Density Functional Theory, we reproduce important qualitative features of [VO(acac)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}_{2}$$\end{document}2] phonon modes, such as the presence of low-lying optical branches. Moreover, we evidence phonon anti-crossings involving acoustic and optical branches, yielding significant transfers of the spin-phonon coupling strength between the different modes.

T he two main pillars of the current research in molecular magnetism are the possibilities of exploiting molecules as classical bits for high-density magnetic memories [1][2][3][4][5] or as quantum bits (qubits) for quantum information processing [6][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 interactions 13 . Molecular vibrations may also play an important role in determining the magnitude and temperature dependence of coherence times in molecular qubits 14 , 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 Γ-point energies have been determined in some compounds through specific heat measurements or Raman and THz spectroscopy techniques 13,[15][16][17] . As also pointed out in a recent critical perspective on the topic 18 , 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 processes 5,[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 model [24][25][26] . Even if they have been successful in understanding nuclear magnetic resonance (NMR) or dynamical susceptibility signals in some systems [26][27][28][29] , the fundamental limitations of the Debye model prevent a quantitative description of relaxation dynamics in the new generations of MNMs 5,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 relaxation 31 . 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 crystals 32 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 spectrometer 34 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 fourdimensional scattering function SðQ; ωÞ in large portions of the reciprocal space possible 8,21,35,36 . This experimental method is largely applied to study phonon dispersions in inorganic systems with extended structures [37][38][39][40][41] , but it has never been applied to complex molecular crystals. This work focuses on VOacetylacetonate ([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 temperatures 14,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 V 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 VO 2+ unit is commonly described by a double bond, a triple bond has also been proposed 42,43 and supported by theoretical calculations 44 . The resulting splitting of the V 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 meV 44 . 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 temperature 14,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 mm in size, required to perform INS experiments, with a high percentage of deuteration (ca. 98%).
Unveiling phonons with 4D-INS. We exploited the cold-neutron LET spectrometer at ISIS 34 to measure the four-dimensional scattering function SðQ; ωÞ of [VO(acac) 2 ] (more details about SðQ; ωÞ 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 positionsensitive detectors covering a wide solid angle, thus a wide 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ðQ; ωÞ was then reconstructed using the software HORACE 46 , which allowed us to slice the 4D datasets into 1D curves to visualise excitations intensities as a function of energy around desired Q values and into 2D slices along specified trajectories in the ð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 Q and the phonon polarisation vector σ d j ðqÞ, thanks to the dot-product in Eq. (4) (see Methods). For instance, longitudinal modes are extracted by inspecting data obtained with Q vectors parallel to the desired σ d j ðqÞ: e.g., longitudinal modes along Γ-Z are obtained by selecting Q vectors of the type (0, 0, Q z ). Conversely, Q-trajectories with sizeable components perpendicular to σ d j ð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., σ d j ðqÞ may not contain displacements only perpendicular or parallel to 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 Q values along Γ-X, Γ-Y and Γ-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.
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 Γ-Z, Γ-X and Γ-Y symmetry directions are reported in Fig. 3. Along Γ-Z, the longitudinal acoustic mode is the only intense one for data along the ½0; 0; ξ longitudinal direction (Fig. 3a). Indeed, it corresponds to the steepest among the three modes along Γ-Z (see below). Highintensity modes in Fig. 3c and d are instead due to the two transverse acoustic branches along Γ-X and Γ-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 Γ-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 Γ-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 (Qq) along Γ-Z, Γ-X and Γ-Y directions ( Fig. 4b-d, respectively).
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 lifetimes [48][49][50] , again possibly affecting the relaxation dynamics 23 . We have therefore performed an analysis of acoustic-optical phonon ACs, considering transverse Γ-Z and Γ-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 Γ-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 Γ-Y, as confirmed by DFT calculations (see the following section).
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) interactions [51][52][53] . The latter are expected to be particularly weak in [VO(acac) 2 ], which is highly volatile as are most β-diketonate neutral complexes.
A periodic 3 × 3 × 3 supercell of [VO(acac) 2 ] comprising 1620 atoms was used to determine all the translational symmetryinequivalent 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 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 lowlying 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 Γ-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 Γ-Z direction.

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 performed 57 . In that contribution the authors showed that in high external magnetic field the dominant mechanism of spin relaxation comes from the modulation of the 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 q: where g α;β are the components of the g tensor and Q j;q are those of the normal modes 57 . Results for the AC along Γ-X are reported in Fig. 8 (see also Supplementary Fig. 2 for results along Γ-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 2 SPh . 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 jqj 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.
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 MNMs 18 . Recently, ab initio calculations have been used to predict the contribution of individual vibrational modes to the spin-phonon coupling 23,31,58 . The agreement with experimental results was limited because only vibrations at the Γ-point were included 31 or due to the gas-phase approximation 5 . 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 branches 48-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 relaxation 23 .
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 arenes 54,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 q directions. At last, the present study represents an important test for vdWcorrected 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 energy 59,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., finitetemperature and zero-point-energy effects 61 ). Conversely, here we have a direct comparison between the practically zero-  (Fig. 3b) and c Γ-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.

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-d 7 ) 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-d 7 ) 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-d 7 ) 2 ] was estimated by electron spray ionisation mass spectrometry, and it was found to be 98 ± 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).  34 . LET is a direct geometry spectrometer with π-steradians positionsensitive 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 a Ã perpendicular to the scattering plane, we were able to explore large portions of the [VO(acac) 2 ] reciprocal space in the Γ-Z and Γ-Y directions and a smaller portion in the Γ-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 software 62 at the level of DFT with the Perdew-Burke-Ernzerhof functional 63 including Grimme's D3 vdW corrections 64,65 . Indeed, the requirement of vdW-corrected functionals for the modelling of molecular crystals is by now confirmed and settled by several results in literature [51][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 ∂g αβ =∂Q jq were calculated as where X l is is the s Cartesian coordinate of the ith atom of N with mass m i , inside the unit-cell replica at position R l and N q is the number of q-points used. Finally, ω jq and L jq is are the jth eigenvalue and eigenvector of the Hessian matrix at the q point. The Cartesian derivatives ∂g αβ =∂X is 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ðQ; ωÞ, i.e. as a function of the transferred energy and of the three components of the transferred momentum Q. The one-phonon coherent scattering function SðQ; ωÞ is defined as 66 SðQ; ωÞ / X j;G;q F j;q ðQÞ 2 1 2ω j ðqÞ n j ðqÞδðω þ ω j ðqÞÞδðQ þ q À GÞ h þ ðn j ðqÞ þ 1Þδðω À ω j ðqÞÞδðQ À q À GÞ i ; where ω j ð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 Q, the reciprocal lattice vector G and the a b Fig. 7 DFT atomic displacements. 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 Γ-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 Γ-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. phonon wave vector q, while the energy conservation depends on the neutron energy transfer ω. The term n j ðqÞ ¼ ðexpðβω j ðqÞÞ À 1Þ À1 is the phonon Bose factor with β ¼ ðk B TÞ À1 .
The structure factor F j;q ðQÞ takes into account interference effects between the different atoms d in the unit cell: where e ÀW d ðQÞ is the Debye-Waller factor, b d is the coherent scattering length, m d the mass and R d the position vector in the real space of each atom d and σ d j ðqÞ are the phonons' polarisation vectors.
Time-of-flight data reduction were performed with Mantid analysis suite 67 and measurements for different rotation angles were combined using HORACE 46 , yielding four-dimensional SðQ; ωÞ datasets for the selected incident energies. HORACE also allows to slice these 4D datasets into 2D slices along specified trajectories in the ð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 Q values. In order to reduce the noise-to-signal ratio, data along each symmetry direction were integrated over the other two within a ±0.1 range (in reciprocal lattice vector units) centred around the Bragg peak (for longitudinal Γ-Z and transverse Γ-X and Γ-Y modes with E i = 7.3 meV it is necessary to integrate within a ±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 ±0.1 interval centred around the desired 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 ω j ðqÞ and normal modes polarisation vectors σ d j ð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 σ = 0.6 meV. For the simulation of 1D cuts in Fig. 2 we also added the contribution of the elastic signal.

Data availability
Raw data from the INS experiment were generated at ISIS Neutron and Muon Source and are available at doi.org/10.5286/ISIS.E.86390979 and doi.org/10.5286/ISIS. E.86390999. Derived data supporting the findings of this study are available from the corresponding authors upon reasonable request.