1 Introduction

Understanding how supercooled liquids become rigid and turn into amorphous solids is still one of the major challenges in condensed matter physics [1,2,3,4]. This so-called glass transition is not a transition in the thermodynamic sense [5], but it is defined by the dramatic increase in viscosity (or relaxation time) upon only a relatively slight change in thermodynamic control parameters, e.g., temperature or density [6, 7]. This sudden and highly nonlinear dynamic response is accompanied by only subtle changes in the microscopic structure, rendering it difficult to identify the main physical mechanisms underlying the glass transition [8,9,10,11,12].

In principle, it is widely accepted that the dynamics of each material is ultimately related to its structure [13], and numerous theories have also aimed to exploit this idea to describe the glass transition [14,15,16,17,18,19,20,21,22]. Among these, mode-coupling theory (MCT) stands out as one of the few theories which is entirely based on first principles [8, 19, 20, 23]. This theory seeks to predict the full microscopic relaxation dynamics of a glass-forming system (as a function of time, temperature, density, and wavenumber k) based solely on knowledge of simple structural material properties, such as the static structure factor S(k). Although the theory is often not fully quantitatively accurate, MCT has enjoyed success in predicting, e.g., multi-step relaxation patterns and universal scaling laws in the dynamics, stretched exponentials, and growing dynamical length scales [24,25,26]. Furthermore, the theory offers a qualitative and physically intuitive account of glass formation in terms of the so-called cage effect [25]. Further extensions of MCT such as the stochastic beta relaxation theory can even explain the crossover from a power law to exponential growth of the relaxation time [21, 22]. The (most severe) limitation of MCT lies, however, in its assumption of Gaussian correlations which causes noticeable discrepancy between the theory and experiments.

So far, promising methodical MCT correction efforts have been put forward for single-component systems—or equivalently systems with a small degree of polydispersity—using higher-order field-theoretic loop expansions [27,28,29,30,31,32]. Results show that such an expansion can indeed be accomplished, producing a novel, hierarchical first-principles theory known as generalized MCT (GMCT). By systematically developing the hierarchical equations, GMCT has already proven to be capable of predicting the microscopic dynamics of glassy materials with near-quantitative accuracy in the low to moderately supercooled regime [29]. Similar to MCT, this generalized framework also requires only static structure as input and has no free parameters. Furthermore, GMCT also provides predictive insights into regimes of previously inaccessible dynamics for single-component glassy systems [29, 33, 34], and notably preserves the celebrated scaling laws of standard MCT [33, 34].

Unfortunately, single-component systems are a poor representation of most studied glasses, which are typically polydisperse and as a result exhibit different overall dynamics compared to monodisperse systems [35, 36]. In experiments, polydispersity is for the most part inevitable, while in computer simulations it is often added to hinder and prevent crystal formation [37]. Noticeably, techniques such as Monte Carlo (MC) swaps capitalize onto polydispersity to achieve faster relaxation dynamics and explore the free energy landscape in uncommon ways [38]. Binary polydisperse systems are the simplest generalization of single-component systems in this direction. They add only a supplementary component to the mix and are able to retain most of the advantages of a polydisperse system while adding the least possible complexity.

In this paper, we set out to extend the GMCT framework to systems with an arbitrary number of species, and seek to apply the newly developed framework to describe arguably the most famous and simple examples of binary glassy systems: the Kob–Andersen binary Lennard-Jones (KABLJ) mixture [37] and its purely repulsive version based on the Weeks–Chandler–Andersen (WCA) potential [39]. These systems have been extensively studied in the past and comparisons with MCT have identified the existence of a region where MCT is already in its non-ergodic phase, thus predicting a glass, while simulations at the same temperature and density indicate a supercooled liquid phase [16, 40,41,42,43,44]. In other words, a discrepancy still exists between the simulations and MCT, even when MCT is extended to multi-component systems [45,46,47,48,49]. Here, after demonstrating the ability of GMCT to systematically tackle this discrepancy, we will also address a fundamental question regarding the simplest ingredients required to describe the dynamics of binary supercooled liquids. Due to the fact that standard MCT—which is based solely on S(k)—fails in predicting the precise location of the glass transition, it could be concluded that higher-order correlations are required [40, 41, 44]. However, GMCT can circumvent this failure. It again uses only S(k) as input, but in a more refined set of equations which can translate structural properties into dynamical ones in a more accurate manner. Applying our multi-component GMCT to both mixtures, we will conclude that each level of the GMCT hierarchy provides a significant improvement in the prediction of the glass transition, finally conjecturing that the infinite hierarchy might be able to accurately predict the glassy dynamics from S(k) only.

2 Multi-component GMCT

Multi-component GMCT is derived starting from the Mori–Zwanzig approach [50, 51] to predict the dynamics of density correlations, similarly to standard MCT. However, while standard MCT amounts to a single integro-differential equation closed by a factorization approximation [8, 19, 20, 23, 47], GMCT is a hierarchy of nested integro-differential equations [29]. Each level of this hierarchy represents an MCT-like dynamical equation for a higher-order, multi-point density correlation function, which we recursively solve and use to predict the dynamics of the correlations at the lower levels. Since a solution of this hierarchy is well defined for any self-consistent closure or truncation of the hierarchy [52], we can formally continue the GMCT scheme up to arbitrary order to include as many higher-order correlations as desired.

In the Supplementary Information, we report the full derivation of multi-component GMCT for an M-component mixture. To summarize it here, we introduce the main objects of the theory, i.e., the species-dependent density modes:

$$\begin{aligned} \rho _{\mathbf {q}}^{\alpha }(t)=\frac{1}{\sqrt{N}}\sum _{i=1}^{N_{\alpha }}e^{-i\mathbf {q}\cdot \mathbf {r}_i^{\alpha }(t)}~. \end{aligned}$$
(1)

Here, \(\mathbf {q}\) is a wavevector of length \(q=|\mathbf {q}|\), t is the time, the index \(\alpha \) represents one of the M species, \(N_{\alpha }\) is the number of particles that belong to the species \(\alpha \), and \(N=\sum _{\alpha =1}^{M}N_{\alpha }\) is the total number of particles in the system. To simplify our equations we introduce the notation that \(\{x_i\}\) is a list \(x_1,\ldots ,x_n\) and \(\{x_i\}/x_j\) is the same ordered list \(\{x_i\}\) where the specific element \(x_j\) has been removed. In solving multi-component GMCT we are interested in determining the dynamical equation of the density correlations of order n. These dynamical correlations are defined as

$$\begin{aligned} F^{(n)}_{\{\alpha _i\};\{\beta _i\}} (\{k_i\},t)=\langle \rho ^{\alpha _1}_{- \mathbf{k} _1}\cdots \rho ^{\alpha _n}_{- \mathbf{k} _n}\rho ^{\beta _1}_\mathbf{k _1} (t)\cdots \rho ^{\beta _n}_\mathbf{k _n}(t)\rangle , \nonumber \\ \end{aligned}$$
(2)

where the angular brackets denote an ensemble average. In particular, when the order \(n=1\), the correlation corresponds to the intermediate scattering function. The only required input of the theory (aside from temperature and density) is the set of wavevector-dependent static correlations \(S^{(n)}\),

$$\begin{aligned} S^{(n)}_{\{\alpha _i\};\{\beta _i\}}(\{k_i\})= & {} F^{(n)}_{\{\alpha _i\};\{\beta _i\}}(\{k_i\},t=0)\nonumber \\\approx & {} \prod _{i=1}^{n}S^{\alpha _i\beta _i}({k_i}), \end{aligned}$$
(3)

which defines the full microstructure of the system at any given temperature and density. Here, we factorize these static correlations as products of the two-point correlation \(S^{\alpha \beta }({k})\equiv S^{(1)}_{\alpha ;\beta }(k)\), which is also known as the static structure factor. This means that all the predictions of the theory are based on the structure factor only. As such, we will conclude later that all the predicted differences in dynamics among the two supercooled liquid systems are already encoded in their static structure factors.

Within our multi-component GMCT hierarchy, each dynamical correlation function \(F^{(n)}_{\{\alpha _i\};\{\beta _i\}}(\{k_i\},t)\) obeys the following equation of motion:

$$\begin{aligned}&{\ddot{F}}^{(n)}_{\{\alpha _i\};\{\beta _i\}}(\{k_i\},t)+\sum _{\{\gamma _i\}} \mu _{\{\alpha _i\};\{\gamma _i\}}{\dot{F}}^{(n)}_{\{\gamma _i\};\{\beta _i\}}(\{k_i\},t)\nonumber \\&\quad +\sum _{\{\gamma _i\}\{\theta _i\}}F^{(n)}_{\{\alpha _i\};\{\gamma _i\}}(\{k_i\},t)\left( S^{(n)}\right) ^{-1}_{\{\gamma _i\};\{\theta _i\}}\nonumber \\&\quad (\{k_i\})J^{(n)}_{\{\theta _i\};\{\beta _i\}}(\{k_i\}) \nonumber \\&\quad +\int _0^t d\tau \!\sum _{\{\gamma _i\}\{\theta _i\}}\!{\dot{F}}^{(n)}_{\{\alpha _i\};\!\{\gamma _i\}}(\{k_i\},t-\tau ) \left( J^{(n)}\right) ^{-1}_{\{\gamma _i\};\{\theta _i\}} (\{k_i\})\nonumber \\&\quad K^{(n)}_{\{\theta _i\};\{\beta _i\}}(\{k_i\},\tau )=0. \end{aligned}$$
(4)

where the effective friction coefficient matrix \(\mu _{\{\alpha _i\};\{\beta _i\}}\) is assumed to be diagonal and species independent, i.e., \(\mu _{\{\alpha _i\};\{\beta _i\}}=\mu \prod _i \delta _{\alpha _i\beta _i}\). The matrices J are static elements defined as

$$\begin{aligned}&J^{(n)}_{\{\alpha _i\};\{\beta _i\}}(\{k_i\})= \left\langle \frac{d}{dt}\left[ \rho ^{\alpha _1}_{k_1}\cdots \rho ^{\alpha _n}_{k_n}\right] |\frac{d}{dt}\left[ \rho ^{\beta _1}_{k_1}\cdots \rho ^{\beta _n}_{k_n}\right] \right\rangle \nonumber \\&\quad \approx \sum _{i=1}^{n}\delta _{\alpha _i,\beta _i}\frac{k_BTx_{\alpha _i}k_i^2}{m_{\alpha _i}S^{\alpha _i\beta _i} ({k_i}) } \prod _{j=1}^{n}S^{\alpha _j\beta _j}(k_j)~. \end{aligned}$$
(5)

with T the temperature, \(m_\alpha \) the particle mass of type \(\alpha \), and \(x_\alpha =N_\alpha /N\) their number ratio. Each level n of the hierarchy is connected to the next via the memory term

$$\begin{aligned}&K^{(n)}_{\{\alpha _i\};\{\beta _i\}}(\{k_i\},\tau )= \frac{\rho }{2}\sum _{\mu '\nu '}\sum _{\mu \nu } \nonumber \\&\quad \int \frac{d\mathbf {q}}{(2\pi )^3}\cdot \sum _{j=1}^{n}\frac{k_BT}{m_{\alpha _j}}\mathcal {V}_{\mu '\nu '\alpha _j}(\mathbf {q,k_j-q,k_j})\nonumber \\&\quad F^{(n+1)}_{\mu ',\nu ',\{\alpha _i\}/\alpha _j;\mu ,\nu ,\{\beta _i\}/\beta _j}(\mathbf {q,k_j-q},\{k_i\}/k_j,\tau ) \nonumber \\&\quad \cdot \mathcal {V}_{\mu \nu \beta _j}(\mathbf {q,k_j-q,k_j})\frac{k_BT}{m_{\beta _j}}~. \end{aligned}$$
(6)

Here, \(\mathcal {V}_{\alpha \beta \gamma }(\mathbf {q,k-q,k})\) is the static vertex function, which remains equal to the standard one of MCT [47] and physically represents the coupling strength among different wavevectors. Explicitly, the vertex function reads

$$\begin{aligned}&\mathcal {V}_{\alpha \beta \gamma }(\mathbf {q,k-q,k})= \delta _{\beta \gamma }\mathbf {q}\cdot \mathbf {k}c_{\alpha \gamma }(q)\nonumber \\&\quad +\delta _{\alpha \gamma }(\mathbf {k}-\mathbf {q})\cdot \mathbf {k}c_{\beta \gamma }(|\mathbf {k}-\mathbf {q}|), \end{aligned}$$
(7)

with the direct correlation function \(c_{\alpha \beta }(q)\). It relates to the static structure factors via the Ornstein–Zernike equation \(c_{\alpha \beta }(q)=\rho ^{-1}(\delta _{\alpha \beta }/x_{\alpha }-(\mathbf {S}^{-1}(q))_{\alpha \beta }\) where \(\rho =N/V\) is the number density of the system [53]. Finally, we note that in this derivation we have assumed the convolution approximation for static three-point correlators [54, 55] and neglected the so-called off-diagonal correlations, similar to what is done in conventional MCT and single-component GMCT [8, 19, 20, 56]. Equation 4 is subject to the initial boundary conditions \({\dot{F}}^{(n)}_{\{\alpha _i\};\{\beta _i\}}(\{k_i\},t=0)=0\) and \(F^{(n)}_{\{\alpha _i\};\{\beta _i\}}(\{k_i\},t=0)=S^{(n)}_{\{\alpha _i\};\{\beta _i\}}(\{k_i\})\) (Eq. 3) for all \(\{\alpha _i\},\{\beta _i\}\), and \(\{k_i\}\).

In principle, the above hierarchical equations can be solved up to arbitrary order n, but in practice we must apply a suitable closure at finite order \(n_{\mathrm{max}}\) to obtain numerically tractable results. We use the following mean-field closure at level \(n_{\mathrm{max}}>2\) [29]

$$\begin{aligned}&K^{(n_{\mathrm{max}}-1)}_{\{\alpha _i\};\{\beta _i\}}(\{k_i\},t)\approx \frac{1}{n_{\mathrm{max}}-2}\nonumber \\&\quad \sum _j K^{(n_{\mathrm{max}}-2)}_{\{\alpha _i\}/\alpha _j;\{\beta _i\}/\beta _j}(\{k_i\}/k_j,t) F^{(1)}_{\alpha _j;\beta _j}(k_j,t). \qquad \end{aligned}$$
(8)

Note that this makes use of the permutation invariance of all wavenumber arguments \(\{k_1,\ldots ,k_n\}\). In terms of the intermediate scattering function the closure is equivalent to the factorization approximation \(F^{(n_\mathrm{max})}(t)\sim F^{(n_\mathrm{max}-1)}(t)\times F^{(1)}(t)\). Hence, at \(n_\mathrm{max}=2\) we obtain the same closure as in standard multi-component MCT [47]. As a reminder, most of the problems with standard MCT come from the fact that such a factorization closure is too strong and unjustified [19, 20, 23]. Multi-component GMCT allows us instead to shift the closure to a larger \(n_\mathrm{max}>2\), meaning that the correlations \(F^{(n')}(t)\) of order \(n'<n_\mathrm{max}\) are not factorized and are more correctly described. This approach has already been shown to be beneficial in single-component glassy systems [27,28,29, 33, 34], and in this paper, we demonstrate that a similar improvement can be gained for binary systems.

3 Numerical details

3.1 Numerical solution of GMCT

Since the system is isotropic and invariant under rotations, we use spherical coordinates to transform the three-dimensional integrals over \(\mathbf {q}\) that appear in any memory function of the hierarchy (Eq. 6) as a double integral over \(q = |\mathbf {q}|\) and \(p = |\mathbf {k}- \mathbf {q}|\). Then, q is discretized over a uniformly spaced grid of \(N_k=100\) points \(q = q_0 + {\hat{q}} \Delta q\) with \({\hat{q}} = 0, 1, \ldots , N_k-1\) and \(\Delta q=40/N_k\). In the Supplementary Information, we show that a grid of \(N_k>70\) wavenumbers is sufficiently converged to predict the MCT critical point. This choice of parameters allows us to replace the double integral by Riemann sums

$$\begin{aligned} \int _0^\infty dk \int _{|q-k|}^{q+k} dp \; \rightarrow \; (\Delta q)^2 \sum _{{\hat{k}}=0}^{N_k-1} \sum _{{\hat{p}} = |{\hat{q}}-{\hat{k}}|}^{\min [N_k-1,{\hat{q}}+{\hat{k}}]}. \end{aligned}$$
(9)

Following Ref. [57], we set \(q_0=\Delta q/2\) in order to prevent any possible divergence for \(q\xrightarrow []{}0\).

To obtain time-dependent solutions for our GMCT equations, we set the effective friction coefficient \(\mu =1\), assume overdamped conditions by dropping the second-order derivative \({\ddot{F}}^{(n)}\) in Eq. 4 (since it only affects the short-time dynamics), and start with a Taylor expansion around \(t=0\) for all dynamical correlation functions \(F^{(n)}\) up to the level \(n_\mathrm{max}-1\); the correlator at the highest level, \(F^{(n_\mathrm{max})}\), follows from our closure relation, Eq. 8. We then integrate Eq. 4 in time using Fuchs’ algorithm [58], where the first \(N_t=64\) time points are calculated with a step size of \(\Delta t = 10^{-6}\), and \(\Delta t\) is subsequently doubled every \(N_t/2\) points. At each point in time, we iteratively update the wavevector-dependent memory kernels (Eq. 6) for all \(n\le n_\mathrm{max}\) until convergence. Note that in these GMCT equations, the partial static structure factors \(S^{\alpha \beta }(k)\) enter both in the initial boundary conditions for \(F^{(n)}\), as well as in the static vertices and the matrices \(J^{(n)}\). In summary, at any given \((T,\rho )\) we only require \(S^{\alpha \beta }(k)\) as input to predict the microscopic relaxation dynamics of the system. While our GMCT framework gives access to all multi-point dynamical density correlations up to order \(n_\mathrm{max}\), in the following we shall restrict the discussion to the intermediate scattering function \(F^{\alpha \beta }(k,t)=F^{(1)}_{\alpha ;\beta }(k,t)\).

3.2 Numerical simulations

We use multi-component GMCT to predict the glassy behavior of two binary mixtures: the Kob–Andersen binary Lennard-Jones (LJ) mixture [37] and its Weeks–Chandler–Anderson truncation (WCA) [39]. Both are three-dimensional 80:20 mixtures of particles A : B which interact with each other via the following potential

$$\begin{aligned}&V_{\alpha \beta }(r)= {\left\{ \begin{array}{ll} 4\epsilon _{\alpha \beta }\left[ \left( \frac{\sigma _{\alpha \beta }}{r}\right) ^{12}-\left( \frac{\sigma _{\alpha \beta }}{r}\right) ^{6} +C_{\alpha \beta }\right] ,&{} r\le r_{\alpha \beta }^c,\\ 0,&{}r>r_{\alpha \beta }^c. \end{array}\right. }\nonumber \\ \end{aligned}$$
(10)

Here, the cutoff radius \(r^c_{\alpha \beta }\) is \(2.5\sigma _{\alpha \beta }\) for LJ, while it corresponds to the potential minimum for WCA [39]. The constant \(C_{\alpha \beta }\) ensures that \(V_{\alpha \beta }(r^c_{\alpha \beta })=0\). We use \(\epsilon _{AA}=1, \epsilon _{AB}=1.5, \epsilon _{BB}=0.5, \sigma _{AA}=1, \sigma _{AB}=0.8, \sigma _{BB}=0.88\) to obtain good glass-forming mixtures [37].

In order to calculate the relevant quantities we need for a comparison with multi-component GMCT, we perform molecular dynamics simulations in the NVE ensemble using HOOMD-blue [59]. We properly equilibrate both systems at different densities \(\rho \) and temperatures T. Periodic boundary conditions are imposed and the box size is set at a length \(L=10\) so that the density is tuned via the number of particles \(N=1200\), 1400 and 1600. All parameters and results are reported in terms of reduced WCA units [39]. From the simulation trajectories, we calculate the partial static structure factors \(S^{\alpha \beta }(k)\) and the collective intermediate scattering functions \(F^{\alpha \beta }(k,t)\). For the multi-component GMCT calculations, we use the simulated \(S^{\alpha \beta }(k)\) as the input of Eq. 4 to predict the theoretical \(F^{\alpha \beta }(k,t)\). In the next section, we compare the output of multi-component GMCT with the \(F^{\alpha \beta }(k,t)\) obtained from simulation, and show that multi-component GMCT becomes progressively closer to the simulated glass transition temperature as we increase the level of the GMCT hierarchy.

4 Results and discussion

4.1 From structure to dynamics

The strength of GMCT is its capability of predicting dynamics from statics. The first result that we show underlines the sensitivity of GMCT to small variations in the static structure. In Fig. 1a, we compare the partial structure factors (normalized for visual purposes) of the binary LJ (yellow) and binary WCA (gray) systems at density \(\rho =1.2\) and temperature \(T=0.65\), which corresponds to low density in the supercooled regime. Notice that all the components of the structure factor are very similar between the two mixtures, consistent with previous simulations [40,41,42, 44]. However, as shown in Fig. 1c, and also in agreement with earlier studies [40,41,42, 44], the simulated relaxation dynamics of the two systems (blue curves) differ significantly. In particular, the structural relaxation of \(F^{AA}(k_\mathrm{peak},t)\), with \(k=k_\mathrm{peak}\) corresponding to the main peak of \(S^{AA}(k)\), is approximately one order of magnitude slower for the LJ mixture. This disparity in dynamics also becomes more pronounced when decreasing T.

Fig. 1
figure 1

Structure and dynamics of supercooled binary LJ and WCA mixtures at \(T=0.65, \rho =1.2\) and \(T=2.50, \rho =1.6\). In panels a and b, we compare the partial static structure factors (where the AA and BB components are normalized for visual purposes) of binary LJ (yellow) and binary WCA (gray). In c and d, we show the dynamics in the supercooled regime of the component \(F^{AA}(k,t)\), which is the term that dominates the dynamics since 80% of the system is type A. The wavenumber \(k=k_\mathrm{peak}\) corresponds to the maximum of \(S^{AA}(k)\). The different curves in panels (c) and (d) show \(F^{AA}(k,t)\) measured from MD simulations (blue), binary MCT (red), binary GMCT with \(n_\mathrm{max}=3\) (orange), and binary GMCT with \(n_\mathrm{max}=4\) (green). When increasing the level of the GMCT hierarchy \(n_\mathrm{max}\), the \(F^{AA}(k,t)\) predicted by multi-component GMCT tends to converge to the simulation results

Full size image

Importantly, binary MCT can only partly account for these dynamical differences based on the input static structure factors, and furthermore, the standard theory cannot reach quantitative accuracy for either system at any given temperature3 [45,46,47,48]. Indeed, it is also demonstrated in Fig. 1c that binary MCT (red curves) fails to predict the correct dynamics at this temperature and density, erroneously predicting a non-ergodic glass phase for both systems. Our multi-component GMCT framework, on the other hand, better approaches the simulated long-time dynamics from the same \(S^{\alpha \beta }(k)\) as input as we increase the level of the hierarchy \(n_\mathrm{max}\). In particular, note that the highest considered GMCT closure level, \(n_\mathrm{max}=4\) (green curves), correctly yields an ergodic phase for both systems, with the LJ mixture having one to two orders of magnitude slower relaxation dynamics than the WCA mixture. This prediction is in good qualitative agreement with simulation.

When the density is high (\(\rho =1.6\)) the attraction that distinguishes LJ from WCA is less significant, since all particles predominantly probe only the short-range repulsive regime. In fact, we see in Fig. 1b that all the components of \(S^{\alpha \beta }(k)\) are virtually identical among the two mixtures. The dynamics in this supercooled regime, reported in Fig. 1d, is also almost the same for the simulated mixtures. Notice that the value of \(T=2.5\) corresponds to approximately \(1.5T_g\), and it is comparable to the value of \(T=0.65\) at \(\rho =1.2\) of Fig. 1c. Similarly, every level of the binary GMCT hierarchy also predicts almost indistinguishable dynamics at this density. However, it is once again noticeable that a higher \(n_\mathrm{max}\) makes multi-component GMCT converge toward the simulations.

Overall, Fig. 1 clearly shows that small differences in the structure are captured by multi-component GMCT and amplified to predict the dynamics in the glassy regime. While on the one hand this sensitivity of the theory means that a high precision is required when measuring the input-S(k), on the other hand this supports the idea that important information about the dynamics is already enclosed in static 2-point density correlations [42, 44].

4.2 The role of polydispersity

By extending the framework of GMCT to include multiple components, we can directly account for polydispersity, i.e., the heterogeneity of sizes of molecules or particles in a mixture. It is ubiquitous in experiments at the colloidal scale because two particles are hardly equal and in the context of glasses it is also useful to avoid crystallization [9]. Furthermore, it has been shown that even in simulations where monodispersity is possible, it can be beneficial to use polydispersity in order to employ algorithms such as Monte Carlo swaps which can significantly improve the performance of computations [38].

If the degree of polydispersity is small it has been shown that single-component GMCT is capable of very accurate predictions [29]. However, for highly polydisperse systems or complex architectures [60,61,62,63] single-component theories require a pre-averaging of the structure. This can severely influence their predictions. In particular, since (G)MCT is very sensitive to the value of the main peak of the static structure factor [60], averaging the AA correlation with the AB and BB components inevitably leads to a decrease in such peak which, in turn, alters the results of (G)MCT.

In Fig. 2, we examine the consequences of approximating a highly polydisperse system, i.e., our binary LJ and WCA mixtures, as being effectively monodisperse. We report the relaxation time \(\tau \) as a function of the inverse temperature at \(\rho =1.2\), comparing simulations to single-component and multi-component GMCT; for single-component GMCT we use the average static structure factors S(k) as input, whereas for multi-component GMCT we explicitly distinguish between all the partial components \(S^{\alpha \beta }(k)\) (see Fig. S2 for a comparison between the pre-averaged monodisperse structure factor and the binary ones). The relaxation time \(\tau \) is defined as

$$\begin{aligned} F^{AA}(k_\mathrm{peak},t=\tau )=\frac{F^{AA}(k_\mathrm{peak},0)}{e}~, \end{aligned}$$
(11)

which grows rapidly during supercooling toward the glass transition temperature \(T_g\) [9]. We find that single-component (G)MCT significantly underestimates the critical glass transition temperature for both systems. In particular, at the supercooled temperatures where the simulated \(\tau \) reaches a value of \(\tau \sim 10^4-10^5\) (i.e., near the simulated \(T_g\)), we find that our single-component (G)MCT approximation yields a relaxation time that is almost three orders of magnitude too low. This underestimation of the glassy dynamics is also consistent with MCT studies of polymeric systems that use pre-averaged static structure factors [60,61,62,63]. Moreover, note that the qualitative shape of the \(\tau (T)\) curves predicted by single-component theory also deviates markedly from the simulation results, and that little improvement is gained by increasing \(n_\mathrm{max}\).

Fig. 2
figure 2

Relaxation time as a function of inverse temperature from simulations (blue), single component (dashed lines), and multi-component (continuous line) (G)MCT, for binary LJ and WCA at \(\rho =1.2\). The data show that single-component GMCT does capture only a very weak slowdown upon supercooling and does not show any trace of a critical point, nor any significant improvement for larger \(n_\mathrm{max}\). The results of multi-component GMCT are discussed more in detail in the next figure; here, we only show that multi-component GMCT reproduces the binary simulations more realistically than single-component GMCT

Full size image

By contrast, when properly taking into account the binary nature of both systems, multi-component (G)MCT yields predictions that more closely resemble the \(\tau (T)\) simulation curves, at least on a qualitative level. We also see that the multi-component theory in fact overestimates the critical temperature, with the highest overestimation found for the lowest \(n_\mathrm{max}\). This general tendency to overestimate the glassiness is also consistent with other multi-component [47, 48] and standard MCT [20] calculations. Overall, these results underline the fact that non-trivial couplings exist in the structure and dynamics of multi-component glassy mixtures, highlighting the need to explicitly account for polydispersity in such systems.

4.3 Relaxation time

We proceed by comparing our numerical simulations with the predictions of multi-component GMCT. The comparison is summarized in Fig. 3 where we report the relaxation time \(\tau \) as a function of the inverse temperature for three different bulk densities. It should be noted that here we solely focus on \(\alpha =\beta =A\) when determining \(\tau \), because particles of type A constitute 80% of the system and therefore dominate the dynamics. Furthermore, as before, we set \(k=k_\mathrm{peak}\), corresponding to the maximum of \(S^{AA}(k)\), thus focusing on the slowest modes in the system [8, 20, 25].

Fig. 3
figure 3

Relaxation time as a function of inverse temperature from simulations (blue) and multi-component (G)MCT, for binary LJ and WCA mixtures at \(\rho =1.2\), 1.4, and 1.6. The relaxation time is evaluated from \(F^{AA}(k_\mathrm{peak},t)\), corresponding to the majority species A and the wavenumber \(k=k_\mathrm{peak}\) where \(S^{AA}(k)\) has its maximum. The data show that by increasing the GMCT closure level, the predictions of multi-component GMCT become increasingly more accurate

Full size image

The results in Fig. 3 show that MCT (red curves, corresponding to GMCT with closure level \(n_\mathrm{max}=2\)) overestimates the value of \(T_g\) obtained from simulations (blue) as expected [40, 41]. However, if we increase the level of the hierarchy to \(n_\mathrm{max}=3\) (orange) and then \(n_\mathrm{max}=4\) (green), the accuracy increases and the critical point of GMCT manifestly converges toward the simulations. This uniform convergence of the multi-component theory with increasing \(n_\mathrm{max}\) is also consistent with earlier findings from single-component GMCT [29, 32,33,34].

From the data in Fig. 3, it is particularly noteworthy that at \(\rho =1.2\), where the difference between the simulated LJ and WCA dynamics is the largest, higher-order multi-component GMCT becomes increasingly better at distinguishing between the two mixtures [41, 64]. Here, it is important to recall that for each temperature and density considered, all our GMCT calculations use the same \(S^{\alpha \beta }(k)\) as input, regardless of the chosen \(n_\mathrm{max}\). The fact that increasing \(n_\mathrm{max}\) leads to better dynamical predictions, and perhaps might even become (near-)exact in the limit of \(n_\mathrm{max}\xrightarrow []{}\infty \) [32], clearly suggests that static 2-point correlations already constitute an important indicator of glassiness—provided that the appropriate dynamical framework is used to translate structure into dynamics. The importance of structural pair correlations has also been verified recently through agnostic machine learning methods [42, 65,66,67], and with the here presented work, we can now place this result on a firmer, first-principles-based theoretical footing.

4.4 Role of attraction

The static structure factor S(k) has also been shown to contain information about higher-order static correlations [44, 68]. However, if we use the relevant \(S^{\alpha \beta }(k)\) as the main input of standard MCT, the theory is not able to efficiently distinguish between LJ and WCA mixtures [40, 41] (also see Fig. 3). This implies that at least on the MCT level, the role of attractive particle interactions in supercooled liquids is not adequately captured.

We show here that higher-order (multi-component) GMCT is more sensitive to small differences in S(k) and thus the theory is able to recursively recognize better the role of attraction. To support this claim, we compare the LJ and WCA dynamics at different temperatures \(T_{\mathrm{LJ}}\) and \(T_{\mathrm{WCA}}\), respectively, where the temperatures are defined such that they yield the same relaxation time \(\tau \). In Fig. 4, we report the measured temperature difference \(T_{\text {WCA}}^{-1}-T_{\text {LJ}}^{-1} \) as a function of the relaxation time \(\tau \). This analysis is based on a power law fitting \(\tau \sim (T-T_0)^{-\gamma } + A_0\), where the parameters \(A_0,T_0\), and \(\gamma \) are fitted to best approximate Fig. 3 for each value of \(\rho \). All the numerical values are reported in the Supplementary Information (Table S I). It can be seen in Fig. 4 that the temperature difference extracted from the simulations (blue) becomes progressively larger as \(\tau \) increases, indicative of the markedly different supercooled LJ and WCA dynamics. In standard MCT, this difference is not properly captured; in fact, binary MCT (red) predicts that the temperature difference is always small and almost constant. However, when we increase the closure level of binary GMCT to \(n_\mathrm{max}=3\) (yellow) and \(n_\mathrm{max}=4\) (green) the difference \(T_{\mathrm{WCA}}-T_{\mathrm{LJ}}\) becomes larger and, similarly to the simulations, it grows approaching the glass transition. We therefore conclude that higher-order GMCT can extract more information from S(k) and hence it is able to better recognize the role of attraction in the emergent supercooled dynamics.

Fig. 4
figure 4

Effect of the attraction in binary LJ mixtures relative to repulsive WCA, measured in terms of the temperature difference \(T_{\text {WCA}}^{-1}-T_{\text {LJ}}^{-1}\) at which both systems yield the same relaxation time \(\tau \). The higher the closure level of binary GMCT, the larger this temperature difference becomes. Hence, we can conclude that higher-order GMCT better recognizes the role of attraction via the corresponding changes in the static structure factor

Full size image

4.5 Conclusions

In this paper, we have derived generalized MCT for multi-component systems, thus extending the earlier version of the theory [27,28,29, 33, 34] to the case of mixtures with an arbitrary number of species. The theory seeks to predict the microscopic relaxation dynamics of glassy mixtures in a fit-parameter-free manner using the static structure factors as its main input. Its hierarchical structure of nested integro-differential equations can be closed and solved self-consistently at any order n. The predictive power of the theory manifestly increases for larger n, providing a promising, and systematically improvable framework to ultimately achieve an accurate description of the elusive structure-dynamics link in glass-forming liquids.

We have used the newly derived multi-component GMCT to describe the glassy dynamics of three-dimensional Kob–Andersen LJ and WCA binary mixtures—systems with almost indistinguishable microstructures but widely different dynamics. We have demonstrated that the theory is able to capture subtle differences in the static structure factors and amplify these to account for the distinct LJ and WCA dynamics. Since the theory only uses \(S^{\alpha \beta }(k)\) as input, all the relevant microstructural information is assumed to be fully encoded in the pair-correlations—a result that is consistent with recent machine learning studies on these systems [42, 67]. Moreover, owing to the improved predictive power of higher-order GMCT compared to standard MCT, we have argued that our theory is also able to better understand the role of attraction in dense supercooled liquids. We have also shown that highly polydisperse systems require a multi-component theory to properly describe the structure-dynamics link in supercooled liquids; this is because the single-component approximation ignores any species-dependent structural correlations in S(k), thus washing away many subtle but important features in the microstructure that subsequently compromise the predictive power of GMCT.

Lastly, we have illustrated that the systematic inclusion of more levels in the multi-component GMCT hierarchy yields quantitatively better predictions for the dynamics, at least based on the first few calculated GMCT levels. This gradual but systematic improvement is also consistent with earlier GMCT studies for single-component systems. In future work, we will aim to push the boundaries of the highest level n we can numerically solve [33, 34], in order to check whether the current GMCT framework might approach the exact scenario in the \(n\xrightarrow []{}\infty \) limit. This effort will also allow us to check if higher-order GMCT can capture the non-trivial glassy phenomenology contained in other successful microscopic theories such as stochastic beta relaxation theory [21, 22]. To conclude, we hope that our multi-component GMCT could be a useful tool to evaluate how static correlations influence the dynamics of supercooled liquids and to make reliable predictions about the dynamics of such liquids from static information only, thereby contributing to the final understanding of the glass transition.