Orientational dynamics in supercooled glycerol computed from MD simulations: self and cross contributions

The orientational dynamics of supercooled glycerol using molecular dynamics simulations for temperatures ranging from 323 K to 253 K, is probed through correlation functions of first and second ranks of Legendre polynomials, pertaining respectively to dielectric spectroscopy (DS) and depolarized dynamic light scattering (DDLS). The self, cross, and total correlation functions are compared with relevant experimental data. The computations reveal the low sensitivity of DDLS to cross-correlations, in agreement with what is found in experimental work, and strengthen the idea of directly comparing DS and DDLS data to evaluate the effect of cross-correlations in polar liquids. The analysis of the net static cross-correlations and their spatial decomposition shows that, although cross-correlations extend over nanometric distances, their net magnitude originates, in the case of glycerol, from the first shell of neighbouring molecules. Accessing the angular dependence of the static correlation allows us to get a microscopic understanding of why the rank-1 correlation function is more sensitive to cross-correlation than its rank-2 counterpart.

Dielectric spectroscopy (DS) is a powerful tool for studying polar supercooled liquid dynamics [1][2][3].The outcome of the measurement, the complex dielectric permittivity ϵ(ω), contains a wealthiness of information regarding the collective orientational motion of the permanent dipoles of the constitutive molecules, and more precisely on the relaxation processes at work in the liquid under scrutiny [4,5].The broad range of available frequencies (10 −5 − 10 13 Hz) for a wide range of temperatures, allows to follow the slow down of the structural α relaxation upon cooling close to the glass transition temperature as well as the emergence of secondary relaxation processes such as JG processes, believed to be an intrinsic characteristic of glassy dynamics [6], or the excess wing, recently associated with dynamical facilitation [7].DS can also be used to characterize the cooperative nature of the α relaxation [8], to determine density scaling of the relaxation time [9,10], or to study physical aging of out-of-equilibrium liquids [11,12].
The complex dielectric permittivity obtained from DS measurements can be linked to the time-dependent equilibrium field free total dipole moment correlation function C 1 (t) through a Fourier-Laplace transform [13][14][15] where k B is Boltzmann's constant, ρ is the liquid density, T its temperature and ϵ ∞ is the permittivity at visible * Corresponding author: marceau.henot@cea.froptical frequencies.The dipole correlation function of rank ℓ is: where N is the number of dipoles in the cavity considered, P ℓ is the Legendre polynomial of rank ℓ and ϑ i,j (t 0 , t 0 +t) corresponds to the angle between molecule i at time t 0 and j at t 0 + t.In DS, this angle is measured between dipole moments and the technique is sensitive to the order ℓ = 1 leading to: The static value can be rewritten in: where g K is the Kirkwood correlation factor that can either be > 1, in which case the dipole-dipole correlation are overall positive, or < 1 meaning that antialignments dominates.The dynamics is also expected to be affected by cross-correlation because there is a priori no reason that the timescales and shapes of the self and cross-correlations function coincide exactly.A striking example of dynamical consequences of intermolecular correlations is the behavior of mono-alcohols which display another relaxation process at low frequencies, called the Debye peak, related to the formation of supramolecular H-bonded structures consisting of chains (g K > 1) or rings (g K < 1) [16].A recent theory from Déjardin et al. [17], showed that the liquid dynamics can be strongly affected by the effect of positive cross-correlations.
Recently, results from DS were compared to other techniques less sensitive to cross-correlations.The fluorescence response of a local probe diluted in a mono-alcohol was shown to be insensitive to the Debye relaxation of the liquid, allowing it to disentangle it from the other relaxation processes [18].Another technique that has been proven useful in that regard is depolarized dynamic light scattering (DDLS) which probes molecular orientations through the anisotropy of the polarizability.The relevant correlation function is given by eq. 2 with ℓ = 2.It follows that the technique does not distinguish between parallel and antiparallel alignments.There is strong experimental evidence that DDLS is insensitive to crosscorrelations.For example, Gabriel et al. [19] showed that in mono-alcohols DDLS displays an α peak but no Debye peak.In addition, in a non-associating liquid, Pabst et al. [20] showed that progressively diluting the system in a non-polar solvent leads the DS spectra to look more and more alike the DDLS spectra.All of this illustrates the importance of cross-correlation effects in DS which can significantly broaden the α peak.Moreover, while the shape of the α peak in DS spectra is system dependent, it was shown in DDLS to follow a generic line shape of slope −1/2 on the high-frequency flank [21].There is still debate, however, on whether this generic response reflects the true structural relaxation better than the dielectric one [22].
When dealing with physical processes taking place at the nanometric scale, molecular dynamics (MD) simulation is an attractive method that can give access to microscopic observables that are otherwise hard, or impossible, to obtain experimentally.This method is however limited to high temperatures or simplified systems, due to its computational cost.To study the generic behavior of liquid glass-formers, model systems can be thought of being made of polydisperse beads interacting through a Lennard-Jones potential.This helped give information on the spatio-temporal nature of relaxations [7,23].Another approach, more suitable for direct comparison with experiments, is to rely on a more precise modelization of specific molecules, taking into account their dipolar nature and electrostatic interactions.This gives access to their dielectric response [24][25][26][27][28]. Recently, MD simulations on a model dipolar system showed that, while the orientational ℓ = 1 correlation function of weakly polar systems is dominated by the self response, strongly polar liquids are much affected by cross-correlations [29].
The wide variety of organic liquids available has led to the choice of some systems, considered as models or representatives.Glycerol, by its apparent simplicity and its low tendency to crystallize, has long been the subject of extensive studies, by various techniques including dielectric spectroscopy [2,11,[30][31][32], neutron spectroscopy [33], nuclear magnetic resonance (NMR) [34], DDLS [35,36] and MD simulations [37][38][39][40][41][42][43][44].Its dynamics is, however, not particularly simple.As a trialcool, it is subject to H-bonds but does not display a Debye peak that would result from linear supramolecu-lar chains.Shear mechanical spectroscopy has shown the existence of a low-frequency mode that is believed to result from the hydrogen-bonded network formed between molecules [45].
In this article, we report a MD study of the orientational dynamics of glycerol, on a large temperature range (from 253 to 323 K) reaching the moderately supercooled regime, simulated from a model already widely used in the literature [37][38][39][40]44] over durations of up to 7 µs.We first compute the self response of the dipolar moment for ranks ℓ = 1 and 2 from which we deduce the loss function χ ′′ ℓ (f ) for frequencies down to 200 kHz.We then analyze the cross-correlation and we exploit the possibility offered by MD to decompose this part of the response as a function of the relative distance and orientation of the dipoles.We compute the total loss function for both ranks as well as the part resulting from crosscorrelations alone which allowed us to verify that crosscorrelations play a major role in the ℓ = 1 response while being almost negligible for ℓ = 2.We compare these data to experimental DS and DDLS spectra and obtain similar temperature dependence for the relaxation time and slope of the high-frequency flank of the α peak.We discuss how the differences in the spectra associated with different ranks can be related to the underlying molecular relaxation mechanisms.Moreover, we show that, for glycerol, the net cross-correlation originates only from the first shell of neighouring molecules.This is the case for both ℓ = 1 and 2 although their different sensitivity to orientational correlations leads, at the end, to significant differences in the importance of contributions coming from cross-correlations.

II. METHODS
The molecular dynamics (MD) simulations were performed using OpenMM [46] on an Nvidia RTX A5000 GPU.Glycerol has been modeled using the reparameterized AMBER force field previously employed in the literature [37][38][39][40][41][42]44] and whose parameters are given in the suppl.mat.Atoms belonging to the same molecule interact through harmonic potentials for bond length and angle and a periodic potential for bond torsion.Non-bounded atoms interact through a Lenard Jones potential with a 1 nm cutoff and a coulomb interaction computed using a Particle Mesh Ewald (PME) algorithm (1 nm cutoff and 0.0005 error tolerance).The simulation does not account for electronic polarizability.Each atom carries a constant partial charge originally derived by Chelli et al. [37] from quantum mechanical calculations.Later, Blieck et al. [39] noticed that this parameterization led, in the temperature range 333 -413 K, to a dynamics 10 times faster than measured experimentally by neutron spectroscopy.They slowed down the dynamics by the right amount by reducing by 5 % the hydroxyl group atomic charges.They also checked that the simulation reproduced fairly well the static structure factor measured by neutron scattering [47].This corresponds to a mean dipole moment of ⟨µ⟩ = 3.2 D which is higher than the µ exp = 2.68 D value measured in a nonpolar solvent [48].This can be seen as a way to compensate for the absence of electronic polarizability which leads, in the real system, through the reaction field, to an effective dipole moment greater than µ exp [4].The same parameters were later used by Egorov et al. [40] (who corrected slightly the charges to ensure the molecule neutrality, and made all bound flexible) to study glycerol-water mixtures and more recently by Becher et al. [44] to reproduce NMR spectra in the 300-540 K range.The parameters used in this work were almost identical with only small modifications intended to reduce the computational cost: the length of bonds involving hydrogen were fixed (as in refs.[37,39]) and hydrogen atom mass was increased by 40% allowing to use an integration time of 4 fs.The simulations were carried out on a system of N = 2160 molecules (30 240 atoms) in a cubic cell of side length a ≈ 65 Å with periodic boundary conditions (PBCs), in the NPT ensemble at eight different temperatures T (from 323 to 253 K) and at pressure P = 1 bar using a Monte Carlo barostat and a Nosé-Hoover thermostat.In order to study the effect of the system size, a simulation at T = 323 K was performed on a system consisting only of N = 540 molecules (a ≈ 41 Å).Random initial states were generated using Packmol [49], equilibrated at 323 K and progressively cooled down to 253 K by 10 K steps by waiting at each step an equilibration time corresponding to 98 to 200 τ α , reaching 7 µs (see details in the suppl.mat.).At each temperature, simulation runs lasted from more than 180 τ α for T ≥ 263 K and 67 τ α at 253 K (corresponding to 4 µs).For all simulation runs, the dipole of each molecule ⃗ µ i (t) (i ∈ [1, N ]) and its position ⃗ r i (t) were determined by computing the barycenter of the positive (q + at ⃗ r + ) and negative charges (q − = −q + at ⃗ r − ) with

III.1. Self correlation
The self dipole correlation function is defined by: It characterizes the molecular relaxation through rotational movement of the permanent dipole.This function is shown in fig.1a at each temperature, for rank ℓ = 1 and 2. Three regimes can be observed: at short time (t < 100 fs) a small decorrelation occurs and a boson peak is visible at t ≈ 70 fs.At long time, there is a complete, non-exponential decorrelation (C self ℓ (t) reaches 0) corresponding to the α relaxation.At intermediate times the correlation is high but slowly decreasing.This regime is almost nonexistent at 323 K but extends over two decades at T = 253 K.While the global shape is the same for ℓ = 1 and 2, the short time decorrelation appears more intense for ℓ = 2.This is simply due to the quicker decreases of P 2 (cos ϑ) compared to P 1 (cos ϑ) for ϑ ≪ 1.The mean self relaxation time is obtained from dt and is shown as a function of 1/T in blue in fig.6a.The relaxation times are shorter for ℓ = 2 (empty markers) than for ℓ = 1 (solid markers) and they both display a super-activated behaviour.
The self loss function was obtained, following eq. 1, by applying the fluctuation-dissipation theorem [50]: where TF is the Fourier transform, computed using the fftlog algorithm adapted to log spaced data [51].The fact that the correlation function was averaged on long times (≈ 100τ α ) leads to a fairly low amount of noise on the spectra, shown in fig.1b.They were all rescaled by superimposing their microscopic peaks at 10 13 Hz.The frequency at which the maximum of the α peak is reached was found to correspond (within the uncertainty) to 1/(2πτ self ℓ ).On the low-frequency side, the spectra follow a power law with slope 1, as expected.On the high-frequency flank of each spectrum, there is a power law regime on one to two decades in frequency with a slope −β self ℓ , interrupted by the fast process [3].The corresponding values of β self ℓ are shown in blue in fig.6b.For ℓ = 1 (solid markers), slope increases with temperature (ranging from 0.36 at 253 K to 0.46 at 323 K) while for ℓ = 2 (empty markers), it is temperature independent and systematically smaller (≈ 0.27).These low β values are associated with the non-exponential nature of the relaxation process.

III.2. Static cross-correlation
As stated in the introduction, experimental methods such as DS and DDLS are sensitive, not only to the self correlation function, but rather to a total correlation made of the self part and of a cross-correlation part.We thus need to get access to the correlation function associated to cross-correlation: However, one has to be careful with the application of this definition directly to the MD simulation box due to the effect of PBCs on the treatment of electrostatic interactions.With the PME method used here, our simulation box can be seen as wrapped in tinfoil, or embedded in an infinite medium in which the macroscopic electric field is null [26,28,52].This effect is responsible for a long-range dipole correlation of significant amplitude, that is an artifact of the simulation, and which cannot be suppressed or diminished by increasing the simulation box size (see fig.S1 of suppl.mat.).This artificial cross-correlation is maximum on average for couples of molecules separated by a distance of the order of the box size a.A way to get around this difficulty is to use a simulation box large enough to decouple the real correlation (occurring at relatively small distances) and the artifact [26,28].We decompose the cross-correlation function C cross ℓ (t) into contributions per unit distance denoted c ℓ (r, t) depending on the distance r = ∥⃗ r j − ⃗ r i ∥ between the reference molecule i and all molecules within [r, r + dr].This quantity, computed for slices of 0.5 Å and averaged over i and t 0 , is plotted in fig.2a for the static case (t = 0) for ℓ = 1 and 2. The dipole density n(r) is plotted in black in fig.2b alongside with the parabola in red that would be obtained for a homogeneous system of the same average density.The difference between these two curves shows a series of maxima that corresponds to the first, second and third neighbour peaks.
For ℓ = 1, it appears that the cross-correlation contribution c 1 (r) (in blue) is maximum for the first neighbours and reaches zero before the second layers of neighbours.Fig. 2c represents, as a function of r, the total contribution r 0 c(r ′ )dr ′ integrated within a sphere of radius r.We see that the cross-correlation reaches a plateau at r ≈ 7 Å while the cross-correlation coming from PCBs starts to be perceptible for r > 20 Å (see suppl.mat.).It is also interesting to study the mean level of static cross-correlation Γ(r) = c(r)/n(r), plotted in fig.2d that shows that the cross-correlation per dipole is positive, is a strictly decreasing function of the distance and is only of the order of 5-15 % in average for the nearest neighbours.From these data, we can deduce the Kirkwood factor g K = 1 + r lim 0 c(r ′ )dr ′ = 1.70 ± 0.02 at 273 K, with r lim = 7.5 Å.This is smaller than the 2.6 ± 0.2 value reported in the literature and deduced from static permittivity measurements [36] but it is compatible with previous numerical results on glycerol for T > 250 K [43].It is also interesting to note that our value of µ 2 g K (which is the quantity accessible from the experiments, see eq. 1) matches exactly the experimental value and displays the same temperature dependence (cf fig.S2 in suppl.mat.).In other words, the simulation gives the expected value of ϵ(0) in the whole temperature range.However, as the dipolar moment of glycerol has been measured with a reasonable accuracy in a non-polar solvent [48], it is likely that the MD simulation underestimate the value of g K .For ℓ = 2, the cross-correlation (in orange) appears to be also due to the first layer of neighbours but is much less intense than for ℓ = 1 and its integrated value saturates for a slightly higher radius.With r lim = 9.5 Å, the quantity analogous to the Kirkwood factor for ℓ = 2 would be only 1.11 ± 0.01.
The dipolar interaction is anisotropic in nature and it makes sense, rather than averaging the cross-correlation over all dipoles situated at a given distance, to distinguish the contribution as a function of the relative orientation.Spherical coordinates (r, θ, ϕ) can be defined with respect to the reference dipole i where θ is the angle between ⃗ µ i and ⃗ r = ⃗ r j − ⃗ r i .By symmetry, the contribution should not depend on the azimuthal angle ϕ.The cross-correlation contribution c ℓ (r, θ) per unit distance and solid angle is shown in the static case in fig.3a  (ℓ = 1) and b (ℓ = 2) as a function of r and θ.For ℓ = 1, similarly to a recent observation in water [28], the spatial distribution of cross-correlation appears strikingly different than the θ averaged curves shown previously.The correlation is positive in an angular sector situated above and below the dipole of reference | cos θ| < cos θ lim and mostly negative on its sides.It is null on lines of constant | cos θ| = cos θ lim shown in grey with θ lim = 52°.The cross-correlation contributions summed over these two angular sectors are shown in fig.3c (top).These positive and negative contributions extend way over what is visible when considering their sum and, while decreasing, are far from negligible at r = 16 Å.However, these contributions cancel each other for r > 7.5 Å.This an-  gular dependence of the cross-correlation is what is expected when considering the energy interaction between two electrostatic dipoles as a function of their relative orientation [53] and drawings of the most favorable dipoles orientations are shown in fig.3a.The net contribution of cross-correlation comes from the first shell of neighbouring molecules.Close molecules situated above and below are strongly positively correlated (even more above than below) with an alignment rate reaching 50 %.This correlation is favored by the dipole-dipole interaction although the top-bottom asymmetry illustrates that at such a close distance, it is not the only interaction playing a role.On the sides, molecules belonging to the first shell but further than 5.2 Å (corresponding to the first neighbor peak) are on average anti-aligned (again as favored by the dipole-dipole interaction), although not enough to compensate for the positive correlation.Finally, side molecules closer than 5.2 Å are positively aligned.This means that there exists other effects (that may be related to constraints constraints on molecule conformation or to the presence of H-bonds) able to compensate for the a priori unfavorable situation of having barycenters of the same sign charges facing each other on average.For ℓ = 2 also, a complex spatial dependence of the cross-correlation is visible in fig.3b.As this quantity is not sensitive to the correlation sign, it shows a very different behavior than ℓ = 1.Some oscillations, that correlate well with the density inhomogeneity, are visible.Those also appear for ℓ = 1 but dominate here.After the first layer of neighbours, the contribution from the two angular sectors are in anti-phase and cancel each other on average.Similarly to the ℓ = 1 case, the net contribution comes from the first neighbour shell but it is interesting to observe, that its different sensitivity to orientational correlations makes this quantity significantly less sensitive to cross-correlations.

III.3. Global dipole dynamics
For the reason mentioned above and related to the effect of PBCs, the cross-correlation function C cross ℓ (t) is computed in the following from eq. 6 by considering only the molecules located within a sphere of radius r lim rather than in all the simulation box.The resulting normalized cross-correlation function is shown in green in fig.4a (ℓ =  but with only a 15 % increase on average between the self and total mean relaxation times for ℓ = 1 and a 25 % increase for ℓ = 2.This means that the Kivelson and Madden relationship [54] (τ tot = g K τ self ) does not seem to be verified in glycerol contrary to other systems such as water [55].
In the same way as for the self part, a loss function χ ′′ (f ) can be computed from the cross part and for the total correlation.The results are plotted for 273 K in fig.4c (ℓ = 1) and d (ℓ = 2).It is well visible that the amplitude of the cross-correlation is of the same order as the self part for ℓ = 1 while it is much smaller for ℓ = 2. On each spectrum, the high frequency flank of the α peak has a power law slope −β such that β self  appears constant close to 0.3.This is already visible on the spectra of fig.5a but it is even clearer in fig.5b where the spectra are shown as a function of a dimensionless frequency ωτ self 1 .The collapse of the high-frequency side is much better for the total spectra with ℓ = 2 (although the α peak slightly broadens upon cooling) than for ℓ = 1 (self or total).

III.4. Comparison with DS and DDLS experiments
The spectra χ ′′ tot 1 (f ) of fig.5a can be compared to experimental glycerol dielectric spectra.The general allure corresponds well to the measurements of Lunkenheimer & Loidl [2] with a clear α peak showing no excess wing in this range of temperature and a small Boson peak around 10 12 Hz.The temperature dependence of the α relaxation time is compared in fig.6a, with DS data as grey solid circles and corresponding MD data τ tot 1 in red.They show the same trend and could be fitted with a WFL law with close parameters (up to a vertical prefactor).However, the MD α relaxation time is systematically shorter than its experimental DS counterpart by a factor 2 to 3. It is important to note the absolute value of the relaxation time is affected in the simulation by the choice of µ which was optimized by Blieck et al. [39] on neutron scattering data for temperatures larger than 333 K (it was already visible that at the smallest temperature simulated by the authors, 313 K, the relaxation time was underestimated).It seems that this is also the case for other studies using the same parameters [44].The temperature evolution of the slope β tot 1 is also comparable in MD and in experiments (for which it comes from a Cole-Davidson fit to the spectra), as shown by grey solid circles in fig.6b but the slopes are systematically underestimated by the simulation.
Contrary to DS, which probes the reorientational dynamics of the molecules by following the permanent dipolar moment, DDLS does it by following the anisotropy of the polarisability tensor.Glycerol molecules are be- lieved to rotate as a rigid entity [44] and we can reasonably assume that both techniques probe the same dynamics (although giving access to different ranks ℓ).It is also worth noting that results from DDLS experiments can also be affected by a scattering mechanism called dipoleinduced-dipole related to the fluctuation of the internal field.Cummins et al. [56] have argued that this effect may be neglected in the experimental spectra of supercooled liquids.Here we follow these authors.The spectra χ ′′ tot 2 (f ), from the simulation, should therefore be comparable with the Fourier transform of the DDLS correlation function reported by Gabriel et al. [36].The authors concentrated mainly on frequencies under a few MHz and thus on temperatures smaller than 260 K (although they also performed a measurement at 323 K).Here also, the MD spectra show a good qualitative agreement with experiments.The DDLS mean relaxation time is shorter than in DS by an amount that appears similar in the experiments and in MD.The inequality β tot 2 < β tot 1 is also verified for both.In the DDLS experiments (see grey empty circles in fig.6d), given the uncertainty, the slope β tot 2 does not appear to depend much on the temperature and it is also the case in the simulation (empty red circles).However, here again, the slopes are underestimated in the simulation.
The differences between the simulation and the real system must be weighed against the relative simplicity of the modeling.Indeed, force field parameters, that control intra and inter-molecular interactions were not adjusted specifically on glycerol but were designed to be applicable  [31,36] for DS and refs.[35,36] for DDLS are shown in grey.
to the widest range possible of organic compounds.The only parameter that was adjusted specifically to glycerol was µ.Moreover, the partial charges in the molecule are considered as fixed point charges at the center of each atom, the electronic polarizability is not taken into account and H-bonds are mimicked only by electrostatic interactions of these fixed charges which can limit their strength [57].Nonetheless, it is already impressive that a classical model can reproduce well some aspects of the real system.With these limitations in mind, MD simulations can be used, as demonstrated recently by Becher et al. [44], to obtain information on microscopic observables or on the relative effect of external parameters such as the temperature.For ℓ = 2, and in contrast to the ℓ = 1 case, the effect of cross-correlation is very weak and the total loss function can reasonably be assimilated to the self loss function alone (see fig. 4d).This is in complete agreement with the observations of Gabriel et al. [36] and previous work by the same team [19,20] on the ability for DDLS to give access to the self orientational dynamics.This is also in agreement with the recent MD observations on a model system of Koperwas et al. [29] who compared the self and total correlation functions for ℓ = 2 at long times (t ≥ τ α ).
The ratio of the self relaxation times τ self 1 /τ self 2 , shown in fig.6c, ranges from 1.4 to 2.0 and is increasing with temperature.This quantity is affected by the details of the molecular relaxation mechanism with two limiting cases leading to values of 3 for isotropic rotational diffusion and 1 for discrete jumps of random amplitude [5].However, a given value of this ratio cannot be associated with a typical orientational jump amplitude.For example, in a mean field model, tuning the interaction parameter allows to change continuously this ratio [58].In the present case, it is more likely that its value reflects a dynamics governed by rare relaxation events, that are less averaged for ℓ = 2 than for ℓ = 1 [59].In this case, its decrease upon cooling could be associated with rarer and rarer events.The slope β < 1 can also be linked with a dynamic governed by relatively rare events, spread on long time scales, that does not average enough to produce an exponential relaxation (β = 1) [59].In this context, as C 2 cancels out for smaller angles than C 1 , it is less averaged over relaxation events leading to β 2 < β 1 , consistently with our observations.
In the DDLS works mentioned above [19,20,36], the authors suggested to use these measurements as a proxy to approach the DS self response, which is not accessible experimentally.This means admitting that the response is reasonably independent of the rank ℓ but this also requires a way to normalize the DDLS data so that they can be quantitatively compared to DS spectra.With a well-chosen, temperature independent, scaling factor, it is possible a get a perfect overlap of the DS and DDLS spectra on the excess wing flank at low temperature (T < 200 K).From this, it was shown that the DS signal could be well fitted by the sum of the DDLS signal and a cross-correlation term (described by a stretched exponential of fixed stretching parameter), with the relative weight of these two terms perfectly matching the Kirkwood correlation factor at all temperatures [36].This demonstrates the usefulness and relevance of the direct comparison between DDLS and DS spectra.
In our numerical work, we were not able to reach the range of low temperatures in which the excess wing appears.This prevented us from using the method described above to determine from scratch a scaling factor between the ℓ = 1 and ℓ = 2 spectra.We instead had to rely on the experimental determination by Gabriel et al. [36].Indeed, while the DS and DDLS spectra collapse on the excess wing for T < 200 K, they cross at a finite frequency f cross for higher temperatures.We determined a (T independent) scaling factor by making sure that the ℓ = 1 and ℓ = 2 total spectra cross at the same f cross /f α than the experiments for T = 323 K and T = 263 K (see vertical dashed lines in fig.5c).The result is shown in fig.5c and it appears that, for all temperatures, the α peak of the total ℓ = 2 and the self ℓ = 1 spectra match fairly well.The difference of slopes β is small enough so that the discrepancy remains low (< 30 %) under the frequency at which the fast process starts to be percept-ible.This allows us to understand the success of the experimental approach consisting of a direct comparison between DS and DDLS data [19,20,36].

IV. CONCLUSION AND PERSPECTIVES.
In this work, we studied using MD simulations, the orientational dynamics of glycerol from which we extracted the self correlation function and the associated loss function for different ranks ℓ of the Legendre polynomial.For ℓ = 1 and 2, we studied the spatial dependence of the cross-correlation and showed that they play a significant role in the ℓ = 1 response while being almost negligible for ℓ = 2.In accordance with recent experimental observations based on a comparison between DDLS and DS spectra, we showed that, although these techniques give access to different ranks ℓ of the correlation function (and consequently do not lead to the exact same spectra, in particular regarding the slope β), the scaling factor that corresponds to a merging of the excess wing of DS and DDLS spectra at low temperatures leads to a fairly good merging of the self part of the ℓ = 1 and the total ℓ = 2 spectra.This strengthens the idea that useful information can arise from a direct comparison between DDLS and DS measurements.
Moreover, we took advantage of the possibility given by MD simulations to access molecular observables to discuss the microscopic origin of the cross-correlation observed in DS.We found out that the net cross-correlation originates from the first shell of neighbouring molecules that tend to positively align independently of their orientation.Investigating in more detail the molecular origin of the cross-correlations, and their link with molecular conformation and H-bonds will be the subject of future work.
Table S4.Lennard-Jones parameters.The potential is of the form VLJ = 4 σ r 12 − σ r 6 .The input parameters of the Amber .frcmodfile is the half atom-atom distance at which the potential reaches its minimum (Rm = 2

Figure 1 .
Figure 1.(a) Dipole self correlation functions for the Legendre polynomial of rank ℓ = 1 (top) and ℓ = 2 (bottom) at different temperature T ranging from 323 K to 253 K with 10 K steps.(b) Dielectric loss function corresponding to the self part of the correlation functions ℓ = 1 and 2. The black dashed lines correspond to a power law fit of slope −β self ℓ on the highfrequency wing.

Figure 2 .
Figure 2. Distance dependence of the dipole static crosscorrelation function (ℓ = 1 and 2) at T = 273 K. (a) Contribution c per unit distance to the cross-correlation of all the molecules situated at r (b) Number of molecules per unit distance at r.The red curve corresponds to a homogeneous medium of the same density.The vertical grey lines on all the plots correspond to the first, second and third neighbour peaks.(c) Contribution to the cross-correlation of all the molecules within a sphere of radius r.(d) Mean level of crosscorrelation Γ (∈ [−1, 1]) for the molecules at r.

Figure 3 .
Figure 3. Relative orientation and distance dependence of the static dipole cross-correlation functions at T = 273 K. (a) Contribution c per unit distance and solid angle to the crosscorrelation (ℓ = 1) of all the molecules situated at r and at angle θ.Red zones are positively correlated while blue zones are anti-correlated as illustrated by the grey dipole drawings.(b) Same plot than a. for ℓ = 2. (c) Same quantity as in fig.2a, for ℓ = 1 (top) and ℓ = 2 (bottom), but with distinguished contribution from the red sector (| cos θ| > 0.62, see grey lines on a,b) and from the complement blue sector.Drawings close to the vertical axes illustrate the physical meaning of the correlation sign for ℓ = 1 (right) and ℓ = 2 (left).
1) and b (ℓ = 2) at T = 273 K.The self part is shown in blue and the total correlation function C tot ℓ (t) = C self ℓ (t) + C cross ℓ (t) is shown in red.The cross part does not display a short time decorrelation but for ℓ = 1, a small amplitude peak is visible at short time (t < 1 ps) (see inset of fig.4a).The mean cross τ cross ℓ and total τ tot ℓ relaxation time, computed by integrating the corresponding normalized correlation functions, are shown for all temperatures in green and red in fig.6a.They follow τ self ℓ < τ tot ℓ < τ cross ℓ

ℓ < β cross ℓ < 1 .
The slope of the total loss function β tot ℓ results from both previous slopes as well as the strength of the crosscorrelation and takes an intermediate value.The total spectra for all temperatures are shown in fig.5a and the values of the slopes β are shown in fig.6b as a function of T : β tot 1 is increasing with T while β tot

Figure 4 .
Figure 4. Correlation and loss functions at T = 273 K. (a, b) Normalized self, cross and total correlation function for ℓ = 1 (a) and 2 (b).The inset of (a) is a zoom at short time with a linear scale.(c, d) Corresponding loss functions for ℓ = 1 (c) and 2 (d).The dashed lines in (a) defines the slopes β cross ℓ

Figure 5 .
Figure 5. Orientational loss function for all temperatures (same color scale than in fig. 1) studied as a function of f (a) and ωτ self 1 (b).The total loss function is shown for ℓ = 1 (top, left scale in a) and ℓ = 2 (bottom, right scale in b).The self part for ℓ = 1 is shown in (a) only for the extreme temperatures and for all temperatures in (b) (middle).(c) For several temperatures, ℓ = 2 total spectra (orange) plotted with a scale factor (chosen to correspond to the one determined experimentally by Gabriel et al. [36]) alongside with the ℓ = 1 self (blue) and total (red) spectra.The vertical dashed lines shows where the crossing between DS and DDLS data occurs on experimental data.

Figure 6 .
Figure 6.(a) Mean relaxation time obtained by integrating the self (blue squares), cross (green triangles) and total (red circles) dipole correlation function for ℓ = 1 (solid markers) and ℓ = 2 (empty markers).Markers are linked to improve readability.Relaxation time measured experimentally by DDLS (empty circles) from refs.[35, 36] and DS (solids circles) from refs.[2, 31, 36] are shown in grey.(b, d) Power law exponent β of the high-frequency wing of the loss function for ℓ = 1 (b) and ℓ = 2 (d).The color code is the same as for (a).(c) Temperature dependence of the ratio of self mean relaxation time for ℓ = 1 and ℓ = 2.In b. and d. experimental data, from refs.[31,36] for DS and refs.[35,36] for DDLS are shown in grey.