Glassy dynamics of sticky hard spheres beyond the mode-coupling regime

Abstract

Sticky hard spheres, i.e., hard particles decorated with a short-ranged attractive interaction potential, constitute a relatively simple model with highly non-trivial glassy dynamics. The mode-coupling theory of the glass transition (MCT) offers a qualitative account of the complex reentrant dynamics of sticky hard spheres, but the predicted glass transition point is notoriously underestimated. Here we apply an improved first-principles-based theory, referred to as generalized mode-coupling theory (GMCT), to sticky hard spheres. This theoretical framework seeks to go beyond MCT by hierarchically expanding the dynamics in higher-order density correlation functions. We predict the phase diagrams from the first few levels of the GMCT hierarchy and the dynamics-related critical exponents, all of which are much closer to the empirical observations than MCT. Notably, the prominent reentrant glassy dynamics, the glass–glass transition, and the higher-order bifurcation singularity classes (A₃ and A₄) of sticky hard spheres are found to be preserved within GMCT at arbitrary order. Moreover, we demonstrate that when the hierarchical order of GMCT increases, the effect of the short-ranged attractive interactions becomes more evident in the dynamics. This implies that GMCT is more sensitive to subtle microstructural differences than MCT, and that the framework provides a promising first-principles approach to systematically go beyond the MCT regime.

---

1 Introduction

The dynamics of supercooled liquids and the nature of the glass transition are still not fully understood from a theoretical perspective.^1,2 It is empirically well established that supercooling of a liquid is accompanied by a spectacular slowdown in dynamics over many orders of magnitude; at the same time, however, the disordered microstructure of the material remains very close to that of a normal liquid.³ One of the major challenges is to explain this emergent slow relaxation dynamics of glass-forming liquids on the basis of the material's static structure. At present, there are many theories that seek to answer this question, ranging from e.g. geometric-frustration-based theories⁴ to spin-glass-inspired approaches such as Random First Order Transition theory.⁵ More recently, machine learning has also been used to search for the link between structure and dynamics in dense disordered systems.^6–9 However, there is still no complete and formally exact theory of glassy dynamics that is founded entirely on first principles.

Mode-coupling theory (MCT), developed in the 1980s by Götze and co-workers,^3,10–13 is an approximate framework that is widely regarded as the only first-principles-based theory of the glass transition. For a review of this theory we refer the reader to ref. 3, 12 and 13. Briefly, among the main successes of MCT is its qualitative prediction of the non-trivial dynamical slowdown – as characterized by the intermediate scattering function F(k, t) for a given wavenumber k and time t – based on only small changes in the static structure factor S(k). For supercooled liquids, MCT typically predicts a two-step decay of F(k, t) via the so-called β- and α-relaxation processes, respectively, with a long plateau emerging at intermediate time scales. This plateau, which becomes longer upon deeper supercooling, is indicative of solid-like behavior and can also be understood in terms of the cage effect. Moreover, MCT predicts several scaling laws for the intermediate scattering function in the β- and α-relaxation regimes that are generally consistent with the behavior seen in experimental glass-forming systems. When the temperature further decreases to a certain value, MCT predicts a liquid–glass transition at which the long-time limit of the intermediate scattering function, also called the non-ergodicity parameter, jumps discontinuously from zero to a positive value. Many of the non-ergodicity parameters of different glass-forming materials are accurately predicted from MCT.^3,14 However, owing to its approximate nature, MCT is also known to become inaccurate or can even break down in certain cases. The main MCT approximation stems from the fact that the intermediate scattering function F(k, t), which in fact is a two-point density correlation function, is determined by a memory function in which four-point density correlation functions are required as the leading terms. In standard MCT, these four-point density correlation functions are approximated by the product of two F(k, t) functions at different wavenumbers.³ This approximation is generally not accurate enough and is even uncontrolled.¹³ Additionally, the four-point density correlation functions may contain information on non-trivial spatiotemporal density fluctuations related to dynamical heterogeneity, but these non-trivial correlations are inherently lost when applying the MCT factorization approximation.^15–18 Hence, a more accurate theory taking multi-point density correlators into account is needed.

Recently, a new theory called generalized mode-coupling theory (GMCT) has been developed based on the framework of MCT.^19–21 Within this theory, the uncontrolled approximation for the four-point density correlation functions is avoided by developing an exact dynamical equation for the four-point density correlators themselves; these are governed by a new memory function that is dominated by six-point density correlations, which in turn are controlled by eight-point correlators, et cetera. One may then apply a generalized MCT closure approximation at an arbitrary level in this hierarchy, thus delaying the uncontrolled factorization to a later stage. In this way, more accurate predictions for the relaxation dynamics can be obtained by systematically including higher-level correlation functions, which in fact establishes a hierarchy of coupled integro-differential equations. We note that, similar to MCT, the only material-dependent input for the theory is still the static structure factor S(k) for any given temperature and density.

GMCT has been well studied for hard-sphere systems with a hard-core repulsive interaction potential.^21–23 It has been found that GMCT indeed systematically improves the predictions from MCT. More specifically, the location of the liquid–glass transition point, which is usually underestimated in MCT, improves in a seemingly convergent manner as higher-order density correlations are included in the first few levels of the GMCT hierarchy.^19–22 Moreover, the qualitative MCT scaling laws near the liquid–glass transition in the β- and α-relaxation regimes are fully preserved in GMCT, but with quantitatively systematically improved exponents.^22,23 These triumphs of GMCT make it a promising first-principles-based theory for the glass transition, although it must be noted that the computational cost of a numerical GMCT calculation heavily increases with the number of included hierarchical levels, both in terms of computing time and required computer memory.²²

Despite the promising GMCT results for hard spheres, many real glass-forming materials are governed by more complex interaction potentials between the particles. Standard MCT has already been widely tested for different glass-forming materials, including model systems interacting through Lennard-Jones potentials,^24,25 amorphous metal alloys,^26,27 colloidal suspensions,^28–30 polymers,^31,32 vitrimers,³³ and systems interacting through many-body potentials.³⁴ In all of these systems, MCT qualitatively captures at least some key features of the glassy dynamics. A natural question is whether GMCT is also applicable for these more complex glass formers. To address this question, we will take the so-called sticky hard sphere model, which is characterized by a hard-core repulsive potential with a short-ranged attractive square well, as our system of study.

The sticky hard sphere model is the simplest model with both repulsive and attractive interactions but with rich dynamical behavior. In the last 20 years, many simulation^35–46 and experimental^47–56 studies have been performed for this model. It is widely accepted that there exists a complex reentrant dynamics when decreasing the temperature of the system at certain fixed values of the packing fraction and potential well width. Moreover, in contrast to the typical power-law decays of the intermediate scattering function near the critical liquid–glass point for hard spheres, an anomalous logarithmic decay has been found for sticky hard spheres. Although still debated recently,⁵⁷ the sticky hard-sphere reentrance is associated with two distinct types of glass: one is the so-called attractive glass, or bonded glass, related to the short deep potential well; the other is the repulsive glass, or non-bonded glass, similar to the hard-sphere glass. In MCT, all the above properties have been successfully predicted. Moreover, MCT predicts a glass–glass transition and even more complex transitions called A₃ and A₄ singularities, which can also explain the logarithmic decay of the intermediate scattering function.⁵⁸ However, we point out that the predictions from MCT are all qualitatively but not quantitatively correct. For example, regarding the phase transition at different well widths, the critical packing fractions are all underestimated by MCT. Furthermore, the A₃ and A₄ singularity points predicted by MCT have to be rescaled to match the simulation or experimental results.^39,49 In this work, we revisit the sticky hard sphere model in the context of GMCT, with the aim to test the applicability of GMCT for systems with both attractive and repulsive interactions, and to establish to what extent first-principles predictions for glassy sticky hard spheres can be improved.

The outline of this paper is as follows. We first introduce the microscopic GMCT framework and the related analytical scaling laws. We also briefly describe the sticky hard sphere model and the numerical details for the application of GMCT to this model. Subsequently, we present the rich phase diagrams of sticky hard spheres predicted from higher-order GMCT and compare them to literature to demonstrate the quantitative improvement attained by GMCT. Next, we discuss the non-ergodicity parameters at the critical points for both the liquid–glass transitions and the glass–glass transition, reflecting the different mechanisms of the two types of glass. The exponent parameters that characterize the dynamics near the corresponding critical points are also provided. Finally, we show the relaxation dynamics near the higher-order A₃ and A₄ singularities, demonstrating that the pronounced logarithmic decay of the dynamical density correlation functions is preserved within higher-order GMCT. We conclude with some critical remarks and perspectives for future work.

2 Theory

2.1 Generalized mode-coupling theory

We first review the microscopic GMCT equations of motion²¹ and several scaling laws derived from them.²³ The microscopic dynamical information of a glass-forming material is assumed to be encoded in the 2n-point density correlation functions F_n(k₁,…,k_n, t), defined as

F_n(k₁,…,k_n, t) = 〈ρ_−k1(0)…ρ_−kn(0)ρ_k1(t)…ρ_kn(t)〉, (1)

where is a collective density mode at wavevector k and time t, N_p is the total number of particles, the angle brackets denote an ensemble average, and the label n (n = 1,…,∞) specifies the level of the GMCT hierarchy. Note that when n = 1, F₁(k, t) is the intermediate scattering function. In the overdamped limit, the GMCT equations readwhere ν_n is an effective friction coefficient,are 2n-point static density correlation functions which we approximate as a product of static structure factors S(k_i), andwith k_B and m denoting the Boltzmann constant and the mass of the particle, respectively. For the memory functions we havewhere ρ is the number density, δ_i,j is the Kronecker delta function, and the V_q,ki−q are so-called static vertices. These vertices are defined as

V_q,k−q = (k·q)c(q) + k·(k − q)c(|k − q|), (6)

where c(q) denotes the direct correlation function.⁵⁹ The direct correlation function is related to the static structure factor as c(q) ≡ [1 − 1/S(q)]/ρ. Therefore, the static structure factor S(q) serves as the only material-dependent input to the GMCT equations. While the above GMCT equations still rely on several approximations, including factorization of all static multi-point correlators [eqn (3)], the neglect of so-called off-diagonal dynamic multi-point correlators,^19,21,60 and the treatment of all dynamics with the normal Liouvillian operator rather than the projected one,^12,61 we emphasize that the theory is still purely first-principles-based.

In principle, one should expand the GMCT equations up to n → ∞ to obtain the correct dynamical density correlations. However, in practice we can only solve the GMCT equations up to a finite level N < ∞, which necessitates the use of a closure approximation for the last included level N. There are two kinds of GMCT closures that have been well-studied before.^21,22 One kind is the so-called mean-field (MF) closure, which approximates the highest-order correlator F_N as a product of lower-order dynamical density correlation functions. In this paper, we only focus on one such example, for N > 2,

where {k_j}^(N−1)_j≠i represents the N − 1 wavenumbers in {k₁,…,k_N} except the k_i, and the permutation invariance of all wavenumber arguments {k₁,…,k_N} in M_N(k₁,…,k_N, t) has been taken into account. This closure, denoted as MF-N[(N − 1)¹1¹], is exactly the same as the one used in the original paper deriving the microscopic dynamical GMCT,²¹ and is qualitatively equivalent to the one used in the paper analyzing the scaling laws for GMCT.²³ For N = 2, the closure is trivially satisfied with F₂(k₁, k₂, t) = F₁(k₁, t)F₁(k₂, t), which is in fact the standard-MCT closure³ that we denote here as MF-2.¹² The other kind of GMCT closure, named exponential (EXP-N) closure, is a simple truncation of the hierarchy such that F_N(k₁,…,k_N, t) = 0, which leads to F_N−1 ∼ exp(−t/τ_N−1) with τ_N−1 = ν_N−1/(S_N−1⁻¹J_N−1). It has been demonstrated numerically that the mean-field and exponential closures provide an upper and lower bound for the relaxation dynamics, respectively.^21,22,62 In the limit N → ∞, the differences of the relaxation dynamics from these two kinds of closures are expected to disappear. From previous numerical results for hard spheres,^21,22 it is known that mean-field closures require relatively few GMCT levels to predict the emergence of strongly glassy behavior; in this work we therefore only focus on the results with mean-field closures to study the glass transition.

Eqn (2), (5), and (7) define a unique time-dependent solution F_n(k₁,…,k_n, t) for n ≤ N.⁶³ Before we go to the specific model and the numerical calculation, let us discuss some universal properties of the GMCT solutions which can be obtained from mathematically analyzing the equations.

Firstly, using the Laplace transform and the final value theorem, the long-time limit of the multi-point density correlation function can be calculated via

F_n(k₁,…,k_n) = S_n(k₁,…,k_n) − [S_n⁻¹(k₁,…,k_n) + J_n⁻²(k₁,…,k_n)M_n(k₁,…,k_n)]⁻¹, (8)

where represents the long-time limit of the memory function. When n < N,and the MF closure eqn (7) becomeswhere for notational convenience we have introduced the symbol [scr M, script letter M] to denote the memory functionals. The solutions of eqn (8)–(10) can exhibit non-trivial singularities that may be analyzed using the Jacobian of these equations. The Jacobian is equivalent to 1 − C, where C is the stability matrix of eqn (8)–(10) which satisfies

H = CH.

Here H is a vector in the form of [{H₁(k₁)}_∀k1, {H₂(k₁, k₂)}_∀k1,∀k2,…,{H_n(k₁,…,k_n)}_{∀k1,…,∀kn}]^T. Explicitly, for n < N,

H_n(k₁,…,k_n) = [S_n(k₁,…,k_n) − F_n(k₁,…,k_n)]²J_n⁻²(k₁,…,k_n) × [scr M, script letter M]_n[H_n+1]. (11)

and

H_N(k₁,…,k_N) = [S_N(k₁,…,k_N) − F_N(k₁,…,k_N)]²J_N⁻²(k₁,…,k_N) × {[scr M, script letter M]_N[M_N−1,H₁] + [scr M, script letter M]_N[[scr M, script letter M]_N−1[H_N], F₁]}. (12)

Since all elements of C are non-negative and the positive elements between different levels ensure that the matrix is irreducible, there is a non-degenerate maximum eigenvalue E of the matrix C according to the Perron–Frobenius theorem.⁶⁴ The singularities, also called bifurcation points,⁶⁵ of the long-time limit of multi-point density correlation functions occur if the Jacobian is a singular matrix, i.e., if the matrix C has an eigenvalue E = 1. In principle, there could be very complex singularities, including a family labeled A_l, l = 2, 3,…, which are topologically equivalent to the bifurcation singularities of the real roots of real polynomials of degree l.⁶⁵ The simplest singularity is the so-called A₂ bifurcation point, which most reported liquid–glass transitions in the literature belong to. In this class, the normalized long-time limits of the dynamic density correlation functions (also usually called non-ergodicity parameters) f_n(k₁,…,k_n) ≡ F_n(k₁,…,k_n)/S_n(k₁,…,k_n) display a simple bifurcation: for liquid states f_n(k₁,…,k_n) = 0, while for glass states f_n(k₁,…,k_n) > 0. Such an A₂ singularity is found in e.g. the Percus–Yevick hard-sphere system,²³ for which the only control parameter of the system is the packing fraction φ; the lowest packing fraction for the glass state, φ^c, is referred to as the critical point of the glass transition. More generally, if there are p > 1 control parameters in the system, f_n(k₁,…,k_n) is a function of these p variables and hence in principle there could exist singularities A_l where l > 2. These are the cases we will pay more attention to in this paper.

At each singularity point, i.e., when E = 1, the left eigenvector ê_n(k₁,…,k_n) and right eigenvector e_n(k₁,…,k_n) of the C matrix can be obtained with conventions

andwhere represents the summation over all possible wavenumbers k₁,…,k_m for level m and the superscript c refers to the critical-point values. The above eigenvectors are used to characterize every transition point by a single number λ, defined asFor bifurcation singularities of the A₂ type, 0 < λ < 1. This λ becomes unity at the end points of the A₂ transition curve (or surface), which could be A₃ or even higher order singularities.^58,65

Let us briefly summarize some leading order scaling laws near A₂ singularities within the framework of GMCT, both for the normalized long-time limits of the density correlators and for the time-dependent ones [i.e., the solutions of eqn (2), (5) and (7)]. Details of these scaling laws can be found in ref. 23. Here the control parameter is the relative distance to a A₂ critical point which, for example, is measured by the packing fraction in hard spheres, ε = (φ − φ^c)/φ^c.

• If ε > 0 and |ε| ≪ 1, the system is in the glass state and the non-ergodicity parameters scale with as

where

h_n(k₁,…,k_n) = e_n(k₁,…,k_n)/S_n(k₁,…,k_n). (16)

• If |ε| ≪ 1, the β-relaxation regime can be characterized by a unique timescale τ_β ∼ |ε|^−1/2a. In this regime, the dynamical normalized density correlation functions f_n(k₁,…,k_n, t) ≡ F_n(k₁,…,k_n, t)/S_n(k₁,…,k_n) obey a time-wavenumber factorization property such that

f_n(k₁,…,k_n, t) = f^c_n(k₁,…,k_n, t) + h_n(k₁,…,k_n, t)G(t). (17)

Here G(t) is the so-called β-correlator and it scales with ε as , with the subscript ± denoting the sign of ε. For the limiting case that ε → 0⁻, i.e., on the liquid side of the transition, the β-correlator satisfiesandwhich describes the decay towards and away from the plateau, respectively. The latter equation for the late β-relaxation regime is also called the von Schweidler law. The exponents a and b can be determined by λ in eqn (14)via

λ = Γ(1 − a)²/Γ(1 − 2a) = Γ(1 + b)²/Γ(1 + 2b). (20)

The t₀ in eqn (18) is a time scale characterizing transient dynamics, while τ in eqn (19) is the α-relaxation time. Within GMCT under MF closures, we have τ ∼ |ε|^−γ with

γ = 1/2a + 1/2b. (21)

• If ε < 0 and |ε| ≪ 1, the system is in the liquid state and the dynamical normalized density correlation functions in the α-relaxation regime obey a time-density superposition principle

f_n(k₁,…,k_n, t) = f̃_n(k₁,…,k_n, t/τ). (22)

All the above scaling laws near A₂ singularities are very similar to those predicted from standard MCT, i.e., GMCT with MF closure level N = 2. However, we point out that in higher-order GMCT (i) all the exponents (a, b, λ, γ, etc.) become closure-level dependent, and (ii) the dynamics of multi-point density correlators can also be predicted. In particular, the latter might contain spatiotemporal information related to dynamical heterogeneity.

For higher-order singularities A_l with l > 2, a major difference with the scaling laws above is that near such singularities, to leading order in , the β-correlators G(t) may follow a logarithmic decay,^58,66i.e.

The technique to obtain this scaling law is similar to the one for the g_±(t/τ) near A₂ singularities (see ref. 23 for more details). The difference is that when λ = 1, the equation of the β-correlator becomes

σ|ε| − s²G²(s) + s[script L][G²(t)] = 0 (24)

where σ is a constant that depends on the system and on the GMCT closure. If σ < 0, the solution is eqn (23), similar to standard MCT.⁶⁷ Therefore, the emergence of logarithmic decay can be used as an indicator for the existence of high-order singularities, both in simulation and experiment.³⁹ Another difference near higher-order singularities is that the superposition principle in the α-relaxation regime need not be applicable anymore, which is related to the breakdown of the von Schweidler law in the late β-relaxation regime. In Section 3, we will discuss more subtle effects for control parameters p = 3 where A_l (l = 2, 3, 4) singularities exist.

2.2 Sticky hard sphere model

We apply our GMCT framework to a system composed of sticky hard spheres. The pair interaction potential of this system includes both a repulsive and attractive part, as shown in Fig. 1. Within the particle diameter d, the value of the potential is u(r) = +∞, which accounts for hard-core repulsion to prevent particle overlap. Outside the hard-core region there is a small square potential well with width Δ ≪ d, i.e., for d < r < d + Δ we have u(r) = −u₀, with u₀ > 0. Thus, if the distance of two particles falls within this range, the particles will attract each other to be sticky. When r > d + Δ, u(r) = 0.

Fig. 1 Pair potential u(r) for sticky hard spheres as a function of interparticle distance r.

There are three independent parameters of this model. One is the packing fraction φ = πρd³/6. The second is the reduced temperature θ, defined as the inverse of the attraction well depth over thermal energy, i.e., θ = k_BT/u₀. The third one is the relative width of the well, δ = Δ/d. Note that when u₀ = 0 (θ → ∞), this sticky hard sphere model reduces to the hard sphere model and the only parameter is the packing fraction φ. We calculate the structure factors S(k) from the pair potential using the Ornstein–Zernike equation together with the mean-spherical approximation.⁵⁹ More specifically, we use the approach motivated by Baxter's theory in ref. 58 to numerically obtain the structure factors in the three-dimension control parameter space (φ, θ, δ); these structure factors subsequently serve as our GMCT input.

For the GMCT calculations, we numerically solve eqn (2)–(7) for the time-dependent dynamics of the multi-point density correlation functions, and eqn (8)–(10) for the corresponding long-time limits, for MF closure levels N = 2, 3, 4. The wavevector-dependent integrals over q in the memory functions are approximated as a double Riemann sum⁶⁸ with an equidistant wavenumber grid of N_k = 500 points that ranges from k_d = 0.2 to k_d = 200.2. Such a fine wavenumber grid is necessary to reach convergence, in particular for low values of θ. For the time-dependent solutions we use the integration algorithm described by Fuchs et al.,⁶⁹ starting with a time step size of Δt = 10⁻⁶ that is subsequently doubled every 32 points. In all our GMCT calculations we set k_B/m = 1 and the effective friction coefficients ν_n = 1 for 1 ≤ n ≤ N.

3 Results and discussion

3.1 Phase diagram

We first report the GMCT phase diagrams in the (φ, θ) plane for several fixed values of δ (i.e., fixed potential well widths), as shown in Fig. 2. Let us start by discussing the N = 2 results, which are also fully consistent with the standard-MCT results of approximately 20 years ago.⁵⁸ For all values of δ considered (δ = 0.09, 0.06, 0.0465, 0.03), a common phenomenon is that for a given reduced temperature θ, there is always a liquid–glass transition point φ^c. That is, for φ > φ^c the non-ergodicity parameters are larger than zero and the system is in the glass state, while for φ < φ^c the non-ergodicity parameters are all 0, representing a liquid state. Thus, for an arbitrary temperature and potential well width, one may always find an ergodicity-breaking transition by going to sufficiently high packing fractions. However, for a given packing fraction, the influence of the reduced temperature is more subtle and depends on the width of the potential well. At relatively large widths, such as δ = 0.09, the picture is fairly trivial: the lower the temperature, the more glassy the system. This is indeed similar to systems with a relatively large attractive potential well, such as the Lennard-Jones potential.²⁴ Conversely, when the width of the attractive well decreases to e.g. δ = 0.06, the scenario changes: at high temperatures the system is in the glass state, at lower temperatures the glass melts into a liquid state, and finally at even lower temperatures the system re-enters the glass phase. This is the famous reentrant phenomenon of systems with a short attractive well, and can be found regardless of the shape of the well.⁷⁰ Physically, the existence of the glass phase at high reduced temperatures (i.e., the repulsive glass) can be understood in terms of the cage effect, similar to the case of hard spheres. That is, the MCT modes with wavelengths close to the average distance of neighbors play the dominant role. By contrast, the glass phase at low temperatures (i.e., the so-called attractive glass) arises from particle bonding, with the MCT modes at wavelengths comparable to the width of the attractive potential well becoming more important. The liquid state in between can be explained by the following two arguments:^49,58 (i) compared to the repulsive glass, the attraction causes the separation of two particles to be smaller so that the average size of the free space increases, facilitating longer-distance particle motion and thus faster structural relaxation; (ii) compared to the attractive glass, the “stickiness” between particles is weaker, hence the particles have a higher probability to break their bonds and move relatively fast away from their initial positions.

Fig. 2 Phase diagrams of sticky hard spheres, obtained from GMCT with MF closure levels N = 2, 3, 4. The glass transition curves are shown as a function of packing fraction φ and reduced temperature θ for several values of the relative attraction well width δ. The A₃ end points are marked by squares and the A₄ end points are marked by stars. The A₃ and A₄ points also exist for GMCT closure level N = 4, but their positions could not be accurately determined due to limiting computing power.

The two-glass picture becomes a glass–glass transition for very small potential well widths. For example, when δ = 0.03, in addition to the reentrant liquid–glass transition, there is another A₂ glass–glass transition predicted from MCT (black solid line from the crossing point near φ = 0.535 to the black square marker in Fig. 2). The existence of this glass–glass transition is still debated, most recently based on swap Monte Carlo simulations,⁵⁷ but from the theoretical point of view, this transition indeed exists since the maximum eigenvalue E of the C matrix introduced in the previous section is unity at the glass–glass transition curve, with λ < 1. Moreover, when δ = 0.03, along the A₂ glass–glass transition curve, λ increases with increasing φ until λ → 1 (black solid line in Fig. 9) indicating an end point with an A₃ singularity (black square in Fig. 2). When increasing the width δ, the length of the glass–glass transition curve shrinks and there will be a δ* at some φ* and θ* when the curve vanishes. The point (φ*, θ*, δ*) is the A₄ singularity (black star in Fig. 2).

While standard MCT already successfully captures the two competing effects, i.e., caging and bonding, underlying vitrification in sticky hard spheres, the MCT-predicted phase diagram is not quantitative accurate. Indeed, as in many other systems, MCT generally tends to overestimate the critical packing fraction of the liquid–glass transition point. Let us therefore turn to the higher-order GMCT results in Fig. 2. We find that overall the shape of the phase diagram under higher-order MF closures remains similar to the MCT case, with all the liquid–glass transitions, glass–glass transitions and A_l (l = 2, 3, 4) singularities preserved. However, when the MF closure level N in GMCT increases, all the phase transition curves for a given δ move to larger packing fraction values. Note also that at high reduced temperatures θ, the critical liquid–glass transition points at different widths δ will converge to the corresponding GMCT critical points for hard spheres. It has been shown that the GMCT-predicted critical packing fraction φ^c for hard spheres shifts toward higher values in a seemingly convergent manner, i.e., the difference φ^c(N) − φ^c(N − 1) becomes increasingly smaller as N increases.²² The phase transition curves in Fig. 2 also imply this uniform convergence for any given δ and θ. This is arguably one of the biggest triumphs of GMCT: by systematically including higher levels of density correlators in the GMCT hierarchy, the critical point can be systematically improved, even in the presence of competing repulsive and attractive particle interactions.

In practice, the phase transition lines from MCT are usually shifted or rescaled to compare with simulation or experimental results. In Fig. 3 we present the shifted GMCT liquid–glass transition curves for N = 2, 3, 4 relative to the corresponding critical points φ^N(θ = 1.4) at each δ. We have chosen θ = 1.4, which is approximately equivalent to the hard-sphere limit, as a reference point for convenience. For each δ value, it is obvious that the shifted phase transition curves at different closure levels N do not collapse to a single curve, especially at low temperatures θ. This indicates that the GMCT-predicted phase diagrams under higher closure levels are non-trivial in the sense that they cannot be obtained by simply shifting the MCT results. At high δ, such as δ = 0.09, increasing N leads to larger shifted critical packing fractions φ^N_shift of the liquid–glass transition (dotted lines in Fig. 3). Moreover, the difference between the φ^N_shift with N = 3 (or N = 4) and the one with N = 2, i.e., Δφ^N=3_shift ≡ φ^N=3_shift − φ^N=2_shift, is larger at low θ. This implies that the higher order density correlations included in GMCT may contribute more to the dynamics when the attractive potential well is deeper. However, when δ decreases to e.g. δ = 0.03 (solid lines in Fig. 3 and 4), the relation between the temperature θ and the shifted liquid–glass transition lines becomes more complex. At low temperatures, when N increases, the shifted critical packing fractions φ^N_shift increase as well, which is similar to the δ = 0.09 case. However, at higher temperatures, φ^N_shift decreases as the level N increases. There thus exists a crossing point between any two shifted phase transition curves with different closure levels N; this crossing point (θ ≈ 0.22, φ^N_shift ≈ 0.015) is close to the ‘convex’ part of the phase diagram, and more specifically at the high-temperature side of the convex region. The Δφ^N_shift shown in the inset of Fig. 4 demonstrates this opposite trend for relatively high and low temperatures.

Fig. 3 Shifted liquid–glass transition lines of sticky hard spheres for different values of δ and different GMCT MF closure levels. The reference critical packing fraction for each closure level N is chosen at the reduced temperature of θ = 1.4, which is high enough to be regarded as the liquid–glass transition point for hard spheres.

Fig. 4 Shifted liquid–glass transition lines of sticky hard spheres for different GMCT MF closure levels at a relative attraction well width of δ = 0.03. The results are the same as in Fig. 3 but enlarged for the low reduced-temperature region. The inset shows the difference between the higher-order GMCT and MCT curves of the main panel, i.e. Δφ^N_shift ≡ φ^N_shift − φ^N_shift^{= 2}. The black curve in the inset represents N = 2 and the values are 0 by definition; the blue and red curves represent N = 3 and N = 4, respectively.

This subtle, complex effect at small δ, which in our theory emerges only when going beyond the standard-MCT level (N > 2), is in fact also seen in earlier empirical studies of sticky hard spheres dating back to at least 20 years ago, but it appears to have been ignored thus far, or was perhaps regarded as a small (numerical) error. For example, in Fig. 1 of the simulation work by Sciortino, Tartaglia, and Zaccarelli,³⁹, the simulation data at high k_BT/u₀ have a smaller packing fraction than the shifted MCT predicted liquid–glass transition curve, while the data points at low k_BT/u₀ are at higher packing fraction values compared to the predicted curve. This trend fully agrees with our predictions from GMCT with higher closure levels N. Another example is from the experimental work by Pham et al.⁴⁹ Although the experimental interaction potential is not the same as the one shown in Fig. 1, both contain a short-ranged attractive well, so that we may expect the effect of higher-order GMCT compared to standard MCT to be similar. Indeed, at high polymer concentrations, which corresponds to low reduced temperatures in our present work, the rescaled phase transition curve in Fig. 1 of ref. 49 has lower packing fractions than the experimental data, while at low polymer concentrations, the critical packing fractions are overestimated by the MCT predictions after rescaling. Notice that the crossing point of the shifted MCT predicted curve and the phase transition curve drawn from data in both simulation and experiment is at the side of high reduced temperatures (high k_BT/u₀ in the simulation and low polymer concentration c_p in the experiment), which further demonstrates that the improvement of the liquid–glass transition phase diagram made by higher-order GMCT is reasonable.

3.2 Non-ergodicity parameters

We now turn to the long-time limits of the normalized multi-point density correlation functions, i.e., the non-ergodicity parameters f_n(k₁,…,k_n). We will focus on a very low attractive potential well width of δ = 0.03, for which the sticky hard sphere model exhibits both glass–liquid–glass reentrance and a glass–glass transition. Let us first consider the liquid–glass transition. Fig. 5 shows the critical non-ergodicity parameters f^c₁(k) and f^c₂(k,k) at two different temperatures, one at the critical repulsive glass state (θ = 0.2) and the other at the critical attractive glass state (θ = 0.13). It can be seen that the repulsive glass state is similar to the hard-sphere case in the sense that f^c₁(k) is strongly modulated by the shape of S(k) and the localization length is around the size of one particle diameter.^22,58,71,72 Conversely, the non-ergodicity parameters of the attractive glass are much higher than those of the repulsive glass, and they decay to zero over a much wider range of wavenumbers k.^58,71,72 This indicates that particle bonding at small distances indeed plays a more important role in the attractive glass phase. Notably, we observe this phenomenology for all considered GMCT closure levels N. Hence, even though higher-order GMCT predicts quantitatively different values of the non-ergodicity parameters, the theory clearly distinguishes between the two types of glasses in a similar manner as standard MCT.

Fig. 5 Non-ergodicity parameters f^c₁(k) and f^c₂(k, k) as a function of wavenumber k at the liquid–glass transition point for different GMCT MF closure levels at δ = 0.03. At θ = 0.2, the critical packing fractions are φ^c = 0.53571, 0.55382, 0.56667 for closure levels N = 2, 3, 4, respectively; at θ = 0.13, φ^c = 0.46714, 0.50959, 0.53817 for N = 2, 3, 4.

Let us further quantify the difference between MCT and higher-order GMCT by careful inspection of Fig. 5. For the lowest-order non-ergodicity parameters f^c₁(k), we find that increasing N generally leads to an increase in f^c₁(k), similar to what has been reported previously for GMCT of hard spheres (see Fig. 1 of ref. 23). The only small exception to this N-dependent trend for sticky hard spheres is observed at low temperature (θ = 0.13) and at small wavenumbers k_d < 7.4, where k₀d = 7.4 corresponds to the first peak of S(k). Overall, at the level of the two-point density correlators, the effect of higher-order GMCT corrections to MCT thus appears to be qualitatively the same for both the attractive and repulsive glass. However, the influence of N on the diagonal four-point density correlators f^c₂(k, k) is more complex. Indeed, for large wavenumbers k (k_d > 15 for θ = 0.2 and k_d > 60 for θ = 0.13), f^c₂(k, k) monotonically decreases with increasing N for both glass types, while for smaller k, the ordering of the f^c₂(k, k) curves changes under different closure levels. Moreover, the number of such crossovers is different for different temperatures. Hence, in general we can conclude that within GMCT f^c₂(k, k) ≠ f^c₁(k) × f^c₁(k), a result that may hint at the presence of dynamical heterogeneities.

The off-diagonal correlators f^c₂(k₁, k₂) at θ = 0.2 and θ = 0.13 are shown in Fig. 6(a) and Fig. 7(a–d), respectively. We find that the pattern of f^c₂(k₁, k₂) is qualitatively rather similar for all closure levels N, with a shape that is modulated by S(k). For MCT this shape is in fact trivial since then f^c₂(k₁, k₂) = f^c₁(k₁) × f^c₁(k₂); for higher-order GMCT, our repulsive-glass results are also consistent with the GMCT results for hard spheres.²¹ To see the effect of N for sticky hard spheres more clearly, we plot the difference between f^c₂(k₁, k₂) for N > 2 and N = 2 in Fig. 6(b, c) and 7(b, c, e, f). It is clear that as N increases, the GMCT deviations of f^c₂(k₁,k₂) with respect to MCT increase both in absolute value and in wavenumber range. However, comparing the results for the repulsive and attractive glass, we see that for the attractive glass phase the deviations extend over a much broader window of wavenumbers. This difference can also be seen if we plot f^c₂(k₁, k₂) − f^c₁(k₁) × f^c₁(k₂) for a given state point, i.e. the absolute difference of the four-point non-ergodicity parameters with respect to Gaussian factorization (see the figures in the appendix). We speculate that the difference may stem from the different physical mechanisms, i.e., caging versus bonding, underlying the vitrification process, and it might also point toward different degrees of dynamical heterogeneity in the two glassy phases. This hypothesis, however, still requires further study and should be investigated in future work.

Fig. 6 (a) Non-ergodicity parameters f^c₂(k₁, k₂) as a function of wavenumbers k₁ and k₂ at the critical packing fraction for closure level N = 2 (i.e., standard MCT) at θ = 0.2 and δ = 0.03. The patterns of f^c₂(k₁, k₂) from N = 3 and N = 4 are similar. (b) The difference of the non-ergodicity parameters f^c₂(k₁, k₂) of closure level N = 3 with respect to N = 2. (c) The difference of the non-ergodicity parameters f^c₂(k₁, k₂) of closure level N = 4 with respect to N = 2. The packing fractions are the same as in Fig. 5.

Fig. 7 (a) Non-ergodicity parameters f^c₂(k₁, k₂) as a function of wavenumbers k₁ and k₂ at the critical packing fraction for closure level N = 2 (i.e., standard MCT) at θ = 0.13 and δ = 0.03. The patterns of f^c₂(k₁, k₂) from N = 3 and N = 4 are similar. (b) The difference of the non-ergodicity parameters f^c₂(k₁, k₂) of closure level N = 3 with respect to N = 2. (c) The difference of the non-ergodicity parameters f^c₂(k₁, k₂) of closure level N = 4 with respect to N = 2. Panels (d–f) represent the zoomed in small-wavenumber regime of panels (a–c), respectively. The packing fractions are the same as in Fig. 5.

To demonstrate the existence of the glass–glass transition within GMCT, we plot the non-ergodicity parameters f₁(k) as a function of θ for different closure levels N. Fig. 8 shows f₁(k) for three selected packing fractions and for three representative wavenumbers k; we mention that the wavenumber dependence of f₁(k) at both sides of the glass–glass transition line is similar to the result of Fig. 5. For all closure levels N and all considered wavenumbers k, we can see that f₁(k) initially decreases as θ increases, until a discontinuous drop occurs; this point corresponds to the A₂ singularity of the attractive-to-repulsive-glass transition. In this process, f₁(k) follows the same scaling law as eqn (15). As θ further increases within the repulsive glass phase, f₁(k) will decrease further for all N (see e.g. the open triangles in Fig. 8). However, after crossing a certain temperature, f₁(k) is seen to increase with θ and finally approaches the hard-sphere limit. This effect is essentially the counterpart of the glass–liquid–glass reentrance within the glassy region. The non-trivial variation of f₁(k) with θ near but above the critical temperature is a precursor phenomenon of the nearby A₃ singularity. Therefore, we conclude that for all GMCT mean-field closure levels N, the glass–glass transition and the A₃ singularity are preserved but with quantitatively changed positions in the three-dimension parameter space of sticky hard spheres.

Fig. 8 Non-ergodicity parameters f₁(k) as a function of the reduced temperature θ for fixed δ = 0.03 and fixed packing fraction φ. The packing fraction values are selected in such a way that the path along θ crosses the glass–glass transition curve for each corresponding GMCT closure level N. For N = 2, the packing fraction is φ = 0.54; for N = 3, φ = 0.556; for N = 4, φ = 0.57. Filled and open symbols represent the attractive and repulsive glass, respectively.

Finally we discuss the exponent parameter λ, which plays a crucial role in the context of A_l singularities. Fig. 9 shows λ as a function of θ for the phase transition lines under MF closure levels N = 2 and N = 3. Let us first focus on the standard MCT (N = 2) results. For relative high values of δ, such as δ = 0.09, 0.06, λ shows a single maximum. As δ decreases, the peak value of λ increases and the corresponding position of θ decreases. However, for very small δ, e.g. δ = 0.03, there is another line, representing the glass–glass transition, above the peak. This line ends at the A₃ singularity with λ = 1. Therefore there must be a δ between 0.06 and 0.03, which we estimate to be δ ≈ 0.0465, for which the glass–glass transition curve shrinks to zero and the peak lies at λ = 1, which defines the A₄ singularity. In fact, the higher-order singularity points shown in Fig. 2 can be found with the guide of the λ curves, as mentioned before. For GMCT with N = 3, the shapes of the λ curves are similar to MCT but quantitatively different. The peak of λ(θ) at δ = 0.0465 for N = 3 is 0.97 ≈ 1, and hence we conjecture the N = 3 A₄ point to exist for a δ value close to or slightly lower than 0.0465. From the values of λ when λ < 1, we can also calculate the corresponding exponents a, b and γ by virtue of eqn (19) and (20), i.e., the scaling laws of the relaxation dynamics near A₂ singularities. Notice that at large θ, λ must go to the value of the hard-sphere limit, i.e., λ = 0.73 for N = 2 (thin black line in Fig. 9) and λ = 0.78 for N = 3 (thin blue line in Fig. 9), for any value of δ. This puts a physical constraint on possible rescaling approaches; in particular, for fixed δ, the curves of λ(θ) for different closure levels N cannot collapse to one curve by simply shifting θ. This is a further manifestation of the non-trivial effects imposed by higher-order GMCT on the glassy dynamics of sticky hard spheres, consistent with the subtle differences in the phase diagrams of Section 3.1.

Fig. 9 Exponent parameter λ for points on the transition lines for GMCT MF closure levels N = 2 and N = 3. The thin horizontal black and blue lines indicate the asymptotic hard-sphere values for N = 2 (λ = 0.73) and N = 3 (λ = 0.78), respectively.

In view of the above findings, let us also make a general remark on the rescaling procedure widely used in standard MCT. For sticky hard spheres, it is in a sense accidental that a linear rescaling of only the packing fraction φ and the reduced temperature θ at constant δ has been so successfully applied for the A₃ and A₄ singularities. The success of this rescaling approach relies on the fact that the shape of φ^c(θ) near the A₃ or A₄ singularities and the shape of the peaks of λ(θ) remain very similar under different closure levels N for a given δ. We expect this similarity to continue as N increases until N → ∞. It is this inherent robustness that renders a simple rescaling of θ and φ in MCT sufficient to predict the A₃ and A₄ singularities within tolerance errors.³⁹ However, regardless of how we rescale φ and θ, it is strictly impossible to merge the full φ^c(θ) and λ(θ) curves without simultaneously rescaling δ. Consequently, we argue that it is a priori somewhat imprudent to test MCT by only rescaling two parameters; ideally one should first solve the first two or three levels of the GMCT hierarchy to establish a robust trend of the transitions before applying such a rescaling approach.

3.3 Relaxation dynamics

We now consider the time dependence of the dynamic density correlators. The effect of higher-order GMCT on the structural relaxation dynamics of monodisperse hard spheres has already been studied in detail.^21–23 As already mentioned in Section 2, one of the conclusions of these recent studies is that for all packing fractions and wavenumbers, GMCT MF closures provide an upper bound to the dynamics while the exponential closures give a lower bound. The two types of closures systematically converge to each other when the closure level N increases. More specifically, if we focus on the MF closures, the predicted relaxation becomes faster in a convergent manner. Here we find that this pattern also applies to sticky hard spheres, i.e., for a given state (fixed packing fraction, fixed reduced temperature, and fixed relative width of the potential well), the relaxation is faster when the MF closure level N increases. We can clearly see this in Fig. 10, where the normalized intermediate scattering functions f₁(k, t) = F₁(k, t)/S(k) at nine different liquid states under different MF closure levels are presented. For example, it can be seen that the N = 3 results (blue lines in Fig. 10) always decay faster than the corresponding N = 2 (black) curves, and more generally the differences between the curves for closure levels N + 1 and N become smaller as N increases. These findings are thus consistent with a uniform convergent trend of the GMCT hierarchy.

Fig. 10 Normalized two-point density correlation functions f₁(k, t) for liquid states at δ = 0.03 and θ = 0.1875 near the A₃ singularities under different GMCT closure levels. The crosses in the inset indicate the considered liquid-state points. (a) Results for packing fractions φ = 0.5, 0.53, 0.5357; these points lie near the A₃ singularity (φ = 0.545, θ = 0.182) predicted from GMCT under closure level N = 2. The solid, dashed, and dash-dotted lines correspond to different packing fractions; different colors correspond to different MF closure levels. The green lines are fits to the expression f₁(k, t) = f^c₁(k, t) − C₁(k) ln(t), where C₁(k) is a wavenumber-dependent constant. (b) Same as in (a), but for packing fractions φ = 0.518, 0.548, 0.554; these points lie near the A₃ singularity ((φ = 0.56, θ = 0.167) predicted from GMCT under closure level N = 3. (c) Same as in (a), but for packing fractions φ = 0.531, 0.561, 0.5666. It takes more than one month to obtain the data for the red dash-dotted line (N = 4); we expect this curve to ultimately decay to 0.

The more important properties of the dynamics are the scaling laws in the β- and α-relaxation regimes, since these are widely used for analyzing experimental and simulation data. Within MCT, the scaling laws of eqn (17)–(19) and (22) are known to be applicable for liquid states that lie near the A₂ singularities but away from the A₃ and A₄ singularities. Note that in this case the exponent parameter λ must be smaller than 1, or more strictly λ < 0.9, to ensure that the liquid state is sufficiently far away from the A₃ and A₄ points. We find that the above scaling laws are also preserved within GMCT under MF closures near the corresponding A₂ singularities, i.e., near both the GMCT-predicted liquid-to-repulsive-glass and liquid-to-attractive-glass transitions. This is natural and reasonable given the similar mathematical form of the MCT and GMCT equations. However, when the liquid states are close to the A₃ and A₄ singularities, the above scaling laws for A₂ break down and new behaviors of the relaxation dynamics will appear. Hence we will focus on the dynamics near these higher-order singularities in the following.

We first study the dynamics near the A₃ singularities for different GMCT closure levels. The A₃ points at δ = 0.03 are (θ ≈ 0.182, φ ≈ 0.545) and (θ ≈ 0.167, φ ≈ 0.56) for N = 2 and N = 3, respectively. Even though we could not accurately determine the location of the A₃ point for N = 4 due to limited computing power, we are certain this point must exist and should lie at a φ slightly above 0.57 and a θ slightly above 0.156 for δ = 0.03, since there is an N = 4 glass–glass transition found at φ = 0.57 and θ = 0.156 for δ = 0.03 (see also the non-ergodicity parameters in Fig. 8). For our A₃ analysis we fix θ = 0.1875 and select three liquid states with different packing fractions for each closure level N; this value of θ is close to the θ of the A₃ singularity for all considered N. The results for the predicted normalized dynamical two-point density correlation functions f₁(k, t) are shown in Fig. 10. As we get closer to the A₃ singularities, we find that the relaxation of f₁(k, t) becomes more stretched compared to the power laws described in eqn (18) and (19), and instead a logarithmic relaxation regime appears. This logarithmic decay of f₁(k, t) for liquid states near A₃ singularities, described viaeqn (23), is applicable for all closure levels N (see the green dashed fitted curves in Fig. 10). The logarithmic law for f₁(k, t) near A₃ singularities within GMCT qualitatively agrees with the one in MCT, and hence our GMCT results also agree with the body of simulation and experimental data that have already been qualitatively compared to MCT.^36,39,49 More generally, we find that higher-order GMCT inherits all the main qualitative features of the MCT dynamics near A₃ singularities. In particular, upon approaching the A₃ point from the liquid side (black, blue, and red dashed lines in Fig. 10(a), (b), and (c), respectively), we can also see that the α-relaxation process does not start with von Schweidler law of eqn (19), and as a consequence the superposition principle of eqn (22) is disobeyed.⁵⁸ The only main difference between MCT and GMCT is that the corresponding packing fractions are improved, since the A₃ point becomes more accurate upon increasing the closure level.

Besides the logarithmic decay for f₁(k, t), we find that this logarithmic decay is also applicable for 2n-point density correlation functions with n > 1. For example, in Fig. 11 we show f₂(k, k, t) near the A₃ singularity under closure level N = 3 for two different wavenumbers. The fit with f₂(k, k, t) = f^c₂(k, k) − C₂(k, k) ln(t), where C₂(k, k) is a wavenumber-dependent constant, is applicable for more than three decades in time. To the best of our knowledge, no empirical results for f₂(k, k, t) or similar higher-order dynamical density correlation functions have been reported in this glassy regime. We argue that the logarithmic decay for higher-order correlators can be another strict test for both GMCT and the A₃ singularities predicted from it for sticky hard spheres.

Fig. 11 Normalized four-point density correlation functions f₂(k, k, t) for liquid states at δ = 0.03, θ = 0.1875 and φ = 0.554 (corresponding to the right cross in the inset of Fig. 10(b)) near the A₃ singularities predicted from GMCT with closure level N = 3. The solid and dotted curves correspond to wavenumber k = 4.2 and k = 7.0, respectively. The green lines are fits to the expression f₂(k, k, t) = f^c₂(k, k) − C₂(k, k) ln(t), where C₂(k, k) is a wavenumber-dependent constant.

Finally let us go to the dynamics of the density correlation functions near the A₄ singularities. The regime of the very stretched and nearly logarithmic decay of f₁(k, t) is found for both closure levels N = 2 and N = 3, as indicated by the solid lines in Fig. 12. The range of this regime can be much wider than the range near A₃ in Fig. 10, although for N = 3 [blue solid line in Fig. 12(b)] the range of the strictly logarithmic decay is not as wide as the one for N = 2 [black solid line in Fig. 12(a)]. We attribute this to the lower numerical accuracy of the A₄ singularity for N = 3 (with a maximum λ = 0.97 at δ = 0.0465 instead of λ = 1). We have also verified that the logarithmic decay is found for f₂(k₁, k₂, t) at several wavenumbers (not shown here). We further test the time-wavenumber factorization property, eqn (17), near the A₄ singularities. If the factorization property holds, the rescaled correlation functions [f_n(k₁,…,k_n,t) − f^c_n(k₁,…,k_n)]/h_n(k₁,…,k_n) should collapse on the common function G(t). In particular, all correlators should cross their plateau value f^c_n(k₁,…,k_n, t) at the same time t₋ with G(t₋) = 0. We show the rescaled correlation function of f₁(k, t) in Fig. 13(a) with the f^c₁(k) and h₁(k) [shown in Fig. 13(b)] from the A₄ singularity point under the corresponding closure level N. We find that for both N = 2 and N = 3, the factorization property is still satisfied (see Fig. 13(a) with the arrows indicating the t₋). However, the validity range is only around two decades in time, which is much narrower than the one near A₂ singularities, for which the factorization property applies over approximately four decades.^22,58

Fig. 12 Normalized two-point density correlation functions f₁(k, t) at different wavenumbers near the A₄ singularities under different GMCT closure levels. (a) f₁(k, t) under closure level N = 2 at φ = 0.524 and θ = 0.234, near the A₄ point (φ = 0.526, θ = 0.234) predicted from MCT. (b) f₁(k, t) under closure level N = 3 at φ = 0.543 and θ = 0.22 near the A₄ point (φ = 0.544, θ = 0.22) predicted from GMCT under closure level N = 3.

Fig. 13 Scaling laws near the A₄ singularity point. (a) Rescaled two-point density correlation functions [f₁(k, t) − f^c₁(k)]/h₁(k) at the same liquid states as in Fig. 12, i.e., near the corresponding A₄ points for GMCT closure levels N = 2 (black lines) and N = 3 (blue lines). The black (blue) arrow indicates the position of G(t₋) = 0 with t₋ = 1.48 × 10³ (t₋ = 1.95 × 10³) for closure level N = 2 (N = 3). (b) The functions f^c₁(k) and h₁(k) obtained from the corresponding A₄ singularities using eqn (8) and (16).

4 Conclusion

In this work, we have presented a detailed numerical analysis of the glassy dynamics of sticky hard spheres using first-principles-based microscopic generalized mode-coupling theory. We have obtained the phase diagrams in the three-dimensional control-parameter space under different mean-field closure levels in GMCT. Upon increasing the closure level N of the GMCT hierarchy, i.e., by expanding the dynamical equations in terms of increasingly higher order density correlations, the phase diagrams become closer to the empirical results. We have clarified that the improvement from GMCT is not merely a simple rescaling of control parameters, as is usually done in standard MCT analyses; rather, higher-order GMCT yields subtle but intrinsically different results compared to MCT, as reflected in the shifted phase diagrams of sticky hard spheres relative to the hard-sphere case, and the exponent parameter λ.

From the phase diagrams and the non-ergodicity parameters, we have demonstrated that the peculiar features of dense sticky hard spheres, i.e., glassy reentrance, a glass–glass transition, and higher-order A₃ and A₄ singularities, are also successfully predicted from GMCT. The positions of these transition points are also more accurate than those obtained from standard MCT. These triumphs indicate that the complex interplay of repulsive and attractive short-ranged interactions can be adequately captured within GMCT. We have also demonstrated that when the attraction is very strong, the difference of the predicted critical packing fractions at any given width of the potential well between GMCT and MCT becomes much larger. This implies that, compared to repulsive interactions, the attractive part of the interaction potential may play a more important role in the higher-order density correlation functions. This non-trivial effect of attractions has also been reflected in the different patterns of the non-ergodicity parameters of the four-point density correlation functions f^c₂(k₁, k₂) for the repulsive-glass and attractive-glass phase. We hypothesize that these higher-order density correlations may also contain new, and non-trivial, information on dynamical heterogeneity,⁴³ but further research is needed to clarify their precise physical interpretation.

Within higher-order GMCT, we have shown that the celebrated scaling laws for A₂ singularities are no longer applicable when the considered state point is in the vicinity of an A₃ or A₄ singularity. In that case, we find a logarithmic decay in the dynamics that is consistent with the MCT form of eqn (23). We emphasize, however, that this logarithmic decay is applicable not only for the conventionally studied intermediate scattering function f₁(k, t), but also for higher-order density correlation functions such as f₂(k, k, t). The latter of course needs to be further verified in simulations or experiments.

From all the results above, we conclude that GMCT can make detailed first-principles-based predictions of the microscopic dynamics of sticky hard spheres well beyond the standard-MCT regime, whilst maintaining many of the qualitative features of MCT. Interestingly, although a significant body of experimental^47–51,56 and simulation^35–44 work has corroborated the qualitative MCT scenario for sticky hard spheres,^35–41 a recent swap Monte Carlo simulation⁵⁷ disputes the existence of the MCT-predicted liquid–glass transition, the glass–glass transition, and all the higher-order A_l singularities, claiming that these MCT-predicted phenomena in fact lie within a unique ergodic region. Our results support this conclusion in the sense that the MCT-predicted glass transitions are shifted toward higher packing fractions within higher-order GMCT, and as N increases they might become close to or even higher than the highest packing fractions accessible in simulation. However, we here argue that the liquid–glass transition, glass–glass transition and higher-order singularities could still exist, possibly at even higher packing fractions than those probed by state-of-the-art simulations.⁵⁷ We also emphasize that our theoretical results are based on a monodisperse sticky hard sphere model, while the swap Monte Carlo simulations of ref. 57 assume a relatively large polydispersity of 23%. It is known in general that polydispersity can significantly affect the glassy dynamics and the corresponding phase diagrams,⁷³ and hence this aspect should be also carefully studied in future work. Finally, the high-order singularities and logarithmic relaxation of sticky hard spheres may also have a deep analogy to the glassy behavior of randomly pinned supercooled liquids and that of the random field Ising model, as discovered within the inhomogeneous MCT.⁷⁴ We believe this analogy would also be preserved in our GMCT framework, although further work on this is needed.

Since the sticky hard sphere model contains both attractive and repulsive interactions, we may expect that GMCT is also a suitable framework for other systems with short-ranged attractive interactions. It remains to be explored, however, whether the current GMCT framework can ultimately provide a universal description of the elusive structure-dynamics link in glass-forming materials. Firstly, there are still some subtle approximations in our current theory, such as the diagonal approximation for the dynamical higher-order density correlation functions in the GMCT memory functions,⁶⁰ and the Gaussian factorization and convolution approximation for all higher-order static density correlation functions.^19,21 Recent studies suggest that higher-order and local structural metrics also encode rich physics and contribute to the dynamics in a non-trivial manner,^6,18,75 but these structural quantities are all simply regarded as being implicitly included in S(k) or neglected in the current GMCT framework. Future work will be done to test how accurate these approximations are and how one might admit more detailed structural information into the theory. Second, until now, GMCT has been tested only for systems with purely repulsive or repulsive with short-ranged attractive isotropic potentials. These systems all belong to the class of fragile glass formers.^76,77 More work needs be done to test GMCT for strong glass formers such as silica¹⁴, or even superstrong systems such as vitrimers.³³

Conflicts of interest

There are no conflicts to declare.

Appendix: breakdown of Gaussian factorization for attractive and repulsive glasses

In Fig. 14 and 15, we plot the off-diagonal non-ergodicity parameters f^c₂(k₁, k₂) and the difference f^c₂(k₁, k₂) − f^c₁(k₁) ×f^c₁(k₂) at the critical liquid–glass transition points for θ = 0.2 and θ = 0.13 at δ = 0.03 under GMCT closure levels N = 2, 3, and 4. The patterns of f^c₂(k₁, k₂) −f^c₁(k₁) × f^c₁(k₂) for the repulsive glass (Fig. 14(b and c)) and the attractive glass (Fig. 15(b, c, e and f)) are very different. Overall these results support the fact that Gaussian factorization can both overestimate and underestimate the four-point non-ergodicity parameters in certain wavenumber regimes, and the deviation pattern depends on the type of glass state.

Fig. 14 (a) The non-ergodicity parameters f^c₂(k₁,k₂) as a function of wavenumbers k₁ and k₂ at the critical packing fraction for closure level N = 2 (i.e., standard MCT) at θ = 0.2. The patterns of f^c₂(k₁, k₂) from N = 3 and N = 4 are similar. (b) The difference of the non-ergodicity parameters f₂^c(k₁,k₂) of closure level N = 3 with respect to the factorization approximation f^c₁(k₁) × f^c₁(k₂) at the same state. (c) Same as in (b) but for N = 4. The packing fractions are the same as in Fig. 5.

Fig. 15 (a) The non-ergodicity parameters f^c₂(k₁, k₂) as a function of wavenumbers k₁ and k₂ at the critical packing fraction for closure level N = 2 (i.e., standard MCT) at θ = 0.13. The patterns of f^c₂(k₁, k₂) from N = 3 and N = 4 are similar. (b) The difference of the non-ergodicity parameters f^c₂(k₁, k₂) of closure level N = 3 with respect to the factorization approximation f^c₁(k₁) × f^c₁(k₂) at the same state. (c) Same as in (b) but for N = 4. Panels (d–f) contain the same data as panels (a–c), respectively, but are enlarged for small wavenumbers. The packing fractions are the same as in Fig. 5.

Acknowledgements

We acknowledge the Netherlands Organisation for Scientific Research (NWO) for financial support through a START-UP grant.

Notes and references

1. P. G. Debenedetti and F. H. Stillinger, Nature, 2001, 410, 259.
2. L. Berthier and G. Biroli, Rev. Mod. Phys., 2011, 83, 587.
3. W. Götze, Complex dynamics of glass-forming liquids: A Mode-Coupling Theory, Oxford University Press, Oxford, 2009.
4. G. Tarjus, S. A. Kivelson, Z. Nussinov and P. Viot, J. Phys.: Condens. Matter, 2005, 17, R1143.
5. T. Kirkpatrick and D. Thirumalai, Rev. Mod. Phys., 2015, 87, 183.
6. H. Liu, Z. Fu, K. Yang, X. Xu and M. Bauchy, J. Non-Cryst. Solids, 2019, 119419.
7. X. Ma, Z. S. Davidson, T. Still, R. J. Ivancic, S. Schoenholz, A. Liu and A. Yodh, Phys. Rev. Lett., 2019, 122, 028001.
8. S. S. Schoenholz, E. D. Cubuk, D. M. Sussman, E. Kaxiras and A. J. Liu, Nat. Phys., 2016, 12, 469–471.
9. V. Bapst, T. Keck, A. Grabska-Barwińska, C. Donner, E. D. Cubuk, S. S. Schoenholz, A. Obika, A. W. Nelson, T. Back and D. Hassabis, et al. , Nat. Phys., 2020, 16, 448–454.
10. U. Bengtzelius, W. Götze and A. Sjolander, J. Phys. C: Solid State Phys., 1984, 17, 5915.
11. E. Leutheusser, Phys. Rev. A: At., Mol., Opt. Phys., 1984, 29, 2765.
12. D. R. Reichman and P. Charbonneau, J. Stat. Mech.: Theor. Exp., 2005, P05013.
13. L. M. C. Janssen, Front. Phys., 2018, 6, 97.
14. F. Sciortino and W. Kob, Phys. Rev. Lett., 2001, 86, 648.
15. W. Kob, C. Donati, S. J. Plimpton, P. H. Poole and S. C. Glotzer, Phys. Rev. Lett., 1997, 79, 2827.
16. S. C. Glotzer, V. N. Novikov and T. B. Schrøder, J. Chem. Phys., 2000, 112, 509–512.
17. N. Lačević, F. W. Starr, T. Schrøder and S. Glotzer, J. Chem. Phys., 2003, 119, 7372–7387.
18. Z. Zhang and W. Kob, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 14032–14037.
19. G. Szamel, Phys. Rev. Lett., 2003, 90, 228301.
20. J. Wu and J. Cao, Phys. Rev. Lett., 2005, 95, 078301.
21. L. M. C. Janssen and D. R. Reichman, Phys. Rev. Lett., 2015, 115, 205701.
22. C. Luo and L. M. C. Janssen, J. Chem. Phys., 2020, 153, 214507.
23. C. Luo and L. M. C. Janssen, J. Chem. Phys., 2020, 153, 214506.
24. W. Kob and H. C. Andersen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 52, 4134.
25. L. Berthier and G. Tarjus, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 031502.
26. V. Zöllmer, K. Rätzke, F. Faupel and A. Meyer, Phys. Rev. Lett., 2003, 90, 195502.
27. B. Nowak, D. Holland-Moritz, F. Yang, T. Voigtmann, T. Kordel, T. Hansen and A. Meyer, Phys. Rev. Mater., 2017, 1, 025603.
28. G. Szamel and H. Löwen, Phys. Rev. A: At., Mol., Opt. Phys., 1991, 44, 8215.
29. M. Siebenbürger, M. Fuchs, H. Winter and M. Ballauff, J. Rheol., 2009, 53, 707–726.
30. L. Schrack and T. Franosch, Philos. Mag., 2020, 100, 1032–1057.
31. J. Farago, A. Semenov, H. Meyer, J. Wittmer, A. Johner and J. Baschnagel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 051806.
32. J. Farago, H. Meyer, J. Baschnagel and A. Semenov, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 051807.
33. S. Ciarella, R. A. Biezemans and L. M. C. Janssen, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 25013–25022.
34. C. Ruscher, S. Ciarella, C. Luo, L. M. C. Janssen, J. Farago and J. Baschnagel, J. Phys.: Condens. Matter, 2020, 33, 064001.
35. A. M. Puertas, M. Fuchs and M. E. Cates, Phys. Rev. Lett., 2002, 88, 098301.
36. E. Zaccarelli, G. Foffi, K. A. Dawson, S. Buldyrev, F. Sciortino and P. Tartaglia, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 66, 041402.
37. G. Foffi, K. A. Dawson, S. V. Buldyrev, F. Sciortino, E. Zaccarelli and P. Tartaglia, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 65, 050802.
38. E. Zaccarelli, G. Foffi, F. Sciortino and P. Tartaglia, Phys. Rev. Lett., 2003, 91, 108301.
39. F. Sciortino, P. Tartaglia and E. Zaccarelli, Phys. Rev. Lett., 2003, 91, 268301.
40. E. Zaccarelli, F. Sciortino and P. Tartaglia, J. Phys.: Condens. Matter, 2004, 16, S4849.
41. G. Foffi, F. Sciortino, E. Zaccarelli and P. Tartaglia, J. Phys.: Condens. Matter, 2004, 16, S3791.
42. I. Saika-Voivod, E. Zaccarelli, F. Sciortino, S. V. Buldyrev and P. Tartaglia, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 70, 041401.
43. D. R. Reichman, E. Rabani and P. L. Geissler, J. Phys. Chem. B, 2005, 109, 14654–14658.
44. A. J. Moreno and J. Colmenero, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 021409.
45. K. González-López, M. Shivam, Y. Zheng, M. P. Ciamarra and E. Lerner, Phys. Rev. E, 2021, 103, 022605.
46. K. González-López, M. Shivam, Y. Zheng, M. P. Ciamarra and E. Lerner, Phys. Rev. E, 2021, 103, 022606.
47. F. Mallamace, P. Gambadauro, N. Micali, P. Tartaglia, C. Liao and S.-H. Chen, Phys. Rev. Lett., 2000, 84, 5431.
48. T. Eckert and E. Bartsch, Phys. Rev. Lett., 2002, 89, 125701.
49. K. N. Pham, A. M. Puertas, J. Bergenholtz, S. U. Egelhaaf, A. Moussad, P. N. Pusey, A. B. Schofield, M. E. Cates, M. Fuchs and W. C. Poon, Science, 2002, 296, 104–106.
50. K. Pham, S. Egelhaaf, P. Pusey and W. C. Poon, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 69, 011503.
51. L. J. Kaufman and D. A. Weitz, J. Chem. Phys., 2006, 125, 074716.
52. S. Buzzaccaro, R. Rusconi and R. Piazza, Phys. Rev. Lett., 2007, 99, 098301.
53. P. J. Lu, E. Zaccarelli, F. Ciulla, A. B. Schofield, F. Sciortino and D. A. Weitz, Nature, 2008, 453, 499–503.
54. Z. Zhang, P. J. Yunker, P. Habdas and A. Yodh, Phys. Rev. Lett., 2011, 107, 208303.
55. Z. Brown, M. J. Iwanicki, M. D. Gratale, X. Ma, A. Yodh and P. Habdas, EPL, 2016, 115, 68003.
56. N. Willenbacher, J. S. Vesaratchanon, O. Thorwarth and E. Bartsch, Soft Matter, 2011, 7(12), 5777–5788.
57. C. J. Fullerton and L. Berthier, Phys. Rev. Lett., 2020, 125, 258004.
58. K. Dawson, G. Foffi, M. Fuchs, W. Götze, F. Sciortino, M. Sperl, P. Tartaglia, T. Voigtmann and E. Zaccarelli, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2000, 63, 011401.
59. J.-P. Hansen and I. R. McDonald, Theory of simple liquids, Elsevier, Amsterdam, 2013.
60. S. Ciarella, PhD thesis, Eindhoven University of Technology, 2021.
61. R. Zwanzig, Nonequilibrium statistical mechanics, Oxford University Press, 2001.
62. L. M. C. Janssen, P. Mayer and D. R. Reichman, J. Stat. Mech.: Theor. Exp., 2016, 054049.
63. R. A. Biezemans, S. Ciarella, O. Çaylak, B. Baumeier and L. M. C. Janssen, J. Stat. Mech.: Theor. Exp., 2020, 103301.
64. C. D. Meyer, Matrix analysis and applied linear algebra, Siam, 2000, vol. 71.
65. V. I. Arnol'd, Catastrophe theory, Springer Science & Business Media, 2003.
66. W. Gotze and L. Sjogren, J. Phys.: Condens. Matter, 1989, 1, 4203.
67. W. Götze and M. Sperl, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 66, 011405.
68. T. Franosch, M. Fuchs, W. Götze, M. R. Mayr and A. Singh, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 55, 7153.
69. M. Fuchs, W. Gotze, I. Hofacker and A. Latz, J. Phys.: Condens. Matter, 1991, 3, 5047.
70. W. Götze and M. Sperl, J. Phys.: Condens. Matter, 2003, 15, S869.
71. L. Fabbian, W. Götze, F. Sciortino, P. Tartaglia and F. Thiery, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 59, R1347.
72. J. Bergenholtz and M. Fuchs, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 59, 5706.
73. S. Ciarella, C. Luo, V. E. Debets and L. Janssen, 2021, arXiv preprint arXiv:2103.16522.
74. S. K. Nandi, G. Biroli, J.-P. Bouchaud, K. Miyazaki and D. R. Reichman, Phys. Rev. Lett., 2014, 113, 245701.
75. J. F. Robinson, F. Turci, R. Roth and C. P. Royall, Phys. Rev. Lett., 2019, 122, 068004.
76. T. Scopigno, G. Ruocco, F. Sette and G. Monaco, Science, 2003, 302, 849–852.
77. L.-M. Martinez and C. Angell, Nature, 2001, 410, 663–667.

---