Critical fluctuations and slowing down of chaos

<jats:title>Abstract</jats:title><jats:p>Fluids cooled to the liquid–vapor critical point develop system-spanning fluctuations in density that transform their visual appearance. Despite a rich phenomenology, however, there is not currently an explanation of the mechanical instability in the molecular motion at this critical point. Here, we couple techniques from nonlinear dynamics and statistical physics to analyze the emergence of this singular state. Numerical simulations and analytical models show how the ordering mechanisms of critical dynamics are measurable through the hierarchy of spatiotemporal Lyapunov vectors. A subset of unstable vectors soften near the critical point, with a marked suppression in their characteristic exponents that reflects a weakened sensitivity to initial conditions. Finite-time fluctuations in these exponents exhibit sharply peaked dynamical timescales and power law signatures of the critical dynamics. Collectively, these results are symptomatic of a critical slowing down of chaos that sits at the root of our statistical understanding of the liquid–vapor critical point.</jats:p>

section_content: F luctuations are sovereign in critical phenomena 1, 2. They rule fluids at the liquid-vapor critical point, a unique instability that punctuates the space of thermodynamic states 3. First established experimentally by Andrews 4, the liquid-gas critical point was given a molecular explanation shortly thereafter by van der Waals 5. In van der Waals' picture, now a paradigm in liquid state theory 6, 7, repulsive forces largely determine the structural arrangements of molecules in non-critical liquids, not the attractive forces. Near the liquid-vapor critical point, however, their roles reverse and the paradigm shifts 8 ; dynamical fluctuations reach macroscopic magnitude and overrule molecular size, shape, and interactions in dictating bulk behavior. These fluctuations are generated by the nonlinear dynamics of classical critical fluids. Yet, the relationship between the microscopic instability of the dynamics and the thermodynamic singularity has never been entirely clear 9. 

While the field of critical phenomena continues to absorb increasingly diverse systems 10, the basic phenomenology is firmly established 1. Its taxonomy is built on scaling and universality, the similar behavior of dissimilar systems. Despite the early discovery of their critical points, fluids were somewhat resistant to classification 5. Simulations 11 and theory 12 were, and continue to be, integral in providing mechanistic insights, the location of the critical point, and estimates of static critical exponents [13] [14] [15] [16]. Through simulations, the classical atomistic dynamics of fluids are known to be chaotic 17, a part of the machinery of nonlinear dynamics. Measuring deterministic chaos has given insights into the physical mechanisms of the jamming transition in granular materials 18, self-organizing systems 19, evaporating collections of nuclei 20, 21, and the phase changes of atomic clusters [22] [23] [24] [25]. In addition, the dynamics of model spatially-extended systems have recently begun to collect into dynamic universality classes 26. These findings are part of efforts to coalesce statistical physics and nonlinear dynamics, and they reinvigorate the question of how fluids, specifically the properties of their molecular dynamics, fit within this phenomenological architecture of critical phenomena. 

At the liquid-vapor critical point, how do correlations in molecular positions overcome the destabilizing force of deterministic chaos in the molecular dynamics? Here, we resolve the instability of the molecular motion that generates this critical phenomenon. Because of the absence of long-range order-and an inability to make a small vibration approximation, as in solids, or a molecular randomness hypothesis, as in gases 27 -the dynamical instability in critical fluids has been largely grounded in purely statistical terms. However, critical correlations imply structural organization that is intrinsically opposed by a chaotic dynamics. Both numerical simulations and analytically tractable models show here that the critical dynamics carry signatures of this internal tension between order and chaos. We find that Lyapunov exponents 28, a measure of chaos and dynamical instability, are minimal at the liquid-vapor critical point. Signatures also appear in the finite-size scaling of fluctuations in these observables. The fluctuations decay as a power law towards the thermodynamic limit with the slowest rate at the critical point. Overall, these results suggest this singular state is a limit of dynamical order that long-range correlations can impose on the dynamics.

section_name: Results

section_content: To analyze the molecular dynamics at the liquid-vapor critical point, we numerically simulate a homogeneous, single-component, non-associated, equilibrium fluid (Fig. 1a ) and apply techniques from nonlinear dynamics. The fluid consists of N molecules interacting pairwise through van der Waals forces, repelling (attracting) at distances of a (few) molecular diameter(s) according to the Lennard-Jones (LJ) potential 29. As order parameters, we use the hierarchy of 6N spatiotemporal (Lyapunov) vectors. Each vector has an associated exponent, λ i indexed i = 1,…,6N 28 and in descending order, that measures the contribution of each vector to the global dynamics. Larger exponents indicate more unstable phase space directions 30. We calculate the full Lyapunov spectrum at the critical density ρ c and over a range of temperatures including the critical temperature, and we choose the energy density to fix the mean kinetic temperature, T. Each λ i is a time-average of a 10 6 time step trajectory. 

Suppression of chaos at the critical point. When approaching the liquid-vapor critical point from the supercritical regime, T > T c, the Lyapunov exponents decrease monotonically (Fig. 1 ). This decline represents a "slowing down" in the divergence of trajectories as T approaches T c and accompanies the transient formation of spatial regions in the system that will separate into coexisting vapor and liquid phases below T c. What drives this decrease in the largest exponent, λ 1, is a suppression of the fastest dynamical events 31. Here, these events are molecular collisions that sample the repulsive part of the intermolecular potential. As the system "practices" phase separating, as Widom put it 5, fewer of these high-energy collisions occur 32. Away from the critical point, the van der Waals picture is the foundation for the statistical-mechanics of liquids 6. Perturbative treatments, for example, assume the structure of a dense, monatomic liquid resembles that of a hard sphere fluid and, to a first approximation, the attractive interactions have little effect on the liquid structure 6. At the critical point, though, repulsive intermolecular forces play a subordinate role compared to critical fluctuations. While the van der Waals picture of fluids was only recently shown to extend to the Lyapunov exponents and their fluctuations 33, it does not hold near the critical point 8. A steady 7% decrease with temperature and the clear minimum near T c in the first Lyapunov exponent are direct evidence that critical conditions inhibit the effect of repulsive forces on the dynamics, consistent with the breakdown of the van der Waals picture. Also evident from Fig. 1 is a minimum in the largest Lyapunov exponent near the known 34 critical point, (T c, ρ c ). The minimum in this exponent, and the maximum in the Lyapunov time, 1/λ 1, near T c means the critical dynamics are predictable over longer timescales because initially similar configurations do not diverge as quickly (Fig. 1c ). For the finite-size systems we simulate, the doubling of the Lyapunov time is in accord with the doubling of the characteristic time of the autocorrelations in the kinetic energy (Supplementary Note 1). A suppression of the dynamical instability in the critical regime aligns with the physical intuition that long-range correlations are at play near T c. However, for the finite-size systems we simulate, the dynamics are not entirely predictable at T c ; the first exponent has a non-zero value through the temperature range including T c, which suggests the continuous liquid-vapor phase change remains chaotic throughout and into the coexistence region for finite-size systems. These chaotic critical dynamics differ from the jamming transition in granular materials, which seems to be a transition from a chaotic to a non-chaotic state 18. 

To more fully resolve dynamical instability across the critical point, we also calculated the full spectrum of Lyapunov exponents. From simulations of two and three-dimensional liquids, the shape of the spectrum depends on the kinetic temperature and density (Supplementary Fig. 2 ) 32, 35, 36. Figure 1b shows the long-time Lyapunov spectrum at the critical density for temperatures spanning the critical temperature, T c. These data for N = 1000 converged to <1% and show all unstable vectors have exponents that decrease when approaching the critical point from above. But, only the most unstable vectors have exponents with minima at the critical point. Critical correlations appear to have the largest impact on the vectors with scaled index up to i/3N ≈ 0. 18 (Supplementary Fig. 3 ). Exponents beyond this point decrease monotonically with temperature. The spectrum is also compressed at the critical point (Fig. 1c inset), meaning there is a weaker preference for trajectories to diverge in the direction of any given vector. 

Many models and molecular simulation techniques have been used to probe the liquid-vapor critical point 37, 38. But, the longtime Lyapunov exponents in Fig. 1 are the first glimpse of its chaotic properties. These exponents are intensive observables in the liquid and supercritical phases, but there is a slight dependence of the spectrum on the N at ρ c that affects the location of the minimum (Supplementary Figs. 4 and 5 ). A similar weak dependence of the leading Lyapunov exponent on N was found for a two-dimensional LJ fluid 39. To support our numerical simulations and pinpoint the location of the minimum, we analyze two analytically tractable approximations to a model system. The Lyapunov exponents and the Lyapunov time τ at the critical temperature T c depend strongly on the nature of the approximation. 

Because the liquid-vapor critical point is generally considered to belong to the Ising universality class 14, 16, we consider a continuous analog of the Ising model. The model is onedimensional system of N interacting, classical particles of mass m in a bistable potential. Two approximations of this model capture the features of the heat capacity at the liquid-vapor critical point: the singularity and the finite-jump discontinuity 8. To analyze the effect of a diverging heat capacity on the Lyapunov exponents, we treat the average effects of the anharmonic term in the Hamiltonian 9 by imposing a temperature-dependent weakening of the anharmonic restoring force on each particle. The Lyapunov exponents are given by

There is a conjugate pair of exponents for each of the N particles. 

Approaching the critical point, the first Lyapunov vector will generate the instability that leads to the phase transition; the Lyapunov exponents become degenerate and their associated Lyapunov vectors become linearly dependent. From the last expression, we see that λ ± → 0 as T → T c. That is, all Lyapunov exponents vanish at the critical temperature according to the power law |λ ± | ∝ |T-T c | 1/2 (Fig. 2 ). The corresponding Lyapunov time diverges τ ± = 1/λ ± ∝ |T-T c | -1/2 at T c and coincides with the divergence in the isochoric heat capacity. Applying another (mean-field) approximation to the coupledparticle system, we can analyze the Lyapunov exponents when the heat capacity has a simple discontinuity and further confirm that the critical point is a limit of dynamical order (Supplementary Fig. 6 ) [40] [41] [42]. As an order parameter, we use the mean particle displacement. This order parameter is nonzero for T < T c and zero

for T ≥ T c (Supplementary Fig. 7 ). While the Lyapunov exponents are nonzero, they have a minimum at T c. The Lyapunov times have a maximum. The rate of change of the Lypaunov exponents and times with temperature have a finite-jump discontinuity at T c mirroring the jump in the isochoric heat capacity, a typical feature of mean-field theories. It is well known that this jump in the heat capacity predicted by mean-field theory is embedded in a logarithmic divergence when fluctuations are included. Also known is that the exact scaling of observables at the critical point is strongly affected by statistical correlations in molecular positions. 

Overall, this model of nonlinearly coupled particles gives evidence that the minima in the Lyapunov exponents are at T c for pure fluids. The minima also support the sensitivity of the most unstable Lyapunov vectors to the critical dynamics and motivate a closer look at their finite-time fluctuations in our simulations. 

Suppression of approach to the thermodynamic limit. Cooling the fluid towards the critical point, there are long wavelength fluctuations in density that cause the correlation length to diverge. In our simulations, though, the correlation length saturates at the size of our simulation cell, L, truncating longer wavelength fluctuations; only replicating periodic fluctuations through our boundary conditions may be a source of the shallowness of the minima in the largest Lyapunov exponents (Fig. 1b ). Fleeting clusters of all sizes up to L, which will eventually become liquid, begin to appear. These clusters' ephemeral existence affects instability in the critical dynamics on short timescales. As the structure evolves, the dynamics continues to temporarily sample phase space domains where trajectories diverge more quickly and more slowly than the average, domains that will determine the Lyapunov vector directions and the finite-time estimates of the Lyapunov exponents. To analyze the fluctuations in finitetime Lyapunov exponents λ i (t), we divided trajectories of 10 6 time steps uniformly into windows of 100 time steps. Distributions of λ 1 (t) are shown in Fig. 3a. 

Non-trivial 43 power laws are apparent in the decay of finitetime Lyapunov exponent fluctuations with system size. At thermodynamic equilibrium, the precise scaling of relative fluctuations is often called self-averaging 44. Loosely, a system is self-averaging with respect to a given property, X, if the value of the thermodynamic observable corresponds to the average over independent subsystems or, in this case, time windows. More precisely, the relative fluctuations of an observable

with ΔX = X-〈X〉, the wandering exponent γ, and a generalized diffusion coefficient D X. If the observable is self-averaging, the relative variance R X of the property X vanishes in the thermodynamic limit: R X → 0 when N → ∞. The wandering exponent γ can have values between 0 and 1-a value of one meaning the observable is strongly selfaveraging. Weakly self-averaging observables have γ < 1 and nonself-averaging observables have γ = 0. 

Because of the statistical independence of spatial domains in the system, equilibrium observables are often strongly selfaveraging with γ = 1. However, these domains become statistically dependent near a critical point because of the divergence of the correlation length 1. Statistical signatures of dynamical observables, like the Lyapunov exponents, are still being elucidated 26, 45. Deep in the liquid state, for example, the fluctuations for all but the first Lyapunov exponent (the "bulk") are strongly self-averaging. The first exponent fluctuations, however, self-average weakly with a rate of decay γ < 1 that depends on the length scale of the interparticle interactions and captures the van der Waals picture of dominant repulsive forces 33. How does this dynamical version of the van der Waals picture change at the critical point? Numerical calculation of the entire Lyapunov spectrum comes at significant computational cost, cost that increases significantly when scaling with system size 46. However, to quantify the scaling behavior of the finite-time Lyapunov exponent fluctuations with system size, we simulated 11 systems ranging from N = 100 to 2000 over the same range of temperatures at fixed density ρ c. Over this range of system sizes, fluctuations in the first finite-time Lyapunov exponent, λ 1 (t), as measured by the diffusion coefficient D λ 1 ðN; T; ρÞ, decay with system size as N -γ. This scaling holds over the temperature range T = 0. 8-1. 4 spanning the critical point (Fig. 3a ). The wandering exponent varies between 0. 9 and 0. 6 (with its smallest value near T c ), showing temperature controls the magnitude of the wandering exponent through the structural changes and spatial correlations it brings about near T c. A higher-order statistical analysis indicates this weak self-averaging is, at least in part, due to the non-Gaussian features of the distributions (Supplementary Note 2). Distributions of the kinetic energy per particle have a strong Gaussian character under all conditions we simulate. 

Most prominent in the temperature dependence of the wandering exponent of the first exponent is the peak at T = Fig. 3 Finite-size scaling of fluctuations in dynamical observables. a Fluctuations in the kinetic energy and Lyapunov exponents decay with increasing number of particles (log-log): the diffusion coefficients follow a power law D X ðN; TÞ $ N Àγ X ðTÞ for 11 system sizes from N = 100 to 2000. Data shown for λ 1 (t). The corresponding probability distributions also concentrate around the long-time value λ 1 (Inset: first Lyapunov exponent at ρ = ρ c = 0. 317 and T = T c = 0. 937). b Wandering exponents, γ X (T), of the first finite-time Lyapunov exponent and the kinetic energy per particle, K/N, peak at the critical temperature, T c, at the critical density ρ c. Error bars indicate the standard deviation (1σ) for parameter estimates from linear fits in a 0. 962 (Fig. 3b ). It shows the critical dynamics self-average most weakly in the direction of the first Lyapunov vector. As a reference, the fluctuations in kinetic energy per particle also decay with system size; the wandering exponent peaks at T = 0. 937, further confirming the location of T c and further quantifying the decay of global fluctuations (Fig. 3b ). This critical temperature agrees with that found from grand canonical Monte Carlo simulations 34. The peak in -γ for the first Lyapunov exponent is just above the critical temperature. However, the wandering exponent -γ and its maximum are affected by the range of system sizes and the statistical error in the linear fits of Fig. 3b. Near the maximum, -γ values appear to converge from below and, so, we take them as a lower bound (Supplementary Fig. 9 ). The order of the thermodynamic limit and the time limit in the definition of the Lyapunov exponents can determine whether a dynamics is chaotic or not 45, 47. Because of the subtleties of the time and thermodynamic limit order, fluctuations in the finitetime Lyapunov exponents of both dissipative and conservative dynamical systems have recently been subject to a finite-size scaling analysis 26, 33, 45. In all cases reported to date, the selfaveraging of the first exponent is distinct from that of the bulk, the set of 3N-1 positive exponents that exclude the first. In spatially-extended dynamical systems where it is known, the scaling of fluctuations is homogeneous across the bulk of the spectrum 26, 48, 49. Liquids show this behavior, for example, and all the bulk exponents are strongly self-averaging 33. However, the critical dynamics break this scaling symmetry-a significant fraction of the bulk exponents self-average weakly as shown in Fig. 4. This scaling feature is so far unique to the liquid-vapor critical point. It is also apparent in the self-averaging of the entire Lyapunov spectrum through the average diffusion coefficient hDðλÞi ¼ P 6N i¼1 Dðλ i Þ=6N, which contrasts that of the largest exponent (Fig. 4a ). The corresponding wandering exponent has an inflection point around T = 0. 962. Increasing the fraction of exponents included in the average shows that a portion of the more unstable vectors have a γ-peak that vanishes with increasing index. 

While fluctuations appear to decay with system size in all unstable phase space directions on our accessible time and length scales, the rate of decay is far from homogeneous. The γ-spectrum quantitatively resolves the rates at which this unique thermodynamic equilibrium state emerges from the molecular dynamics (Fig. 4a ). The good data collapse in Fig. 4c reveals clear scaling functions for the diffusion coefficient spectrum, DðT; ρ ¼ ρ c Þ; only three representative temperatures are shown. This scaling function is a system-size independent measure of the finite-time fluctuations in the Lyapunov spectrum, fluctuations caused by the local heterogeneities in phase space sampled by our simulated trajectories. The basic form of this scaling function is similar to that seen for simple liquids 33 and Hamiltonian lattices 45. Our calculation of D λ 1, a dynamical invariant in the long time limit, leads directly to a dynamical timescale for λ 1 (t) fluctuations, τ D ¼ 1=D λ 1 N γ. This timescale peaks just above the critical temperature 34. Although these vectors are highly active on short timescales, the fluctuations destructively interfere on longer timescales (showing a net suppression, Fig. 1 ). The peak in the correlation time of the kinetic energy per particle is evidence of critical slowing down (Fig. 4b ), and confirms the location of the critical temperature.
section_name: Dynamics of a simple fluid.

section_content: Based on early foundational work, there are known connections between hydrodynamics and Lyapunov vectors with small, but finite, exponents 50. The dynamics of these so-called hydrodynamic modes can characterize macroscopic transport 36, 51. Here, above the liquid-vapor critical point, the Lyapunov exponents of these modes depend strongly on temperature but show no clear signs of the critical dynamics. Recently, though, it was found that hydrodynamics can shape the most unstable directions as well; numerical simulations of two prototypical Hamiltonian lattices, the FPU-β and the Φ 4 models, demonstrated that longrange correlations can have a dramatic effect on the fluctuations of the leading finite-time Lyapunov exponent and even cause them to diverge with system size 45. The minima we observe in the long-time exponents of fluids near T c show the most unstable Lyapunov vectors are also sensitive to the correlations induced by critical phenomena. Dynamical systems have begun to collect into universality classes through the finite-size scaling of Lyapunov exponent fluctuations. Numerical estimates of the wandering exponent γ reported 26, 48, 49 so far in the literature suggest the Lyapunov exponents of dissipative systems are weakly self-averaging, γ < 1, with the dynamics of the first Lyapunov vector belonging to the Kardar-Parisi-Zhang (KPZ) universality class 52. Fluctuations of the maximum Lyapunov exponent in Hamiltonian lattice models, FPU-β and the Φ 4, however, are not self-averaging and each belong to their own universality class 45, 53. There is mounting evidence that their non-KPZ behavior is a consequence of longrange spatiotemporal correlations and not a more universal feature of Hamiltonian dynamics 33. The weakening of the selfaveraging behavior and peak in -γ near the liquid-vapor critical point is further support for this hypothesis. Because the molecular dynamics at the liquid-vapor critical point and in the liquid state 33 are Hamiltonian, however, the dynamic universality class of the Lyapunov vector dynamics remains an open question. 

Though distinct from the fluids here, models for non-additive systems with long-range (gravitational or electrostatic) forces have shed some light on the relationship between microscopic chaos at phase transitions. The dynamics of the Hamiltonian mean-field model, for example, are chaotic for finite-size systems but not for an infinite size system 47. The dynamics and thermodynamics that have been studied both numerically [54] [55] [56] [57] and analytically 58. And, in contrast to the results here, near the second-order phase transition, both the largest Lyapunov exponent and kinetic energy fluctuations have a maximum near the critical temperature 55. 

Here, by treating the nonlinear dynamics directly with numerical simulations and analytical models, we have resolved the phase space directions responsible for the thermodynamic instability and the breakdown of the van der Waals picture at the liquid-vapor critical point. The long-time Lyapunov spectra show that critical dynamics of finite-size fluidic systems are less sensitive to the detailed features of intermolecular forces but also initial conditions. Through numerical simulations and analytical models, critical conditions appear to constrain the dynamics so that different phase space directions have a relatively homogeneous degree of instability and scaling features in finite-time Lyapunov exponent fluctuations that are so far unique to critical dynamics. 

Correlations in molecular positions at the critical point span the length scales of intermolecular forces to the entire system. They imply structural organization that is intrinsically opposed by the chaotic dynamics. The mechanisms balancing this internal tension between order and dynamical instability, however, are subtle. As a result, theoretical explanations for the mechanical origins of critical phenomena are uncommon. Continuous transitions in crystals are a notable exception, where structural changes arise through the instability of a lattice vibration 9. There, the mode responsible for the phase transition is a collective excitation whose frequency decreases anomalously during an approach to the transition point. For example, in SrTiO 3 the frequency of a soft phonon mode decreases substantially and, ultimately, freezes at the transition temperature when approached from below 59. 

Unlike continuous crystal-crystal phase transitions, there is not one unique unstable vector with a vanishing frequency at the liquid-vapor critical point. Instead, the whole spectrum softens with a subset that have extrema near the critical point; both the long-time Lyapunov exponents and their fluctuations on short times reflect their high sensitivity to long-range correlations of molecular positions. These vectors do not appear to completely freeze, at least not in finite-size systems, but do exhibit large fluctuations with a peak in dynamical timescales indicating the critical slowing down of chaos, the stabilization of unstable vectors, and a longer memory of initial conditions. In short, the relative mechanical stability of molecular motion underlies the bulk behavior of fluids at this thermodynamic instability.
section_name: Discussion

section_content: Model non-associated fluid. Our system is the three-dimensional, periodic Lennard-Jones (LJ) fluid 29. The Hamiltonian of this collection of N particles is

The pairwise interaction potential between particles i and j a distance rij apart is given by

The parameter ε ij corresponds to the well depth at the equilibrium distance and measures the strength of the interaction between particles i and j. For the singlecomponent fluid we consider, all interactions are identical (ε ij = ε). The first term of the potential energy function takes into account short-range repulsive interactions. The second term corresponds to the attractive part of the interaction, which acts over a comparatively longer range. The parameter σ stands for the distance at which the attractive and the repulsive forces are equal and can serve as a measure of the particle size. We work in LJ reduced units with the distance r ¼ r=σ, density

, and Lyapunov exponents λ ¼ λ ffiffiffiffiffiffiffiffiffiffiffiffiffi mσ 2 =ε p. All particles have unit mass m = 1, ε = 1, and σ = 1. 

Numerical simulation methods. All numerical data are from molecular dynamics simulations of constant energy trajectories. The initial conditions were sampled from a constant temperature trajectory using the Berendsen thermostat with a time constant of 0. 5. Subsequent NVE trajectories are updated with the velocity Verlet algorithm using a time step of 10 -3. Their total energy is well conserved and the kinetic temperature fluctuates around the temperature of the sampled NVT trajectory. We shift and truncate the interparticle potential 60 and use a Verlet list with a 2. 5σ cutoff and a "skin" of 0. 1σ. All data are from single-precision calculations. 

At fixed number of particles N, volume V, and energy E, we simulate deterministic trajectories of this equilibrium fluid and the dynamics of Lyapunov vectors in tangent space 46. The ith Lyapunov vector components are first variations in position and momentum (δq ij, δp ij ) T with i, j = 1, …, 6N and evolve according to their own linearized Hamiltonian equation of motion. We numerically solve this equation of motion with the linearized form of the velocity Verlet algorithm used to evolve trajectories and orthonormalization at every time step 46. The initial basis sets are random and orthonormal. During a transient, that we discard, the first vector orients itself parallel to the maximally changing tangent space direction. Regular orthonormalization restricts the collapse of the remaining vectors onto the most expanding tangent space direction. The algorithm requires the second derivatives (Hessian) of the interaction potential at every time steps. We use forward differences of the analytical gradients with a displacement of 10 -4. 

Within the linearized limit, the expansion or the compression factor along the phase-space direction of the ith Lyapunov vector over time t is e Γ i ðtÞ. The corresponding finite-time Lyapunov exponent is λ i (t) = Γ i (t)/t. The complete finitetime Lyapunov spectrum, {λ i (t)} is calculated from the set of Gram-Schmidt vectors with standard methods 61, 62. We evaluate the full Lyapunov spectrum at each time step using the norm

The ith finite-time exponent over a time interval t = t′-t 0 has the form:

At each thermodynamic condition, simulations of long trajectories are 10 6 time steps and span a time of t final = 1000 in reduced time units. For each long trajectory we obtain estimates of the long-time exponents {λ i }. 

Each trajectory is also used to calculate the finite-time Lyapunov exponents by dividing the long trajectory into constant time intervals of size of 0. 1. Every time segment has a unique initial condition, so the partitioning of the long-time trajectory produces an ensemble of 10,000 finite-time trajectories. Finite-time Lyapunov exponents are fluctuating variables. We estimate the magnitude of their fluctuations over a time interval, t, with the diffusion coefficients {D(λ i )} 26, 33, 45 and the variance,

Averages 〈⋅〉 are over an ensemble of 10 4 trajectories, each having time span t = 0. 1 in reduced time units. The average 〈λ i 〉 = λ i is the average of the ith Lyapunov exponent over an entire long trajectory. While the finite-time Lyapunov exponents depend on the chosen norm, the diffusion coefficient is a dynamical invariant 45. 

To probe the self-averaging property of finite-time Lyapunov exponent fluctuations, we ran trajectories over a range of system sizes. From N = 100 up to N = 2000, we calculated the diffusion coefficient of the finite-time {λ i (t)} over the trajectory ensemble for each system size. We scaled the number of molecules N and the volume V to ensure the thermodynamic limit of the microcanonical ensemble: N, V → ∞ keeping the number density ρ = N/V and energy density e = E/V constant. According to the equipartition theorem, the kinetic temperature of the system is given by T = 2〈E kin 〉/3Nk B. We analyze the scaling of fluctuations in dynamical variables with system-size N for a range of T with ρ fixed at ρ c = 0. 317. From the scaling of {D(λ i )} with N, we estimate the values of the wandering exponents γ. Estimates of the wandering exponent with up to N = 1000 or N = 2000 particles does not affect the location of maximum in -γ or qualitative features of its temperature dependence, suggesting our estimates are well converged. The power-law decay of the diffusion coefficient of the finite-time Lyapunov exponents with system size was also found to be robust to changes in the chosen time interval. 

Analytical model of nonlinearly coupled oscillators. We consider a onedimensional system of N interacting, classical particles of mass m in a bistable potential

i. The single-particle Hamiltonian is

where 〈ij〉 indicates all pairwise particle interactions with interaction strength J. 

Here, we use A = κ 4 /4 and B = -κ 2 /2 with positive constants κ 2, κ 4 > 0. We drop the particle index in what follows. The minima of the potential wells are at x = ±x 0 with x 0 = κ 2 /κ 4 with potential energy V 0 ¼ À 1 4 κ 2 x 2 0. The interaction energy for particles at their minima is Jx 2 0. The Ising limit is when the barrier separating the two wells is high relative to the interaction energy, jV 0 j=Jx 2 0 ) 1, so that each particle is localized in the left (0 state) or right (1 state) well. Particles may randomly occupy these 0 and 1 states. 

Two approximations of this model capture the singularity and the finite-jump discontinuity of the isochoric heat capacity at the liquid-vapor critical point. Both approximations decouple the oscillators. Let δx i = (δq i, δp i ) be a variation of a single-particle phase point. Its equation of motion is

where Ω is the symplectic matrix, D 2 H 1 is the second derivative of the relevant single-particle Hamiltonian, and

The submatrix T is diag(1/m) and D 2 x H 1 is the Hessian matrix. The eigenvalues of

In both approximations below, each particle contributes one Lyapunov exponent. As a result, the Lyapunov spectrum and the Kolmogorov-Sinai entropy are system-size extensive, h KS = Nλ +. 

Soft-mode approximation: To analyze the effect of a diverging heat capacity on the Lyapunov exponents, we treat the average effects of the anharmonic term in Eq. ( 5 ) by imposing a temperature-dependent weakening of the anharmonic restoring force on each particle 9. We set J = 0 and x 4 → x 2 〈x 2 〉 0, using the canonical average of the displacement squared for the harmonic oscillator 〈x 2 〉 0 = k B T/κ 2. The Hamiltonian of a single particle becomes

The critical temperature is k B T c ¼ 2κ 2 2 =κ 4, which we can use to express the free energy of a single particle

for T > T c with the free energy of a single harmonic oscillator F 0,1. Using the thermodynamic relation C V = -T∂ 2 F/∂T 2, the isochoric heat capacity of each particle is

and C V,1 → ∞ as T → T c from above. The canonical ensemble average Lyapunov exponents are

If the space and time averages are equivalent, then λ ± ¼ λ ± h i ¼ lim t!1 λ ± ðtÞ. The right eigenvectors of ΩD 2 H ′ 1 are (1/m, λ + ) T and (1/m, -λ + ) T. At the critical point, the Lyapunov exponents are degenerate and these eigenvectors are linearly dependent. Notice that because λ ± = λ ± (x), there are no spatial fluctuations:

The precise form of the critical exponent will depend on the nature of the correlations near the transition temperature. While this approximation neglects these correlations, the heat capacity does diverge at T c, which is also a feature of the liquid-vapor critical point, with a concomitant divergence in the Lyapunov time and degeneracy of the Lyapunov spectrum. 

Mean-field approximation: Applying a mean-field approximation to the original Hamiltonian [40] [41] [42], leads to an isochoric heat capacity with a finite-jump discontinuity (Supplementary Fig. 7 ). We replace the bilinear interaction between particles by -J 〈x〉 x, where 〈x〉 is the ensemble average displacement. The addition of J 〈x〉 2 /2 avoids double counting interactions. The Hamiltonian of a single particle is then

For our choice of A = κ 4 /4 and B = -κ 2 /2, the potential V(x) is a double well. We neglect the kinetic energy in what follows. The equilibrium free energy of each particle is

where ρ MF (x) = exp{-βV(x) + βJ〈x〉x} and Tr[⋅] is the configuration integral. The free energy is bistable for T < T c and monostable for T ≥ T c. Taking 〈x〉 to be a free parameter, the minimum of the free energy at

gives a constraint on the possible microstates consistent with this macrostate. Using the mean displacement as an order parameter, subject to the constraint above, the free energy is a minimum only when 〈x〉 satisfies the self-consistency condition

The minimum value M is the mean displacement of a particle at equilibrium. This order parameter is nonzero when T < T c (the ordered phase) and is zero when T ≥ T c (the disordered phase). The critical temperature is found by making contact with the Landau free energy for continuous phase transitions. At equilibrium, the free energy becomes F MF ðT; MÞ ¼ Àk B TlnTr½ρ MF ðxÞ þ 1 2 JM 2. We consider symmetric potential functions V(x) = V(-x) so that only even terms survive in the cumulant expansion ln e x h i 0 ¼

where 〈O〉 0 represents Tr[O(x)ρ 0 (x)]/Tr[ρ 0 (x)] and ρ 0 (x) = e -βV(x). Up to second order, the expanded free energy is

The second term vanishes at the critical temperature k B T c = J 〈x 2 〉 0 (Supplementary Fig. 6 ). 

Taylor expanding the coefficient of the second term in the free energy when T is near T c gives an expression that agrees with the Landau free energy to second order

This expression does not explicitly use the form of the potential V(x), only that it is symmetric 42. At a given β = 1/k B T, we self-consistently solve for M to find the optimal approximation of the original Hamiltonian. The canonical ensemble average of each observable O is Tr[Oρ MF (x)]/Tr[ρ MF (x)] with 〈x〉 = M. Assuming ergodicity, the ensemble average Lyapunov exponent λ ± = 〈λ ± 〉 is equal to the time average lim t!1 λ ± ðtÞ. 

The susceptibility diverges at the critical point (Supplementary Fig. 7 ). Following a standard procedure in Landau-Ginzburg theory 63, the order parameter is coupled to an external field, ξ. In this case, the minimization of the free energy FMF ðT; M; ξÞ ¼ F MF ðT; MÞ À ξM: ð20Þ determines the equilibrium value M. We see that ∂ FMF =∂M ¼ 0 when JM 1 À βJ x 2 h i 0 À Á ¼ ξ. From this expression, the inverse susceptibility is

and χ -1 → 0 as T → T c.
section_name: Methods

section_name: Acknowledgements

section_content: We acknowledge the donors of The American Chemical Society Petroleum Research Fund ( ACS PRF #55195-DNI6 ) for support of this research. We also acknowledge the use of the supercomputing facilities managed by the Research Computing Group at the University of Massachusetts Boston. We thank Bala Sundaram for comments on the manuscript and Anthony B. Costa for helpful discussions.
section_name: Acknowledgements

section_name: Data availability

section_content: The data are available from the corresponding author upon request. 

Received: 29 October 2018 Accepted: 11 April 2019
section_name: Data availability

section_name: Author contributions

section_content: M. D. and J. R. G. contributed to the design and implementation of the research, the analysis of the results, and the writing of the manuscript.
section_name: Author contributions

section_content: Supplementary Information accompanies this paper at https://doi. org/10. 1038/s41467-019-10040-3.
section_name: Additional information

section_content: The authors declare no competing interests. 

Reprints and permission information is available online at http://npg. nature. com/ reprintsandpermissions/ Journal peer review information: Nature Communications thanks Juan López and Andrea Rapisarda for their contribution to the peer review of this work. 

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
section_name: Competing interests: