In quantum physics, an exponentially decaying state is characterized by a complex value of the frequency, with the imaginary part giving the lifetime of the particle and its decay probability (1). Well-known examples of unstable states with overwhelming imaginary part arise in the nuclear $$ \alpha $$-decay (Gamow states) and, in particle physics—the $$ W$$ and $$ Z^0$$ bosons (2)—where they are usually called “resonances.” Moreover, within hydrodynamics and gravitational theories, these excitations are labeled “quasinormal modes” (3), and they are the responsible for the black holes’ ringdown recently observed by LIGO Scientific Collaboration and Virgo Collaboration (4). In condensed matter physics, states with imaginary frequency arise frequently in disordered dissipative systems, in the form of purely relaxational overdamped modes, and even in heated crystals (5). They play a crucial role in liquids and glasses, where they are often called instantaneous normal modes (INMs) (67–8) and correspond to saddle points, with negative eigenvalues, in the energy landscape (910–11).

In all cases, defining the spectrum, or density of states (DOS), of systems that contain saddle points is challenging (12). The major obstacle resides in the fact that the Plemelj identity, that formally provides the connection between the propagator and the DOS, is defined on the real axis; hence only states with real frequencies $$ \omega $$, and positive eigenvalues $$ \lambda ={\omega }^2$$ of the Hessian, are allowed. As a consequence, so far, it has not been possible to analytically derive the DOS of systems with imaginary modes from a physical model, and the DOS of systems with INMs has been computed analytically only in one dimension where negative eigenvalues are absent (13).

In liquids [broadly defined to include plasma (14)], due to the presence of a large quantity of overdamped, exponentially decaying modes with purely imaginary frequency (the INMs), this problem has hampered the derivation of a universal law for the DOS of vibrational excitations, $$ g\left(\omega \right)$$. This is in stark contrast with the case of low-temperature crystals, where the vibrational modes are all real, and the Debye law $$ g\left(\omega \right)\approx {\omega }^2$$ has served as a fundamental law for the DOS of solids since 1912. In solid glasses, the Debye law is still valid at the lowest frequencies, although hybridized with $$ {\omega }^4$$ modes due to anharmonicity (15, 16). In liquids, we know from numerical simulations (7) that, as highlighted in ref. 17, $$ g\left(\omega \right)\approx \omega $$ at low frequency, but none of the theoretical models have been able to reproduce this law analytically. This is due to the impossibility of including modes with negative eigenvalues which dominate the low-energy part of the spectrum.

Here we provide a solution to this long-standing problem, by developing the fundamental law for the DOS of liquids, which takes the imaginary modes into account analytically. The key step in our derivation is the analytic continuation of the Plemelj identity to the complex plane. This leads to the possibility of defining a DOS for systems with imaginary frequency/energy. Our analytical model is successfully tested against numerical simulations from the literature on model (Lennard-Jones [LJ]) liquids.

The vibrational DOS of condensed matter systems is defined as

with the index $$ j$$ labeling different normal modes. Here $$ \mathcal{N}$$ denotes the total number of atoms in the medium. We consider decaying modes which are purely relaxational, as they arise, for example, from an overdamped Langevin dynamics in a liquid,where $$ \mathbf{v}$$ is the particle velocity, and $$ \tau $$ is the relaxation time. The Langevin differential operator, associated with Eq. 2, is given by $$ \mathfrak{L}\hspace{0.17em}≔\hspace{0.17em}d/dt+\mathrm{\Gamma }$$, and the corresponding Green’s function is obtained by solving $$ \mathfrak{L}\hspace{0.17em}G\left(t,x\right)=\delta \left(t,x\right)$$. After Fourier transforming the previous identity, we obtain the following form for the Green’s function:In order to relate the Green function for the single INM Eq. 3 with the total DOS Eq. 1, we need to use the generalized form of Plemelj identity (see refs. 18 and 19) valid for integration paths that do not necessarily lie on the real line,where $$ z$$ and $$ z^{\prime }$$ are arbitrary points in the complex plane region $$ {\mathcal{D}}^{+}$$, which coincides with the complex plane $$ \mathcal{C}$$ minus the wedge $$ -\pi /4<\arg \left(z\right)<5\pi /4$$, and $$ P$$ indicates the Cauchy principal value. This limitation arises upon considering the derivation of the Plemelj identity via distribution theory. The starting point of the derivation is usually given by an integral,where, clearly, as long as $$ z$$ is real, no fundamental problem appears, and the integral converges everywhere. Typically (20), one then proceeds by showing that this integral is equal to what one obtains upon applying the Fourier transform to the Heaviside function in $$ \mathcal{S}$$ (the space of Schwartz distributions) (21), which then leads to Eq. 4 above, with $$ z,z^{\prime }$$ as all real variables.

Now, the analytic continuation of this integral, that is, promoting $$ z$$ to be a complex variable, leads to a divergence for $$ \text{Im}\hspace{0.17em}z<0$$. Following Zeldovich (1), the Gaussian regularization

can be used (upon subsequently taking the limit $$ \lambda \to 0^{+}$$). Following ref. 18, one can show that Eq. 5 converges everywhere in the complex plane apart from the wedge defined by $$ -\pi /4<\arg \left(z\right)<5\pi /4$$. Hence, one retrieves the generalized Plemelj identity Eq. 4 valid for arbitrary pathways in the allowed complex plane region $$ {\mathcal{D}}^{+}$$.

We therefore write

thus identifying $$ z\equiv i\hspace{0.17em}\omega +\xi $$, and $$ z^{\prime }\equiv -\mathrm{\Gamma }+\xi $$, where $$ \xi $$ is real. Next, we apply the generalized Plemelj identity (Eq. 4) to $$ G\left(z\right)$$ to evaluate the DOS,and, upon taking the imaginary part of the left-hand side, we thus obtainThe above Eq. 9 is our main result and provides a simple yet universal law for the DOS of liquids. Eq. 9 is the density of states for a single INM with Green function Eq. 3. Given the linearity of the problem, we can generalize the result to a set of $$ j$$ INMs with different relaxation times aswhere $$ {\mathrm{\Gamma }}_j$$ is the relaxation rate of the $$ j$$th INM.

Expanding Eq. 10 in the limit of low frequency, $$ \omega \ll {\mathrm{\Gamma }}_j$$, we obtain

which recovers a common trend observed in many molecular simulation studies of liquids, where $$ \alpha $$ is treated as a fitting parameter (17, 22).

In Eq. 11, the various relaxation rates sum up in parallel,

implying that, in the presence of a separation of scales $$ {\mathrm{\Gamma }}^{*}\ll {\mathrm{\Gamma }}_2,{\mathrm{\Gamma }}_3,\dots \hspace{0.17em}$$, the average relaxation rate $$ \mathrm{\Gamma }$$ would be given by the smallest of them, that is, $$ {\mathrm{\Gamma }}^{*}$$. This is tantamount to saying that the low-frequency dynamics is governed by the longest living imaginary mode—a well-known result in the realm of hydrodynamics and effective field theory around equilibrium. In the rest of the paper, whenever we will discuss the relaxation rate, we will mean the total one defined in Eq. 12.

We now present predictions of the main result of our paper, Eq. 11, in comparison with numerical simulations data of simple liquids, that is, the LJ system. It is worth recalling that $$ g\left(\omega \right)$$ is obtained numerically (from diagonalization of the Hessian matrix of instantaneous snapshots of particle positions) by retaining also the imaginary frequencies, because $$ g\left(\omega \right)\equiv g\left(\left|\omega \right|\right)$$ (6). Hence, in the following, $$ \omega $$ stands for the absolute value of the excitation frequency.

In Fig. 1, we show the comparison between the model predictions and original MD simulations data of LJ liquids from ref. 17. The model Eq. 11 uses an effective relaxation time $$ \mathrm{\Gamma }$$ as a fitting parameter, besides the normalization prefactor.

The simulation data can be nicely fitted using Eq. 11 with just two parameters up to the maximum of the DOS, after which a Gaussian cutoff $$ \sim e^{-{\left(\omega /{\omega }_D\right)}^2}$$ is needed to capture the fall-off [note that, in low-$$ T$$ solid glasses, a simple exponential cutoff is instead used (24)]. In the dataset at the highest density, a small peak becomes visible, which is a relic of a pseudo-van Hove singularity upon approaching the solid state.

In Fig. 1, we compare predictions of Eq. 11 with the more recent simulations data from ref. 23, also for the LJ liquid. Also in this case, excellent agreement between Eq. 11 and simulations is found. As a further confirmation of the validity of our Eq. 11, we have also checked that the mean relaxation time $$ \mathrm{\Gamma }$$ follows the Arrhenius law as a function of $$ T$$, as expected for equilibrium liquids (17), as shown in Fig. 1 A, Inset. We found that the activation energy for these relaxations is $$ \sim 0.15ϵ$$, where $$ ϵ$$ is the depth of the LJ attractive well, whereas the attempt frequency is $$ \sim 10ϵ/\hslash $$, consistently about 10 times the escape rate from the LJ well. Since $$ \mathrm{\Gamma }$$ corresponds to the frequency of the maximum of the DOS (obtainable by setting the derivative of the summand in Eq. 10 equal to zero), this estimate gives an upper bound for the thermally activated relaxation processes from anharmonic saddle points in the energy landscape.

This comparison shows that, especially in the low-frequency region of the DOS up to the maximum, the spectrum is dominated by the relaxational modes or INMs. At higher frequencies, also, phonons are present in the DOS (25), which are not described by our minimal model but can be included in future work. At present, the mean relaxation time $$ \mathrm{\Gamma }$$ acts as an effective parameter, which effectively takes into account, also, other vibrational excitations that are not explicitly implemented in the model and may become important around the frequency of the maximum and above.

In sum, the proposed theoretical model provides a universal law $$ g\left(\omega \right)\sim \omega $$ for the low-frequency DOS of liquids, in agreement with observations from many simulations studies in the literature (6, 7, 17, 23). This law plays, for liquids, the same pivotal role that the Debye law $$ g\left(\omega \right)\sim {\omega }^2$$ plays for solids, and it explains, for example, that the liquids are mechanically unstable [by leading to diverging negative nonaffine contributions to the shear modulus (26)] and that INMs may play a role, also, for the thermal properties of solids. Furthermore, it could explain the observation of DOS scalings $$ \sim {\omega }^{\alpha }$$ with $$ \alpha \in \left[\mathrm{1,2}\right]$$ in glasses (27) and other complex systems. This law has not been derived before, due to the difficulty of dealing with the imaginary frequencies [associated with saddle points in the energy landscape (9, 10)] that dominate the low-frequency spectrum of liquids. In this work, we have solved this problem by analytically continuing the Plemelj identity to the complex plane using input from recent developments (18, 19). This methodology is general and extends far beyond the case of liquids or glasses. Further applications of our result will lead to the possibility of analytically describing the DOS of unstable states in quantum mechanics and the spectrum of black holes where imaginary modes and quasinormal modes play an important role.