Rapid onset of molecular friction in liquids bridging between the atomistic and hydrodynamic pictures

Friction in liquids arises from conservative forces between molecules and atoms. Although the hydrodynamics at the nanoscale is subject of intense research and despite the enormous interest in the non-Markovian dynamics of single molecules and solutes, the onset of friction from the atomistic scale so far could not be demonstrated. Here, we fill this gap based on frequency-resolved friction data from high-precision simulations of three prototypical liquids, including water. Combining with theory, we show that friction in liquids emerges abruptly at a characteristic frequency, beyond which viscous liquids appear as non-dissipative, elastic solids. Concomitantly, the molecules experience Brownian forces that display persistent correlations. A critical test of the generalised Stokes–Einstein relation, mapping the friction of single molecules to the visco-elastic response of the macroscopic sample, disproves the relation for Newtonian fluids, but substantiates it exemplarily for water and a moderately supercooled liquid. The employed approach is suitable to yield insights into vitrification mechanisms and the intriguing mechanical properties of soft materials. How friction in liquids emerges from conservative forces between atoms is currently not well-understood, but it is a crucial parameter for dynamic processes in liquid matter. Here, the authors combine frequency-resolved simulation data with theory to show that the friction felt by a single molecule occurs abruptly below a certain frequency.

M olecular friction is a key ingredient for dynamic processes in fluids: it limits diffusion, governs dissipation, and enables the relaxation towards equilibrium. In a liquid environment, the friction experienced by solvated molecules and nanoprobes exhibits a delayed response to external stimuli, indicating non-Markovian dynamics [1][2][3][4] . Such memory is found on sub-picosecond up to microsecond scales; it has repercussions on macromolecular transition rates [5][6][7] and is manifest in the visco-elastic behaviour of soft materials [8][9][10][11] . However, the origin of friction from conservative forces between molecules and atoms remains as one of the grand challenges of the physics of fluids [12][13][14] . Stokes's friction law, describing the resistance to a steadily dragged immersed sphere of radius a, links the friction ζ 0 to the (macroscopic) shear viscosity η 0 , and the relation ζ 0 = 6πη 0 a scales down remarkably well to single molecules 15,16 . Stokes's hydrodynamic treatment 17 from 1851 was actually more general and addressed slow oscillatory motions in viscous fluids, motivated by inaccuracies of pendulum clocks caused by air flow (Fig. 1a). These predictions of a dynamic friction ζ(ω) that depends on the oscillation frequency ω have been interpreted in terms of a delayed response, also referred to as hydrodynamic memory, and have only recently been shown to be quantitative also for micron-sized particles 1,2 . In the domain of microrheology, measurements of ζ(ω) are used to infer the mechanical properties of complex fluids [18][19][20][21][22] .
From the perspective of individual molecules or atoms, fluids are governed by conservative interactions and obey Newton's equations of motion, yielding smooth and time-reversible trajectories (Fig. 1b, c). In particular, a single molecule is not subject to friction in this picture, and the mechanism of the required entropy production is far from obvious. Macroscopic friction and other transport coefficients have been linked to microscopic chaos and the Lyapunov spectrum of the liquid [23][24][25] , yet the connection of the latter to the corresponding Green-Kubo integrands, or equivalently, to the dynamic friction ζ(ω), is an open issue 24 . First-principle theories to friction are hampered by the fact that liquids are strongly interacting systems. An insightful, formal relation between the many-particle Liouville operator and dissipation spectra was derived in the seminal works by Zwanzig, Mori, and others 26,27 , but the analytic evaluation of the resulting expressions hinges on uncontrolled approximations. For example, a systematic short-time expansion of the motion at all orders would yield zero friction (see "Methods"). Early work on dissipation spectra recognised the importance of exact sum rules 28,29 , the proposed ad hoc models, however, violate the sum rules at higher orders. Theoretical progress was made for hardsphere fluids, where billiard-like collisions generate an instantaneous, Markovian contribution to friction 30 , thereby rendering the frictionless regime inaccessible.
To gradually bridge between the atomistic and hydrodynamic regimes, one would ideally like to have a magnifying glass that allows for zooming from the slowest to the fastest processes, thus obtaining an increasingly sharper view of the molecular details (Fig. 1a). Spectral quantities such as ζ(ω) can serve this purpose with the frequency ω as the control knob. Implementing this idea in simulations of liquids and numerically tracing the friction of molecules over wide frequency windows from fully developed dissipation all the way down to the non-dissipative regime would reveal sizable variations of ζ(ω). Such deviations of ζ(ω) from a constant friction ζ 0 signify non-Markovian motion that is widely cast in the generalised Langevin equation ["Methods", Eq. (1)], parametrised by an associated memory function γ(t). The quantities ζ(ω) and γ(t) are related to each other by a cosine transform, but the determination of either of them from data is a formidable challenge, with substantial progress in the past years [31][32][33][34][35][36][37][38] . Reaching the high-frequency regime was precluded so far by practical limitations, e.g., statistical noise and insufficient dynamic windows.
Here, we have overcome these limitations by high-precision simulations and an advanced data analysis that utilises physical principles. For three distinct liquids-water, a dense Lennard-Jones (LJ) fluid representing a simple, mono-atomic liquid, and a supercooled binary mixture serving as a model glass former-we obtained low-noise, cross-validated data for the dynamic friction ζ(ω) and the frequency-dependent viscosityηðωÞ over windows that span from the fastest to the slowest processes in each liquid. Corroborated by these data, we give rigorous arguments that dissipative processes in molecular liquids are exponentially fast suppressed as frequency is increased. As a consequence, the liquid's response is purely elastic above a characteristic frequency ω c , a feature that goes beyond popular models of friction and viscoelasticity. Furthermore, the magnitude of friction cannot be inferred from the instantaneous behaviour of molecular trajectories: its origin is non-local in time. Finally, having data available for both ζ(ω) andηðωÞ, we tested the connection between microscopic friction and the macroscopic mechanical properties of complex fluids, postulated by the generalised Stokes-Einstein Fig. 1 Dynamic friction bridges between the hydrodynamic and atomistic pictures of liquids. a A pendulum that oscillates in a viscous fluid with frequency ω experiences a dynamic friction ζ(ω) as calculated by Stokes (1851) for slow motion. The associated flow pattern (stream lines) leads to hydrodynamic memory of the motion. Magnifying the microscopic scale, the fluid consists of molecules that obey Newton's equations, which are nondissipative and generate smooth trajectories. Stokes's result for the pendulum scales down to single molecules, and the function ζ(ω) provides the bridge between the frictionless (microscopic) and the hydrodynamic (macroscopic) descriptions. b, c In liquids, the short-time trajectories of molecules (b) are smooth, but chaotic curves due to molecules rattling in transient cages formed by their neighbours (c brown spheres). The onset of friction is driven by the momentum transfer to the cage, as supported by control simulations of a single particle (red sphere) in a matrix of pinned particles (brown and yellow). Illustration for a mono-atomic fluid in two dimensions.
relation. The latter is found to either fail completely (monoatomic liquid), serve as a qualitative description (water), or being a nearly quantitative relation (supercooled liquid).

Results
Molecular friction in liquids emerges rapidly. Instead of observing the response to an oscillatory force, we recorded the Brownian position fluctuations in equilibrium and link them to the friction, taking advantage of a fluctuation-dissipation relation. For the three liquids under investigation, we carried out molecular dynamics simulations to compute the mean-square displacement (MSD). Using the MSD as sole input, we estimated both the dynamic friction ζ(ω) and the memory function γ(t) in an ansatz-free approach, following two complementary routes that allow for cross-validation (see Fig. 2 and "Methods"). The first route invokes complex analysis and is based on the Fourier-Laplace transform of correlation functions, sampled on a sparse time grid ("adapted Filon algorithm"). Second and independently, we computed the antiderivative of γ(t) using a stable deconvolution technique for uniform time grids.
Although all three liquids display rather different dynamics, leaving distinct fingerprints in their friction spectra, their highfrequency behaviours of ζ(ω) share significant similarities (Fig. 3a-c). Most remarkably, the data demonstrate that beyond a liquid-characteristic frequency, ω ≳ ω c , the friction ζ(ω) goes exponentially fast to zero. Such a rapid spectral variation has to be contrasted to the typical algebraic peaks, i.e., the Lorentz-Debye shape, and we argue in the following that our finding is generic for molecular fluids. Upon decreasing frequency, the onset of friction is followed by liquid-specific behaviour over several decades in time until the hydrodynamic value ζ 0 is established: in water, our results for the friction ζ(ω) exhibit a local maximum at ω/2π ≈ 5 THz, followed by a slow increase towards the limiting value ζ 0 , which is reached near frequencies of 0.1 THz. For the LJ fluid, ζ(ω) varies more smoothly with a global maximum at an intermediate frequency, and ζ 0 is approached slowly from above in accord with the hydrodynamic square-root singularity [Eq. (18)], essentially calculated by Stokes already 17 . On theoretical grounds, this feature of ζ(ω) is generic for all liquids, yet it is suppressed in our data for the other two liquids due to a small prefactor. In the supercooled liquid, we observe a scale separation by three dynamic decades of (i) the rapid onset of friction and (ii) the slow emergence of the hydrodynamic limit. The second process is associated with cage relaxation, strongly delayed in the glassy state, which suggests that the driving mechanism behind  22)); the dynamic friction ζ(ω) follows via Eq. (7). A Fourier backtransform then yields the memory function γ(t). The latter can be obtained more directly along the lower route, which is based on the velocity autocorrelation function (VACF) Z(t) and employs a deconvolution in time domain (Eq. (23)).  Fig. 2, with the mean-square displacement (MSD, not shown) as initial input. a-c The dynamic friction (red) is the real part of the memory kernelγðωÞ, as computed from Eq. (6). It interpolates between the hydrodynamic value ζ 0 /m (horizontal lines) and a rapid decrease to zero at high frequencies (ω ≫ ω c ); thin red lines mark exponential decays, $ e Àω=ω c , to guide the eye. In case of the Lennard-Jones fluid (a), ζ(ω) is consistent with Stokes's small-ω asymptote [grey line, Eq. (18)], with parameters taken from a fit to the long-time tail of Z(t) (inset of Fig. 4a). The elastic response ImγðωÞ (turquoise) exhibits a local maximum at high frequencies, defining the characteristic frequency ω c (vertical lines), and follows the high-frequency asymptote [Eq. (15), grey line], with all parameters fixed by short-time fits to the MSD data. The frequency ω c is close to, but different from the Einstein frequency ω 0 . d-f Numerical results for the generalised mobilityŶðωÞ, which describes the response to an applied small, oscillatory force (see "Methods"). The real part ReŶðωÞ (green) approaches the hydrodynamic mobility 1/ζ 0 (horizontal lines) as ω → 0 and vanishes rapidly at large frequencies. The imaginary part (yellow) encodes the elastic response, which has a resonance near ω c (vertical lines); at larger frequencies, the data match with theoretical predictions [grey lines, Eq. (14)].
the onset of friction is not in the structural relaxation of the fluid. Close to the glass transition, the small-frequency friction ζ 0 is governed by self-similar relaxation processes and obeys asymptotic scaling laws 39,40 ; the magnitude (prefactor) of these laws, however, is set at high frequencies.
Friction depends on the coupling of fast and slow processes. The obtained data of ζ(ω) cover the full range of the dynamic response, thereby connecting physics at different scales. Before discussing the remaining panels (d-f) of Fig. 3, we rationalise key features of ζ(ω) by tracing their origins to the dynamics of the fluid particles at short and long times (going backwards in Fig. 2). The relevant properties are prominently visible in the second derivative of the MSD, the velocity autocorrelation function (VACF), ZðtÞ :¼ ∂ 2 t MSDðtÞ=6, after numerical differentiation with respect to the lag time t (Fig. 4a-c). The following should be contrasted to Ornstein's model for Brownian motion, employing a single exponential decay of velocity correlations, ZðtÞ % v 2 th expðÀζ 0 t=mÞ, which implies a constant (Markovian) friction, ζ(ω) ≈ ζ 0 ; by v th we denote the thermal velocity, and m is the molecular mass. As a distinct feature of molecular fluids, obeying Newton's equations, the VACF's true short-time decay is parabola-shaped 16 À Á , introducing the Einstein frequency ω 0 . For dense liquids, the VACF, after a sign change, generically displays a regime of anti-correlations, which reflect the transient caging by neighbouring molecules (Fig. 1c). For water and the supercooled liquid, these anti-correlations relax slowly with an intermediate power-law decay, Z(t)~−t −5/2 (insets of Fig. 4b, c); such a decay was observed earlier in supercooled liquids 41,42 and it is a well-established long-time tail in colloidal suspensions [43][44][45] and for diffusion in an arrested, disordered environment 46,47 .
The famous long-time tail encoding hydrodynamic memory 16,48,49 , Z(t → ∞)~t −3/2 , is clearly developed in our data for the LJ fluid after another sign change (inset of Fig. 4a), and in this situation, Stokes's hydrodynamics describes the slow motion of single molecules (Fig. 3a). For the other two liquids studied, the tail is not visible in our data due to a small prefactor, which following mode-coupling arguments decreases as either viscosity or diffusivity increases 16 .
The dynamic friction is closely linked to the complex-valued, generalised mobilityŶðωÞ via ζðωÞ ¼ Re½ŶðωÞ À1 . This mobility encodes the response to a small, oscillatory force and is accessible in, e.g., scattering experiments 50 . Here, we computedŶðωÞ from the VACF upon employing a fluctuation-dissipation relation (see "Methods" and Fig. 3d-f). Our data reveal a generic, rapid decrease of the dissipative part,Ŷ 0 ðωÞ :¼ ReŶðωÞ, upon increasing frequency towards the microscopic regime, ω ≫ ω c ; concomitantly, the elastic response,Ŷ 00 ðωÞ :¼ ImŶðωÞ, has a resonance near ω c due to vibrational motion of molecules in their cages. In the low-frequency limit, the reciprocal of the macroscopic friction is recovered,Ŷðω ! 0Þ ¼ 1=ζ 0 . In the examples studied, both regimes are separated by at least two decades in time, which show material-specific features: the mobility of water molecules is 2.5-fold enhanced over its macroscopic value near ω/2π ≈ 1.3 THz; similarly, a factor of 30 is observed for the supercooled liquid. At variance, the hydrodynamic long-time tail, as for the LJ liquid, demandŝ Y 0 ðω ! 0Þ to be approached from below. The interplay of slow and fast processes enters the responseŶðωÞ, and thus the dynamic friction ζ(ω), at all frequencies, which is borne out by the Kramers-Kronig relations. In particular, non-analytic behaviour ofŶ 00 ðωÞ as ω → 0 influences the detailed onset of friction at large frequencies.
The origin of friction is non-local in time. The rapid decrease of Y 0 ðωÞ is mathematically justified from the short-time properties of the VACF. Physical molecular trajectories, being solutions to Newton's equations, are smooth and yield a smooth function Z(t); in particular, all derivatives Z (n) (t) exist at t = 0 and are finite. Thus, invoking exact sum rules [Eq. (10)], all moments of the spectrumŶ 0 ðωÞ are finite, which requires an exponentially fast decay as ω → ∞. (This is a special case of a more general characterisation of exponentially decaying probability measures 51 .) Combining with the large-ω asymptote of the imaginary part, Y 00 ðωÞ ' 1=mω )Ŷ 0 ðωÞ (see "Methods") and using ζðωÞ 1 Y 0 ðωÞ=Ŷ 0 ðωÞ 2 þŶ 00 ðωÞ 2 Â Ã proves that ζðωÞ ' ðmωÞ 2Ŷ0 ðωÞ as ω → ∞ and thus an exponentially fast suppression of the friction. We stress further that such behaviour is not contained in representations of ζ(ω) as a truncated continued fraction 16,52 .
It is tempting to use a systematic short-time expansion of Z(t) to predict the large-frequency behaviour of the friction. However, Z(t) being an even function due to time-reversal symmetry in equilibrium renders the large-ω asymptotes zero,Ŷ 0 ðωÞ 0 and thus ζ(ω) ≡ 0, even if the complete Taylor series of Z(t) in t = 0 was known (see "Methods"). Note thatŶ 00 ðωÞ and the elastic counterpart of ζ(ω) are well captured by such an expansion ( Fig. 3d-f). This observation underlines that, on all scales, friction emerges as a phenomenon that is non-local in time, i.e., it cannot be anticipated from the local behaviour of the molecular trajectories.
Our numerical and theoretical findings are supported by an analytically solvable example. The choice combines the smoothness and time-reversal symmetry of molecular autocorrelation functions with a long-time tail. It yields an exponential decay of the mobility,Ŷ 0 ðωÞ ¼ ðπτ=2mÞ e Àjωτj (dissipative part), and hence a rapid suppression of friction, ζðω ! 1Þ$ðωτÞ 2 e Àωτ , as demanded above (see "Methods" and Fig. 5).
Irreversible momentum transfer drives the onset of friction. A pressing question is about the physical mechanism that generates the onset of friction. Motivated by our results for the supercooled liquid (Fig. 3c), we performed a control simulation for the LJ fluid with the structural relaxation switched off by pinning all particles but one. The rattling motion in such a frozen-in cage (Fig. 2) experiences a dynamic friction that closely resembles our generic findings for ζ(ω) at high frequencies, ω ≳ ω c (Fig. 6). Upon decreasing frequency from ω c to zero, the two dynamics deviate strongly as is most evident in the elastic response: whereas ImγðωÞ decreases for the unconstrained fluid and exhibits a zero crossing enforced by hydrodynamics [Eq. (17) and Fig. 3a], it remains positive for the pinned case and grows as Imγðω ! 0Þ$1=ω, reminiscent of what one obtains for an Ornstein-Uhlenbeck (OU) particle in a harmonic trap 1 . For even smaller frequencies, ω ≲ τ À1 LJ , the dissipation diverges too, approximately as ζ(ω → 0)~ω 1/2 , which we attribute to the irregular shape of the confining cages; for the OU model with harmonic confinement, ζ(ω) ≃ const. At high frequencies, however, the confinement is not relevant, and we conclude that it is the fast, yet irreversible momentum transfer to neighbouring molecules that drives the onset of friction. This is corroborated by the observation that instantaneous momentum exchange implies a non-zero limit, ζ(ω → ∞) > 0, as in the case of hard spheres 30 . The fact that dissipation is linked to trajectories for which the time-reversed path is extremely improbable leads us to speculate that the onset frequency ω c is intimately related to the largest Lyapunov exponent λ max of the fluid, which is close to, but different from the Einstein frequency ω 0 (Posch & Hoover 25 ).
Dynamic friction implies intricate memory of Brownian motion. Within the framework of the generalised Langevin equation [Eq. (1)], momentum relaxation is governed by a memory function γ(t) that is fully determined by ζ(ω) [Eqs. (7) and (9)]. At the same time, γ(t) is also the autocorrelator of Brownian, random forces on the molecules, up to a constant  (20), combining the smoothness and time-reversal symmetry of molecular autocorrelation functions with a long-time tail. As a sensitive test of our numerical procedures, symbols show numerical results from the mean-square displacement as input, sparsely sampled on a geometrically spaced grid. The panels show a the complex-valued, generalised mobilityŶðωÞ, b the dynamic friction ζ(ω) and its elastic counterpart ImγðωÞ, and c the memory function γ(t) in time domain. The latter inherits the long-time tail $t À2 from the velocity autocorrelation function, but of opposite sign (inset). The tail induces nonanalytic behaviour of ζ(ω) at small frequencies, which crosses over to $ω 2 due to the smooth, exponential cutoff of the Fourier integrals (inset of panel b). Same colour code as in Figs. 3 and 4. Fig. 6 Friction emerges due to rattling motion in immobile cages. Dynamic friction ζ(ω) (red symbols) obtained in a control simulation of a single particle moving in a frozen-in cage formed by neighbouring particles (Fig. 1c). The setup was created by pinning the particles of the Lennard-Jones fluid except for one; results correspond to an ensemble average over 10 6 typical cages. The imaginary part of the memory kernelγðωÞ (turquoise) ties in with the high-frequency prediction (grey). Dashed lines show the results for the unconstrained fluid at the same conditions for comparison (Fig. 3a).
COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-020-0389-0 ARTICLE COMMUNICATIONS PHYSICS | (2020) 3:126 | https://doi.org/10.1038/s42005-020-0389-0 | www.nature.com/commsphys prefactor, and ζ(ω) encodes the corresponding spectrum [Eq. (2)]. Within Ornstein's idealised model of Brownian motion, one assumes independent Brownian forces, implying a flat, white spectrum, ζ(ω) ≈ ζ 0 , and a delta-peaked memory function γ(t). For molecular liquids, however, the memory functions display a universal parabola-shaped short-time decay [ Fig. 4d-f and Eq. (16)]. For water and the LJ fluid, the short-time regime of γ(t) is followed by oscillatory behaviour including sign changes and, finally, different power-law decays for the two liquids encoding different physics (insets of Fig. 4d, e also see Fig. 5c); generally, power-law tails of the memory function are directly inherited from the VACF without a change of exponent [Eq. (19)]. For the supercooled liquid, γ(t) remains positive and exhibits the onset of a plateau (near t ≈ 0.3τ LJ ), which decays logarithmically slowly over 2 decades in time (Fig. 4f). From the modelling perspective, it is desirable to approximate the memory functions such that the initial decay, the typical persistence time, and the integral of γ(t) are reproduced, the latter yielding ζ 0 [Eq. (13)]. For all three liquids, the complexity and the long-lived nature of the memory, however, preclude simple models of γ(t) such as the superposition of few exponential decays. The quantitative knowledge of ζ(ω), as obtained here, paves the way for more favourable approximations of memory in the frequency domain, which can yield mathematically and physically consistent interpolations of Brownian motion from the fastest to the slowest scales.

Discussion
Based on theoretical arguments and corroborated by simulation data that cover up to three orders in magnitude and 4 decades in frequency, we have shown that the friction ζ(ω) in molecular liquids emerges rapidly at a large, liquid-specific frequency ω c from the time-reversible motion of the atoms. This onset of friction is driven by the irreversible momentum transfer to neighbouring molecules, but also influenced by the slowest processes at hydrodynamic scales. The exponentially fast decay of dissipation spectra such as ζ(ω) is easily shadowed by approximations, and as a modelling constraint it is under-investigated in the existing microscopic theories of liquids. The observation that the high-frequency behaviour of ζ(ω) is reproduced by the motion of a single particle in an immobile cage refines the question on the quantitative link between friction and microscopic chaos. Concretely: whether and how does the frequency dependence of transport coefficients relate to the Lyapunov spectrum of the liquid 23,24 ?
Going beyond the dynamics of single molecules and their friction, the relation to the visco-elastic properties of complex fluids is of ongoing interest for the physics of polymers, living cells, and the glass transition. The potentially tight coupling between single-particle and collective responses was phrased as an ad hoc extension of Stokes's friction law to the frequency domain, ζðωÞ ¼ 6πRe½ηðωÞa, referred to as generalised Stokes-Einstein relation (GSER), which has found wide applications in the context of microrheology experiments [18][19][20][21] . It links the dynamic friction ζ(ω) of a probe particle to the dynamic shear modulus, GðωÞ ¼ ÀiωηðωÞ, a complex-valued function encoding the stress response of the macroscopic fluid sample to a small, oscillatory shear strain. The generalised viscosityηðωÞ tends to the hydrodynamic shear viscosity η 0 as ω → 0, with ReηðωÞ representing the spectrum of shear stress fluctuations (up to a constant factor) by a fluctuation-dissipation relation. Thus, if the GSER holds the single-particle memory γ(t) is proportional to the autocorrelator of shear stresses, which means that the Brownian force on the particle and the fluctuation of the shear stress are statistically equivalent variables.
A critical assessment of the validity of the GSER is permitted by comparing our data for ζ(ω) to results forηðωÞ, calculated within the same simulations (see Fig. 7 and "Methods"). Generically, the dissipative part, ReηðωÞ, decays exponentially fast as ω → ∞, which is required by analogous arguments as given for ζðωÞ andŶ 0 ðωÞ. For high frequencies, only the imaginary part  (Fig. 3a,  c), and correspondingly for the imaginary counterparts of the elastic responses (orange symbols and turquoise dashed lines). The effective particle radius a for each liquid is chosen such that the viscosity and friction curves coincide at ω = 0. The pink solid line in a is an empirical fit of a compressed exponential, ReηðωÞ ' η 0 expðÀðω=ω η Þ β Þ with β = 1.29 and ω η = 1.19 ω 0 . In b, c the data for the elastic responses are shifted by the indicated factors for clarity. d-f The GSER is tested by plotting the ratios ζðωÞ=6π½ReηðωÞa (violet) and m½ImγðωÞ=6π½ImηðωÞa (orange); deviations from unity quantify the GSER violation.
survives due to its slow decay,ηðωÞ ' G 1 =ðÀiωÞ, inducing a non-zero and real-valued modulus,Ĝðω ) ω c Þ % G 1 > 0. Therefore, our data clearly demonstrate that, indeed, liquids respond to high-frequency shear like a non-dissipative, elastic solid as put forward by the classical work of Frenkel 53 . However, the attempt to predict the elastic modulus G ∞ from the vibrational motion of molecules in their cages, by virtue of the GSER, would considerably overestimate the modulus by factors of ≈2 for the three liquids studied (Fig. 7d-f). In passing, we note that Maxwell's model for viscoelasticity 16,52,53 ,ηðωÞ ¼ G 1 τ=ð1 À iωτ M Þ with some relaxation time τ M , breaks down at high frequencies as it implies a slowly decaying real part, Reηðω ! 1Þ $ ω À2 , in sharp contrast to exact sum rules 28 and to the exponentially fast decay for molecular liquids. Therefore, treatments of sound-like, elastic waves on the footing of this and similar models appear incomplete.
The passage from the elastic to the viscous limit occurs upon decreasing frequency, leading in case of the LJ fluid to a monotonic increase of ReηðωÞ, which, empirically, follows a compressed exponential over almost the full frequency domain (Fig. 7a). In particular,ηðωÞ % η 0 is constant for ω≲2τ À1 LJ , which defines the hydrodynamic regime. For these frequencies, the single-molecule response is very well described by Stokes's dynamic friction [Eq. (18)], making the GSER violation apparent for Newtonian fluids, for whichηðωÞ ¼ η 0 . It is evidenced further by the dissimilarity of the elastic parts, ImγðωÞ and ImηðωÞ.
For water, the viscosity and friction spectra share similar features and coincide fairly well (Fig. 7b), including the elastic parts. Thus, the GSER serves as a reasonable qualitative description, in particular for frequencies below ≈2 THz, i.e., slower than the vibrations of the first hydration shell 54 . In supercooled liquids, the Stokes-Einstein relation for molecules (i.e., the GSER for ω → 0) holds in the presence of a huge variation of the viscosity. In particular, the ratio ζ 0 /η 0 is observed to be constant over a wide temperature range-in line with the mode-coupling theory of the idealised glass transition 40 . At very low temperatures, however, marked deviations from the Stokes-Einstein relation (mostly studied for ω → 0) have received considerable attention as they signify the opening of additional relaxation channels not included in the standard theory [55][56][57][58][59] . For the moderately supercooled liquid studied here exemplarily, both viscous and elastic responses satisfy the GSER over a wide frequency window (Fig. 7c). Notably, the elastic components collapse almost perfectly in this case, ImγðωÞ $ ImηðωÞ, which we attribute to the same (apparent) power-law scaling, ≈ω −0.75 , at intermediate frequencies. Yet, the collective response lacks the elastic peak of ImγðωÞ at ω c , causing the breakdown of the GSER at large frequencies. This suggests that a future frequency-resolved study of systematic deviations from the GSER upon further supercooling can clarify the separate contributions of fast and slow processes to the decoupling of diffusion and viscosity ("Stokes-Einstein violation") close to the glass transition temperature.
We conclude that molecular friction in liquids arises from a complex interplay of processes on very disparate time scales. The large variability of ζ(ω) observed over orders of magnitude reveals the strongly non-Markovian nature of Brownian motion in any liquid environment, with far-reaching implications for nanoscale processes. Examples are as diverse as reaction rates and barrier crossings in macromolecular dynamics [5][6][7] and flows near liquid-solid interfaces 44,[60][61][62] ; the ability to quantify the corresponding memory is vital for their realistic descriptions. From a numerical point of view, our ansatz-free approach has immediate applications to and extends current methods 37,38 for the analysis of high-resolution microrheology data [1][2][3]21 , which involves deducing frequency-dependent moduli from the displacement statistics along the same lines as done here for the dynamic friction. Relying merely on the existence of a steady state [Eq. (12)], the developed methodology is not limited to friction, but can be transferred to the analysis of non-Markovian time series from simulation and experiment. It finds novel uses in, e.g., the anomalous diffusion within living cells 11 and the kinetics of chemical reactions 63 . It opens a promising avenue for research on the migration of malignant cells in tissue 64 and on predictive stochastic models of financial market 65 and geographic 66 data.

Methods
Generalised Langevin equation. A labelled fluid particle of mass m, position r(t), and momentum pðtÞ ¼ m_ rðtÞ obeys the generalised Langevin equation (GLE) 67 : where the Brownian force ξ(t) is a stochastic process with zero mean and covariance to satisfy the fluctuation-dissipation theorem. Rewriting Eq. (1) for the the velocity autocorrelation function (VACF), Z(t) = 〈p(t) ⋅ p(0)〉/(3m 2 ), describing momentum relaxation, yields Its Fourier-Laplace transform [Eq. (8)] provides the link to and the definition of the (complex-valued) memory kernelγðωÞ, Linear response. For a mass m driven by a periodic force FðtÞ ¼ F ω cosðωtÞ with frequency ω and amplitude F ω , the stationary response vðtÞ, averaged over many realisations of the experiment, is given by 67 corresponding to Eq. (1) after shifting the lower integral boundary to −∞ to ensure relaxation of transients. The upper boundary can be shifted to +∞ with the convention that the response function γ(t < 0) = 0. By linearity of the equation, the solution is of the form vðtÞ ¼ Re v ω e Àiωt ½ with complex amplitude v ω , and inserting into Eq. (5) yields v ω ¼ŶðωÞ F ω in terms of the generalised mobility, also referred to as complex-valued admittance. Its central ingredient is the onesided Fourier transform of the response function,γðωÞ :¼ R 1 0 e iωt γðtÞ dt. Comparing to Eq. (4), which describes equilibrium correlations, yields the fluctuation-dissipation relation:ẐðωÞ ¼ k B TŶðωÞ.
Friction describes the resistance to a prescribed velocity, as in Stokes's pendulum experiments 17 . Thus, inverting the above argument, the force response to an oscillatory velocity vðtÞ ¼ v ω cosðωtÞ has complex amplitude F ω ¼ŶðωÞ À1 v ω , and we identifyŶðωÞ À1 as a generalised friction. However, merely the real part describes dissipation and deserves to be called a friction, which is seen from the mean dissipated power: , averaged over a full cycle of length T p = 2π/ω. Thus, we set the dynamic friction as in particular, ζ(ω) ≥ 0. This is in line with the conventional (Markovian) Langevin equation, _ pðtÞ ¼ Àðζ 0 =mÞ pðtÞ þ ξðtÞ. There, the response is governed bŷ YðωÞ ¼ ½Àiωm þ ζ 0 À1 , implying a static friction, ζ(ω) = ζ 0 .
Equations (6) and (7) (and variants thereof) are the basis of (passive) microrheology experiments [18][19][20][21] , which use observations of Brownian motion to infer the friction ζ(ω) and ImγðωÞ of a probe particle in a complex medium and relate it via the GSER to the local visco-elastic properties. The macroscopic shear viscosity, η 0 ¼ ðk B TÞ À1 R 1 0 C Π ðtÞ dt, is the Green-Kubo integral of the autocorrelation, C Π (t) = 〈δΠ ⊥ (t) δΠ ⊥ (0)〉/V, of shear stress fluctuations δΠ ⊥ (t), given as an off-diagonal element of the stress tensor 16 ; V denotes the sample volume. Similarly by a fluctuation-dissipation relation, the frequency-dependent response coefficientηðωÞ to oscillatory shear is the Fourier-Laplace transform [Eq.
Mathematical framework. For the harmonic analysis of the autocorrelation function C(t) of a stationary time series, we use the Fourier-Laplace transform which is well-defined for all frequencies z in the upper complex plane, C þ ¼ fz j Im z > 0g. Along the imaginary axis, z = iy, it recovers the conventional Laplace transform. For real frequencies ω, the real and imaginary parts ofĈðωÞ describe physically accessible spectra, which are related to each other by Kramers-Kronig integrals 52,67 ; for example, ReĈðωÞ for fixed ω is determined by the full function ImĈðωÞ. The real part is positive, ReĈðωÞ ≥ 0, and most importantly, we have the unique Fourier backtransform: If C(t) is n-times continuously differentiable at t = 0, this implies sum rules for the spectrum (k = 0, 1, …, n): In equilibrium, only positive frequencies are needed as ReĈðωÞ ¼ ReĈðÀωÞ, and the integrals are real-valued. Next, we introduce a memory function of C(t) solely by invoking results from complex analysis 68 .ĈðzÞ as above is a holomorphic function with ReĈðzÞ ≥ 0, i.e., iĈðzÞ is of Herglotz-Nevanlinna type, and ðIm zÞ jĈðzÞj is bounded in C þ . Suppose that C(t) has a regular short-time expansion, Cðt ! 0Þ ' C 0 1 À νt À 1 2 at 2 Â Ã , which implieŝ CðzÞ ' C 0 Àiz ð Þ À1 À νC 0 Àiz ð Þ À2 À aC 0 Àiz ð Þ À3 ð11Þ for large frequencies, |z| → ∞ with j arg zj > δ for some δ > 0. Under these mild requirements, one shows 68 : For givenĈðzÞ there is a unique memory kernelMðzÞ such thatĈ with iMðzÞ of Herglotz-Nevanlinna type andMðzÞ ' ν þ a=ðÀizÞ as |z| → ∞. In particular,MðzÞ corresponds to the autocorrelation function of another (a priori unknown) observable. Iterating the argument yields the continued-fraction representation ofĈðzÞ, well-known from the Zwanzig-Mori projection formalism 16 .
By means of Eq. (9),γðzÞ specifies the memory function γ(t), which has a physical interpretation as the autocorrelator of the fluctuating acceleration ξ(t)/m in Eq. (1), divided by v 2 th . At low frequencies, mγðz ! 0Þ ¼ ζ 0 implies a Green-Kubo relation for the hydrodynamic friction: Short-time expansion. The smoothness of physical molecular trajectories, being solutions to Newton's equations, allows for a short-time expansion of the VACF. Combining with the time-reversal symmetry in equilibrium, Z(t) = Z(−t), only even powers in t contribute and one obtains Zðt ! 0Þ ' k B T P 1 k¼0 c k t 2k =ð2kÞ! with Taylor coefficients c k given from equilibrium matrix elements of powers of the underlying Liouville operator 52 . To connect with the notation of the main text, c 0 = 1/m, c 1 =c 0 ¼ Àω 2 0 , and we put c 2 /c 0 =: Ω 4 . Fourier-Laplace transforming term by term yields the high-frequency expansion ofẐðωÞ and thus ofŶðωÞ ¼ ðk B TÞ À1Ẑ ðωÞ, which is purely imaginary:Ŷðω ! 1Þ ' P 1 k¼0 c k ðÀiωÞ À1À2k ¼ c 0 =ðÀiωÞ þ c 1 =ðÀiωÞ 3 þ . Using Eq. (7), we have ζðωÞ ¼ jŶðωÞj À2 ReŶðωÞ, which implies that for high frequencies the friction vanishes, ζ(ω) ≡ 0, at all orders in ω → ∞. A similar situation is familiar from calculus text books: f(x) = e −1/x has a Taylor series f(x) ≡ 0 at x = 0; the radius of convergence is 0.
The expansion ofŶðωÞ can be represented as a continued fraction that has the same large-ω asymptotics up to terms of order ω −5 :

YðωÞ ' 1=m
Àiω þ introducing ω 2 1 :¼ Ω 4 À ω 4 0 À Á =ω 2 0 for brevity. This truncation is an excellent description of our data forŶ 00 ðωÞ at high frequencies, with ω 0 and Ω obtained from fits to Z(t), see Fig. 3d-f. For the memory kernelγðωÞ, one reads off using Eq. (6), implying for the memory function in time domain: Long-time tails. In an unbounded fluid, momentum conservation leads to persistent velocity correlations, Zðt ! 1Þ ' v 2 th ðt=τ 1 Þ À3=2 , which was explained in terms of hydrodynamic backflow and diffusion of a momentum vortex 16,48,49 . The tail induces a small-ω singularity in the frequency domain 69 , which reads for the memory kernel: mγðω ! 0Þ ' ζ 0 1 þ ðτ 1 ζ 0 =mÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi À4πiωτ 1 p Þ ½ , using Eq. (4) and the hydrodynamic friction ζ 0 :¼ mγð0Þ ¼ k B T= R 1 0 ZðsÞ ds. In the framework of the creeping flow equations, Stokes found 17 in terms of the vorticity diffusion time τ fl and an effective particle mass m fl ; matching with the previous expression for the asymptote ofγðωÞ, one identifies τ fl ¼ 4πðζ 0 =mÞ 2 τ 3 1 . The real part yields the dynamic friction, showing that its macroscopic limit, ζ 0 = 6πη 0 a, is approached from above as ω → 0 (see Fig. 3b). For water and the supercooled liquid, a different type of power-law decay, Z(t)~−t −5/2 , was found (Fig. 3e, f). For a general long-time tail of the VACF, Z(t)~t −σ with σ > 1, the memory function γ(t) asymptotically inherits a power-law decay with the same exponent, but of opposite sign 70 : which follows from Eq. (3) and by invoking a Tauber theorem 69 . Without any adjustable parameter, the prediction is in excellent agreement with our data for γ(t) in case of the LJ fluid. Inspection of a few examples (Figs. 4d, e and 5c) suggests that, in order to accommodate the sign change of the tail, the number of zero crossings (knots) in γ(t) is increased by one relative to Z(t).
Adapted Filon algorithm. The computation of the frequency-dependent friction requires a robust numerical Fourier transform, for which we developed a physicsenriched version of Filon's quadrature formula. The goal is to evaluatef ðωÞ ¼ R 1 0 e iωt f ðtÞ dt for a function f(t) sparsely sampled on an irregular grid t 0 = 0, t 1 , …, t n for an arbitrary set of frequencies (ω j ). The idea of Filon's algorithm is to interpolate f(t) by elementary functions between the grid points (usually parabolas), thereby reducing the Fourier integral to a finite sum of integrals, for which analytic expressions exist. Anticipating that the normal physical decay of correlation functions is exponential, we approximate f ðtÞ % a k e Àγ k t in the interval [t k , t k+1 ] with a k and γ k fixed by f(t k ) and f(t k+1 ). Then, t n a n e ðiωÀγ n Þt dt : Spurious low-frequency oscillations of the transform are removed by smoothly truncating the integral at t n , here by assuming a terminal exponential decay of f(t), which leads to the last term on the r.h.s. of Eq. (22). In order to preserve the shorttime properties of f(t) we fit a polynomial in t 2 to the first few data points and solve the integral analytically; this improves the high-frequency behaviour off ðωÞ.
The dynamic friction ζ(ω) and the memory function γ(t) are obtained from MSD data as follows (Fig. 2): The timescale-dependent diffusion coefficient, D(t): = ∂ t MSD(t)/6, is obtained from numerical differentiation. In all cases studied, it grows out from zero, has a maximum, and converges slowly towards the longtime diffusion constant D ∞ = D(t → ∞), see Supplementary Fig. 1. Using the above algorithm, we compute 71Ẑ ðωÞ ¼ D 1 À iω R 1 0 dt e iωt ½DðtÞ À D 1 . Then, ζ(ω) is given by Eq. (4) and is transformed back to the time domain with the same algorithm [Eq. (9)]; in particular, we use again a smooth, exponential cutoff. In Fig. 5, the numerical procedure is successfully tested against the analytical model with high accuracy.
Deconvolution in time domain. Inversion of the convolution in Eq. (3) yields the memory function γ(t) directly 72 , without resorting to the frequency domain. Numerically, it is not easy to obtain accurate and robust results, and a variety of algorithms have been developed, see ref. 36 for a comparative study. The presence of _ ZðtÞ in Eq. (3) is removed by integration, yielding ZðtÞ ¼ Zð0Þ À R t 0 ds Gðt À sÞ ZðsÞ with the integrated memory GðtÞ :¼ R t 0 ds γðsÞ. Discretising on a uniform time grid, t i = iΔt (i = 0, 1, …), and employing the trapezoidal rule for the integral, a recursion relation for G i :¼ Gðt i Þ with the initial value G 0 = 0 follows 36 : Going beyond Kowalik et al. 36 , we introduce a predictor-corrector scheme to stabilise the numerical solutions: In the predictor step, one evaluates G Ã i and G Ã iþ1 from Eq. (23). Afterwards, the weighted average G i :¼ ðG iÀ1 þ 3G Ã i þ G Ã iþ1 Þ=5 manifests itself as the corrector step. Results for G(t) can be found in the Supplementary Fig. 2. Finally, the memory function γ(t) = ∂ t G(t) is obtained by central differences. If one starts from MSD data on a sparse (e.g., geometrically spaced) time grid, a cubic spline interpolation of the MSD in the variable t 2 is suitable to sample Z(t) on a uniform grid of up to 10 5 points.
Molecular dynamics simulations. Simulations of liquid water were performed with the GROMACS 5.1 package 73 using the SPC/E water model, which was shown to accurately reproduce the linear absorption spectra of water from experiments and ab-initio MD simulations up to frequencies of about 30 THz 74 . The system of 3007 molecules in a cubic box of linear size 4.49 nm was equilibrated at 300 K and 1 bar following standard procedures 36 . Correlation functions were obtained from an constant-particle-number, constant-volume, constant-energy (NVE) simulation over 275 ps with the velocity-Verlet integrator and a time step of 1 fs, using double floating-point precision to achieve good energy conservation. The frequencydependent viscosity was computed from additional NVE runs, totalling to 49 ns.
For the other two liquids, we used the massively parallel software HAL's MD package 75 (version 1.0α6), permitting the sampling of dynamic correlations on a sparse time grid and featuring smoothly truncated potentials to virtually eliminate any energy drift. The mono-atomic fluid consists of 10 5 particles interacting pairwise via the LJ potential, UðrÞ ¼ 4ε ðr=σÞ À12 À ðr=σÞ À6 Â Ã , truncated for r ≥ 2.5σ; a unit of time is defined by τ LJ :¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi mσ 2 =ε p . Equilibration in the NVE ensemble at number density ϱ = 0.8σ 3 and thermal energy k B T = 1.3ε followed the protocol given by Roy et al. 76 . The supercooled liquid was realised by a Kob-Andersen 80:20 binary mixture 77 of 64,000 LJ beads at ϱ = 1.2σ −3 and T * : = k B T/ε = 0.6, equilibrated over a time span of 9,000 τ LJ , and we traced the species of the larger particles. The chosen temperature is well below the melting point 78 , T * ≈ 1.03, and is on the onset of universal scaling behaviour according to mode-coupling theory (MCT) 40 ; here, the value of the critical temperature is T Ã MCT % 0:43. The simulations generate trajectories r(t) of an ensemble of labelled particles for each fluid; our main observable is the mean-square displacement MSDðtÞ :¼ jrðtÞ À rð0Þj 2 for lag time t. For both liquids, single-particle MSDs were averaged from ten production runs, each over 10 8 integration steps of length 0.001 τ LJ .
Control simulations of a single particle in its pinned cage are based on an equilibrated sample of the LJ fluid with 10 6 particles. MSDs were recorded after equilibration of the mobile particle in its static environment over 100 τ LJ and were averaged over 10 6 different cages, computed in parallel, to remove spurious oscillations. Technically, the setup was realised by making two initially identical copies of the fluid interact with each other: the first copy contains the immobile matrix, the second one the tracers (not interacting with each other).

Data availability
The data sets generated and analysed during the current study are available from the corresponding author upon reasonable request.

Code availability
Primary data were generated with open source software as indicated in the "Methods" section. The source code used to analyse the data for the current study is available from the corresponding author upon reasonable request.