Finite-size scaling and thermodynamics of model supercooled liquids: long-range concentration fluctuations and the role of attractive interactions

Abstract

We compute partial structure factors, Kirkwood–Buff integrals (KBIs) and chemical potentials of model supercooled liquids with and without attractive interactions. We aim at investigating whether relatively small differences in the tail of the radial distribution functions result in contrasting thermodynamic properties. Our results suggest that the attractive potential favours the nucleation of long-range structures. Indeed, upon decreasing temperature, Bathia-Thornton structure factors display anomalous behaviour in the k→0 limit. KBIs extrapolated to the thermodynamic limit confirm this picture, and excess coordination numbers identify the anomaly with long-range concentration fluctuations. By contrast, the purely repulsive system remains perfectly miscible for the same temperature interval and only reveals qualitatively similar concentration fluctuations in the crystalline state. Furthermore, differences in both isothermal compressibilities and chemical potentials show that thermodynamics is not entirely governed by the short-range repulsive part of the interaction potential, emphasising the nonperturbative role of attractive interactions. Finally, at higher density, where both systems display nearly identical dynamical properties and repulsive interactions become dominant, the anomaly disappears, and both systems also exhibit similar thermodynamic properties.

---

1 Introduction

The supercooled state challenges our understanding of the theory of liquids. In particular, the connection between dynamics, which varies considerably upon supercooling, and structure, which appears to remain essentially unchanged, is the subject of intense research.^1–7 Model systems with reduced complexity, still retaining essential physical features, provide a direct route to investigate this problem. For example, Kob–Andersen mixtures⁸ with purely repulsive Weeks-Chandler-Andersen interactions (KAWCA)⁹ exhibit substantially different dynamics compared to their Lennard-Jones counterpart (KALJ).⁸ By contrast, their structure, investigated from the point of view of radial distribution functions, is somewhat similar.^7,10

The connection between pair correlations and dynamical properties has been extensively investigated.^10,11 On the one hand, a variety of studies conclude that two-body contributions are not enough to account for the difference in dynamics between the KAWCA and KALJ systems. Perhaps the most well-known example is mode-coupling theory, based on pair correlation functions, which underestimates these dynamical differences.¹⁰ Additionally, deviations in many-body structural descriptors such as triplet¹² and point-to-set correlations,¹³ as well as bond-order distributions¹⁴ and the packing capabilities of local particle arrangements,² have been observed between the KALJ and KAWCA systems. These results indicate that higher-order features may be necessary to resolve the difference in their dynamical properties.¹⁵

On the other hand, several studies indicate that two-body structure is enough to describe particular aspects of the dynamics of model supercooled liquids. For example, features based on the pair structure have been used to predict diffusion constants from short-time trajectories of the KALJ model.^16–18 Concerning the comparison between models, Bhattacharyya and coworkers^19,20 directly explored structure-dynamics relationships in KALJ and KAWCA systems. In particular, they used the Adam–Gibbs relation,²¹ to connect relaxation time to the configurational entropy. Their results demonstrated that the two-body contribution to the entropy plays a significant role in distinguishing the dynamics of the two systems.

To further contribute to the discussion, recent research efforts have focused on the detailed characterization of the liquid's two-body structure. In particular, softness parameters, defined via weighted integrals of pair-correlation functions^4,22 or multi-dimensional integrals of partial structure factors,²³ respond to minor structural changes and can accurately describe dynamical differences. However, either non-trivial reweighting procedures or combinations of local and nonlocal terms prevents an unambiguous identification of the dominant, short- versus long-range, contributions to the resulting structure-dynamics relationship.

The potentially dominant role of short-range pair correlations brings with it yet another dilemma. According to perturbation theory, short-range repulsive interactions mostly dominate the liquid's structure.⁹ By contrast, based on Kirkwood–Buff theory,²⁴ long-range fluctuations in the tail of the pair correlation function have a significant effect on the system's solvation thermodynamics.^25–29 The studies mentioned above investigating KALJ and KAWCA dynamics have mainly focused on short-range contributions. Nevertheless, evidence for the nucleation of long-range structures in glassy systems at low temperatures^30–33 highlights the necessity to carefully address this point. Finite-size effects present in computer simulations dramatically affect the tail of the pair correlation function and the k→0 limit of the structure factor, i.e., the long-range structure properties, which in turn sensitively impact thermodynamic quantities. Consequently, a careful evaluation of finite-size effects becomes critical for investigating these properties in the supercooled regime.

In this paper, we investigate various thermodynamic properties of KALJ and KAWCA a–b mixtures in the supercooled liquid state. We calculate structure factors of density, S_ρρ(k), and concentration, S_cc(k), while highlighting the k→0 limit. The KALJ liquid exhibits anomalous behaviour reflected in a major increase in concentration fluctuations. This anomaly closely resembles the nucleation of nanometric clusters reported by Fischer in low-temperature ortho-terphenyl,^30,31 and it has been recently identified as a general feature present in polydisperse colloidal models.⁵ By contrast, the purely repulsive KAWCA system remains perfectly miscible in the supercooled state. A finite-size Kirkwood–Buff analysis confirms this picture by enabling the precise identification of the k→0 limit. Furthermore, we show that the isothermal compressibility and chemical potential of the two models exhibit similar trends with temperature, apart from constant shifts. These differences highlight the nonperturbative role of attractive interactions in the system. To sum up, we demonstrate that seemingly small differences in the tail of the radial distribution function result in significantly different structural and thermodynamic properties for supercooled systems with and without attractive interactions.

The paper is organised as follows: we provide the computational details in Section 2, present the results in Section 3 and conclude in Section 4.

2 Computational details

We have simulated the Kob–Andersen model, which is a binary mixture (80 : 20) of Lennard-Jones (KALJ) particles.⁸ The inter-atomic pair potential between species α and β, U_αβ(r), with α,β = a,b is described by a shifted and truncated Lennard–Jones potential,³⁴ as given by:where U_αβ^(LJ)(r;σ_αβ,ε_αβ) = 4ε_αβ[(σ_αβ/r)¹² − (σ_αβ/r)⁶] and r_αβ^(c) is equal to 2.5σ_αβ for LJ system and r_αβ^(c) is equal to the position of the minimum of U_αβ^(LJ) for the WCA systems (KAWCA).⁹ We have added a linear correction so that both the potential and the force go to zero continuously at the cutoff distance.³⁴ We have used LJ natural units, such that length, temperature and time are measured in σ, k_BT/ε and τ = √(mσ²/ε), respectively. For all the simulations, we have used the following interaction parameters σ_aa = 1.0σ, σ_ab = 0.8σ, σ_bb =0.88σ, ε_aa = 1.0ε, ε_ab = 1.5ε, ε_bb = 0.5ε, m_a = m_b = 1.0m.

We have performed two different sets of simulations: the first for the calculation of structural properties, and the second for the calculation of chemical potentials, which employed a different box geometry and number of particles. All simulations have been carried out using the LAMMPS molecular dynamics software.³⁵ We have performed the first set of simulations in a cubic box with periodic boundary conditions in the canonical ensemble (NVT), using the Nosé–Hoover thermostat³⁶ with an integration timestep of 0.005τ and a time constant of 100 timesteps. The system is composed of N = 23 328 particles, with N_a = 18 664 particles of type a. We have simulated this system at two different densities, ρ = 1.2/σ³ and 1.6/σ³ for different temperatures. Starting from the high temperature case, the final configuration of the simulation has been used as an initial configuration for the simulation one (temperature) step below. The same procedure has been followed for the KALJ and KAWCA systems. For all state points, three to five independent simulations with run lengths >100τ_α (τ_α is the α-relaxation time estimated from ref. 20) have been performed.

To calculate the excess chemical potential, we have used the LAMMPS³⁵ implementation of SPARTIAN already described in ref. 37. The SPARTIAN method, a variant of the adaptive resolution method,^38–43 simulates the coexistence of an atomistic system to its ideal gas representation at a constant density and temperature. We have computed the excess chemical potential of the system as the external potential required to balance the density across the simulation box. To guarantee enough statistics, we have used a slab geometry (An anisotropic box with L_x = 36σ, L_y = 578σ and L_z = 10σ), also with periodic boundary conditions and at density ρ = 1.2/σ³, resulting in a system with N = 2 50 000 and N_a = 2 00 000. The same protocol as described above has been used to quench the system before performing the SPARTIAN calculation. For the SPARTIAN method calculation, we have considered an slab geometry with atomistic region of length of 10σ and hybrid regions of linear size 10σ aligned along the x direction.

After equilibration, we have performed the SPARTIAN calculations in the canonical ensemble (NVT), using a Langevin thermostat with dt = 0.001τ and damping parameter of 10τ. In order to get the correct density profiles and therefore, chemical potential, we have simulated for 3 × 10⁶ simulation steps.

3 Results and discussions

3.1 Kirkwood Buff analysis

We consider temperatures in the range 0.45ε/k_B ≤ T ≤ 6ε/k_B for KALJ system and 0.3ε/k_B ≤ T ≤ 6ε/k_B for KAWCA system (see Section 2). Visual inspection of the radial distribution functions (RDFs) for both systems reveals that they are almost indistinguishable (Fig. S1, ESI†), and only the RDF g_bb(r) for the minor component shows relatively small differences, visible at r <3σ (Fig. 1(a)).^20,44 However, this direct comparison is misleading: a few thermodynamic quantities are quite sensitive to small fluctuations in the tail of the RDFs.

Fig. 1 Differences between KALJ and KAWCA systems in terms of the g_bb(r) component and the KBIs at T = 0.45ε/k_B. (a) Differences between the RDF for the low-concentration b-component of the mixture seem to be small and mostly coming from the local structure of the fluid. (b) G^R_bb as obtained from eqn (3) shows a different short-range behaviour and, more importantly, the tails do not converge due to finite-size effects. (c) KBIs obtained using the method described in ref. 29 (eqn (4)). The KBIs in the thermodynamic limit G^∞_αβ are obtained from the slope of a linear fitting of the region 0 < λ < 0.3. This straight line is indicated for the bb case. Horizontal, dark lines correspond to the asymptotic limit −δ_αβ/ρ_α with δ_αβ the Kronecker delta and ρ_α the number density of the α-species. The G^∞_αβ values obtained in this way are plotted in panel (b) as horizontal lines.

One such quantities are the Kirkwood–Buff integrals (KBIs),²⁴ which relate the microscopic structure of a liquid mixture to its solvation thermodynamics. For a multi-component system of species α and β, in equilibrium at temperature T, the KBIs in the thermodynamic limit (TL) take the form

where g_αβ is the radial distribution for an infinite, open system. Here, it is obvious from eqn (2) that small deviations for large r result in important contributions to G_αβ. In computer simulations, usually far from the thermodynamic limit, eqn (2) is often approximated aswhere g^c_αβ(r) is the RDF of the closed, finite, system and R is a truncation radius. It is essential to choose R larger than the correlation length of the system. Nevertheless, this expression seldom converges due to different finite-size effects. Here, it is already clear that G^R_bb for the KALJ and KAWCA systems displays different behaviour (see Fig. 1(b)).

By explicitly including finite-size effects due to the thermodynamic ensemble and the finite integration domains, we compute the KBIs as²⁹

with , G^∞_αβ being the value of the KBIs in the thermodynamic limit, c_αβ a constant with units of length, δ_αβ the Kronecker delta and ρ_α the number density of the α-species. We can compute G_αβ(λ), the KBIs for a subdomain of volume V inside a simulation box of volume V₀, in terms of fluctuations of the number of particles^{25–29,45,46}where G_αβ(λ) ≡ G_αβ(V;V₀) and the average number of α-particles, 〈N_α〉′ ≡ 〈N_α〉_V,V0, depends on both the subdomain and simulation box volumes. Fig. 1(c) shows the results obtained from eqn (4) and (5) for the KALJ and KAWCA systems at T = 0.45ε/k_B. These curves are rather similar in both cases, with a major difference appearing for the G_bb^∞ component, which can be obtained as the slope of a linear fit of G_bb(λ) within the region λ <0.3. The resulting values of G_bb^∞ are plotted as dashed lines in Fig. 1(b) to indicate the value at which the KBIs should converge.

3.2 Density and concentration structure factors

As anticipated, fluctuations in the tail of the radial distribution function affect the long-range structure of the fluid. Hence, to investigate these effects, we compute partial structure factorswhere k is the norm of a reciprocal-lattice vector, ρ = ρ_a + ρ_b is the total number density and x_α = N_α/N is the mole fraction of the α-species. To avoid numerical instabilities at the low k limit,⁴⁷ we compute the structure factor directly from the simulated trajectory¹⁰ using the following expression as,where α and β denote the species, and the indexes i and j run over particles belonging to α and β, respectively. The average runs over the values of k such that |k| = k and over the ensemble. We follow the standard procedure in computer simulations to compute the structure factor. Namely, for a cubic simulation box of linear size L with periodic boundary conditions, we discretise the reciprocal space by considering wavevectors k = 2π(n_x,n_y,n_z)/L with n_x, n_y and n_z integer values. We use a maximum number n_x^max = n_y^max = n_z^max = 76.

Partial structure factors are difficult to interpret for liquid mixtures. Hence, we focus on density, S_ρρ(k), and concentration, S_cc(k), structure factors⁴⁸ which carry a direct physical meaning.⁴⁹S_ρρ(k) and S_cc(k) describe the correlation of density and concentration fluctuations in the liquid mixture. They are defined as

For large k-values, the behaviour of S_ρρ and S_cc is rather similar for both systems (see Fig. 2). This includes a first peak at k₀ ≈ 7.13/σ, followed by a second peak at approximately 1.7k₀ that develops at low temperatures. This second peak is associated with the nucleation of structural motifs that precede the complete crystallisation of the system. As it has been reported for various metallic glasses, the splitting of this second peak⁵⁰ results from the optimal facet-sharing configurations of such structural (icosahedral and tetrahedral) motifs that grow upon decreasing temperature.^51,52 In our particular case, we do not observe this feature down to T = 0.45ε/k_B, thus confirming that both systems remain liquid-like. Concerning the difference between the KALJ and the KAWCA systems, the first and second peaks in the S_cc show slightly more structure for the KALJ system at T = 0.45ε/k_B, as expected from the RDF (see Fig. 1(a)).

Fig. 2 Density, S_ρρ (top), and concentration, S_cc (bottom), structure factors for both KALJ (blue) and KAWCA (red) systems for the temperatures considered here.

Perhaps more interesting, it is apparent from the inset in Fig. 2 that the KALJ and KAWCA systems show substantially different behaviour in the region of small k (large r). On the one hand, the KAWCA liquid behaves like a normal liquid with monotonically decreasing density fluctuations upon decreasing temperature. On the other hand, the KALJ system exhibits anomalous behaviour, similar to SAXS curves obtained for ortho-terphenyl³¹ and supercooled water,⁵³ with clear density fluctuations starting around k ∼ 2/σ (r ∼ 3σ) appearing at temperatures lower than the onset temperature of glassy dynamics T = 1ε/k_B (see Fig. S2, ESI†).⁵⁴ These results indicate that the two systems display stark structural differences in the supercooled regime, with clear long-range density domains (r >3σ) induced by the presence of attractive interactions in the KALJ mixture.

The extrapolation to the k→0 limit by using eqn (6) or (7) is not trivial because finite-size effects in the simulation affect the precision in computing structure factors as we approach the linear size of the simulation box. In the next subsection, we use the relation between the structure factor in the limit k→0 and the KBIs to investigate this limiting case in more detail.

3.3 KBIs and the k→0 limit

Similar to the single component case, the extrapolation to the k→0 limit provides useful physical information.⁵⁰ Here, we use the relation between the limit k→0 in the structure factor, evaluated in eqn (6), and the KBIs in the thermodynamic limit, eqn (2). We obtainthus

The last relation in eqn (10) gives two contributions that allows us to connect long-range density fluctuations to both the isothermal compressibility κ_T of the system and to concentration fluctuations modulated by the difference in partial molar volumes v_a − v_b, with δ = ρ(v_a − v_b).⁵⁰ The isothermal compressibility and the partial molar volumes can also be written in terms of the KBIs, namely:

andwhere η = ρ_a + ρ_b + ρ_aρ_b(G^∞_aa + G^∞_bb − 2G^∞_ab).

We use the definition in eqn (7) and (8) to compute S_ρρ(k), and compare with the obtained from the KBIs (eqn (10)). The results, presented in Fig. 3 (top panel), confirm the information given by the partial structure factors (see Fig. S2, ESI† for a comparison between the values obtained from the structure factor and the KBIs). Namely, in contrast to the KAWCA system, the KALJ system exhibits increasingly large density fluctuations upon decreasing temperature. To investigate the origin of the anomaly, we investigate the contributions to S_ρρ separately as given by the r.h.s. of eqn (10). The middle and lower panels of Fig. 3 splits S_ρρ into isothermal compressibility and concentration fluctuation terms, respectively. There, it is apparent that the anomalous behaviour exhibited by the KALJ system at low k values is due to the formation of long-range concentration domains (red and blue triangles). By contrast, the isothermal compressibility contribution remains nearly the same for both systems (red and blue circles).

Fig. 3 Density–density correlation function S_ρρ(k). (Upper panel) obtained from the KBIs (eqn (10)). At high temperature, both systems present a similar monotonically decreasing behaviour upon decreasing temperature. At the onset temperature of glassy dynamics (T = 1ε/k_B),⁵⁴ the data for the KALJ system shows an inflexion point which signals density–density correlations visible for distances longer than r = 2.5σ. Individual components of : (middle panel) κ_Tρk_BT and (lower panel) with δ = ρ(v_a − v_b) the product of the total density with the difference in partial molar volumes. It is apparent that the contrast in S_ρρ originates from major concentration fluctuations present in the KALJ system, as indicated by S_cc(k).

The anomaly observed in the limit k→0 in Fig. 2 has also been reported in ref. 55. The authors suggest a plausible explanation involving the gas–liquid phase separation of the KALJ system. However, we note that the gas–liquid coexistence region for the KALJ system is still far from the point ρ = 1.2/σ³, T = 0.45ε/k_B^56,57 (see also the discussion in Section IV in ref. 55). Moreover, we observe the non-monotonic behaviour starting just below the onset temperature of glassy dynamics T = 1ε/k_B, which is even farther away from the coexistence region. Moreover, the virial part of the pressure remains positive even at the lowest temperature considered for the KALJ model (see Fig. 1 in ref. 58 and Fig. S3, ESI†). For the KAWCA model, as expected, the system's pressure is systematically higher than the KALJ pressure due to the absence of attractive interactions. In general, the positive pressure of our simulated state points already suggests that the anomalous behaviour is not due to gas–liquid coexistence.

Plots of the excess coordination number (N_αβ = ρ_βG^∞_αβ), which gives the change in the number of α particles when one β particle is removed from the system, as a function of temperature provide a clear insight (see Fig. 4). As expected from the model, the effective interaction between a and b particles is favoured in both systems at all temperatures: excess coordination numbers are close to zero. Below the onset temperature of glassy dynamics, the excess coordination number shows a collective tendency for the KALJ mixtures to increase b–b effective interactions upon cooling. This tendency is not apparent for the a particles because they are the majority component in the mixture. More importantly, this demixing propensity is not observed in the KAWCA case. We underline here that these concentration domains for the KALJ system resemble the behaviour discovered by Fischer³⁰ for supercooled ortho-terphenyl. Namely, anomalies in the structure factor at low k-values, which are not commensurate with the isothermal compressibility, are connected to the nucleation of nanometric structures.⁵⁹ Furthermore, our results agree with recent theoretical efforts demonstrating that the low k portion of the structure factor for polydisperse colloidal systems can be separated into a compressibility contribution and a term related to composition fluctuations.⁵

Fig. 4 Excess coordination number (N_αβ = ρ_βG^∞_αβ) as a function of temperature for both KALJ and KAWCA systems. N_ab close to zero corresponds to a preferential a–b effective interaction. Below the onset temperature of glassy dynamics upon cooling, N_bb gets close to zero for the KALJ system, indicating a growing preferential b–b effective interaction, ultimately leading to phase segregation.

We now focus on the isothermal compressibility (eqn (11)). In Fig. 5, we present a log–log plot of κ_Tvs. T, where it is apparent that the KALJ system is systematically more compressible than the KAWCA system at all temperatures considered here. Hence, it is again clear that small differences in the tail of the RDFs result in sizeable differences in their thermodynamic properties. Furthermore, a power-law behaviour κ_T = κ_T⁰T^−γ is apparent with γ = 0.46 ± 0.01 for the KALJ system and γ = 0.45 ± 0.01 for the KAWCA system. Below the onset temperature of glassy dynamics, both systems deviate from this power law and become comparatively less compressible in the deeply supercooled region. One would expect that, for a system undergoing a gas–liquid separation, compressibility substantially increases upon approaching the gas–liquid region. By contrast, the behaviour observed in Fig. 5 suggests that this is not the case. Finally, the existence of this power law, including the low-temperature deviations,⁵³ is somewhat similar (γ = 0.40 ± 0.01)⁶⁰ to the one observed experimentally in supercooled water.

Fig. 5 Bulk isothermal compressibility κ_T, calculated from eqn (11), as a function of temperature for KALJ and KAWCA systems (log–log representation). We observe that a power law relationship holds as κ_T = κ_T⁰T^−γ with γ = 0.46 ± 0.01 and 0.45 ± 0.01 for KALJ and KAWCA, respectively. The dashed lines are the corresponding power law fitting.

3.4 Chemical potential

Finally, we compute the excess chemical potential for both systems (Fig. 6) using the SPARTIAN method.³⁷ Recent calculations of the chemical potential for the KALJ system in the range of temperature 0.5ε/k_B < T < 1.0ε/k_B are in excellent agreement with our results.⁶¹ At the onset temperature of glassy dynamics, there is a transition between two regimes, reflecting the tendency for the system to minimise its free energy. The fact that the curves for the KALJ and the KAWCA systems are identical up to a constant factor is a consequence of writing the LJ potential energy as U_LJ ≈ U_WCA + U_Attractive. This expression lies at the foundation of perturbation theory that assumes that U_Attractive is very small compared to U_WCA. However, the sizeable difference in chemical potential (≈5ε) indicates that this approximation does not hold in this case.

Fig. 6 Difference of excess chemical potentials between species a and b for both, KALJ and KAWCA systems. The KAWCA system results were shifted by a constant in order to mach the lowest temperature T = 0.5ε/k_B, indicating that the potential energy can be approximated to U_LJ ≈ U_WCA + U_Attractive. A change in the behaviour with T, indicated by the dashed-grey lines, is apparent at the onset temperature of glassy dynamics. The inner plot shows the difference of chemical potential for KALJ and KAWCA (without shifting). Results for the KALJ system in the temperature range 0.5–1.0ε/k_B well compare with results available in the literature.⁶¹

Similarly to other thermodynamic properties like excess²⁰ and configurational entropy,¹⁹ isothermal compressibility and chemical potential results confirm that perturbation theory is not valid for the KALJ and KAWCA systems at ρ = 1.2/σ³ since attractive interactions induce non-perturbative structural effects. In the following section, we investigate these systems at higher density, namely ρ = 1.6/σ³, where we expect that repulsive interactions play an increasingly dominant role.^7,10

3.5 KALJ and KAWCA mixtures at ρ = 1.6/σ³

We perform a similar thermodynamic analysis for KALJ and KAWCA systems at ρ = 1.6/σ³. Our results show that density and concentration structure factors (Fig. 7) are nearly identical for both systems in the range of temperature considered. Structure factors in the limit k→0, in particular, show no evidence for the nucleation of long-range structures. Dynamical properties for these mixtures available in the literature^20,58,62 reveal that both systems exhibit similar structural and dynamical properties at this density. Therefore, we conclude that long-range concentration fluctuations might be closely connected to the significant mismatch between dynamical properties of the two systems at ρ = 1.2/σ³.

Fig. 7 Density and concentration structure factors, S_ρρ and S_cc, for the KALJ and KAWCA systems at a higher density (ρ = 1.6), in the range of temperature considered here.

Concerning the isothermal compressibility (Fig. 8), the two systems are essentially indistinguishable in the whole temperature range. As a reference, the onset temperature of glassy dynamics for KALJ and KAWCA systems at this density is close to 2.80ε/k_B.⁵⁴ The two systems hence behave similarly well below the onset temperature of glassy dynamics. This result highlights the dominant role played by attractive interactions in determining thermodynamic properties of high-density liquids.

Fig. 8 Bulk isothermal compressibility κ_T for the KALJ and KAWCA systems at a higher density (ρ = 1.6/σ³).

3.6 Crystallisation of the KAWCA system

In the last section, we investigate the crystallisation of the KAWCA system at ρ = 1.2/σ³. We further decrease the temperature down to T = 0.35ε/k_B. Fig. 9 shows snapshots of the system at T = 0.45ε/k_B (top) and T = 0.35ε/k_B (bottom). It is apparent from the figure that the system at T = 0.45ε/k_B appears like a miscible liquid. Conversely, the system at T = 0.35ε/k_B shows crystalline domains with a clear tendency for phase-segregation.

Fig. 9 Snapshot of the KAWCA system at T = 0.45ε/k_B (top) and T = 0.35ε/k_B (bottom).

Density and concentration structure factors (Fig. 10) enable us to validate this crystallisation-demixing scenario. In particular, we observe the splitting of the second peak of S_ρρ(k) at 1.7k₀ with k₀ ≈ 7.13/σ that indicates the presence of facet-sharing domains between neighbouring crystalline regions. Perhaps more interesting, the structure factors at k < 2/σ show a marked formation of long-range domains qualitatively similar to the ones present for the KALJ case below onset temperature.

Fig. 10 Density and concentration structure factors, S_ρρ and S_cc, for the KAWCA system in the temperature range 0.3ε/k_B < T < 0.45ε/k_B.

These results indicate that the crystallisation of the KAWCA system is a process driven by phase segregation. We note that similar behaviour to that observed for the low-temperature density and concentration structure factors for KAWCA has also been observed in polydisperse glass-forming systems.⁶³ Furthermore, recent GPU simulations report the crystallisation of the KALJ system⁵⁵ due to composition fluctuations similar to the ones investigated in this work. Indeed, our results also support previous claims pointing out demixing as a precursor for crystalisation in the modified KA model⁶⁴ and in liquid metals.⁶⁵

We conclude here that the presence of long-range concentration fluctuations is a qualitatively common feature for both KALJ and KAWCA systems. Perhaps more important, it is not directly related to gas–liquid coexistence present in the KALJ system. Indeed, in the crystalline state (T = 0.35ε/k_B), the KAWCA system exhibits a significant growth of concentration fluctuations, apparent in the lower panel in Fig. 10. The attractive interactions present in the KALJ system favour the nucleation of concentration fluctuations starting at relatively high temperatures (T = 1.00ε/k_B). Although the KAWCA system crystallises in our simulation timescale as its dynamics are considerably faster than for the KALJ system, our results suggest that both systems crystallise upon demixing by following a similar pathway.

4 Conclusions

We compute various thermodynamic properties of model supercooled liquids, with (KALJ) and without (KAWCA) attractive interactions at density ρ = 1.2/σ³. We aim at studying whether fluctuations in the tail of the two-body correlation function induce significant thermodynamic differences between the two systems. Density and concentration structure factors in the limit k→0 indicate that the KALJ system exhibits anomalous structural behaviour that we identify as the nucleation of long-range concentration domains. Conversely, the KAWCA system behaves like a normal liquid, with density and concentration structure factors decreasing monotonically. A finite-size Kirkwood–Buff analysis used to extrapolate to the k→0 limit confirms this picture. Differences in isothermal compressibilities and chemical potentials highlight the non-perturbative role of attractive interactions. Results of the crystallisation of the KAWCA system suggest that the anomaly, enhanced by the presence of attractive interactions, is a common feature of both models. All our results indicate that these long-range concentration fluctuations are not connected to gas–liquid coexistence, implying that demixing precedes crystallisation in both systems. Finally, upon increasing density (ρ = 1.6/σ³), where KALJ and KAWCA systems show similar dynamical properties, the KALJ anomaly disappears, and both systems exhibit nearly identical thermodynamic properties. Hence, we speculate that there might be a connection between large-scale concentration fluctuations and the significant dynamical slow down of the KALJ system in the deeply supercooled regime.

Author contributions

Atreyee Banerjee: conceptualization, methodology, formal analysis, writing – original draft Mauricio Sevilla: methodology, software, formal analysis, visualisation Joseph F. Rudzinski: investigation, writing – review & editing Robinson Cortes-Huerto: conceptualisation, methodology, writing – original draft, review & editing.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thank Kurt Kremer for his insightful discussions and his critical reading of the manuscript. They are also grateful to Pietro Ballone, Burkhard Dünweg, Smarajit Karmakar and Werner Steffen for their valuable feedback and suggestions. R. C.-H. thankfully acknowledge funding from SFB-TRR146 of the German Research Foundation (DFG). Open Access funding provided by the Max Planck Society.

References

1. E. Boattini, S. Marn-Aguilar, S. Mitra, G. Foffi, F. Smallenburg and L. Filion, Nat. Commun., 2020, 11, 5479.
2. H. Tong and H. Tanaka, Phys. Rev. Lett., 2020, 124, 225501.
3. J. Chattoraj and M. P. Ciamarra, Phys. Rev. Lett., 2020, 124, 028001.
4. F. P. Landes, G. Biroli, O. Dauchot, A. J. Liu and D. R. Reichman, Phys. Rev. E, 2020, 101, 010602.
5. L. Klochko, J. Baschnagel, J. P. Wittmer, O. Benzerara, C. Ruscher and A. N. Semenov, Phys. Rev. E, 2020, 102, 042611.
6. M. Leocmach and H. Tanaka, Nat. Commun., 2012, 3, 974.
7. L. Berthier and G. Tarjus, Phys. Rev. Lett., 2009, 103, 170601.
8. W. Kob and H. C. Andersen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 52, 4134.
9. J. D. Weeks, D. Chandler and H. C. Andersen, J. Chem. Phys., 1971, 54, 5237–5247.
10. L. Berthier and G. Tarjus, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 031502.
11. W. Götze and L. Sjögren, Z. Phys. B: Condens. Matter, 1987, 65, 415–427.
12. D. Coslovich, J. Chem. Phys., 2013, 138, 12A539.
13. G. M. Hocky, T. E. Markland and D. R. Reichman, Phys. Rev. Lett., 2012, 108, 225506.
14. S. Toxvaerd, Phys. Rev. E, 2021, 103, 022611.
15. F. Sciortino and W. Kob, Phys. Rev. Lett., 2001, 86, 648.
16. V. K. de Souza and D. J. Wales, J. Chem. Phys., 2008, 129, 164507.
17. M. P. Ciamarra, R. Pastore and A. Coniglio, Soft Matter, 2016, 12, 358–366.
18. J. F. Rudzinski, M. Radu and T. Bereau, J. Chem. Phys., 2019, 150, 024102.
19. A. Banerjee, S. Sengupta, S. Sastry and S. M. Bhattacharyya, Phys. Rev. Lett., 2014, 113, 225701.
20. A. Banerjee, M. K. Nandi, S. Sastry and S. M. Bhattacharyya, J. Chem. Phys., 2016, 145, 034502.
21. G. Adam and J. H. Gibbs, J. Chem. Phys., 1965, 43, 139–146.
22. E. D. Cubuk, S. S. Schoenholz, J. M. Rieser, B. D. Malone, J. Rottler, D. J. Durian, E. Kaxiras and A. J. Liu, Phys. Rev. Lett., 2015, 114, 108001.
23. M. K. Nandi and S. M. Bhattacharyya, Phys. Rev. Lett., 2021, 126, 208001.
24. J. G. Kirkwood and F. P. Buff, J. Chem. Phys., 1951, 19, 774–777.
25. M. Rovere, D. W. Heermann and K. Binder, J. Phys.: Condens. Matter, 1990, 2, 7009.
26. F. L. Román, J. A. White and S. Velasco, J. Chem. Phys., 1997, 107, 4635.
27. S. K. Schnell, T. J. Vlugt, J.-M. Simon, D. Bedeaux and S. Kjelstrup, Chem. Phys. Lett., 2011, 504, 199–201.
28. P. Krüger, S. K. Schnell, D. Bedeaux, S. Kjelstrup, T. J. Vlugt and J.-M. Simon, J. Phys. Chem. Lett., 2013, 4, 235–238.
29. R. Cortes-Huerto, K. Kremer and R. Potestio, J. Chem. Phys., 2016, 145, 141103.
30. E. Fischer, Phys. A, 1993, 201, 183–206.
31. A. Patkowski, T. Thurn-Albrecht, E. Banachowicz, W. Steffen, P. Bösecke, T. Narayanan and E. W. Fischer, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2000, 61, 6909–6913.
32. P. S. Salmon, R. A. Martin, P. E. Mason and G. J. Cuello, Nature, 2005, 435, 75–78.
33. Z. Zhang and W. Kob, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 14032–14037.
34. S. Toxvaerd and J. C. Dyre, J. Chem. Phys., 2011, 134, 081102.
35. S. J. Plimpton, J. Comput. Phys., 1995, 117, 1.
36. D. J. Evans and B. L. Holian, J. Chem. Phys., 1985, 83, 4069–4074.
37. M. Heidari, K. Kremer, R. Cortes-Huerto and R. Potestio, J. Chem. Theory Comput., 2018, 14, 3409–3417.
38. M. Praprotnik, L. Delle Site and K. Kremer, J. Chem. Phys., 2005, 123, 224106–224114.
39. M. Praprotnik, L. Delle Site and K. Kremer, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 73, 066701.
40. M. Praprotnik, L. Delle Site and K. Kremer, J. Chem. Phys., 2007, 126, 134902.
41. M. Praprotnik, L. Delle Site and K. Kremer, Annu. Rev. Phys. Chem., 2008, 59, 545–571.
42. R. Potestio, S. Fritsch, P. Español, R. Delgado-Buscalioni, K. Kremer, R. Everaers and D. Donadio, Phys. Rev. Lett., 2013, 110, 108301.
43. R. Potestio, P. Español, R. Delgado-Buscalioni, R. Everaers, K. Kremer and D. Donadio, Phys. Rev. Lett., 2013, 111, 060601.
44. U. R. Pedersen, T. B. Schrøder and J. C. Dyre, Phys. Rev. Lett., 2010, 105, 157801.
45. M. Heidari, K. Kremer, R. Potestio and R. Cortes-Huerto, Mol. Phys., 2018, 116, 3301–3310.
46. M. Heidari, K. Kremer, R. Potestio and R. Cortes-Huerto, Entropy, 2018, 20, 222.
47. F. Sedlmeier, D. Horinek and R. R. Netz, J. Am. Chem. Soc., 2011, 133, 1391–1398.
48. A. B. Bhatia and D. E. Thornton, Phys. Rev. B: Solid State, 1970, 2, 3004–3012.
49. P. Kumari, V. V. S. Pillai, D. Gobbo, P. Ballone and A. Benedetto, Phys. Chem. Chem. Phys., 2021, 23, 944–959.
50. W. Knoll and S. Steeb, Z. Naturforsch., 1978, 33a, 472–479.
51. B. W. van de Waal, J. Non-Cryst. Solids, 1995, 189, 118–128.
52. C. Desgranges and J. Delhommelle, Phys. Rev. Lett., 2018, 120, 115701.
53. K. H. Kim, A. Späh, H. Pathak, F. Perakis, D. Mariedahl, K. Amann-Winkel, J. A. Sellberg, J. H. Lee, S. Kim and J. Park, et al. , Science, 2017, 358, 1589–1593.
54. A. Banerjee, M. K. Nandi, S. Sastry and S. Maitra Bhattacharyya, J. Chem. Phys., 2017, 147, 024504.
55. T. S. Ingebrigtsen, J. C. Dyre, T. B. Schrøder and C. P. Royall, Phys. Rev. X, 2019, 9, 031016.
56. S. Sastry, Phys. Rev. Lett., 2000, 85, 590.
57. V. Testard, L. Berthier and W. Kob, Phys. Rev. Lett., 2011, 106, 125702.
58. L. Berthier and G. Tarjus, J. Chem. Phys., 2011, 134, 214503.
59. J. D. Stevenson and P. G. Wolynes, J. Phys. Chem. A, 2011, 115, 3713–3719.
60. A. Späh, H. Pathak, K. H. Kim, F. Perakis, D. Mariedahl, K. Amann-Winkel, J. A. Sellberg, J. H. Lee, S. Kim, J. Park, K. H. Nam, T. Katayama and A. Nilsson, Phys. Chem. Chem. Phys., 2019, 21, 26–31.
61. H. Vinutha and D. Frenkel, J. Chem. Phys., 2021, 154, 124502.
62. A. Banerjee and D. J. Wales, J. Chem. Phys., 2020, 153, 124501.
63. A. Ninarello, L. Berthier and D. Coslovich, Phys. Rev. X, 2017, 7, 021039.
64. U. K. Nandi, A. Banerjee, S. Chakrabarty and S. M. Bhattacharyya, J. Chem. Phys., 2016, 145, 034503.
65. C. Desgranges and J. Delhommelle, J. Am. Chem. Soc., 2014, 136, 8145–8148.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d2sm00089j

‡ These authors contributed equally to this work.

---