Excitation spectra in fluids: How to analyze them properly

Although the understanding of excitation spectra in fluids is of great importance, it is still unclear how different methods of spectral analysis agree with each other and which of them is suitable in a wide range of parameters. Here, we show that the problem can be solved using a two-oscillator model to analyze total velocity current spectra, while other considered methods, including analysis of the spectral maxima and single mode analysis, yield rough results and become unsuitable at high temperatures and wavenumbers. To prove this, we perform molecular dynamics (MD) simulations and calculate excitation spectra in Lennard-Jones and inverse-power-law fluids at different temperatures, both in 3D and 2D cases. Then, we analyze relations between thermodynamic and dynamic features of fluids at (Frenkel) crossover from a liquid- to gas-like state and find that they agree with each other in the 3D case and strongly disagree in 2D systems due to enhanced anharmonicity effects. The results provide a significant advance in methods for detail analysis of collective fluid dynamics spanning fields from soft condensed matter to strongly coupled plasmas.

calculation of excitation spectra in fluids, we performed MD simulations whose details are described in Methods. We considered several model systems at the same density and different temperatures, including a Lennard-Jones fluid, and fluids of particles interacting via inverse power law repulsion ∝1/r 12 (IPL12), and a more soft variant ∝1/r 8 (IPL8), to analyze features arising in systems with different kinds of interactions. We used the Lennard-Jones system as a representative and well-studied model describing noble gases 5,36,39,77 , while the inverse-powerlaw fluids [5][6][7]21,26 were considered to compare the results with the Lennard-Jones system and, hence, to reveal features associated with attraction and repulsion softness.
For the systems under consideration, after particle velocities were obtained from MD simulations, we calculated the velocity current spectra 9 is the velocity of the s-th particle; the summation is performed over all N particles in the system; = ⋅ q j q j q ( )/ L 2 and = ⋅ ⊥ ⊥ j j e e ( ) T are the longitudinal and transverse components of the current, respectively; ⊥ e is a unit vector normal to q; the brackets 〈…〉 denote the canonical ensemble average.
Since simple fluids are isotropic, velocity current spectra depend only on the frequency ω and wavenumber = | | q q . Therefore, we can average ω C q ( , ) L T , over all directions of the wavevector q to suppress noise caused by a limited number N and finite simulation time, where N q is the number of directions used for averaging. Note that this approach is not suitable at  π q 2 /L, where L is the size of the simulation box.
To obtain excitation spectra using ω C q ( , ) L T , , we assumed that the velocity current spectra corresponding to damped oscillating excitations are Fitting (separately) velocity current spectra (2) (obtained from MD) by the Lorentzian fit (4), we calculated the frequencies and damping rates of longitudinal and transverse excitations (single mode analysis).
Another approach to the same problem (which is, however, rarely used) is related to the joint mode analysis (two-oscillator model), at which the total velocity current spectrum is fitted by the sum of high-and low-frequency double-Lorentzian terms where D is the space dimension.
In the both methods, using single mode analysis and two-oscillator model, we obtained dispersion relations ω q ( ) L T , and damping rates Γ q ( ) . To compare the methods, we present in Fig. 1 our results for the Lennard-Jones fluid that were calculated at temperatures (a)-(d) 1.83 and (e)-(h) 36.41 (in Lennard-Jones units). Figure 1(a) presents the total velocity current spectra ω C q ( , ) in a color-coded format (normalized to the maximum at each q). Blue rhombs correspond to the dispersion relations ω q ( ) L T , calculated using single mode analysis of longitudinal (LM) and transverse (TM) modes. Results obtained using the two-oscillator model are shown by red circles for ω q ( ) L T , and by triangles for ω ± Γ q ( ) L T L T , , to illustrate a relation between damping rates and oscillation frequencies. The damping rates obtained from single mode analysis are almost the same as those obtained www.nature.com/scientificreports www.nature.com/scientificreports/ using the two-oscillator model and, thus, we do not show them in Fig. 1(a,e). Numerically-obtained profiles of ω C q ( , ), as well as their longitudinal and transverse components ω C q ( , ) L T , , at different wavenumbers are shown in Fig. 1(b-d) by symbols, while the solid lines are fits by Eqs. (4) and (5). The description of the results shown in Fig. 1(e-h) is the same as in panels (a)-(d).
It is noteworthy that the both approaches yield close results, if the total current spectrum ω C q ( , ) has two pronounced maxima (at lower and higher frequencies, which are typically associated with longitudinal and transverse excitations, respectively). This corresponds to lower temperatures and wavenumbers corresponding to the first Brillouin pseudo-zone. Discrepancy between dispersion relations obtained using different methods increases with the temperature, as seen in Fig. 1(a,e). One can see in Fig. 1(b-d,f-h) that the current spectra obtained from MD simulations are described well by theoretical fits (4) and (5). The worst results are observed for excitations at high temperatures and short wavelengths, as indicated in Fig. 1(h). This is related to change in the shape of velocity current spectra in the high-q limit (see corresponding discussion below). Results of analysis using maxima of ω C q ( , ) T is presented in the next section. Note that the single mode analysis becomes unsuitable if longitudinal and transverse modes cross, as can be seen at π . −  qn / 19 1/3 in Fig. 1(a). This is because structural disorder in a fluid mixes longitudinal and transverse modes 66,67 , leading to an effective interaction between the modes and resulting in their anticrossing accompanied by the strong redistribution and hybridization of spectra 78 . Strictly speaking, due to the mode anticrossing, highand low-frequency hybridized modes with mixed polarizations appear instead of longitudinal and transverse modes. Mode anticrossing should be properly taken into account in the q-region where the modes cross each other, and the detailed analysis of this problem stands beyond the scope of this work. However, we should stress that the mode mixing cannot be taken into account properly using single mode analysis in principle, comparing to the two-oscillator model.  T 36 41 and = n 1, description is the same as in panels (a-d), respectively. www.nature.com/scientificreports www.nature.com/scientificreports/ tions in thermodynamic equilibrium. The corresponding wavenumber range forms a gap in the q-space, which is referred to as the q-gap and is already clearly seen in Fig. 1(a,e). The long-standing problem of the calculation of the q-gap has recently attracted considerable interest in view of crossover between liquid-and gas-like collective dynamics in fluids with change in temperature and pressure 8,26,32,[34][35][36][37][38][39][40][41][42][43][44][45][47][48][49][50][51][52][53][54][55][56][57][58] . Here, we present the calculated excitation spectra near the q-gap and compare the dispersion relations obtained using the two-oscillator model to the results based on the analysis of transverse current spectra.
Analysis of ω C q ( , ) L T , maxima at given q-values is widely used to obtain dispersion relations ω q ( ) L T , in crystals (see, e.g., refs 8,[32][33][34][35][36][38][39][40][41]52,53,58 ). Such an approach is suitable in this case, since the current spectra C L,T in crystals typically consist of narrow bands on the ω q ( , ) plane, while the damping rates are negligible and, thus, are assumed to be equal to zero. On the contrary, the high-and low-frequency branches of the current spectra ω C q ( , ) L T , in fluids 21,26,66,67 or crystals with developed anharmonicity 7 are significantly broaden and overlapped with each other. For this physical reason, we should use the two-oscillator model to calculate excitation spectra, since other approaches become unsuitable.
As seen, the two-oscillator model enables in-depth analysis of excitations spectra. For instance, using model (4) for the transverse part of velocity current spectra, we readily derive the frequency ω m corresponding to the maximum of ω C q ( , ) T at a given q value (the dispersion relation based on the analysis of where overdamped excitations become non-overdamped, one can expand ω m 2 in a q-series and obtain where q * is the wavenumber at which Γ T /ω = 3 T and * means the values at = ⁎ q q . However, a posteriori analysis of the excitation spectra in fluids justifies that it is sufficient to use the first-order terms in the expansions for ω T and Γ T , since the second q derivatives are small at = ⁎ q q and can be neglected. Figure 2 presents the low-frequency branches of dispersion curves in the 3D Lennard-Jones fluid at the density = n 1 and temperatures (a) 2.85 and (b) 6.03. The curves were obtained using the two-oscillator model (orange circles), maxima of the transverse part of current spectra ω C q ( , ) T in the two-oscillator model (black pentagons), and maxima of transverse current spectra ω C q ( , ) T (blue rhombs). The black solid and dashed lines correspond to Eq. (7) and its linearized variant with β ≡ ⁎ 0, respectively, while the frequencies, damping rates, and q-derivatives were provided by the two-oscillator model (5). Figure 2 demonstrates the following three domains in the q-space: , the range of missing oscillating excitations, and domains of (ii) overdamped excitations The (transition) region of overdamped oscillating excitations is shown in Fig. 2 by the gradient blue area.
At low temperatures, ω q ( ) T obtained using maxima of the ω C q ( , ) T part of the two-oscillator model is close to that obtained from the analysis of transverse current spectra, as shown in Fig. 2(a). At high temperatures, the two-oscillator model still yields dispersion curves ω q ( ) T with a clear q-gap, while the method based on ω C q ( , ) T maxima provides irregular points and becomes inappropriate. The observed irregular behavior is related to the flat form of ω C q ( , ) T T values, at which even weak data noise (due to thermal fluctuations and limited statistics) can strongly affect a position of the detected maximum of ω C q ( , ) T (e.g., it is clearly seen in ref. 53 ). Second, the high-frequency branch of excitations is significantly broaden and, therefore, also affects the positions of the ω C q ( , ) T maxima due to the mode mixing. Figure 2 proves that approaches based on maxima are unsuitable for the proper analysis of the true q-gap where ω = 0 T . Comparing results represented by orange-and black-colored symbols, one can see that the discrepancy between ω T and ω m is maximal near the q-gap and exhibits a systematic difference of about 50% between q g and q * . Therefore, it is crucially important to properly include damping in the analysis of excitation spectra when the frequency and damping rate are comparable with each other. To accurately calculate q g , using discrete values of ω q ( ) T obtained within the two-oscillator model, we fitted ω q ( ) T in the region near the gap boundary by the linear dependence: (2019) 9:10483 | https://doi.org/10.1038/s41598-019-46979-y www.nature.com/scientificreports www.nature.com/scientificreports/ where θ q ( ) is the Heaviside step function and c T characterizes the velocity of the overdamped low-frequency excitations. Results of q-gap calculations in the Lennard-Jones, IPL12, and IPL8 fluids at different temperatures are presented and discussed below.
Concluding this section, note once again that the two-oscillator model can be used to obtain accurately the parameters of the square root relation ω q ( ) m (7) near the transition between overdamped and non-overdamped domains of excitations. Namely, the square root relation has been typically discussed in the context of collective excitations in fluids (see, e.g. 27,53 and references therein).
Excitation spectra at high temperatures. To elaborate the results of previous sections, we studied in detail the dynamics and specific heat of the 3D Lennard-Jones fluid at the density = n 1 in a wide range of temperatures. Using different approaches to analyze excitations spectra, we aimed to analyze changes in thermodynamics of the system, as well as collective and individual dynamics of particle motion. The results are summarized in Fig. 3. Figure 3(a-c) present the dispersion relations of the excitation spectra at temperatures = .
T 17 2, 20.0, and 23.2, respectively. The results are shown in the same manner as in Fig. 2, together with the linear asymptotic curve for the high-frequency (acoustic) branch: where P is the pressure and C V and C P are the specific heats at constant volume and pressure, respectively. The specific heat at constant pressure can be calculated from MD data as: According to Fig. 3(a-c), Eq. (9) accurately describes the long-wavelength behavior of the high-frequency branch of excitations, typically associated with sound propagation. Calculated adiabatic indices γ = C C / P V used in Eq. (9) are shown in Table 1. Trends in discrepancies between different approaches to the analysis of excitation spectra are the same as in the previous section.
The results in Fig. 3(a-c) are shown for temperatures chosen to discuss crossover from liquid-to gas-like dynamics of fluids with increasing temperature (which is referred to as "Frenkel line 8 "). There is a set of features www.nature.com/scientificreports www.nature.com/scientificreports/ used in different previous works and attributed to the transition from the liquid-to gas-like state of fluids 8 : (i) individual dynamics feature: the velocity autocorrelation function (VAF) becomes monotonically decreasing rather than oscillatory, (ii) collective dynamics features: positive sound dispersion vanishes, transverse excitations extend beyond the first Brillouin pseudo-zone in the fluid, and (iii) thermodynamic feature: due to the absence of transverse excitations, the specific heat per particle at constant volume (in the quasi-harmonic approximation) becomes To discuss relations between these features, we presented in Fig. 3(d-f) results for the specific heat, velocity autocorrelation functions, and calculated q-gaps at different temperatures, respectively. Note that, since an unambiguous definition of the first Brillouin pseudo-zone in fluids is absent, its boundary can be estimated in the range π π < < − qn 1 / (6/ ) 1/3 1 /3 in the 3D case shown in Fig. 3(f) by dashed black lines. Analysis of Fig. 3(a-c)  and ω ≠ 0 m ) extend beyond the first Brillouin pseudo-zone with a farther increase in the temperature. Positive sound dispersion vanishes at the same temperatures, since orange symbols at = .
T 17 2 lie slightly above the red line (9), while this is not the case at ≥ . T 20 0. Moreover, these results agree with the analysis of velocity autocorrelation functions shown in Fig. 3(e).  www.nature.com/scientificreports www.nature.com/scientificreports/ However, the situation with all excitations, including overdamped ones, is different, since Fig. 3(a-c)  To explain this statement, we showed in Fig. 3(f) the temperature dependencies of q * (black pentagons) and q g (orange circles). Here, the vertical dashed blue line corresponds to the temperature at which = C 2 V , while the blue region marks the temperature range where the form of velocity autocorrelation function becomes monotonic instead of oscillating. The main conclusion here is that different features attributed to the Frenkel line agree with each other in general, but the collective dynamic features should be reformulated only for non-overdamped excitations, since stable transverse excitations from the first-pseudo zone do not necessarily vanish.
We found that the discussed trends are fairly generic and qualitatively the same results, as in Fig. 3, are observed for fluids with other interactions between particles. In particular, we considered fluids with inverse-power-law (IPL) repulsions ∝1/r 12 (IPL12) and ∝1/r 8 (IPL8). The IPL12 fluid is a Lennard-Jones fluid without the attractive branch ∝1/r 6 , which makes it possible to understand the role of the attractive branch of the potential in excitations of fluids. Then, comparison of results for the IPL12 and IPL8 fluids allows us to analyze trends when repulsive interactions become softer. Figure 4 shows our results for the (a) Lennard-Jones, (b) IPL12, and (c) IPL8 fluids, presented in the same manner as in Fig. 3(f) (details of the description are the same). Figure 4(a) coincides with 3(f) and is reproduced here to highlight the trends appearing under the variation of the interaction. The comparison of Fig. 4(a,b) indicates that the absence of the attractive branch of interactions leads to discrepancy between thermodynamic and individual dynamic features at the Frenkel line ( = C 2 V and change in the form of the velocity autocorrelation function). When the repulsion becomes softer, agreement is recovered, but all low-frequency excitations leave the first Brillouin pseudo-zone much before the velocity autocorrelation function becomes monotonic, as seen in Fig. 4(b,c). It is remarkable that the form of temperature dependencies of the q-gap, q T ( ) g , obtained from the analysis of excitations within the two-oscillator model, are weakly sensitive to the kind of interaction. Correspondence between different approaches to the Frenkel line identification requires additional detailed study.

Excitations in 2D fluids: Enhanced anharmonicity at work. Discrepancy between different dynamic
and thermodynamic features attributed to the Frenkel line increases dramatically in 2D systems due to the enhanced role of anharmonicity. To illustrate this, we simulated the 2D Lennard-Jones fluid in a wide range of temperatures (see details in Methods). Then, the velocity currents, excitation spectra, specific heat, and velocity autocorrelation functions at different temperatures were analyzed in the same manner as in the case of the 3D Lennard-Jones system. The results color-coded in the same manner as in Fig. 3, are shown in Fig. 5.
We focused on the temperature region around , which is observed at .  T 12 15, as seen in Fig. 5(d). However, the excitation spectra shown in Fig. 5(a-c) clearly prove that a lot of non-overdamped excitations still exist in the first pseudo-zone at = T 12, while positive sound dispersion indeed disappears near the temperature corresponding to The same conclusions are confirmed by the temperature dependencies of q * and q g shown in Fig. 5(f). Figure 5(e) justifies an oscillatory behavior of the velocity autocorrelation functions at = T 12 and 15, while it becomes monotonic only at  T 40 (not shown here). Thus, in contrast to the 3D Lennard-Jones system, where thermodynamic, collective, and individual dynamic features of Frenkel crossover agree with each other, the situation in 2D fluids is essentially different.
This discrepancy has clear physical reasons related to anharmonicity and can be explained as follows. Fluctuations in 2D systems play a much more significant role than in 3D ones. In particular, these fluctuations lead to the independent behavior of translational and orientation order parameters and rich scenarios of melting 79 . Enhanced fluctuations enlarge the amplitudes of particle motions in 2D systems, increasing the anharmonicity of direct interactions between the nearest particles. For this reason, the potential part of the interaction energy (per degree of freedom) becomes significantly smaller than 1/2. As a result, the thermodynamic feature = C 3/2 V , which was derived in the quasi-harmonic approximation for the analysis of 3D systems 8 , becomes unsuitable for 2D fluids. www.nature.com/scientificreports www.nature.com/scientificreports/ In some sense, 2D fluids of softly interacting particles behave similarly to 3D fluids of hard-sphere particles, where they move mainly beyond the interparticle interaction potential. Developed anharmonicity leads to a rapid decrease in the specific heat with increasing temperature, which is similar to that observed in strongly-anharmonic crystals 7 . However, the analysis of the behavior of the q-gap in 2D systems shows that it exhibits an opposite trend, which is similar to 3D fluids with soft interactions between particles. The comparison of different features associated with the Frenkel line is beyond the scope of this work and should be considered in future.

Excitations at large q values: Getting back to individual dynamics. To complement our results
previously obtained for long-wavelength collective fluctuations in 3D Lennard-Jones fluids shown in Figs. 1 and  3, we considered velocity currents in a regime associated with individual particle dynamics, i.e., at high q values 65 . To illustrate change in velocity currents at the transition from collective to individual dynamics, we considered excitations at the temperature = .
T 36 41 and density = n 1. Figure 6 presents the sections of longitudinal (LM) and transverse (TM) modes of velocity currents ω C q ( , ) L T , . The dispersion relation for a free particle is ω = vq, where v is the particle velocity. In thermodynamic equilibrium, the velocity distribution has a Maxwellian form and, hence, for the velocity squared (corresponding to free particle dynamics), we have where m is the particle mass. Note that the maximum of the distribution (11) at each fixed q value is reached at is the most probable velocity 26,58 . Transverse current spectra in the regime of individual particle dynamics correspond to particles that can be treated as independent oscillators. As a result,  (11) and (12). Remarkable agreement is seen between MD and theoretical results for longitudinal and transverse velocity current spectra, justifying our assumption about the individual regime of particle motion in the high-q limit.
The comparison of Figs 1 and 6 indicates that an increase in q is accompanied by the transition from collective to individual dynamics whose fingerprints are clearly seen in the form of velocity current spectra. At high q values, the Lorentzian profiles (4) and (5) become inaccurate and profiles (11) and (12) should be used instead of them.

Conclusions
In conclusion, we have compared different methods for calculating excitation spectra in fluids, including the analysis of maxima of current spectra, its longitudinal and transverse components, and the joint mode analysis using the two-oscillator model to fit the total velocity current spectra. Performing MD simulations, we have studied the excitations at low-and high-temperatures in fluids of particles interacting with each other via the Lennard-Jones, IPL12, and IPL8 potentials, to understand trends arising due to potential softness and in the presence of attraction. Then, we have analyzed excitations at high temperatures corresponding to crossover from liquid-to gas-like dynamics in fluids (Frenkel line) and have studied the form of velocity current spectra in the regime of individual dynamics corresponding to excitations with high q-values.
We have found that the method based on the two-oscillator model provides the most accurate and comprehensive results for the frequencies and damping rates of both high-and low-frequency branches of the excitation spectra. Other considered approaches, including those based on the analysis of the transverse part of the total velocity current spectrum and on the analysis of maxima of the transverse current spectrum, can be used only for rough estimation of dispersion relations at low temperatures (not far from the melting line) and become unsuitable at higher temperatures.
The comparison of approaches based on the two-oscillator model and maxima of its transverse part has allowed us to distinguish three regions in the q-space corresponding to q-gap, as well as ranges of overdamped and non-overdamped low-frequency excitations. We have demonstrated that the results based on the analysis of the maxima of transverse current spectra can be obtained with the parameters of spectra calculated using the two-oscillator model. We have derived corresponding expansions near the transition between overdamped and non-overdamped regimes of low-frequency excitations.
We have analyzed the excitation spectra, specific heat, velocity autocorrelation functions, and temperature evolution of the q-gap in the 3D Lennard-Jones fluid at high temperatures and considered different features attributed to crossover from the liquid-to gas-like state of the fluid (Frenkel line). Considering this crossover, we observed that only non-overdamped low-frequency excitations rather than all of them leave the first pseudo-zone. However, the dynamic and thermodynamic features of the Frenkel crossover agree with each other in terms of low-frequency non-overdamped excitations instead of excitations in a broad sense (which include overdamped ones). In the case of 2D fluids, the dynamic and thermodynamic features strongly disagree with each other in general due to the increased role of anharmonicity.
We have found that the method based on the two-oscillator model is suitable for analysis of excitations in the regime of collective dynamics in 3D and 2D fluids with different interaction potentials. However, in the limit of fluctuations with high q values, we observed the return to the regime of individual particle dynamics, at which longitudinal and transverse current spectra are determined by Maxwellian distributions of particle velocities.
In addition to MD simulations, a powerful approach to experimental studies of fluid dynamics is provided by complex (dusty) plasmas, i.e., microparticles in a plasma discharge that acquire a negative charge and exhibit weakly damped individual particle dynamics, which can be imaged using video recording 80,81 . Particle-resolved experiments with complex (dusty) plasmas make it possible to analyze collective dynamics in strongly coupled many-body systems and, in particular, in a fluid state [82][83][84][85][86] . Although we considered only simple fluids in this work, we expect that the results can be applied to analyze excitations in fluid complex (dusty) plasmas.
The results of this work pave the way for detailed analysis of excitations in fluids of different natures, from simple fluids and noble gas fluids to liquid metals, molecular and complex fluids, complex (dusty) plasmas and www.nature.com/scientificreports www.nature.com/scientificreports/ related systems. We believe that our results will be of interest in context of various problems in condensed matter physics, chemical physics, physical chemistry, plasma physics, materials science, and soft matter.

Methods: Details of MD Simulations
We considered systems of particles interacting via the Lennard-Jones (LJ) and inverse power law (IPLk) pair potentials 1 : where r is the distance between a pair of particles; ε and σ are the energy and length scale of the interaction, respectively; and k is the exponent. We have studied 3D and 2D fluids consisting of = N 10 4 particles in an NVT ensemble with a Nose-Hoover thermostat. The cutoff radius of interaction was chosen as = . − r n 7 5 c D 1/ , where = n N V / is the numerical density and D is the space dimension. The numerical time step was chosen as σ ε ∆ = × − t T m T 5 10 / ( ) 3 0 2 and ε = . T / 001 0 . All simulations were performed in a cube (a square in 2D case) with a size = . L 21 5 ( = L 100 in 2D case) with periodic boundary conditions for × 2 10 6 time steps, where the first × 5 10 5 steps were used to equilibrate the system and the rest, to calculate the properties. Simulations were performed within the LAMMPS 87 and HOOMD-blue 88,89 packages.
In liquids, due to the isotropy, absence of translational order, and strong damping of excitations (especially with short wavelengths), a small part of the fluid inside the total volume (simulation box) can be considered as an independent subsystem. Due to this, the spectra can be calculated with fine q-resolution and in different directions (providing statistically the same results). We used this for spatial averaging of the spectra and for obtaining fine q-resolution at short wavelengths, while, in the long-wavelength region, the q-resolution was determined by 2π/L.