Many successful theories of liquids near the melting temperature assume that small length scale density fluctuations follow Gaussian statistics. This paper presents a numerical investigation of density fluctuations in the supercooled viscous regime using an enhanced sampling method. Five model systems are investigated: the single component Lennard-Jones liquid, the Kob-Andersen binary mixture, the Wahnström binary mixture, the Lewis-Wahnström model of ortho-terphenyl, and the TIP4P/Ice model of water. The results show that the Gaussian approximation persists to a good degree into the supercooled viscous regime; however, it is less accurate at low temperatures. The analysis suggests that non-Gaussian fluctuations are related to crystalline configurations. Implications for theories of the glass transition are discussed.

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:10 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.16 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 time35–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)40 model, the binary Kob-Andersen mixture (KABLJ),41 the Wahnström binary mixture (WABLJ),42 a coarse grained model of o-terphenyl (LWoTP),43 and the TIP4P/Ice44 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 = {r1, r2, …, rN} in a volume V with periodic boundaries so that the average density is ρ = N/V. Let the real-space density field be ρ(r)=nNδ(rnr), 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(ikr)/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

ρk=n=1Nexp(ikrn)/N,
(1)

where k = (2πnx/Lx, 2πny/Ly, 2πnz/Lz), nx, ny, and nz are integers, and Lx, Ly, and Lz are the length of the volume that confines the liquid (V = LxLyLz).

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 = Lx = Ly = Lz, it is possible to consider k vectors of lengths k = 2πn/L, where n is an integer (nx = n and ny = nz = 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 Sk = ⟨|ρk|2⟩ is the structure factor routinely measured in scattering experiments. If the statistics of density fluctuations follows a Gaussian (G) distribution, then

PG(|ρk|)=2|ρk|exp(|ρk|2/Sk)/Sk.
(2)

The central limit theorem states that in the thermodynamic limit, the density fluctuations become Gaussian (see also discussion in Sec. V): P(|ρk|) → PG(|ρ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.45 

The fourth moment ⟨|ρk|4⟩ of the Gaussian distribution PG(|ρk|) is 2|ρk|22. Thus, a non-Gaussian parameter can be defined as

αρk=|ρk|4/2|ρk|221.
(3)

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)=nNexp(ikrn)h(rn)/Vh, where Vh is the size of the subvolume. k = 0 relates to the average density ρh = Nh/Vh in the subvolume. Here, Nh=nNh(rn) is the number of particles in the subvolume. The Gaussian approximation for the ρh density fluctuations is

PG(ρh)=exp([ρhρh]2/2m2)/2πm2,
(4)

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

αρh=m4/3m221.
(5)

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 method46). 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:47 

Hκa(R,Ṙ)=H(R,Ṙ)+κ[|ρk(R)|a]2/2,
(6)

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

P(|ρk|)=NκaPκa(|ρk|)exp(κ[|ρk|a]2/2kBT),
(7)

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.48 Alternatively, the distribution function Pκa can be investigated directly. The prediction for the Gaussian approximation is found by combining Eqs. (2) and (7)

PG(κa)(|ρk|)=2|ρk|exp(|ρk|2/Skκ[|ρk|a]2)/NG  (κa)Sk,
(8)

where the normalization constant is

NG  (κa)=exp(A2Skκ)1+Skκexp(A2)+Aπ12[1+erf(A)],
(9)

A=aκ/Sk1+κ, κ′ = κ/2kBT, and erf(x)=2π120xexp(s2)ds is the error function. The first moment of the biased distribution is

|ρk|PG(κa)=exp(a2κ)N2aκ[1+Skκ]2[2Sk[aκ]2+B[1+erf(A)]],
(10)

where B=Aπ12eA2(1+Sκ(1+2a2κ)). Predictions from the Gaussian approximation can be used as initial guesses for the iterative MBAR method.48 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 is47 

Fj(κa)=Fj(0)κ(|ρk|a)j|ρk|,
(11)

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

j|ρk|=kR[iρk*exp(ikrj)]/|ρk|.
(12)

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 Nh 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.48 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.40 This study investigates a truncated version where the potential energy surface is U(R)=n>mNu(|rmrn|), with u(r) = 4ϵ[(σ/r)12 − (σ/r)6 − 2.5−12 + 2.5−6] 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.41 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 Tm = 1.027(3).52 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 PuBr3.52 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%:42ϵ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 MgZn2 crystal structure.27 In the supercooled regime, the liquid exhibits icosahedral25,31,53 and Frank-Kasper order.27 The former structure is an A particle with a first shell consisting of six A particles and six B particles forming a distorted icosahedron.53 The latter structure is a geometric arrangement where two touching B particles have six smaller A particles as common neighbors.27 The amount of these structures increases with the decrease in the temperature since they pack space well.27,54

  • LWoTP:

    Lewis and Wahnström43 suggested a coarse grained model of ortho-terphenyl (C18H14) 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: ε = kB(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.44 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,56 RUMD,57 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.47)

FIG. 1.

Natural logarithm of the probability distribution P(|ρk|) for k = (2πn/5.0273, 0, 0) from simulations of the single component LJ model in the normal liquid regime. The solid black lines are the distribution function computed from reweighed biased simulations. The red dashed lines are the Gaussian approximations [Eq. (2)]. The distributions have been shifted vertically for clarity. The inset shows the structure factor, where dots indicate the investigated k-vectors.

FIG. 1.

Natural logarithm of the probability distribution P(|ρk|) for k = (2πn/5.0273, 0, 0) from simulations of the single component LJ model in the normal liquid regime. The solid black lines are the distribution function computed from reweighed biased simulations. The red dashed lines are the Gaussian approximations [Eq. (2)]. The distributions have been shifted vertically for clarity. The inset shows the structure factor, where dots indicate the investigated k-vectors.

Close modal
FIG. 2.

Configurations representing the system in the non-Gaussian tails of the distributions shown in Fig. 1. The left panels show configuration of the single component LJ model, and the right panels show the corresponding scattering spectra Sk in the xy-plane (average over several configurations). From top to bottom, the panels show data from simulations with wave vectors with n = 6, n = 8, n = 10, and n = 1, respectively.

FIG. 2.

Configurations representing the system in the non-Gaussian tails of the distributions shown in Fig. 1. The left panels show configuration of the single component LJ model, and the right panels show the corresponding scattering spectra Sk in the xy-plane (average over several configurations). From top to bottom, the panels show data from simulations with wave vectors with n = 6, n = 8, n = 10, and n = 1, respectively.

Close modal

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 Tm = 1.027(3)52 where the structural relaxation time is about 103 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.

FIG. 3.

The distribution function Pκa(|ρk|) in a series of |ρk| biased simulations (κ = 4 and a = {0, 1, 2, 3, 4, 5}) of the KABLJ liquid in the supercooled regime. The red dashed lines show the Gaussian prediction [Eq. (8); no fitting parameters]. The agreement is good. This type of data has been reweighed to get the P(|ρk|) distributions in Fig. 5.

FIG. 3.

The distribution function Pκa(|ρk|) in a series of |ρk| biased simulations (κ = 4 and a = {0, 1, 2, 3, 4, 5}) of the KABLJ liquid in the supercooled regime. The red dashed lines show the Gaussian prediction [Eq. (8); no fitting parameters]. The agreement is good. This type of data has been reweighed to get the P(|ρk|) distributions in Fig. 5.

Close modal

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.

FIG. 4.

Examples of crystallizing trajectories: (a) the |ρk| trajectory in a simulation with the bias field 5(|ρk|6.5)2 added to the Hamiltonian of the KABLJ mixture. A crystallite is formed in the last third of the simulation. The crystallite consists of pure A particles, and thus, crystallization is accompanied by a phase separation. (b) The number of A’s that have 12 A’s in the first neighbor shell. This order-parameter is an indicator of the crystallization. (c) The configuration in the last step of the simulation. The A particles are green colored, while the B’s are red colored. (d) The crystalline qubatic order-parameter50 and (e) the potential energy in biased simulations of the LWoTP model. (f) The final configuration of a trajectory where a crystal is formed. Molecules are given separate colors.

FIG. 4.

Examples of crystallizing trajectories: (a) the |ρk| trajectory in a simulation with the bias field 5(|ρk|6.5)2 added to the Hamiltonian of the KABLJ mixture. A crystallite is formed in the last third of the simulation. The crystallite consists of pure A particles, and thus, crystallization is accompanied by a phase separation. (b) The number of A’s that have 12 A’s in the first neighbor shell. This order-parameter is an indicator of the crystallization. (c) The configuration in the last step of the simulation. The A particles are green colored, while the B’s are red colored. (d) The crystalline qubatic order-parameter50 and (e) the potential energy in biased simulations of the LWoTP model. (f) The final configuration of a trajectory where a crystal is formed. Molecules are given separate colors.

Close modal

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.

FIG. 5.

Natural logarithm of probability distribution P(|ρk|) for (a) the KABLJ mixture (T = 0.45, ρ = 1.2, and N = 1000), (b) the WaBLJ mixture (T = 0.6, ρ = 0.75, and N = 1024), (c) LWoTP trimer molecules (T = 350 K, ρ = 1.09 g/ml, and N = 324), and (d) tip4p/ice water model (T = 280 K, ρ = 1 g/ml, and N = 432). The red dashed lines are Gaussian predictions. The insets show scattering functions. The agreement with the Gaussian prediction is good; however, deviations are apparent in the tails of the distributions. These are likely attributed to crystalline structures.

FIG. 5.

Natural logarithm of probability distribution P(|ρk|) for (a) the KABLJ mixture (T = 0.45, ρ = 1.2, and N = 1000), (b) the WaBLJ mixture (T = 0.6, ρ = 0.75, and N = 1024), (c) LWoTP trimer molecules (T = 350 K, ρ = 1.09 g/ml, and N = 324), and (d) tip4p/ice water model (T = 280 K, ρ = 1 g/ml, and N = 432). The red dashed lines are Gaussian predictions. The insets show scattering functions. The agreement with the Gaussian prediction is good; however, deviations are apparent in the tails of the distributions. These are likely attributed to crystalline structures.

Close modal
FIG. 6.

The natural logarithm of the distribution P(|ρk|) with k = 1.4 nm−1 of the LWoTP model for system sizes of N = 324 (blue) and N = 2592 (green) molecules, respectively. The Gaussian approximation (red dashed) is more accurate for the larger system as expected from the central limit theorem. The non-Gaussian parameters αρk [Eq. (3)] are 0.029 and 0.0011 for N = 324 and N = 2592, respectively.

FIG. 6.

The natural logarithm of the distribution P(|ρk|) with k = 1.4 nm−1 of the LWoTP model for system sizes of N = 324 (blue) and N = 2592 (green) molecules, respectively. The Gaussian approximation (red dashed) is more accurate for the larger system as expected from the central limit theorem. The non-Gaussian parameters αρk [Eq. (3)] are 0.029 and 0.0011 for N = 324 and N = 2592, respectively.

Close modal

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 −kBT ln(P) required to form such a bubble is the sum of bulk and surface contributions.

FIG. 7.

Density fluctuations in a 3 × 3 × 3 subvolume, h, of the LJ model in the normal liquid regime with a gas-liquid interface (T = 0.7, N = 3000, Lx = Ly = 10). The inset shows a typical configuration of the system with a gas-liquid interface, and the subvolume h located in the bulk liquid part.

FIG. 7.

Density fluctuations in a 3 × 3 × 3 subvolume, h, of the LJ model in the normal liquid regime with a gas-liquid interface (T = 0.7, N = 3000, Lx = Ly = 10). The inset shows a typical configuration of the system with a gas-liquid interface, and the subvolume h located in the bulk liquid part.

Close modal

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(NA, NB). 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.52 This suggests that non-Gaussian features in the tails are related to the first order-transition from a liquid to a crystal.

FIG. 8.

Structural relaxation time τα as a function of temperature of the KABLJ mixture with a gas-liquid interface (N = 3000, Lx = Ly = 10). The structural relaxation time is here defined as Fs(k = 2π, t = τα) = 1/e, where Fs is the self intermediate scatter function of A’s located inside the slap (−5 < z < 5). The inset shows the density of the liquid slap as a function of temperature.

FIG. 8.

Structural relaxation time τα as a function of temperature of the KABLJ mixture with a gas-liquid interface (N = 3000, Lx = Ly = 10). The structural relaxation time is here defined as Fs(k = 2π, t = τα) = 1/e, where Fs is the self intermediate scatter function of A’s located inside the slap (−5 < z < 5). The inset shows the density of the liquid slap as a function of temperature.

Close modal
FIG. 9.

Density fluctuations at different temperatures in a 3 × 3 × 3 subvolume of the KABLJ mixture with a gas-liquid interface. The Gaussian approximation becomes better at higher temperature. The non-Gaussian parameter [Eq. (5)] for the high temperature T = 0.55 is αρh=3×104 while it is significantly larger, αρh=4×102, for the two lower temperatures, T = 0.40, 0.35.

FIG. 9.

Density fluctuations at different temperatures in a 3 × 3 × 3 subvolume of the KABLJ mixture with a gas-liquid interface. The Gaussian approximation becomes better at higher temperature. The non-Gaussian parameter [Eq. (5)] for the high temperature T = 0.55 is αρh=3×104 while it is significantly larger, αρh=4×102, for the two lower temperatures, T = 0.40, 0.35.

Close modal
FIG. 10.

(a) The ln(P(NA, NB)) distribution in a 3 × 3 × 3 subvolume of the KABLJ mixture at T = 0.325 (the white squares were not computed due to bad statistics). Non-gaussian features are seen as the contour lines that deviate slightly from being oval. (b) A configuration from the tail of the distribution with equimolar composition. The upper half of particles has been made invisible to reveal the arrangement of particles in the 3 × 3 × 3 subvolume h. The structure corresponds to a cubic CsCl crystallite. This is one of the thermodynamically stable phases of the mixture.52 

FIG. 10.

(a) The ln(P(NA, NB)) distribution in a 3 × 3 × 3 subvolume of the KABLJ mixture at T = 0.325 (the white squares were not computed due to bad statistics). Non-gaussian features are seen as the contour lines that deviate slightly from being oval. (b) A configuration from the tail of the distribution with equimolar composition. The upper half of particles has been made invisible to reveal the arrangement of particles in the 3 × 3 × 3 subvolume h. The structure corresponds to a cubic CsCl crystallite. This is one of the thermodynamically stable phases of the mixture.52 

Close modal

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

F(x)=kBTln(P(x)),
(13)

where P(x) is the probability distribution. The function F(x) can be expressed as an infinite polynomial expansion around the minimum at x0. 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 energy5 

FG(x)=a2[xx0]2+const.
(14)

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)

PG(x)=exp(a2[xx0]2/kBT).
(15)

(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 Hamiltonian59,60 includes higher-order terms to give a description of the density fluctuations near the gas-liquid critical point: FL(x)=a2[xx0]2+a4[xx0]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-workers22,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 (MgZn2).27 This contrasts the KABLJ mixture where the locally preferred structure25,32 is different from the crystalline ground state.52 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 landscape45,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.

1.
G.
Hummer
,
S.
Garde
,
A. E.
García
,
A.
Pohorille
, and
L. R.
Pratt
, “
An information theory model of hydrophobic interactions
,”
Proc. Natl. Acad. Sci. U. S. A.
93
,
8951
8955
(
1996
).
2.
G. E.
Crooks
and
D.
Chandler
, “
Gaussian statistics of the hard-sphere fluid
,”
Phys. Rev. E
56
,
4217
4221
(
1997
).
3.
L. R.
Pratt
and
D.
Chandler
, “
Theory of the hydrophobic effect
,”
J. Chem. Phys.
67
,
3683
3704
(
1977
).
4.
D.
Chandler
, “
Gaussian field model of fluids with an application to polymeric fluids
,”
Phys. Rev. E
48
,
2898
2905
(
1993
).
5.
M.
Kardar
,
Statistical Physics of Fields
(
Cambridge University Press
,
2007
).
6.
C. A.
Angell
, “
Formation of glasses from liquids and biopolymers
,”
Science
267
,
1924
1935
(
1995
).
7.
E.
Zaccarelli
,
G.
Foffi
,
F.
Sciortino
,
P.
Tartaglia
, and
K. A.
Dawson
, “
Gaussian density fluctuations and mode coupling theory for supercooled liquids
,”
Europhys. Lett.
55
,
157
(
2001
).
8.
P. G.
Debenedetti
and
F. H.
Stillinger
, “
Supercooled liquids and the glass transition
,”
Nature
410
,
259
(
2001
).
9.
J. C.
Dyre
, “
Solidity of viscous liquids. IV. Density fluctuations
,”
Phys. Rev. E
74
,
021502
(
2006
).
10.
J. C.
Dyre
, “
Colloquium: The glass transition and elastic models of glass-forming liquids
,”
Rev. Mod. Phys.
78
,
953
972
(
2006
).
11.
M. D.
Ediger
and
P.
Harrowell
, “
Perspective: Supercooled liquids and glasses
,”
J. Chem. Phys.
137
,
080901
(
2012
).
12.
G.
Biroli
and
J. P.
Garrahan
, “
Perspective: The glass transition
,”
J. Chem. Phys.
138
,
12A301
(
2013
).
13.
L.
Berthier
and
G.
Biroli
, “
Theoretical perspective on the glass transition and amorphous materials
,”
Rev. Mod. Phys.
83
,
587
645
(
2011
).
14.
L.
Berthier
,
G.
Biroli
,
J.-P.
Bouchaud
,
L.
Cipelletti
,
D. E.
Masri
,
D.
L’Hôte
,
F.
Ladieu
, and
M.
Pierno
, “
Direct experimental evidence of a growing length scale accompanying the glass transition
,”
Science
310
,
1797
1800
(
2005
).
15.
J.
Schmelzer
and
T.
Tropin
, “
Glass transition, crystallization of glass-forming melts, and entropy
,”
Entropy
20
,
103
(
2018
).
16.
R.
Becker
and
W.
Döring
, “
Kinetische behandlung der keimbildung in übersättigten dämpfen
,”
Ann. Phys.
416
,
719
752
(
1935
).
17.
E.
Donth
, “
The size of cooperatively rearranging regions at the glass transition
,”
J. Non-Cryst. Solids
53
,
325
330
(
1982
).
18.
K.
Schmidt-Rohr
and
H. W.
Spiess
, “
Nature of nonexponential loss of correlation above the glass transition investigated by multidimensional NMR
,”
Phys. Rev. Lett.
66
,
3020
3023
(
1991
).
19.
W.
Kob
,
C.
Donati
,
S. J.
Plimpton
,
P. H.
Poole
, and
S. C.
Glotzer
, “
Dynamical heterogeneities in a supercooled Lennard-Jones liquid
,”
Phys. Rev. Lett.
79
,
2827
2830
(
1997
).
20.
A.
Widmer-Cooper
and
P.
Harrowell
, “
On the relationship between structure and dynamics in a supercooled liquid
,”
J. Phys.: Condens. Matter
17
,
S4025
(
2005
).
21.
J. E.
Hallett
,
F.
Turci
, and
C. P.
Royall
, “
Local structure in deeply supercooled liquids exhibits growing lengthscales and dynamical correlations
,”
Nat. Commun.
9
,
3272
(
2018
).
22.
H.
Shintani
and
H.
Tanaka
, “
Frustration on the way to crystallization in glass
,”
Nat. Phys.
2
,
200
(
2006
).
23.
H.
Tanaka
, “
A simple physical model of liquid-glass transition: Intrinsic fluctuating interactions and random fields hidden in glass-forming liquids
,”
J. Phys.: Condens. Matter
10
,
L207
L214
(
1998
).
24.
M. D.
Ediger
, “
Spatially heterogeneous dynamics in supercooled liquids
,”
Annu. Rev. Phys. Chem.
51
,
99
128
(
2000
).
25.
D.
Coslovich
and
G.
Pastore
, “
Understanding fragility in supercooled Lennard-Jones mixtures. I. Locally preferred structures
,”
J. Chem. Phys.
127
,
124504
(
2007
).
26.
C. P.
Royall
,
R. W.
Stephen
,
O.
Takehiro
, and
T.
Hajime
, “
Direct observation of a local structural mechanism for dynamic arrest
,”
Nat. Mater.
7
,
556
(
2008
).
27.
U. R.
Pedersen
,
T. B.
Schrøder
,
J. C.
Dyre
, and
P.
Harrowell
, “
Geometry of slow structural fluctuations in a supercooled binary alloy
,”
Phys. Rev. Lett.
104
,
105701
(
2010
).
28.
I.
Gallino
,
J.
Schroers
, and
R.
Busch
, “
Kinetic and thermodynamic studies of the fragility of bulk metallic glass forming liquids
,”
J. Appl. Phys.
108
,
063501
(
2010
).
29.
H.
Tanaka
, “
Bond orientational order in liquids: Towards a unified description of water-like anomalies, liquid-liquid transition, glass transition, and crystallization–bond orientational order in liquids
,”
Eur. Phys. J. E
35
,
113
(
2012
).
30.
R.
Pastore
,
A.
Coniglio
, and
M. P.
Ciamarra
, “
Dynamic phase coexistence in glass–forming liquids
,”
Sci. Rep.
5
,
11770
(
2015
).
31.
F.
Turci
,
T.
Speck
, and
C. P.
Royall
, “
Structural-dynamical transition in the Wahnström mixture
,”
Eur. Phys. J. E
41
,
54
(
2018
).
32.
F.
Turci
,
C. P.
Royall
, and
T.
Speck
, “
Devitrification of the Kob-Andersen glass former: Competition with the locally favored structure
,” e-print arXiv:1808.01791 (
2018
).
33.
J.
Russo
and
H.
Tanaka
, “
Crystal nucleation as the ordering of multiple order parameters
,”
J. Chem. Phys.
145
,
211801
(
2016
).
34.
D.
Turnbull
and
M. H.
Cohen
, “
Free-volume model of the amorphous phase: Glass transition
,”
J. Chem. Phys.
34
,
120
125
(
1961
).
35.
L. O.
Hedges
,
R. L.
Jack
,
J. P.
Garrahan
, and
D.
Chandler
, “
Dynamic order-disorder in atomistic models of structural glass formers
,”
Science
323
,
1309
1313
(
2009
).
36.
D.
Chandler
and
J. P.
Garrahan
, “
Dynamics on the way to forming glass: Bubbles in space-time
,”
Annu. Rev. Phys. Chem.
61
,
191
217
(
2010
).
37.
J. P.
Garrahan
and
D.
Chandler
, “
Geometrical explanation and scaling of dynamical heterogeneities in glass forming systems
,”
Phys. Rev. Lett.
89
,
035704
(
2002
).
38.
F.
Ritort
and
P.
Sollich
, “
Glassy dynamics of kinetically constrained models
,”
Adv. Phys.
52
,
219
342
(
2003
).
39.
G. H.
Fredrickson
and
H. C.
Andersen
, “
Kinetic Ising model of the glass transition
,”
Phys. Rev. Lett.
53
,
1244
1247
(
1984
).
40.
J. E.
Lennard-Jones
, “
On the determination of molecular fields. I. From the variation of the viscosity of a gas with temperature
,”
Proc. R. Soc. London, Ser. A
106
,
441
462
(
1924
).
41.
W.
Kob
and
H. C.
Andersen
, “
Scaling behavior in the β-relaxation regime of a supercooled Lennard-Jones mixture
,”
Phys. Rev. Lett.
73
,
1376
1379
(
1994
).
42.
G.
Wahnström
, “
Molecular-dynamics study of a supercooled two-component Lennard-Jones system
,”
Phys. Rev. A
44
,
3752
3764
(
1991
).
43.
L. J.
Lewis
and
G.
Wahnström
, “
Relaxation of a molecular glass at intermediate times
,”
Solid State Commun.
86
,
295
299
(
1993
).
44.
J. L. F.
Abascal
,
E.
Sanz
,
R. G.
Fernández
, and
C.
Vega
, “
A potential model for the study of ices and amorphous water: TIP4P/Ice
,”
J. Chem. Phys.
122
,
234511
(
2005
).
45.
A.
Heuer
, “
Exploring the potential energy landscape of glass-forming systems: From inherent structures via metabasins to macroscopic transport
,”
J. Phys.: Condens. Matter
20
,
373101
(
2008
).
46.
D.
Frenkel
and
B.
Smit
, in
Understanding Molecular Simulation: From Algorithms to Applications
, Volume 1 of Computational Science Series, 2nd ed., edited by
D.
Frenkel
,
M.
Klein
,
M.
Parrinello
, and
B.
Smit
(
Academic Press
,
2002
).
47.
U. R.
Pedersen
, “
Direct calculation of the solid-liquid Gibbs free energy difference in a single equilibrium simulation
,”
J. Chem. Phys.
139
,
104102
(
2013
).
48.
M. R.
Shirts
and
J. D.
Chodera
, “
Statistically optimal analysis of samples from multiple equilibrium states
,”
J. Chem. Phys.
129
,
124105
(
2008
).
49.
A. J.
Patel
,
P.
Varilly
, and
D.
Chandler
, “
Fluctuations of water near extended hydrophobic and hydrophilic surfaces
,”
J. Phys. Chem. B
114
,
1632
1637
(
2010
).
50.
U. R.
Pedersen
and
P.
Harrowell
, “
Factors contributing to the glass-forming ability of a simulated molecular liquid
,”
J. Phys. Chem. B
115
,
14205
14209
(
2011
).
51.
J.-P.
Hansen
and
L.
Verlet
, “
Phase transitions of the Lennard-Jones system
,”
Phys. Rev.
184
,
151
161
(
1969
).
52.
U. R.
Pedersen
,
T. B.
Schøder
, and
J. C.
Dyre
, “
Phase diagram of Kob-Andersen type binary Lennard-Jones mixtures
,”
Phys. Rev. Lett.
120
,
165501
(
2018
).
53.
U. R.
Pedersen
, “
Long-time simulations of viscous liquids–from strong correlations to crystallization
,” Ph.D. thesis,
Roskilde University (Denmark)
,
2009
.
54.
J. K.
Kummerfeld
,
T. S.
Hudson
, and
P.
Harrowell
, “
The densest packing of AB binary hard-sphere homogeneous compounds across all size ratios
,”
J. Phys. Chem. B
112
,
10773
10776
(
2008
).
55.
U. R.
Pedersen
,
T. S.
Hudson
, and
P.
Harrowell
, “
Crystallization of the Lewis-Wahnström ortho-terphenyl model
,”
J. Chem. Phys.
134
,
114501
(
2011
).
56.
S. J.
Plimpton
, “
Fast parallel algorithms for short-range molecular dynamics
,”
J. Comput. Phys.
117
,
1
19
(
1995
).
57.
N.
Bailey
,
T.
Ingebrigtsen
,
J. S.
Hansen
,
A.
Veldhorst
,
L.
Bøhling
,
C.
Lemarchand
,
A.
Olsen
,
A.
Bacher
,
L.
Costigliola
,
U.
Pedersen
,
H.
Larsen
,
J.
Dyre
, and
T.
Schrøder
, “
RUMD: A general purpose molecular dynamics package optimized to utilize GPU hardware down to a few thousand particles
,”
SciPost Phys.
3
,
038
(
2017
).
58.
U. R.
Pedersen
,
T. B.
Schrøder
, and
J. C.
Dyre
, “
Repulsive reference potential reproducing the dynamics of a liquid with attractions
,”
Phys. Rev. Lett.
105
,
157801
(
2010
).
59.
L. D.
Landau
, “
On the theory of phase transitions
,”
Zh. Eksp. Teor. Fiz.
7
,
19
32
(
1937
).
60.
A.
Tröster
,
C.
Dellago
, and
W.
Schranz
, “
Free energies of the ϕ4 model from Wang-Landau simulations
,”
Phys. Rev. B
72
,
094103
(
2005
).
61.
M.
Goldstein
, “
Viscous liquids and the glass transition: A potential energy barrier picture
,”
J. Chem. Phys.
51
,
3728
3739
(
1969
).
62.
S.
Sastry
,
P. G.
Debenedetti
, and
F. H.
Stillinger
, “
Signatures of distinct dynamical regimes in the energy landscape of a glass-forming liquid
,”
Nature
393
,
554
557
(
1998
).
63.
B.
Vorselaars
,
A. V.
Lyulin
,
K.
Karatasos
, and
M. A. J.
Michels
, “
Non-Gaussian nature of glassy dynamics by cage to cage motion
,”
Phys. Rev. E
75
,
011504
(
2007
).
64.
H.
Kramers
, “
Brownian motion in a field of force and the diffusion model of chemical reactions
,”
Physica
7
,
284
304
(
1940
).
65.
J. C.
Dyre
,
N. B.
Olsen
, and
T.
Christensen
, “
Local elastic expansion model for viscous-flow activation energies of glass-forming molecular liquids
,”
Phys. Rev. B
53
,
2171
2174
(
1996
).