Competing Relaxation Channels in Continuously Polydisperse Fluids: A Mode-Coupling Study

Systems with a high degree of size polydispersity are becoming standard in the computational study of deeply supercooled liquids. In this work we perform a systematic analysis of continuously polydisperse fluids as a function of the degree of polydispersity within the framework of the Mode-Coupling Theory of the glass transition (MCT). Our results show that a high degree of polydispersity tends to stabilize the liquid phase against vitrification, the magnitude of which depends on the shape of the polydispersity distribution. Further, we report on a separation between the localization lengths of the smallest and largest particles. A diameter-resolved analysis of the intermediate scattering functions reveals that this separation significantly stretches the relaxation patterns, which we quantitatively study by an analysis of the dynamical exponents predicted by the theory. Our observations have strong implications for our understanding of the nature of dynamical heterogeneities and localization lengths in continuously polydisperse systems. These results suggest that the dynamics of the smallest particles is of central importance to understand structural relaxation of continuously size polydisperse fluids, already in the mildly supercooled regime where MCT is usually applicable.


I. INTRODUCTION
Size polydispersity, i.e. a heterogeneity in particle sizes, is often inevitable and even necessary to produce stable supercooled liquids. Indeed, without such dispersity, most liquids and colloidal suspensions soon crystallize upon supercooling or compression [1,2]. In fact, size-polydisperse systems form the bulk of the modern literature on computational studies of the glass transition, in particular with the advent of enhanced sampling procedures such as swap Monte Carlo, which require extremely high degrees of size polydispersity for high efficiency [2,3]. While size polydispersity has been shown to stabilize the liquid phase in the supercooled regime, for both simple colloidal systems [2][3][4][5] as well as for polymeric ones [6], the introduction of different particle species also inherently complicates the structural and dynamical aspects of glass-forming systems [7][8][9]. In particular, the behavior of the smallest particles appears to be of crucial importance for the system's global structural relaxation [10,11].
From a theoretical perspective, the Mode-Coupling Theory of the glass transition (MCT) is the only microscopic framework capable of making quantitative predictions on the dynamics of dense, supercooled liquids [12][13][14] from purely structural considerations (i.e. the liquid's pair correlation function). The theory provides a set of equations of motion for several dynamical observables: the coherent and incoherent intermediate scattering functions, as well as the mean-square displacement if appropriate limits are taken [14]. In the case of simple glass formers, MCT predicts a phase transition to an ideal glass state which can be rigorously studied using techniques from the bifurcation of iterated maps [15]. Beyond this, MCT's principal successes lie in the prediction of accurate form factors (e.g. Debye-Waller and Lamb-Mössbauer (LM) factors) as well as the hallmarked two-step, and stretched exponential relaxation of intermediate scattering functions. Furthermore, the theory predicts non-trivial associ- * Contact: l.m.c.janssen@tue.nl ated scaling laws and dynamical exponents (and algebraic relations between them) which many liquids in the supercooled regime have been found to obey both in computer simulations and experimentally [16][17][18][19]. Moreover, the theory can describe more intricate dynamics such as re-entrant behavior in e.g. mildly asymmetric binary mixtures [20], systems with short-range attractions [21,22] or systems in strong geometric confinement [23], as well as more exotic phase transitions, in particular glass-glass transitions [21,24,25] and partially arrested states [25,26]. Recent work has also demonstrated that the theory is capable of making predictions on dynamical heterogeneities, something which was initially deemed impossible due to the inherent absence of many-body correlations in the theory [14]. Indeed, a careful treatment reveals that MCT predicts a diverging (dynamical) correlation length associated with its ideal glass transition [27].
Within the context of MCT, polydispersity-induced effects have been widely studied in binary and ternary mixtures, including for large size ratios (i.e. high polydispersity degrees) [28,29]. Furthermore, Weysser et al. [30] have tested MCT for a system of mildly continuously polydisperse (∼ 5%) quasi-hard spheres where the diameters are drawn from a uniform distribution. By performing an appropriate multicomponent analysis of the mixture, they reported good agreement for the collective and single-particle observables predicted by the theory, except at low wave-numbers for the latter (a typical shortcoming of the theory, especially for singleparticle observables). Their work has clearly highlighted the importance of treating polydispersity accurately by comparing predictions from effective monodisperse MCT and multicomponent MCT for up to five components. Beyond this, little is known on how the degree of size polydispersity quantitatively affects MCT's predictions in continuously polydisperse mixtures.
This work aims to extend these prior studies by providing a systematic and quantitative study of the glassy dynamics predicted by MCT as the degree of polydispersity of continuously polydisperse mixtures is varied. The focus of this work is on MCT predictions solely, without resorting to explicit comparisons with experimental or simulation results.
There are various reasons for this choice; firstly, this removes the need for computationally intensive simulations since the MCT equations of motion can be solved within reasonable computational time, at least up to ten different particle components using the latest open-source code implementation of MCT [31]. Furthermore, there are already many experimental and computational studies that have demonstrated that continuously size polydisperse mixtures display non-trivial dynamical features [10,11,32] which are by-products of the continuous aspect of polydispersity. Secondly, and perhaps more importantly, it is now universally accepted that MCT does not account for the experimentally observed dynamical crossover, which is a widely seen phenomenon that signals a change of dynamics upon entering the more deeply supercooled regime [33,34]. The physical mechanisms that drive the crossover are known as activated events, and are believed to be associated with increasingly cooperative relaxation [35][36][37][38][39]. This in fact forms one of the main criticisms against the theory, and in practice, limits MCT's range of applicability to only the mildly supercooled regime. Here we turn this deficiency into an advantage, as it allows us to probe how the degree of polydispersity affects one well-known, and dominant relaxation mechanism in the mildly supercooled regime in a way that is completely independent from the other mechanisms that are beyond the scope of the theory. This enables us to disentangle potential polydispersity-induced effects which have remained unstudied in recent simulation works on very deeply supercooled liquids [38,39], yet are present already at mildly supercooled states where MCT would be applicable.

II. THEORY
We consider an equilibrium, underdamped mixture of N particles which can be divided into n distinct components, each made up of N σ particles of type σ, where σ = 1, 2, . . . , n. We define the incoherent intermediate scattering functions (ISF) per species at wave-vector k as F where r σ i (t) denotes the position of particle i of species σ at time t. We note that the ISF only depends on the wave-number k = |k| since we consider a translationally and rotationally invariant system. The angular brackets denote an ensemble average. Using standard projection operator methods, the incoherent ISF for species σ is found to satisfy the following integro-differential equation with Ω (s) where m σ is the mass of particle type σ, k B Boltzmann's constant and T denotes the temperature. The MCT memory kernel K (s) σ (k, t) reads where ρ is the bulk density, and c αβ (k) is the partial twobody direct correlation function, which is related to the partial static structure factor S αβ (k) by [40]. Throughout this work, Greek indices refer to particle species labels. The superscript (s) refers to a 'self' (tagged) quantity. Furthermore, Einstein summation convention is used throughout the manuscript. Within MCT, The incoherent ISF depends on the coherent one, which is defined as The coherent ISF satisfies the following equation of motion: in which Ω αβ (k) 2 ≡ k 2 k B T x α /m α · (S −1 (k)) αβ and where the collective memory kernel K αβ (k, t) is approximated by a bi-linear functional of the coherent ISF. This is the celebrated mode-coupling approximation, which reads We denote by x σ ≡ N σ /N the fraction of particle species σ. The coupling constants V αβγ (k, q) can be expressed in terms of the partial two-body direct correlation function [13,14,20]. Overall MCT thus culminates in a closed set of dynamical equations, which can be solved self-consistently once the system-specific static inputs (in the form of the static structure factors, bulk density, temperature, and particle mass) are known. Since our goal is to investigate how the degree of size polydispersity in continuously polydisperse mixtures modifies the theory's principal predictions [14,41], let us briefly recall them. There exists a critical packing fraction φ c beyond which an 'ideal glass' forms (i.e. the point where the density correlation functions no longer decay to zero at any finite time). This marks an ergodicity breaking point. When approached from the liquid side, the transition point is accompanied by a strict divergence of the structural relaxation time τ α ∼ |φ − φ c | −γ with a critical exponent γ. The structural relaxation time τ α is typically defined as the point where the correlation functions have decayed to a small threshold value (e.g. 1/e or 0.1), although more formal definitions also exist [41].
Close to φ c , both the coherent and incoherent ISF develop a distinct two-step relaxation pattern with a long-lived plateau at intermediate times. Around this plateau, the ISFs admit the following asymptotic expansions to leading order: and where f αβ (k) denotes the plateau height of the partial coherent ISF [also known as the partial Debye-Waller (DW) factor], and f (s) σ (k) that of the incoherent contributions to the ISF [also known as the Lamb-Mössbauer (LM) factor] respectively. One should consider the (+) solutions on the approach towards the plateau, and the (−) solutions upon departure from it. The quantity t 0 is some reference timescale which determines the range of validity of the asymptotic expansion, usually defined as the point at which the coherent (incoherent) ISF equals the value of the associated DW (LM) factor at the critical point [41]. The factorization theorem [14] indicates that every term in the asymptotic series separates into a purely wave-number-dependent contribution K (also known as the critical amplitude), and a solely time-dependent function which, to leading order behaves as a power-law with exponents a, b depending on whether one is above or below the critical plateau. The exponent b is commonly known as the von-Schweidler exponent [42]. Furthermore, the three dynamical exponents are related by the following two relations where Γ(x) denotes the gamma function. We note that the exponents are system dependent, but for a given system they are universal for the correlation functions discussed above. Since all species share the same exponents, the species averaged quantity F (s) (k, t) ≡ x σ F (s) σ (k, t) also admits the same asymptotic expansions, albeit with different critical amplitudes.
In addition to accounting for the highly non-trivial scaling laws discussed above, MCT also provides a physically intuitive picture for glassy dynamics in terms of the cage effect. Briefly, this feedback effect stems from the self-consistent, non-linear form of the memory kernel of Eq. (5), which renders MCT incredibly sensitive to subtle changes in the microstructure of the liquid (as quantified by the static structure factor). In particular at the wave-vector corresponding to the first peak of the static structure factor, which typically changes the most upon supercooling, the relaxation dynamics is predicted to slow down significantly, which in turn also drives the slowdown of all surrounding wave-vectors via the self-consistent mode-coupling. Since this process is initiated at the length scale of the first solvation shell, the microscopic origin of the macroscopic slowdown is attributed to caging, whereby particles become trapped in local cages formed by their nearest neighbors. This ultimately culminates in the collective freezing of density fluctuations at the ideal glass transition.

III. NUMERICAL METHODS
In the remainder of this work, we examine in detail how the degree of size polydispersity affects the dynamical behavior of a three-dimensional hard-sphere fluid near the ideal MCT transition. We study results for three different continuous distributions P (D) of particle diameters D: where D max and D min denote the maximal and minimal diameters of the distribution respectively.
• A Gaussian distribution P Gauss (D) defined as with mean D and standard-deviation ∆.
• An inverse cubic distribution P inv. (D) defined as where D max and D min denote the maximal and minimal diameters of the distribution respectively, and the normalization reads A = 1/2 + 1/(σ − 1) where σ ≡ D max /D min is the size ratio of the largest to the smallest particle [3].
The three distributions are shown on the right-most panels of Fig. 1. We focus on different probability distributions as recent work has shown that not only the degree of polydispersity, but also the shape of the distribution matters near dynamical arrest [10]. We note that the inverse cubic distribution has recently become very popular in computational studies of deeply supercooled liquids, owing to its exceptional equilibration efficiency when combined with the enhanced sampling swap Monte Carlo method [3].
To numerically solve Eqs. (1)-(4), we approximate the size distributions by discretized versions. The particles' effective diameters are determined from a quantization of the probability distributions into n = 10 quantiles (i.e. the distribution of particle diameters is sampled uniformly and thus a randomly selected particle has an equal probability of belonging to any particular bin). The effective diameter of particles of each bin is then set equal to the center of mass of the given bin. We show in Appendix A that an effective 10-component description is sufficient for the regimes of polydispersity investigated in this work. We also rescale the quantized diameter list such that it has unit mean. The width of the distributions (and thus the degree of size polydispersity) can be tuned by varying the range of the allowed diameters (D min. , D max. for the Uniform and Inverse Cubic distributions) or the half-width ∆ for the Gaussian distribution. In order to allow for a fair comparison between different distributions, we characterize the degree of polydispersity by a polydispersity index denoted δ, defined as the standard deviation of the (quantized) diameter list. Since the mean is fixed, the value of δ implicitly sets the upper and lower bounds of the particle diameters for a given distribution P (D).
To permit a completely fit-parameter free, first-principles study, we use analytic partial static structure factors for hard spheres obtained from the multicomponent version of the Percus-Yevick (PY) approximation [43]. These serve as our structural input to the MCT equations at a given density. We note that the multicomponent PY approximation has been used in polydisperse MCT studies in the past, leading to convincing predictions compared to simulation results [44]. Furthermore Frenkel et al. [45] have demonstrated that the PY approximation was capable of accurately capturing the structure of (highly) continously polydisperse hard-sphere mixtures of log-normally distributed diameters.
The equations of motion (1)-(4) are solved for the effective 10-component mixture over a logarithmically coarse-grained time grid [28,46], which has been recently been implemented in an open-source solver for integro-differential equations of the mode-coupling type [31]. The wave-vector integrals were performed (in spherical coordinates) over a grid of 100 equidistant points between k min = 0.2 and k max = 39.8. We set the mass of all particles as well as the thermal energy to unity: m σ = 1 ∀ σ and k B T = 1. More details on the numerical routines, in particular with respect to numerical stability, can be found in Appendix B.

IV. STATE DIAGRAM
We first study the location of the predicted glass transition for continuously polydisperse systems as a function of the degree of polydispersity. In order to determine the state diagram, we solve Eq. (4) for the long time limit (t → ∞) by using a Picard iteration procedure and performing a bisection search for the critical point φ c . In Fig. 1 we show the state diagram as determined by MCT for the three different distributions defined above. At low polydispersities, we find re-entrant behavior for all three cases, which has also been observed in the past in mildly size asymmetric binary hard-spheres in the context of MCT [20,[47][48][49]. As the polydispersity index is increased further, we see that the critical points are shifted to higher packing fractions for all distributions considered. This indicates that, according to MCT, a higher degree of size polydispersity stabilizes the liquid phase. This trend has also been observed in recent simulation studies of continuously polydisperse systems [10], and in studies of strongly size-asymmetric binary mixtures, both in simulation [48] and within the context of MCT [20]. Physically, the shift of φ c toward higher values with increasing δ can be attributed to entropically induced effective attractions between particles with large size differences [50], which are known to stabilize the liquid phase [21,51].
For all three distributions studied here, considering a line at fixed packing fraction φ could therefore lead to non-monotonic dependence of the relaxation time. This is particularly true near packing fractions close to that of the single component MCT transition φ = 0.519. For instance, at this point in the low to moderate polydispersity regime (δ < 0.15), one initially expects an increase of relaxation time as the degree of polydispersity increases, followed by a glass region which eventually re-enters the liquid regime at high degrees of polydispersities, for δ ≥ 0.2. In the latter regime we observe that the dynamics become faster as the degree of polydispersity is increased. This last observation is in line with prior studies on polydisperse Lennard-Jones fluids, where increasing the degree of polydispersity also leads to a speed-up of the dynamics [32].
It can be seen that the uniform and Gaussian distribution yield nearly identical critical points, except at very high polydispersities. The stabilization of the liquid phase is largest for the inverse cubic distribution, for which the critical packing fraction can be shifted beyond φ c = 0.540. Note that this particular distribution also permits supreme equilibration efficiency in the supercooled liquid phase with swap Monte Carlo [3]; whether these two observations bear the same physical origin deserves further investigation. We also recall that the nature of the inverse cubic distribution is such that there is an abundance of small particles compared to the large ones. Our findings therefore indicate, on purely theoretical grounds, that not only the relative size of the particles is crucial for the system's stability against vitrification but also the number of particles of each size. This also explains the deviation of the Gaussian distribution from the uniform one for δ > 0.4, as the former has more small particle outliers at fixed δ, which would fluidize the system.
In the limit of very high polydispersity indices δ, we find that for the uniform distribution, the critical point saturates to a constant value. This suggests that, at least at the level of MCT, there is a 'maximal' polydispersity index beyond which the long-time dynamics are no longer affected. This behavior has also been observed in highly-asymmetric binary mixtures in the past [20]. The saturation at large δ is more subtle in the case of the Gaussian and inverse cubic distributions, but we anticipate it should also be reached at sufficiently high polydispersities (δ > 0.5). We note however that at such high degrees of polydispersity, we cannot rule out the existence of exotic partially arrested states, as is predicted by MCT for binary mixtures where the size ratio of the smallest to the largest particles is below D min /D max ≤ 0.35 [25]. In the present study, the size ratio D min /D max ≤ 0.25 and 0.20 for the uniform and Gaussian distributions of diameters respectively, if δ ≥ 0.40. It is therefore not unthinkable that a partially arrested glass exists in this portion of parameter space. The size ratio for the inverse cubic distribution remains smaller for equivalent values of δ, although the same scenario is anticipated for higher degrees of polydispersity. The study of these more exotic glass states in continuously polydisperse mixtures is however beyond the scope of this work and is left for future study.
Lastly, we have verified that the same MCT state diagram is found for a finer wave-vector grid (see Appendix Fig. 12), as well as for other quantizations for the particle diameters (fewer bins), which indicates that all the trends observed in Fig. 1 are genuine and not an artifact of our numerical routines.
In the remainder of this work, we focus on the dynamics of polydisperse systems where the particle diameters follow a uniform or inverse cubic distribution. This choice is motivated by the similarity between the MCT predictions for the Gaussian distribution and the uniform one. Furthermore, all the results that follow have been evaluated at the wavenumber k corresponding to the peak of the species-averaged static structure factor S(k) ≡ αβ S αβ (k). This peak position is weakly polydispersity dependent near the critical point, as shown in Fig. 13 (see Appendix C). MCT state diagram in the polydispersity-packing fraction plane for polydisperse mixtures of hard-sphere liquids whose diameters are sampled from the three distributions described in the main text. The panels on the right-hand side illustrate the distributions of particle sizes P (D).

A. Incoherent Intermediate Scattering Functions and Form Factors
We now move on to study the dynamical features of structural relaxation in polydisperse mixtures. In order to isolate the effects of polydispersity, we compare systems at a fixed relaxation time τ α = 10 10 [a.u.], where τ α is defined as the point where the diameter-averaged incoherent ISF has decayed to some small value ε, i.e. F (s) (k, τ α ) = ε with ε = 10 −5 . We choose τ α = 10 10 [a.u.] to probe the dynamics fairly close to the MCT critical point, and ε = 10 −5 to ensure that all particle species have fully relaxed. One could also perform a study at fixed packing fraction φ c instead of fixed relaxation time τ α , while varying the polydispersity, but this would render the comparisons more difficult as τ α is highly polydispersity dependent.
We show in Fig. 2 the species-averaged (full line) and the partial incoherent ISF (dashed lines), for three different levels of polydispersity δ = 0.1, 0.24, 0.38 for a uniform size distribution. In particular, looking at the partial incoherent ISF, we see that for δ = 0.24 there is a decade of difference in the relaxation times for the smallest (lowest dashed line) and the largest particles (highest dashed line), and nearly four decades for δ = 0.38. The difference in particle-size resolved relaxation times becomes even more extreme if one uses a more conventional definition of the structural relaxation time, e.g. the point where F (s) (k, τ α ) = 1/e. This also implies that in the case of highly polydisperse mixtures, some caution is warranted when studying structural relaxation: on general grounds one expects the relaxation time to be defined as some point beyond which the correlator has decayed from its plateau by a certain amount, yet for the smallest particles the plateau height can become incredibly low. Hence, the meaning of 'structural relaxation' may be qualitatively affected when using the definition F (s) σ (k, τ α ) = 1/e (indicated by the horizontal black line in Fig. 2), since for small particles this would probe only the short-time dynamics.
As can already be inferred from Fig. 2, the Lamb-Mössbauer factor, which is essentially the height of the plateau of the ISF, is highly species dependent. This indicates that local mobility is strongly dependent on particle diameter, since the inverse Fourier transform of the LM factor is related to the Van Hove function, i.e. the probability of observing a given displacement during an infinite observation window [14]. More precisely, the half-width of the LM factor is inversely proportional to the root mean squared displacement of a given particle [13], which we refer to as a localization length. Let us first consider how the diameter-averaged form factors are affected by the polydispersity. In Fig. 3 we show the species-averaged LM factors f (s) (k) ≡ x σ f (s) σ (k) evaluated at the MCT critical point for various values of the polydispersity index. For both distributions and for all wavenumbers, the LM factors decrease monotonically with increasing δ. Additionally, they decay over longer wave-lengths when the polydispersity index is increased, suggesting that the average localization length of all particles is increased by increasing size polydispersity. Furthermore, we find that the averaged LM factor is progressively less Gaussian with increasing δ, which indicates that the local energy landscape around a given particle is progressively anharmonic [14].
To understand this, we argue that as the polydispersity index is increased, the notion of a local trapping cage becomes less relevant for the smallest particles. We show in Fig. 4 the diameter-resolved LM factors on a semi-logarithmic scale. For low degrees of polydispersity [panels (a)-(b)], we find the Gaussian result f (s) σ (k) ∝ exp(−(kD) 2 ), which suggests that all particles effectively behave as if they were in a harmonic trap [14]. As the polydispersity index is increased [panels (c)-(d)], we find that, while the largest particles retain a Gaussian LM factor, the smallest ones strongly deviate from one. Furthermore, the latter tend to decay faster as a function of k, signaling a much larger localization length. This supports the idea that the smallest particles are able to navigate heterogeneously (while still being localized) through narrow log(t)  channels within the matrix of larger, more strongly localized ones. The red curves in Fig. 4 which represent the diameteraveraged LM factor demonstrate that at high polydispersity indices [panels (c), (d)], the non-Gaussian character of the diameter-averaged LM factor can be imputed to the smallest particles. These findings are consistent with earlier results on bidisperse mixtures [26], as well as experimental and computational work on continuously polydisperse systems [10,52]. We find however that MCT is incapable of resolving another known observation in strongly polydisperse mixtures, namely that small particles escape their cage much earlier than large ones. Indeed, such reports have been made for simulated binary mixtures in the past [28], and more recently for the same inverse cubic distribution used in the first part of this work with δ ≈ 0.23 [11]. In these two studies, it was found that there is a dynamical separation between the smallest and the largest particles. Yet, MCT's asymptotic predictions are such that near the critical point considered here, species-resolved intermediate scattering functions decay from their respective plateaus at the same time, thus failing to cap-ture the dynamical separation. This failure of MCT suggests that the theory effectively overestimates the coupling between different species, ultimately leading to an ideal glass transition scenario. This limitation is attributed to the general absence of activated events in the theory, which, if included, would manifest as additional relaxation channels, restoring ergodicity. The way in which these couple to the degree of polydispersity remains however an open question.

B. Critical Dynamical Exponents
Let us now focus on the critical exponents associated with the dynamical MCT scaling laws. We recall that these exponents are generally system-dependent, and therefore also depend on the chosen size distribution. However, given the qualitatively equivalent behavior of the uniform and inverse cubic size distributions above, we here consider only the case where the particle diameters follow a uniform distribution. Similar qualitative trends are expected to hold for other continuously polydisperse mixtures.
Although the critical scenario of MCT is unaffected by polydispersity (at least at reasonable degrees of polydispersity), we nonetheless observe important quantitative effects of polydispersity on the critical exponent γ and the dynamical exponents a, b. We determine γ by fitting the power law to the relaxation time τ α ∼ |φ − φ c | −γ for a set of numerical solutions in the range 10 −2 ≤ |φ − φ c | ≤ 10 −5 . The dynamical exponents a, b are then determined by solving the system Eq. (8) using standard root finding techniques.
We find that the exponent γ increases monotonically with the level of polydispersity δ, as shown in Fig. 5(a). This increase implies a stronger divergence of the relaxation time, and thus, a higher degree of fragility with increasing polydispersity. Note here that we quantify the degree of fragility of a supercooled liquids as the deviation from an Arrhenius scaling of the relaxation time in terms of a power law, whereas fragility is generally measured from fits to the Vogel-Fulcher-Tammann equation [53]. Concomitantly, the dynamical exponents a and b are found to decrease monotonically with δ [see Fig. 5(b)-(c)]. The quantitative values of the critical exponents that we report here are also in line with earlier studies of polydisperse quasi-hard sphere mixtures [30]. Moreover, we find that all exponents saturate around a polydispersity index of δ = 0.4, which coincides with the point at which the predicted value of φ c becomes constant for this size distribution (see Fig. 1). These results further confirm that, within MCT, there is an upper bound for δ beyond which the dynamics are no longer affected by polydispersity for this distribution.
On general grounds, the critical exponents can be related to a dynamical length scale within the framework of inhomogeneous MCT [27]. Briefly, the growth of the (non-linear) dynamical susceptibility (i.e. the response of the coherent ISF to an infinitesimal localized perturbation of the density field near criticality) is governed by the values of the two dynamical exponents a and b. The dynamical length scale, which measures the size of a strongly dynamically correlated region, grows as ξ(t) ∝ t a/2 . Since a decreases with δ, this implies that increasing the polydispersity index hampers the growth of the dynamical length scale. This observation is in agreement with previous experimental findings [32], where it was found that lower polydispersities imply stronger dynamical heterogeneities. The above result is also consistent with the intuitive picture that we propose here: there exists a sub-population of particles with large localization lengths (the smallest ones), which fluidizes the system and thus makes it more difficult for the correlation length to grow in time, due to the inherent scrambling of the small particles in the voids surrounding the largest ones.

C. Stretched Exponential Relaxation
It is well known that the coupling of different relaxation modes can lead to stretched exponential behavior [54]. In standard MCT, such stretching is usually attributed to the coupling of different wave-vectors, but the additional coupling of different particle species can lead to further stretching of the dynamics. Indeed, our species-averaged dynamics for the intermediate scattering functions (full lines in Fig. 2), as well as the averaged von Schweidler exponent b (Fig. 5), indicate that additional relaxation channels provided by size polydispersity yield a stretching of the relaxation.
In order to more quantitatively analyze the influence of polydispersity on the structural relaxation, we perform a detailed analysis of the long-time tail of the incoherent ISFs. Specifically, we fit a Kohlrausch function to the long-time de-cay of the species-averaged correlator F (s) (k, t), in order to extract the Kohlrausch-Williams-Watts (KWW) stretching exponent β KWW (see Appendix D for details). Note that, within the context of MCT, this asymptotic form is strictly valid only in the infinite wave-number limit k → ∞ [46,55], but it is also frequently employed to fit MCT predictions at finite wave-numbers [17,30] and it is a standard quantifier of relaxation in complex media [54]. Our results for β KWW for a uniform size distribution are presented in Fig. 6(a). As the polydispersity degree is increased, we find that β KWW decreases monotonically until it saturates beyond δ ≈ 0.4 at a low value around β KWW ≈ 0.60. This saturation, which we recall is also observed in the state diagram and in the critical exponents, again strengthens the idea there is a degree of polydispersity beyond which MCT's prediction are no longer affected. We subsequently perform a diameter-resolved analysis of the relaxation by computing the KWW exponent from the partial incoherent ISFs, which we denote β KWW . The results are shown in Fig. 6(b) and (c) for δ = 0.24, 0.38 respectively. The KWW exponent monotonically increases with the particle diameter in a continuous mixture. Furthermore, we find that the species-averaged exponent [dashed black lines in Fig. 6 (b)-(c)] is very close to that of the value for the smallest particles.
Since the stretching exponent is related to the degree of heterogeneity in relaxation timescales, it is not surprising that the diameter-averaged exponent is in fact governed by that of the more dynamically heterogeneous particles (i.e. the smallest ones). This interpretation is coherent with the notion that the growth rate of correlated regions is restrained by polydispersity, since the smallest particles are much more mobile than the largest ones.

D. Competing Relaxation Channels
Lastly, we study the superposition of relaxation mechanisms induced by size polydispersity in a mixture with a uniform size distribution. To this end it is useful to translate our MCT results into frequency space. We follow standard procedure [56] and define the imaginary part of the response spectra Im[χ(k, ω)] ≡ χ ′′ (k, ω) associated with F (s) (k, t) as and analogously for the species-specific quantities, which we denote by χ ′′ σ (k, ω). Figure 7(b) shows the response spectra for different polydispersity degrees at a fixed relaxation time τ α , rescaled by the average LM factor such that the low-frequency peaks (associated with structural α-relaxation) all collapse at the same height. Note that minor deviations from a perfect collapse are noticeable, which are attributed to our numerical accuracy in determining solutions at fixed relaxation time. It can be seen that the α-peak broadens toward higher frequencies as the polydispersity index δ increases. This fanning is associated with an increased stretching of the long-time structural relaxation of F (s) (k, t), in line with the results of Fig. 6(a) and explicitly shown in Fig. 7(a). Compared to the monodisperse case (full line in Fig. 7), we infer that polydispersity induces additional relaxation channels that now compete with each other, thus resulting in a significantly broader spectrum [see Fig. 7(b)]. We also find that as the degree of polydispersity increases, the relative height of the boson peak (i.e. the high-frequency peak of the spectrum) increases with respect to that of the α-peak. Since the two spectral peaks are separated by the caging regime (i.e. the minimum of the spectrum), this further corroborates that the degree of polydispersity influences the ratio between structural relaxation before and after the caging regime. In particular, for the highest polydispersity index considered in this work, the magnitude of the boson peak exceeds that of the principal α-relaxation, which implies in the time domain [ Fig. 7(a)] that more than 50% of the decay in density correlations occurs before reaching the caging plateau.
Let us now inspect the α-peak of the spectrum more closely. The dashed lines in Fig. 8(a), (b) represent the spectrum of the fitted stretched exponentials at polydispersities δ = 0.24, 0.38 for varying packing fractions towards the critical point. For all cases, we find that the signal is well captured up to the vicinity of the spectral minimum, where a smooth transition occurs between stretched exponential relaxation and the von Schweidler excess: χ ′′ (ω) ∼ ω −b [shown for clarity in panels (c) and (d)]. Note however that this tran-sition is hard to see, since the von Schweidler exponent and the slope of the stretched exponential are very similar. In order to understand the microscopic origin of the excess in the relaxation spectrum, we study the diameter-resolved susceptibility spectra, shown in Fig. 8(c) and (d) by the colored lines. Here it is important to recall that the asymptotic properties of MCT around the minimum of the spectra are universal. In particular, the theory predicts that the von Schweidler exponent b for the species-resolved susceptibili- [29,41] [although it is system dependent, as previously shown in Fig. 5(c)]. The results of Fig. 8(c) and (d) reveal that the particle-resolved susceptibilities (full colored lines) indeed all share the same power laws, but the location of the low-frequency α-peak is however diameter dependent. To the left of the α-peak, the highest curve corresponds to that of the largest diameter while the lowest curve to the smallest diameter. At the minimum ω 0 the reverse situation is found, as indicated by the black arrows pointing towards increasing diameters. Hence, the curves must cross at some point beyond which the smallest particles contribute 'in excess' to the spectrum.
In fact, MCT guarantees that the distance between the αpeaks of particles with different diameters is fixed as one approaches the transition, since they all share the same critical exponent γ. Hence, as we get closer to the critical point the spectral contributions from the smallest particles at frequencies past the caging regime do not increase in magnitude, and the entire effect remains subtle. However in molecular dynamics simulations, we know that this is not valid, and that in fact, beyond the spurious transition predicted by MCT, the distance between the α-peaks of particles with different diameters must increase, since the time interval between which the smallest and the largest particles relax increases with deeper supercooled systems [11]. This raises an important question regarding the exact nature and the true microscopic origins of relaxational excess observed in the spectrum of simple, yet polydisperse glass formers. In particular it would be of interest to compare the value of the excess wing exponent reported in [38] with that of the von-Schweidler exponent for the same system, and to perform a similar diameter-resolved analysis of the relaxation spectra in order to check whether or not there is a homogeneous spectral contribution with respect to particle diameters.

VI. CONCLUSION
In this work, we have analyzed the dynamics of continuously polydisperse glass-forming liquids through the lens of Mode-Coupling Theory. We have demonstrated that increasing the degree of size polydispersity first slows down the dynamics, before pushing the critical point φ c to higher packing fractions, thus stabilizing the liquid phase. Of the three size distributions that we have considered, the liquid stabilization effect is largest for the inverse cubic distribution. This generalizes earlier MCT results on simpler, binary mixtures. A detailed study of the Lamb-Mössbauer factors has  revealed that the smaller particles have much longer localization lengths than the larger ones, and strongly deviate from the standard Gaussian behavior. This leads to visible deviations from Gaussian behavior also in diameter-averaged quantities. Our analysis of the particle-averaged dynamics near the critical point shows that increasing the degree of polydispersity induces a significant stretching of the KWW exponent, a lowering of the critical dynamical exponents a, b, and a net growth of the critical exponent γ that governs the divergence of the relaxation time. Our results thus indicate that increasing the polydispersity index affects the fragility of the liquid, as well as dynamical heterogeneities whose growth is dictated by the dynamical exponents a and b [27,57].
By performing a diameter-resolved analysis, we have shown that the Kohlrausch-Williams-Watts stretching exponent of the diameter-averaged dynamics is dominated by the dynamics of the smaller particles, which have significantly lower Kohlrausch-Williams-Watts exponents compared to their larger counterparts. This leads us to infer that temporal heterogeneities in the dynamics are also mainly due to the smaller particles. This conclusion also is in partial agreement with recent observations in molecular dynamics simulations [11], where it has been shown that smaller particles have strongly heterogeneous dynamics compared to their larger counterparts. We have further studied the manifestation of this effect in the susceptibility spectra of the relaxation, which exhibit a broadening and von Schweidler excess that becomes more pronounced upon increased polydispersity. We can attribute this broadening to two factors: firstly, to the increasing temporal heterogeneities of the small particles as they exhibit low stretching exponents (compared to the larger ones), and secondly to the shift in spectral peak positions of the smallest particles. It is important to note, however, that the von Schweidler excess remains strictly governed by the critical exponent b, but since b itself changes with polydispersity, so does the spectrum. A detailed diameter-resolved investigation using computer simulations could verify if these observations persist beyond the MCT regime.
Overall, our results demonstrate that size polydispersity imposes important and non-trivial effects on glassy dynamics. This is particularly relevant in the context of deeply supercooled liquids that typically require a large degree of polydispersity in order to permit equilibration at temperatures below the MCT crossover temperature. Our work reveals, on purely theoretical grounds, that even in the mildly supercooled regime where MCT is usually deemed applicable, the presence of polydispersity can induce complex dynamical features. It is possible that some of these polydispersity-specific features might carry over into the more deeply supercooled regime.
Furthermore, we have shown that strongly continuously polydisperse mixtures of hard spheres display features of both collective glassy behavior for the largest particles, and localization in narrow channels for the smallest ones. In this sense, such systems provide an interesting intermediate step between effectively monodisperse glasses and the more minimal model of the random Lorentz gas used to study transport in heterogeneous environments [58], which is governed by percolation in physically relevant dimensions.  We first discuss the validity of the discrete representation of a continuous distribution of particle sizes. We show in Fig. 9 the structure factors for a uniform continuous distribution of particle diameters, as determined from multi-component Percus-Yevick hard-spheres, for two degrees of polydispersity δ = 0.2, 0.4 as a function of an effective n-component representation at fixed packing fraction φ = 0.521. The top panels show the convergence of the peak of the averaged structure factor S(k * ) as a function of the number of resolved components. We judge that an effective 10-component description is sufficient to accurately resolve the total structure of the system, motivating the choice in the main text. Similar results are found for other degrees of polydispersity, packing fractions, as well as for other type of distributions. To further illustrate this point we plot in Fig. 10 the convergence of the resulting MCT solutions as a function of the number of resolved species for the diameter-averaged incoherent intermediate scattering function. We find that the curves for n = 6, 8 and 10 effective components have converged to the same value. To make the convergence manifest we have dashed the n = 10 curve. Similar results are found for other log(t) degrees of polydispersity and other ranges of packing fractions. Lastly we show in Fig. 11 that the qualitative picture of the phase diagram is independent of the number of resolved diameters for a uniform distribution of diameters. Similar results are found for the other distributions.

Appendix B: Numerical Details
The equations of motion (1)-(4) were solved using a timedoubling algorithm [28,46], which has been recently been implemented in an open-source solver for integro-differential equations of the mode-coupling type [31]. The wave-vector integrals were performed (in bi-polar coordinates) over a grid of N k = 100 equidistant points between k min = 0.2 and k max = 39.8. We have checked the stability of our results against a finer grid with k min = 0.2 and k max = 60.0, with N k = 300 equidistant points in between (which doubles the resolution compared to the one presented in the main text). Results are shown in Fig. 12; we observe no qualitative differences between the two wave-vector grids.

Appendix C: Static Structure Factors
We show in Fig. 13 the averaged static structure factor S(k), as determined from multi-component Percus-Yevick hard-spheres, for various polydispersity indices, evaluated at the critical point φ c for an effective 10-component continuously polydisperse mixture. For both distributions, as the degree of polydispersity is increased, we see that the magnitude of the first peak of the averaged structure factor drastically diminishes. In the case of the uniform distribution, we find that the peak is modestly shifted to lower wave-numbers, signaling a slight increase in the average cage size. We see however that in the case of the inverse cubic distribution, the structure is completely washed out as the polydispersity index is increased, and in fact the qualitative aspect of the structure factor is significantly changed by the abundance of very small particles.