Statistics of small length scale density fluctuations in supercooled viscous liquids

Many successful theories of liquids near the melting temperature assume that small length scale density fluctuations follow Gaussian statistics. In this paper I present numerical investigations of fluctuations in the supercooled viscous regime using an enhance sampling method. I present results for the single component Lennard-Jones liquid, the Kob-Andersen binary mixture, the Wahnstr\"om binary mixture, the Lewis-Wahnstr\"om model of ortho-terphenyl and for the TIP4P/Ice model of water. Results show that the Gaussian approximation persist to a good degree into the supercooled viscous regime, however, the approximation is less accurate at low temperatures. I relate the non-Gaussian fluctuations to crystalline configurations. Implications to theories of the glass transition are discussed.


I. INTRODUCTION
Small length scale density fluctuations in normal homogeneous liquids above the melting temperature obey Gaussian statistics over many orders of magnitude. 1,2][8][9][10][11][12] In this study I directly investigate to what extend Gaussian statistics of small length scale density fluctuations persist into the supercooled viscous regime near the glass-transition.Viscous liquids are highly nontrivial as characterized 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 relaxation.Thus, it is not obvious that Gaussian statistics will persist into the supercooled viscous regime.
][15] In general, liquids can be cooled below melting temperature due to the existence of a free-energy barrier in the form of a critical nucleus. 16The dynamics of a supercooled liquid near the glass-transition is dramatically slower than near the melting temperature.If dynamics were governed by a fixed free energy barrier the slowdown would follow an Arrhenius law.However, many liquid have super-Arrhenius behavior (the first "non").8][19][20] It is an appealing idea that the dynamical heterogeneity is linked to geometric arrangements of locally preferred structures of well-packed particles.Several studies have identified accumulation of a) Electronic mail: ulf@urp.dk1][22][23][24][25][26][27][28][29][30] This gives a picture of a less homogeneous structure with non-Gaussian small length scale density fluctuations.A disadvantages of "locally preferred structure" approach is that it is system specific.In this paper I propose to study statistics in the collective density field.This is a generic approach that can be applied to widely different systems.This is demonstrated by investigating systems belonging to chemically different classes.
2][33] The motivation is the experimental observation that the structure factor is similar in the normal-and in the supercooled liquid regime.The quadratic scaling law of the temperature dependency of the relaxation time [32][33][34][35] originates from generic kinetic contained models 35,36 .These models have trivial thermodynamic statistics, but nontrivial slow dynamics leading to a glass-transition.This picture suggests that statistics of small length scale density fluctuations near the glass transition inherit the Gaussian statistics of the normal liquid regime.
In this study I examine the statistics of small length scale density fluctuations for the single component Lennard-Jones (LJ) 37 model, the binary Kob-Andersen mixture (KABLJ) 38 , the Wahnström binary mixture (WABLJ) 39 , a coarse grained model o-terphenyl (LWoTP) 40 and the TIP4P/Ice 41 water model.Enhanced sampling molecular dynamics methods are used to sample statistics into the wings of the distributions.Results show that the Gaussian hypothesis is fair in the supercooled regime, however, deviations are more significant at low temperatures.The analysis suggest that non-Gaussian features are related to first-order transitions to crystals.
The remainder of the paper is organized as follows.Section II introduce the formalism used the describe density fluctuations and some theory of the Gaussian hy-pothesis.Section III describe numerical methods used for enhanced sampling of the density field and give descriptions of the investigated models (e.i.energy surfaces).Section IV present the results, and implications are discussed in Section V.

II. FORMALISM AND THEORY
In experiments (X-ray or neutron scattering) and theories of the liquid states it is often convenient to work in reciprocal space.In this section I will give the formalism used to describe density fluctuations in both reciprocal k-space and a subvolume.
A. The collective density field in k-space Consider a liquid of N particles located at R = {r 1 , r 2 , . . ., r N } in a volume V with periodic boundaries so that the thermodynamic density is ρ = N/V .Let the real-space density field be ρ(r) = N n δ(r n − r) where δ is Dirac's delta function.The collective density field in reciprocal space is then defined as the Fourier transform of the real-space density field: where k = k k is the scattering-or wave vector (sometimes the letter "q" or "Q" is used).The 1/ √ N factor ensure system size scale invariance of amorphous configurations (liquids).The scaling is √ N for configurations with long-range translational order (crystals) along the k direction.For a system of N point particles in a periodic orthorhombic cell the 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 x L y L z ).Due to the anisotropy of liquids the investigation can be limited to k-vectors along the x-direction without loss of information.For a given cubic box with a size of L = L x = L y = L z we consider vectors of lengths k = 2πn/L where n is an integer (n x = n and n y = n z = 0).We note that the anisotropy of the liquid is in principle destroyed by the constraint of the anisotropic periodic boundaries.However, such effects are expected to be small and are ignored in this study.
In the following we consider the probability distribution P (|ρ k |) where |ρ k | is the norm of the collective density field.We do not need to consider the tedious twodimensional complex plane since the P (ρ k ) distribution is radial symmetric for a liquid.The second moment S k = |ρ k | 2 is the structure factor routinely measured in scattering experiments.If statistics of density fluctuations follow Gaussian (G) statistics then probability ( The central limit theorem dictates that in the thermodynamic limit density fluctuations become Gaussian (see also discussion in Section V): for N → ∞.Thus, we limit our analysis to small length scale fluctuations by studying systems of about 100 to 1000 particles (unless otherwise stated).It have been shown that a small system can represent viscous dynamics of larger system. 42he fourth moment . Thus, we define a non-Gaussian parameter for the |ρ k | fluctuations as ( This parameter quantifies deviations from Gaussian statistics near the center of the distribution, |ρ k | 0. Deviations in the tails of the distribution cannot be expected to be represented by this parameter.For this, higher order moments are relevant.

B. The density fluctuations in a subvolume
The central limit theorem dictates that non-Gaussian feature are more pronounced in smaller systems.Thus, it can be illustrative to investigate density fluctuations in small subvolumes of a larger system.We define a subvolume though 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 ρ

√
V h where V h is the size of the subvolume.The k = 0 relates the actual density ρ h = N h /V h in the subvolume.Here N h = N n h(r n ) is the number of particles in the subvolume.The Gaussian approximation of the ρ h density fluctuations is where m 2 = (ρ h − ρ h ) 2 is the variance.For this distribution the fourth central moment m 4 = (ρ h − ρ h ) 4 equals 3m 2 2 .Thus, we define a non-Gaussian parameter as

III. METHODS AND MODELS
To highlight non-Gaussian features, I suggest to apply a potential that bias the system towards rare configurations that would not be sampled otherwise.Specifically, a harmonic potential is added to the Hamiltonian that will bias the system towards large ρ k values.The Gaussian hypothesis can then either be investigated directly on statistics of the biased simulations, or by re-weighting statistics of a series of simulations (referred to as the "umbrellas sampling method" 43 ).Below is a description of the suggested method: A. Sampling rare fluctuations of the collective density field Let H(R, Ṙ) be the Hamiltonian of a given system.To sample rare ρ k fluctuations at some density and temperature we simulations a Hamiltonian with added harmonic bias field 44 : where κ is a spring constant and a is an anchor point of the bias field.By reweighing we can get that the |ρ k | probability distribution of the unperturbed system: where P κa (|ρ k |) is the distribution of the Hamiltonian with the harmonic bias field.For a series of overlapping distributions (with different a's and κ's) the normalization constants N κa can be determined numerically using the iterative multistate Bennett acceptance ratio (MBAR) method. 45Alternatively, the distribution function P κa can be investigated directly.By combining equations 2 and 7 we get that the Gaussian approximation predicts where ment of biased distribution is Results of the Gaussian approximation can be used as initial guesses for the iterative MBAR method. 45In practice this lead to fewer iterations before reaching convergence.
To perform molecular dynamics simulations forces on particles from the bias field needs to be evaluated.The total force acting on particle j is 44 where is the force of the unbiased Hamiltonian, and Although the force on particle j depends on the positions of all particles it is possible to design an N -scaling algorithm: 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 parallelize to several processes since both the computation of ρ k and the F (κa) j forces involve sums of independent contributions (assuming the same for F (0) j ).Computational efficiency of the algorithm is crucial, since we wish to conduct long-time simulations in the viscous regime where dynamics are slow.

B. Sampling rare density fluctuations in a subvolume
The overall idea of the above mentioned method for computing rare fluctuations of collective density field can be reused to sample rare density fluctuations in a subvolume.In order to perform molecular dynamics with a bias field we define a continuous quantity Ñh that is strongly correlated with the number of particles N h inside the volume h.In practice this is done by using a switching functions on the boarders of the volumes as described in Ref. 46 .The unbiased P (ρ h ) distribution is obtained by reweighing biased P bias ( Ñh ) distribution using the MBAR method 45 .For binary mixtures a bias potential can be applied to both kinds of particles.Again statistical information of the unbiased system can then be determined with reweighing.

C. Energy surfaces
We investigate statistics of density fluctuations for several models defined as a 3N dimensional energy surface.The 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 repulsive term representing Pauli repulsions and an attractive term representing London forces. 37In this study we investigate a truncated version: the potential energy surface is a sum over pair energies and zero otherwise.The LJ model is not a good glass-former, however, the it is included in this study since it is a standard system in computational condense matter physics.Temperature T = 0.8 and density ρ = 0.85 (L = 5.0273) is used as a representation of the "normal liquid regime".This state is close to the freezing temperature. 44,48BLJ: Kob and Andersen suggested a binary LJ mixture as model of a good glass former. 38It is an 80/20 mixture of particles that have a strong affinity towards unlike atoms.This parametrization makes a good glass-former on time-scales and system sizes typically investigated in silico.This model is the standard model for computational studies of low temperature liquid dynamics.It is custom to study the system at the density ρ = 1.2 where the melting temperature is T m = 1.027(3). 49elow this temperature the particles will eventually phase separate in long-time simulation.The major constituent, the A's, will form a face centered cubic crystal.If crystallization is avoided, however, the low-temperature liquid accumulate locally preferred structures where one of the small particle is surrounded by ten larger particles forming a twisted bicapped square prism 24,30 .
WaBLJ: Wahnstöm suggested a 50/50 binary LJ mixture with a size ration of 80% 39 .Unlike the KABLJ mixture the interactions 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 can form a MgZn 2 crystal structure 26 .In the supercooled regime the liquid collect local structures of icosahedral order and Frank-Kasper order. 24,26,29The latter is a geometric arrangement where two touching larger particles have six smaller particles as common neighbors.These structures are favored by the low-temperature liquid since they pack space well, and are also a part of the crystalline structures. 26,50oTP: Lewis and Wahnström 40 suggested a coarse grained model of ortho-terphenyl (C 18 H 14 ) where molecules are constructed from three LJ particles placed in the corners of an isosceles triangle.Each LJ particle represent a benzene ring.To avoid that LJ particles crystallize into a close-packed structure the molecule have an inner angle of 75 • (that is in between 60 • and 90 • degree -the angles found between neighbor 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 molecule. 47,51We study a system of N = 324 molecules (unless else stated) at the temperature T = 350 K at density ρ = 1.09 g/ml (L = 4.84 nm).
Water: Abascal et al. 41 suggested the TIP4P/Ice atomistic water model.This four site model reproduce the complicated phase-diagram of real water, suggesting that it also gives a good representation of hydrogen-bonds in the liquid state.The model is studied at temperature T = 280K and density ρ = 1 g/ml.There are no signs of spontaneous crystallization.
Numerical computations are performed using the software packages LAMMPS 52 , RUMD 53 , and home-made code available at the website http://urp.dk/tools.Implementation of the ρ k bias field is available in the official LAMMPS and RUMD packages.The P (|ρ k |) distribution on a logarithmic scale (k = 1.4 nm −1 ) of the LWoTP model for system sizes of N = 324 and N = 2592 molecules, respectively.The Gaussian is better for the larger system as expected from the central limit theorem.The non-Gaussian parameters αρ k are 0.029 and 0.0011 for N = 324 and N = 2592, respectively.are studied in systems with a gas-liquid interface.This is done 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.
Results for the LJ, KABLJ and WaBLJ models are reported in reduced Lennard-Jones units, while physical units units are used for the LWoTP model and the TIP4P/ice water model.

A. Fluctuations of the collective density field
Before investigating the glass forming liquids, we first focus on the LJ model near the melting temperature.Figure 1 shows probability distributions P (|ρ k |) on a logarithmic scale.The figure includes k-vectors from lengths of k = 1.25 (n = 1) up to k = 15.0 (n = 12).The solid black lines are the reweighed distributions using a series of biased simulations and the red dashed lines are the Gaussian predictions (shifted vertically for clarity).The first impression is that the Gaussian hypotheses give a good description.This confirm the consensus that small length scales fluctuations are Gaussian in the normal liquid regime. 1,2The 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 configura- tion from the tail of the distribution for n = 6 is shown on Fig. 2(a).A crystalline structure is apparent in both the real-space configuration and the scattering spectrum shown on Fig. 2(b).Consistent with this, 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 (Fig. 2(c) and Fig. 2(d)) also have a fat tail.This can also be attributed to a crystalline configuration, but with another orientation.As an aside, bias simulation similar to the ones presented here, can be used to compute the melting point of crystals.This is refereed to as the "interface pinning" method. 44ther k-vectors have thin tails.As an example Figs.2(e) show a configuration from the tail of the n = 10 wave vector.This structure it not crystalline but have disordered.Figs.2(g) and 2(h) show a configuration from the longest wave vector of the investigated system size (n = 1).The liquid have responded to a strong bias field by forming a vapor slap and a crystalline slap.
Next we investigate the glass forming models.First, we consider the KABLJ mixture at T = 0.45 (ρ = 1.2) of a system size of N = 1000 particles.This state-point is well below the melting temperature of T m = 1.027(3) 49 .The structural relaxation is about 10 3 times larger than in the normal liquid regime. 38,54The solid black lines on Fig. 3  trajectories of the KABLJ mixture and LWoTP trimer are shown on Fig. 3. Crystallizing trajectories are discarded in the analysis.Figure 5 show P (|ρ k |) distributions of the glass forming liquids KABLJ, WaBLJ, LWoTP and water for several kvectors.The red dashed lines indicate the Gaussian approximation.The agreement is good, but the tails of the distributions deviates 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 increased from N = 324 to N = 2592.

B. Density fluctuations in a subvolume
Fluctuations in small subvolumes of a larger system can give a further insight to the structural origin of non-Gaussian small length scale density fluctuations.Figure 7 show the distribution function of the ρ h density in a 3 × 3 × 3 subvolume, h.The points are reweighed statistics from simulations with a bias potentials that push the system towards configurations with a certain amount of particles in the subvolume.The red dashed line is the prediction from the Gaussian approximation.The agreement is good, 1,2 however, some deviations are seen in the tails of the distributions.The low density limit correspond to the formation of a cubic vapor bubble in the liquid.As described by classical nucleation theory the free energy −k B T ln(P ) for forming such a bubble is associated with both a bulk-and a surface contribution.
To investigate the supercooled regime we setup a gasliquid coexistence simulation of the KABLJ mixture (in the same way as we did for the single component LJ model).Figure 8 shows the structural relaxation time in liquid slap.The relaxation time is non-Arrhenius in the investigate temperature regime, 0.28 < T < 0.60.The points on Fig. 9 shows the distribution of ρ h fluctuations in a 3 × 3 × 3 subvolume h for the temperatures T = {0.35,0.40, 0.55}.Gaussian statistics are shown as red dashed lines.The conclusion from the analysis of the |ρ k | fluctuations remains -the Gaussian approximation gives a fair description, but becomes increasingly less accurate at lower temperatures.Deviations are both seen in the tails of the distribution, and near the mean as shown by the non-Gaussian parameter α ρ h .Figure 10(a) shows the P (N A , N B ) distribution.We would expect to see elliptical shaped contour lines for Gaussian approximation, but see some deviations from this. Figure 10(b) show a configuration from a fat-tail part of the distribution at equimolar composition in the subvolume.The configuration is a cubic CsCl crystallite.This structure is one of the thermodynamically stable crystal structures of the KABLJ model. 49Thus, I conclude that non-Gaussian feature are related to the first order-transition to a crystal.

V. DISCUSSION
Let us summarize the results before moving on with a discussion of the implications.I have presented an investigation of a range of chemically different glass formers, and the overall conclusion is that the Gaussian hypothesis gives a fair description of the small length scale density fluctuations, however, as temperature is lowered the Gaussian approximation is less accurate.
Gaussian statistics usually comes about in two ways: (i) from the central limit theorem, or (ii) from an harmonic approximation: (i) The central limit theorem dictates that if random variables from any underlying distribution are added the resulting distribution will follow Gaussian statistics.For a non-flowing equilibrium liquid of a sufficient size it can be assumed that subvolumes fluctuate independently.As an example think of the ideal gas model of non interacting particles.For the ideal gas the number of particles in a given subvolume h follow the Poisson distribution.If the size of the subvolume is increased, then Raikovs theorem dictates that fluctuation follow another Poisson distribution with a higher average.This distribution will be closer to the Gauss function.For a liquid with interactions the underlying distribution differs from Poisson statistics, however, by studying small length scale density fluctuations we can gain insight on the nontrivial underlying distribution.This brings us the other less trivial way of arriving at Gaussian statistics, (ii) i.e. by an harmonic expansion around a local minimum of the free energy function F : If x is a order parameter, like ρ k or ρ h , then the free energy along this coordinate is where P (x) is the probability distribution.The function F (x) can be expressed as a polynomial expansion around the minimum at x 0 .It is often convenient to assume that only the second order term is of relevance, thus giving a harmonic approximation for the free energy 5 The truncation of the expansion series is non-trivial, as discussed below.By equating Eqs. 13 and 14 and isolating P we arrive at Gaussian statistics for the probability distribution (I remind the reader that the reason we find near-Gaussian statistics is not due to the harmonic bias field).
To understand first order transitions, e.g. the gas-liquid transition, higher order terms are important.As a classic example, Landau's (L) effective Hamiltonian 55,56 includes higher-order terms to give a description of the density fluctuations near the gas-liquid critical point: . With the Landau theory in mind, we expect deviations from Gaussian statistics for liquid density fluctuations when other phases (gas or crystals) interfere with the liquid state.In the normal liquid regime deviation from Gaussian statistics can be due to the formation of a vapor bubble as exemplified on Figs. 1 and 2. In the supercooled regime the crystal basin of the free energy becomes large suggesting that statistics becomes less Gaussian due to the presence of a crystal.The microscopic picture 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 as exemplified on Fig. 4. I conjecture that it is possible that non-crystalline structures could be important for the statistics of small length scale fluctuations.To address this more refined methods are needed, and I leave this to future studies.Some structural candidates are the "locally favored structures" identified for some of the models [20][21][22][23][24][25][26][27][28][29][30] .These structures have been suggested as an important component to understand the dynamics of the highly viscous liquid near the glass transition.It would be a valuable insight to show that statistics of small length scale density fluctuations is related to dynamics.
Another angle is to view small length scale density fluctuations from is the "energy landscape" perspective 42,[57][58][59] .In this picture, the 3N dimensional energy surface of the liquid is partitioned into basins identified by local minimums.Below a certain onset temperature the system explores confrontational space by 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 timescales the system will explore different basing (activated relaxation).From this perspective the non-Gaussian feature at low temperatures are related to density fluctuations between basins (the inherent states).However, this needs to be investigated.Some theories directly or indirectly assume Gaussian statistics of small length scale density fluctuations.3][34][35][36] As mentioned in the introduction, the picture is that thermodynamics and structural details of glass forming materials are note crucial for understanding the dynamics of a highly viscous liquid.Finally, the Gaussian approximation of density fluctuations play a role in some elastic models. 10,60,61The idea behind these approaches to understand viscous dynamics is though elastic deformation that allow a small length scale subvolume to rearrange.Here the harmonic approximations enters in some theories to make predictions related to experiments.

VI. ACKNOWLEDGMENTS
This study was initiated by discussion with the late David Chandler.I joined his group in 2010-2012, and I 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.For the preparation of this manuscript I also received comments and suggestions from Andreas Tröster, Jeppe C. Dyre and Thomas B. Schrøder.This work was supported by the VILLUM Foundations Matter Grant No. 16515.

FIG. 1 .
FIG. 1. Probability distribution P (|ρ k |) on a logarithmic scale with k = (2πn/5.0273,0, 0) for the single component LJ model in the normal liquid regime.The solid black lines are the distribution function computed from reweighed biased simulations and the red dashed lines are the corresponding Gaussian predictions (Eq.2).The distributions have been shifted vertically for clarity.The insert show the structure factor where dots indicate the investigated k-vectors.

FIG. 2 .
FIG. 2. Configurations taken from the non-Gaussian tails of distributions shown on Fig. 1.The left panels show a representative configuration, and the right panels show the scattering spectra S k in the xy-plane (average over several configurations).From top down the panels shows representations taken from biased simulations with wave vectors with n = 6, n = 8, n = 10 and n = 1.

FIG. 4 .
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 consist of pure A's and crystallization event 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 a indicator of the crystallization.(c) The configuration in the last step of the simulation.The larger A's are colored green, while the smaller B's are red.(d) The qubatic orderparameter 47 and (e) the potential energy in biased simulations of the LWoTP model.(f) The final configuration of a trajectory where a crystal is formed.Each molecule have been give an individual color to tell them apart.
FIG. 6.The P (|ρ k |) distribution on a logarithmic scale (k = 1.4 nm −1 ) of the LWoTP model for system sizes of N = 324 and N = 2592 molecules, respectively.The Gaussian is better for the larger system as expected from the central limit theorem.The non-Gaussian parameters αρ k are 0.029 and 0.0011 for N = 324 and N = 2592, respectively.

FIG. 7 .
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 show a typical configuration of the system with an gas-liquid interface, and the subvolume h located in the bulk liquid part.

ρ
FIG. 8. Structural relaxation time τα as a function of temperature of KABLJ mixture with 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 show the density of the liquid slap as a function of temperature.

FIG. 10 .
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).Nongaussian features are seen as the contour lines that deviates slightly from being oval.(b) A configuration from the tail of the distribution with equimolar composition .The upper half of particles have been made invisible to reveal the arrangement of particles in the 3 × 3 × 3 subvolume h.The structure correspond to a cubic CsCl crystallite.This is one of the known stable structures of the mixture.49