Statistics of small length scale density fluctuations in supercooled viscous liquids

Small length scale density fluctuations in normal homogeneous liquids above the melting temperature obey Gaussian statistics over many orders of magnitude.^1,2 This dependence underlies many successful theories of the liquid state.^1,3–5 The Gaussian approximation is sometimes assumed to hold in the supercooled regime when the liquid approaches a glass transition.^6–12 This study investigates to what extent Gaussian statistics of small length scale density fluctuations persists in the supercooled viscous regime near the glass-transition. Viscous liquids are highly non-trivial as emphasized by the three non’s:¹⁰ non-exponential relaxation of equilibrium fluctuations, non-Arrhenius temperature dependence of structural relaxation time, and nonlinear out-of-equilibrium aging. Thus, it is not obvious that Gaussian statistics will persist in the supercooled viscous regime.

Statistics of density fluctuations add insight into the ongoing debate about the origin of slow dynamics in structural glass-formers.^11,13–15 In general, liquids can be cooled below their melting temperature due to the existence of a free-energy barrier related to the formation of a critical nucleus.¹⁶ The dynamics near the glass-transition is by definition dramatically slower than that near the melting temperature. When the dynamics is determined by a single constant free energy barrier, then the slowdown would follow the Arrhenius law. However, many liquids exhibit a super-Arrhenius slowdown (the first “non”). Experiments and computer studies of model liquids have shown that the dynamics becomes spatially heterogeneous upon cooling.^17–21 It is an appealing idea that the dynamical heterogeneity could be related to geometric arrangements of locally preferred structures of well-packed molecules or atoms. Several studies have identified accumulation of locally preferred structures in specific systems.^20–33 The resulting representation is a heterogeneous structure with non-Gaussian small length scale density fluctuations. A disadvantage of the “locally preferred structure” approach is that it is system specific. In this paper, it is proposed to study statistics in the collective density field. This general approach can be applied to distinct systems, which is shown by investigating the model belonging to chemically different classes.

In contrast to the aforementioned structural origin of the glass-transition, some explanations regard structure as less important.^10,34–36 The motivation lies on the experimental observation that the structure factor is similar in the normal and in the supercooled liquid regime. As an example, the quadratic scaling law of the temperature dependency of the relaxation time^35–38 originates from generic kinetic constraint models.^38,39 These models have trivial thermodynamics but non-trivial dynamics leading to a glass-transition. This concept suggests that statistics of small length scale density fluctuations near the glass transition inherit the Gaussian statistics of the normal liquid regime.

This study examines the statistics of small length scale density fluctuations for the single component Lennard-Jones (LJ)⁴⁰ model, the binary Kob-Andersen mixture (KABLJ),⁴¹ the Wahnström binary mixture (WABLJ),⁴² a coarse grained model of o-terphenyl (LWoTP),⁴³ and the TIP4P/Ice⁴⁴ water model. Enhanced sampling molecular dynamics methods are used to investigate the distributions over many orders of magnitude. The results show that the Gaussian approximation persists to a good degree in the supercooled viscous regime; however, deviations are larger at low temperatures. The analysis suggests that non-Gaussian features are related to first-order transitions to crystals (other possibilities are discussed in Sec. V).

This paper is organized as follows: Sec. II introduces the formalism used to describe density fluctuations and theory of the Gaussian approximation. Section III presents the numerical methods used for the enhanced sampling of the density field and provides descriptions of the investigated models. Section IV presents the results, and the implications are discussed in Sec. V.

This section presents the formalisms employed to characterize density fluctuations.

In experiments (X-ray or neutron scattering) and theories of the liquid states, it is often convenient to work in the reciprocal space. Consider a liquid of N particles located at R = {r₁, r₂, …, r_N} in a volume V with periodic boundaries so that the average density is ρ = N/V. Let the real-space density field be $ρ(r)=∑nNδ(rn−r)$⁠, where δ is Dirac’s delta function. The collective density field in the reciprocal space is then defined as the Fourier transform of the real-space density field: $ρk=∫Vdrρ(r)exp(−ik⋅r)/N$⁠, where $k=kk^$ is the scattering or wave vector (sometimes the letter “q” or “Q” is used). The $1/N$ factor ensures the system size invariance of amorphous configurations (liquids). The scaling is $N$ for configurations with long-range translational order (crystals) along the $k^$ direction. For a system confined in a periodic orthorhombic box, the density field can be written as

where k = (2πn_x/L_x, 2πn_y/L_y, 2πn_z/L_z), n_x, n_y, and n_z are integers, and L_x, L_y, and L_z are the length of the volume that confines the liquid (V = L_xL_yL_z).

The investigation can be limited to k-vectors along the x-direction without the loss of information due to liquid’s isotropy. For a cubic box with sides of length L = L_x = L_y = L_z, it is possible to consider k vectors of lengths k = 2πn/L, where n is an integer (n_x = n and n_y = n_z = 0). The isotropy is, in principle, destroyed by the box’s anisotropic periodic boundaries. However, such effects are expected to be small and are ignored in this study.

Here, we consider the probability distribution P(|ρ_k|), where |ρ_k| is the norm of the collective density field. It is not needed to consider the full two-dimensional complex plane of P(ρ_k) due to the symmetry of the liquid. The second moment of the distribution S_k = ⟨|ρ_k|²⟩ is the structure factor routinely measured in scattering experiments. If the statistics of density fluctuations follows a Gaussian (G) distribution, then

The central limit theorem states that in the thermodynamic limit, the density fluctuations become Gaussian (see also discussion in Sec. V): P(|ρ_k|) → P_G(|ρ_k|) for N → ∞. Thus, the analysis is limited to non-trivial small length scale fluctuations by studying systems of 100–1000 particles (unless otherwise stated). It has been demonstrated that a small system of that size represents viscous dynamics of a large system.⁴⁵ 

The fourth moment ⟨|ρ_k|⁴⟩ of the Gaussian distribution P_G(|ρ_k|) is $2⟨|ρk|2⟩2$⁠. Thus, a non-Gaussian parameter can be defined as

A consequence of the central limit theorem is that non-Gaussian features are more prominent in smaller systems. Thus, it can be illustrative to investigate density fluctuations in small subvolumes of a larger system. A subvolume can be defined through the function h(r) so that h is unity inside the volume and zero outside. The collective density field in this subvolume can then be written as $ρk(h)=∑nN⁡exp(−ik⋅rn)h(rn)/Vh,$ where V_h is the size of the subvolume. k = 0 relates to the average density ρ_h = N_h/V_h in the subvolume. Here, $Nh=∑nNh(rn)$ is the number of particles in the subvolume. The Gaussian approximation for the ρ_h density fluctuations is

where $m2=⟨(ρh−⟨ρh⟩)2⟩$ is the variance. For this distribution, the fourth central moment $m4=⟨(ρh−⟨ρh⟩)4⟩$ equals $3m22$⁠. Thus, a non-Gaussian parameter is defined as

A harmonic potential is added to the Hamiltonian in order to bias the system towards large |ρ_k| values. The Gaussian approximation can then either be investigated directly on statistics of the biased simulations or by re-weighting statistics of a series of simulations (a variation of the umbrellas sampling method⁴⁶). The description of the suggested method is as follows:

Let $H(R,Ṙ)$ be the Hamiltonian of a given system. To sample rare ρ_k fluctuations at some density and temperature, the following Hamiltonian is investigated with molecular dynamics:⁴⁷ 

where κ is the spring constant and a is the anchor point of the bias field. The |ρ_k| probability distribution of the non-biased system is given by

where P_κa(|ρ_k|) is the distribution of the Hamiltonian with the harmonic bias field applied and $Nκa$ needs to be determined numerically or derived from the Gaussian approximation.

For a series of overlapping distributions (with different a’s and κ’s), $Nκa$ can be determined numerically using the iterative Multistate Bennett Acceptance Ratio (MBAR) method.⁴⁸ Alternatively, the distribution function P_κa can be investigated directly. The prediction for the Gaussian approximation is found by combining Eqs. (2) and (7)

where the normalization constant is

$A=aκ′/Sk−1+κ′$⁠, κ′ = κ/2k_BT, and $erf(x)=2π−12∫0x⁡exp(−s2)ds$ is the error function. The first moment of the biased distribution is

where $B=Aπ12eA2(1+Sκ′(1+2a2κ′))$⁠. Predictions from the Gaussian approximation can be used as initial guesses for the iterative MBAR method.⁴⁸ In practice, this leads to fewer iterations before reaching convergence. (The unique solution that the MBAR method converges to can, in principle, be far away from the initial Gaussian guess.)

To perform molecular dynamic simulations, forces on particles from the bias field need to be evaluated. The total force acting on particle j is⁴⁷ 

where $Fj(0)$ is the force of the unbiased Hamiltonian $H(R,Ṙ)$ and

It is possible to design an N-scaling algorithm although the force acting on the particle j depends on the position of all the particles: First loop over all particles to compute ρ_k and then loop over all particles again to get particle forces using Eqs. (11) and (12). The algorithm can be parallelized to several processes since both the computation of ρ_k and the $Fj(κa)$ forces involve sums of independent contributions (assuming the same for $Fj(0)$⁠). Computational efficiency of the algorithm is crucial since long-time simulations are needed in the supercooled viscous regime.

The overall idea of the above mentioned method for computing rare fluctuations of the collective density field can also be used to sample rare density fluctuations in a subvolume. First define a continuous quantity $Ñh$ that is strongly correlated with the number of particles N_h inside the volume h. This is done by applying a switching function to the boundaries of the subvolume as described in Ref. 49 by Patel et al. The unbiased distribution P( ρ_h) is obtained by reweighing biased distribution P(⁠$Ñh$⁠) using the MBAR method.⁴⁸ For binary mixtures, a bias potential can be applied to both types of particles. Statistical information of the unbiased system is determined using the MBAR method.

This study investigates density fluctuations in several atomistic and models. The following examples have been chosen to represent different chemical classes of liquids:

• LJ:

  In 1924, Lennard-Jones suggested a simple model of interaction between atoms by summing a repulsive term representing Pauli repulsions and an attractive term representing London forces.⁴⁰ This study investigates a truncated version where the potential energy surface is $U(R)=∑n>mNu(|rm−rn|)$⁠, with u(r) = 4ϵ[(σ/r)¹² − (σ/r)⁶ − 2.5⁻¹² + 2.5⁻⁶] if r/σ < 2.5 and zero otherwise. Molecular dynamics is conducted for N = 108 particles with σ = ε = m = 1 in a periodic simulation cell at density ρ = 0.85 (L = 5.0273) and temperature T = 0.8. This state point is in the normal liquid regime and close to the freezing temperature.^47,51 The LJ model is not a good glass-former; however, it is included in the analysis since it is a standard system in computational condensed matter physics.

• KABLJ:

  Kob and Andersen suggested a binary LJ mixture as a model of a good glass former.⁴¹ This 80/20 mixture of AB particles has a strong affinity towards unlike atoms: ϵ_AA = σ_AA = 1, ϵ_AB = 1.5, σ_AA = 0.8, ϵ_BB = 0.5, and σ_BB = 0.88. This parametrization makes a good glass-former on time-scales, and the system sizes are typically investigated in silico. The mixture is the standard model for computational studies of low temperature liquid dynamics. It is custom to study the density ρ = 1.2, where the melting temperature is T_m = 1.027(3).⁵² Below this temperature, the particles will eventually phase separate and solidify in long-time simulations. The thermodynamically stable solid consists of two crystals: one, where A particles are in a face centered cubic structure and the other where A and B particles form the crystal structure of PuBr₃.⁵² If crystallization is avoided, however, the low-temperature liquid accumulates locally preferred structures where one of the small B particles is surrounded by ten larger A particles, thus forming a twisted bicapped square prism.^25,32

• WaBLJ:

  Wahnström suggested a 50/50 binary LJ mixture with a size ratio of 80%:⁴² ϵ_AA = ϵ_BB = σ_BB = 1, σ_AB = 1.1, and σ_BB = 1.2. Unlike the KABLJ mixture, the interaction parameters (ε’s and σ’s) follow the Lorentz-Berthelot rule of mixing. The system is a good glass former (in silico); however, in long-time simulations, the mixture may spontaneously crystallize into a MgZn₂ crystal structure.²⁷ In the supercooled regime, the liquid exhibits icosahedral^25,31,53 and Frank-Kasper order.²⁷ The former structure is an A particle with a first shell consisting of six A particles and six B particles forming a distorted icosahedron.⁵³ The latter structure is a geometric arrangement where two touching B particles have six smaller A particles as common neighbors.²⁷ The amount of these structures increases with the decrease in the temperature since they pack space well.^27,54

• LWoTP:

  Lewis and Wahnström⁴³ suggested a coarse grained model of ortho-terphenyl (C₁₈H₁₄) where molecules are constructed from three LJ particles placed at the corners of an isosceles triangle with two sides of length σ. The LJ particles are parametrized to represent benzene rings: ε = k_B(600 K) and σ = 0.483 nm. To avoid crystallization in a close-packed structure, the molecule has an inner angle of 75° (that is, between 60° and 90° degrees—the angles found between neighboring triplets in close packed structures). In long-time simulations, however, the system can crystallize into a structure where the LJ particles form a base centered cubic lattice with random orientations of molecules.^50,55 This study investigates a system of N = 324 molecules (unless otherwise stated) at temperature T = 350 K and density ρ = 1.09 g/ml (L = 4.84 nm).

• Water:

  Abascal et al.⁴⁴ introduced the TIP4P/Ice atomistic water model. This four site model reproduces the complicated phase-diagram of real water, suggesting that it could give a good representation of the hydrogen-bonds in the liquid state. The model is studied at temperature T = 280 K and density ρ = 1 g/ml. There are no signs of spontaneous crystallization.

Numerical computations are performed using the software packages LAMMPS,⁵⁶ RUMD,⁵⁷ and a home-made code available at http://urp.dk/tools. The implementation of the ρ_k bias field is available in the official LAMMPS and RUMD packages. Statistics of fluctuations of the ρ_k density field are investigated in the constant NVT ensemble. Statistics of the ρ_h density are studied in systems with a gas-liquid interface. This is obtained by constructing an elongated simulation cell with periodic boundaries in the x and y directions and walls at the boundaries of the longer z direction.

The results for the LJ, KABLJ, and WaBLJ models are reported in reduced Lennard-Jones units, while the physical units are used for the LWoTP model and the TIP4P/ice water model.

Figure 1 shows the natural logarithm of probability distributions P(|ρ_k|) of the LJ model in the region of the phase diagram just above the melting temperature (the normal liquid regime). Figure 1 includes k-vectors from lengths of k = 1.25 (n = 1) up to k = 15.0 (n = 12, shifted vertically for clarity). The solid black lines are the reweighed distributions obtained using a series of biased simulations, and the red dashed lines are Gaussian approximations. Overall the Gaussian approximation gives a good description of the data. This supports the common result that small length scales fluctuations are Gaussian in the normal liquid regime.^1,2 The tails of the distributions, however, show non-Gaussian features. As an example, the k-vectors with n = 6 and n = 8 show fat tails compared to the Gaussian reference. A representative configuration from the tail of the distribution for n = 6 is shown in Fig. 2(a). A crystalline structure can be identified in both the real-space configuration and the scattering spectrum as shown in Fig. 2(b). A cubic box with 108 LJ particles have an ideal crystal structure with 3 × 3 × 3 fcc unit cells giving a Bragg peak at n = 6. The distribution of the n = 8 vector [Figs. 2(c) and 2(d)] also has a fat tail. This is also attributed to a crystalline structure but with another orientation. (Bias simulation similar to the ones presented here can be used to compute the melting point of crystals. This is referred to as the “interface pinning” method.⁴⁷)

Other k-vectors have thin tails. For example, Fig. 2(e) shows a configuration from the tail of the n = 10 wave vector. This structure is not crystalline but amorphous. Figures 2(g) and 2(h) show a structure from the longest wave vector investigated (n = 1). The liquid responds to a strong bias field by forming a vapor slap and a crystalline slap.

Next the glass forming models are investigated. Figure 3 shows the data for the KABLJ model at temperature T = 0.45 ( ρ = 1.2) for a system size of N = 1000 particles. This state-point is well below the melting temperature T_m = 1.027(3)⁵² where the structural relaxation time is about 10³ times slower than the normal liquid regime.^41,58 The solid black lines in Fig. 3 show biased distributions P_κa(ρ_k) for κ = 4 and several anchor points: a = {0, 1, 2, 3, 4, 5}. The red dashed lines are Gaussian approximation. The agreement is good.

Interestingly, the systems are more prone to crystallization when a |ρ_k| bias field is applied. Figure 4 shows examples of crystallizing trajectories of the KABLJ mixture and LWoTP trimers. Crystallizing trajectories are discarded in the analysis of density fluctuations of the liquid state.

Figure 5 shows the distributions P(|ρ_k|) for the glass forming liquids KABLJ, WaBLJ, LWoTP, and water for several k-vectors. The red dashed lines indicate the Gaussian approximation. The agreement is good, but the tails of the distributions deviate from the Gaussian prediction. The deviations are system-size dependent as expected from the central limit theorem. Figure 6 shows that the non-Gaussian fat tail for k = 1.4 nm of the LWoTP systems is greatly diminished when the system size increases from N = 324 to N = 2592.

Fluctuations in small subvolumes of a larger system can give a further insight into the structural origin of the non-Gaussian small length scale density fluctuations. Figure 7 shows the distribution function of the ρ_h density in a 3 × 3 × 3 subvolume, h. The points are reweighed statistics from simulations with bias potentials. The red dashed line is the prediction from the Gaussian approximation [Eq. (4)]. The agreement is good and is consistent with the findings in Refs. 1 and 2. Some deviations from the Gaussian approximation can be observed in the tails of the distributions. The low density limit corresponds to the formation of a cubic vapor bubble in the liquid. As described by the classical nucleation theory, the free energy −k_BT ln(P) required to form such a bubble is the sum of bulk and surface contributions.

Gas-liquid coexistence simulations of the KABLJ mixture are performed to investigate density fluctuations in the supercooled viscous regime. Figure 8 shows the structural relaxation time of particles in the liquid slap. The relaxation time is non-Arrhenius in the investigated temperature regime, 0.28 < T < 0.60. The points on Fig. 9 shows the distribution of the density ρ_h fluctuations in a 3 × 3 × 3 subvolume h for temperatures T = {0.35, 0.40, 0.55}. The Gaussian approximations [Eq. (4)] are shown as red dashed lines. The conclusion from the analysis of the |ρ_k| fluctuations remains—the Gaussian approximation provides a fair description but becomes less accurate when temperature is lowered. Deviations are seen in both the tails of the distributions and near the mean as quantified by the non-Gaussian parameter $αρh$ [Eq. (5)]. Figure 10(a) shows the joint distribution of the two types of particles, P(N_A, N_B). The Gaussian approximation predicts elliptical shaped contour lines. However, the distribution shows deviations in the tails. Figure 10(b) shows a configuration from a fat-tail region of the distribution at equimolar composition. The configuration is a cubic CsCl crystallite. This structure is one of the thermodynamically stable crystal structures of the KABLJ model.⁵² This suggests that non-Gaussian features in the tails are related to the first order-transition from a liquid to a crystal.

The results for chemically different glass formers have been presented. The overall conclusion is that the Gaussian approximation provides a fair description of the small length scale density fluctuations; however, as temperature is lowered, the Gaussian approximation becomes less accurate.

Gaussian statistics is usually described in two ways: (i) from the central limit theorem, or (ii) from a harmonic approximation. (i) The central limit theorem states that if random variables from any underlying distribution are added together, the resulting distribution will tend to follow Gaussian statistics. For a non-flowing equilibrium liquid of sufficient size, it can be assumed that subvolumes fluctuate independently. As an example, consider the ideal gas model. The number of particles in a given subvolume h follows the Poisson distribution. If the size of the subvolume is increased, then Raikov’s theorem states that the fluctuations again follow a Poisson distribution with a larger average value. This distribution in the thermodynamic limit becomes the Gauss distribution. For a liquid with interactions, the underlying distribution differs from the Poisson distribution; however, by studying small length scale density fluctuations, one can get insight into the nontrivial underlying distribution. This brings us to the other, less trivial, way of arriving at Gaussian statistics. That is, (ii) by a harmonic expansion around a local minimum of the free energy function F: if x is an order parameter, such as ρ_k or ρ_h, then the free energy is

where P(x) is the probability distribution. The function F(x) can be expressed as an infinite polynomial expansion around the minimum at x₀. It is often convenient to approximate F and assume that only the second order term is of relevance. Thus, if we ignore higher order terms, we get a harmonic approximation for the free energy⁵ 

The truncation of the expansion series is non-trivial, as discussed below. The probability distribution of the Gaussian approximation is found by equating Eqs. (13) and (14)

(The reason that statistics is found to be near-Gaussian is not due to the harmonic bias field added to the Hamiltonian since the distribution functions are reweighed using the MBAR method.) To understand the first order transitions, e.g., the gas-liquid transition, higher order terms are relevant. As a classic example, Landau’s (L) effective Hamiltonian^59,60 includes higher-order terms to give a description of the density fluctuations near the gas-liquid critical point: $FL(x)=a2[x−x0]2+a4[x−x0]4+const.$⁠. With Landau’s theory in mind, one expects deviations from the Gaussian approximation when other phases (gas or crystals) interfere with the liquid state. Deviations from the Gaussian statistics can be attributed to the formation of a vapor bubble, as shown in Fig. 2. In the supercooled regime, the crystal basin in the free energy landscape becomes large, suggesting that statistics becomes less Gaussian due to the presence of a crystal. The microscopic image is the formation of subcritical crystallites (Fig. 6). In agreement with this, the systems are more prone to crystallization when a strong field biasing |ρ_k| is applied (Fig. 4).

Beside subcritical crystallites, some supercooled viscous liquid may accumulate crystal-like structures that share the local arrangement of the crystal, i.e., rotational order in the first shell. Tanaka and co-workers^22,23,29,33 have demonstrated that for some model liquids, there is a link between slow dynamics and crystal-like order. Interestingly, the locally preferred structures of the WaBLJ supercooled viscous liquid (the distorted icosahedron and the Frank-Kasper bond) are shared with the crystal (MgZn₂).²⁷ This contrasts the KABLJ mixture where the locally preferred structure^25,32 is different from the crystalline ground state.⁵² In general, non-crystalline structures could also be important for the statistics of small length scale fluctuations. Non-crystalline structures have been suggested as an important component for understanding the dynamics of the highly viscous liquid near the glass transition.^{21,25,26,30–32} More refined methods are needed to clarify the coupling between a specific local structure and small length scale density fluctuations. The joint distribution function of density ( ρ_k or ρ_h) and other order-parameters can give insight into this (see, e.g., Fig. 3 in Ref. 33). Such investigations are left to future studies.

Small length scale density fluctuations can also be analyzed from an energy landscape^45,61–63 perspective. In this perspective, the 3N dimensional energy surface of the liquid is partitioned into basins identified by the local minima. Below a certain onset temperature, the system explores the configuration space by means of two mechanisms. At short times, the system vibrates in a basin that is, to a good approximation, harmonic. Thus, it is expected that these vibrations will give rise to Gaussian statistics of density fluctuations. On longer time-scales, the system will explore different basing (activated relaxation). From this perspective, the non-Gaussian feature at low temperatures is related to density fluctuations between basins (the inherent states).

Some theories directly or indirectly assume Gaussian statistics of small length scale density fluctuations. As an example, the fact that small length scale fluctuations persist to be near-Gaussian in the supercooled regime is in line with the picture behind the generic kinetic constraint models.^35–39 As mentioned in the Introduction, this concept suggest that thermodynamics and structural details of glass forming materials are not crucial for understanding the dynamics of highly viscous liquids. Finally, the Gaussian approximation of density fluctuations plays a role in some elastic models.^10,64,65 The idea behind these approaches is that the viscous dynamics is connected with elastic deformations which allow particles in a subvolume to rearrange. Here, the Gaussian approximation enters in some theories to make predictions that can be related to experiments.

This study was initiated by discussions with David Chandler. I joined his group in 2010, and greatly benefited from interactions with Patrick Varilly, Amish J. Patel, Thomas Speck, David Limmer, Lester O. Hedges, Yael Elmatad, and Adam P. Willard. David had a remarkable talent of challenging conventional beliefs and thereby moving the scientific field forward. I moreover received comments and suggestions from Lorenzo Costigliola, Andreas Tröster, Jeppe C. Dyre, and Thomas B. Schrøder. This work was supported by the VILLUM Foundation’s Matter under Grant No. 16515.