Theory for Swap Acceleration near the Glass and Jamming Transitions

I Introduction

Understanding the mechanisms underlying the slowing down of the dynamics near the glass transition is a long-standing challenge in condensed-matter Ediger et al. (1996); Angell (1991). Unexpectedly, swap algorithms Grigera and Parisi (2001); Fernández et al. (2006) (in which particles of different radii can swap in addition to the usual moves of particle positions) were recently shown to allow for equilibration of liquids far below the glass transition temperature $T_{g}$ Berthier et al. (2016a); Gutiérrez et al. (2015); Ninarello et al. (2017); Berthier et al. (2017). For judicious choice of poly-dispersity, one finds that: $(i)$ the glass transition is shifted to lower temperatures: with swaps the $\alpha$-relaxation time at $T_{g}$ is only two or three order of magnitudes slower that in the liquid, instead of 15 orders of magnitude for regular dynamics. The slowing down of the dynamics occurs at a lower temperature, which we refer to as $T_{0}^{\mbox{\scriptsize swap}}$. $(ii)$ The spatial extent of dynamical correlations, which are significant near $T_{g}$, are greatly reduced with swap and only occur at $T_{0}^{\mbox{\scriptsize swap}}$. $(iii)$ The mean square displacement of the particles on vibrational time scales is increased significantly in this temperature range Ninarello et al. (2017). These observations constrain theories of the glass transition. In particular, current formulation of theories based on a growing thermodynamic length scale appear inconsistent with these observations Wyart and Cates (2017). A theory of the glass transition should explain both swap and non-swap dynamics. Goldstein Goldstein (1969) proposed that the glass transition is initiated by a transition in the free energy landscape: at high temperature, the system resides near saddles, whereas below some temperature $T_{0}$ the dynamics can only occur by activation (whose nature is debated), and is thus much slower. In mean-field models of structural glasses such a transition in the landscape is predicted Lubchenko and Wolynes (2007); Broderix et al. (2000); Grigera et al. (2002); Kurchan et al. (2016) and corresponds to a Mode Coupling Transition (MCT) where the relaxation time diverges. It was suggested that the MCT transition would be shifted to lower temperature with swap dynamics in Wyart and Cates (2017), as proven and confirmed numerically in a mean-field model of glasses Ikeda et al. (2017). Yet, understanding the real-space mechanisms underlying the speed up induced by swap in finite dimensions (where the relaxation time cannot diverge) as well as the nature of the very stable glassy configurations swap can reach remains a challenge.

In this work, we tackle these questions by first reviewing the equilibrium statistical mechanics theory of polydisperse systems Briano and Glandt (1984); Kofke and Glandt (1987), to show that they are equivalent to a system of identical particles that can individually deform according to a chemical potential $\mu(R)$, where $R$ is the particle radius. In the (practically important) case where poly-dispersity is continuous, $\mu(R)$ is smooth, allowing us to define normal modes of the generalised Hessian that includes radii as degrees of freedom. We prove that requiring its stability is strictly more demanding than for the usual Hessian. Second, we show that these results stringently constrain the glassy states generated by swap algorithms. We illustrate this point by studying the jamming transition in soft repulsive particles, which we prove must be profoundly altered: hyper-staticity is found with an excess number of contacts $\delta z$ with respect of the Maxwell bound $\delta z\!\sim\!\alpha^{1/2}\!>\!0$, where $\alpha$ characterises the width of the radii distribution $\rho(R)$. Although we find that the vibrational spectrum of the generalised Hessian is marginally stable with respect to soft extended modes near jamming, these modes are gapped in the regular Hessian, unlike for packings obtained with regular dynamics Wyart et al. (2005); DeGiuli et al. (2014a); Franz et al. (2015). These results are verified numerically by introducing a novel algorithm performing a steepest descent in the generalised potential energy that includes $\mu(R)$, which can generate extremely stable glasses without any activation. Third, we investigate the glass transition. We show that the inherent structures obtained after a rapid quench with the regular dynamics are unstable with respect to this new algorithm, which reaches significantly smaller energies. This result indicates that metastable states appear at lower energies with swap, and therefore at lower temperatures when the liquid is equilibrated. Thus the Goldstein transition must be shifted to a lower temperature with swap dynamics, suggesting a natural explanation for its speed up which specifies the collective modes facilitating the dynamics for $T_{0}^{\mbox{\scriptsize swap}}<T<T_{g}$. We predict this shift to be proportional to $\alpha$ in general, and to be inversely proportional to the distance to jamming for sufficiently compressed soft spheres. Lastly, we argue that these results apply to hard spheres as well, if the energy is replaced by a coarse-grained free energy landscape as previously studied in Brito and Wyart (2006, 2009); DeGiuli et al. (2014a); Altieri et al. (2016). We use this approach to provide a simple phase diagram where the Goldstein transition and the emergence of marginality Wyart et al. (2005) (referred to as a Gardner transition in infinite dimension Altieri et al. (2016)) can be related to structure for both swap and non-swap dynamics.

II Grand-canonical description of poly-disperse systems

We now show that a poly-disperse system can be described by an effective potential that includes the radius as degree of freedom, an idea already present in the early works of Briano and Glandt (1984); Kofke and Glandt (1987). We consider a system of $N$ particles with continuous polydispersity $\rho(R)$, of width $\alpha\!=\!\langle(\langle R^{2}\rangle\!-\!\langle R\rangle^{2})^{1/2}\rangle/\langle R\rangle$. Here $\{R\}$ indicates the set of particle radii and $\{r\}$ their positions. In what follows we denote $\langle R\rangle\equiv R_{0}$, and ${\cal U}(\{r\},\{R\})$ the total potential energy in the system. We define the partition function $Z(\{r\})$ annealed over the particle radii:

where the sum is on all the permutations ${\cal P}(\{R\})$ of the particle radii. In the thermodynamic limit, a grand-canonical formulation is equivalent, in which particles of different radii correspond to different species. The associated partition function writes:

where $\mu(R)$ is the chemical potential at radius $R$. It is chosen such that in the thermodynamic limit, the distribution of radii that follows from Eq. (S.2) is

A key remark is that once Eq. (S.2) is integrated on particle positions $\{r\}$, one obtains the partition function for the coupled degrees of freedom $\{r\}$ and $\{R\}$ with an effective energy functional:

II.2 Mechanical stability under swap

Let us consider first an athermal system. In the thermodynamic limit, mechanical stability under swap dynamics requires ${\cal V}$ to be at a minimum. Beyond the usual force balance condition, it implies:

where $f_{ij}$ are the contact forces between particle $i$ and $j$ (positive in our notations for repulsive forces), and ($\{r^{*}\}$, $\{R^{*}\}$) the particle positions and radii at the minimum. For unimodal distribution $\rho(R)$, one expects $\mu(R)$ to be unimodal too. In an amorphous solid the fluctuations of the left hand side of Eq. (S.5) are of order $pR_{0}^{d-1}$ where $p$ is the pressure and $d$ the spatial dimension. To achieve a distribution of radii of width $\alpha$, the stiffness $k_{R}$ acting on each particle radius must thus be of order:

where the average is made on all particles i.

II.3 Generalised vibrational modes

Stability also requires the Hessian $H_{\mbox{\scriptsize swap}}$ (the matrix of second derivatives of ${\cal V}$) to be positive definite. Since they are now $N(d+1)$ degrees of freedom, $H_{\mbox{\scriptsize swap}}$ is a $N(d+1)\times N(d+1)$ symmetric matrix, of eigenvalues $\omega_{\mbox{\scriptsize swap}}^{2}$. It contains a block of size $Nd\!\times\!Nd$ which is the regular Hessian $H_{ij}\!=\!\partial^{2}U/\partial r_{i}\partial r_{j}$. We denote by $\omega^{2}$ its eigenvalues. Because hybridisation with additional degrees of freedom can only lower the minimal eigenvalues of the Hessian, $H_{\mbox{\scriptsize swap}}$ has lower eigenvalues than $H$ (as quantified below), implying that mechanical stability is more stringent with swap dynamics.

Let us illustrate this result perturbatively when $k_{R}\!\gg\!k$, where $k$ is the characteristic stiffness of the interaction potential ${\cal U}$. In general, the eigenvalues of $H$ are functions of the set of stiffnesses $\{k_{ij}\}$, but also of the interaction forces $\{f_{ij}\}$ Alexander (1998). We first ignore the effects induced by such pre-stress. Moving along a normal mode of $H$ by a distance $x$ (while leaving the radii fixed) leads to an elastic energy $\sim\!\omega^{2}x^{2}$ and change forces by a characteristic amount $\delta f$ satisfying $\delta f^{2}/k\!\sim\!x^{2}\omega^{2}$. Because of such change, Eq. (S.5) is not satisfied anymore. Thus the potential ${\cal V}$ can be reduced further by an amount of order $\delta f^{2}/k_{R}\!\sim\!\omega^{2}x^{2}k/k_{R}$ if the radii are allowed to adapt. This reduced energy can be approximatively written as $x^{2}\omega_{\mbox{\scriptsize swap}}^{2}$, where $\omega\mbox{\scriptsize swap}^{2}$ is the eigenvalue associated to that mode in the effective Hessian. We thus obtain:

III Soft sphere systems

To illustrate these ideas, we consider soft spheres with half-sided harmonic interactions, so that:

where $r_{ij}$ is the distance between particles $i$ and $j$, and $\Theta(x)$ is the Heaviside step function.

When materials with such finite-range interactions are quenched to zero temperature, they can jam into a solid or not depending on their packing fraction. At the jamming transition separating these two regimes, vibrational properties are singular Liu et al. (2010); Wyart (2005), and the effects of swap are expected to be important, as we now show. The vibrational spectrum of the regular Hessian is strongly affected by excess number of constraints $\delta z$ with respect to the Maxwell threshold where the numbers of degrees of freedom and constraints match. Effective medium Wyart (2010) or a variational argument Yan et al. (2016) imply that in the absence of pre-stress, soft normal modes in the Hessian must be present with eigenvalues:

For swap with a small poly-dispersity $\alpha<<\Delta$, where we introduced the dimensionless particle overlap $\Delta\equiv pR_{0}^{d-2}/k$, then from Eq.S.6 $k_{R}>>k$ and Eqs. (S.7) applies. It implies that soft normal modes will be present at lower eigenvalues $\omega_{\mbox{\scriptsize swap}}^{*}\!{}^{2}\!\sim\!\delta z^{2}k(1-C_{0}\alpha/\Delta)$ where $C_{0}$ is a numerical constant. Pre-stress can be shown to shift eigenvalues of the Hessian by some amount $\approx\!-C_{1}k\Delta$ Wyart et al. (2005); DeGiuli et al. (2014b), leading to eigenvalues satisfying $\omega_{\mbox{\scriptsize swap}}^{0}\!{}^{2}\!\sim\!\delta z^{2}k(1-C_{0}\alpha/\Delta)-C_{1}k\Delta$. Mechanical stability requires positive eigenvalues and we obtain:

Eq.S.10 indicates that away from jamming, the relative effects of swap on the structure are proportional to $\alpha/\Delta$. Certain materials are marginal stable, corresponds the saturation of inequalities of the kind of Eq. (S.10). As we shall see below, we provide numerical evidence that it is also the situation if swap is used, at least near jamming. Here this assumption gives an expression for $\delta z$ which is above (but very close to in the limit $\alpha<<\Delta$) the bound for non-swap dynamics of Wyart et al. (2005), recovered by setting $\alpha\!=\!0$. Thus in this limit we expect very small change of structure in the glass phase between swap and non-swap dynamics.

For swap with a large poly-dispersity $\alpha>>\Delta$, the situation is completely different. We then have $k_{R}<<k$: in this regime the strong interactions correspond to interactions between particles in contact. As far as the low-frequency end of the spectrum is concerned, these interactions can be considered to be hard constraints (i.e. $k=\infty$), whose number is $Nz/2$. The dimension of the vector space satisfying such hard constraints is $N(d+1)\!-\!Nz/2\!=\!N\!(1-\!\delta z/2)$. These modes gain a finite frequency due to the presence of the weaker interactions of strength $k_{R}$ associated with the change of radius, of strength $k_{R}$. Importantly, the number of these weaker constraints left is simply the number of particles $N$. If $\delta z$ is small the number of degrees of freedom $N(1-\delta z/2)$ is just below the number of constraints $N$: for this vector space we are close to the ”isostatic” or Maxwell condition where the number of constraints and degrees of freedom match. Thus we can use the same results for the spectrum valid near the jamming transition introduced above. They also apply in that situation, with the only difference that the stiffness scale $k$ is replaced by $k_{R}$. In particular if pre-stress is not accounted for, a plateau of soft modes must appear above some frequency given by Eq.(S.9):

This plateau survives up to the characteristic frequency $\omega_{i}\!\sim\!\sqrt{k_{R}}$. When pre-stress is accounted for, eigenvalues of the Hessian are again shifted by $\sim\!-k\Delta$. Mechanical stability then implies $\Delta/\alpha\delta z^{2}\!>\!C_{2}\Delta$ and:

In this regime, marginal stability (the saturation of the stability bound of Eq. S.12) corresponds to a coordination independent of pressure, with $\delta z\sim\sqrt{\alpha}$ and $\omega^{*}_{\mbox{\scriptsize swap}}{}^{2}\!\sim\!k\Delta$. We thus predict that swap dynamics destroys isostacity, and significantly affect structure and vibrations. For sufficiently large $\alpha$, this regime will include the entire glass phase, and vibrational properties and stability will be affected in the vicinity of the glass transition (which sits at a finite distance from the jamming transition Ikeda et al. (2012)) as well.

Note that these predictions apply to algorithms that allow for swap moves up to the jamming threshold. This is not the case e.g. in Coslovich et al. (2017), where swaps are used to generate dense equilibrated liquids, that are then quenched without swap toward jamming. We also expect isostaticity to be restored in algorithms for which the set of particle radii is strictly fixed, but only below some pressure $p_{N}$ that vanishes as $N\!\rightarrow\!\infty$, above which our predictions should apply.

III.2 Numerical Model

As shown in Eq. (S.4), swap dynamics is equivalent to a system of interacting particles which can individually deform. To test our predictions, we consider soft spheres as defined in Eq.S.8, whose radius follow the internal potential:

where ${\bar{k}}_{R}$ is a characteristic stiffness. We considered a potential diverging as $R_{i}\!\rightarrow\!0$ to avoid particles shrinking to zero size. To avoid crystallisation we further considered that particles are of two types: for 50% of them, $R_{i}^{(0)}\!=\!0.5$ while for the others $R_{i}^{(0)}\!=\!0.7$. This choice leads to a bimodal distribution of size $\rho(R)$, as shown in Fig. S1. Our model corresponds to a swap dynamics where swap is allowed only between particles of the same type. Note that broad mono-modal distributions can be optimised to make swap more efficient while avoiding crystallisation Ninarello et al. (2017), which would be similar to having a very large $\alpha$ in our theoretical description. The spatial dimension is $d=2$ in our simulations and $k=1$ is our unit stiffness, leading to a simple relation $\Delta=p$.

To study the jamming transition, we consider a pressure-controlled protocol at zero temperature described in the S.I. The chemical potential of Eq. (S.13) must evolve with pressure to maintain a fixed polydispersity. As shown in Fig. S1.B, it can be achieved within great accuracy simply by imposing that ${\bar{k}}_{R}\!=\!p/{\bar{\alpha}}$, where ${\bar{\alpha}}$ is a parameter that controls the width $\alpha$ as shown in the inset of Fig. S1. For this bimodal distribution, $\alpha$ is defined as $\alpha=(\alpha_{1}+\alpha_{2})/2$, where $\alpha_{1}$, $\alpha_{2}$ are the relative width of each peak in $\rho(R)$. In the limit where the non-swap dynamics is recovered – which happens when ${\bar{\alpha}}\rightarrow 0$ – $\alpha$ and ${\bar{\alpha}}$ are proportional.

III.3 Structure and stability

Our central prediction is that for swap dynamics materials must display a larger coordination to enforce stability. This prediction is verified in Fig. S2.A, which shows $\delta z$ v.s. $p$ for various values of ${\bar{c}}$. Isostaticity is indeed lost and the coordination converges to a plateau as $\Delta$ decreases. Strikingly, we find for the plateau value $\delta z\!\sim\!\sqrt{\alpha}$, consistent with a saturation of the stability bound of Eq. (S.12). This scaling behavior is implied by the scaling collapse in Fig. S2.B which also confirms that the characteristic overlap below which swaps affects the dynamics scale as $\Delta\!\sim\!\alpha$. Overall, these results supports that the numerical curves $\delta z(\Delta)$ in Fig. S2.A correspond to the marginal stability lines under swap dynamics, shown for different polydispersity (see more on that below).

III.4 Packing fraction

For traditional dynamics, polydispersity tends to have very limited effects on the value of jamming packing fraction $\phi_{J}$. We have confirmed this result in the S.I., by showing that although our model can generate very different distributions $\rho(R)$, the values we obtain for $\phi_{J}$ cannot be distinguished if jamming is investigated using non-swap dynamics. However, for swap dynamics we expect the situation to change dramatically: since stability requires much more coordinated packings, they presumably need to be denser too. We denote the jamming packing fraction for swap $\phi_{c}\!\equiv\!\lim_{\Delta\rightarrow 0}\phi(\Delta)$. The inset of Fig. S3.A confirms that $\phi_{c}$ increases significantly as $\rho(R)$ broadens. To quantify this effect we consider $\phi(\Delta,\bar{\alpha})$, as shown in the main panel. Assuming a scaling form for this quantity, and requiring that it satisfies the known results for the jamming transition for $\Delta\!\gg\!\bar{\alpha}$ implies $\phi(\Delta,\bar{\alpha})\!-\!\phi_{J}\!=\!f(\Delta/\bar{\alpha})\bar{\alpha}^{\beta}$ where $f(x)$ is some scaling function and $\beta\!=\!1$. Since the coordination does not change for $\Delta\!\ll\!\bar{\alpha}$, we expect that it is true for the structure overall and for $\phi$, implying that $f(x)\!\sim\!x^{0}$ as $x\!\rightarrow\!0$. These predictions are essentially confirmed in Fig. S3.B. Note however that the best scaling collapse is found for $\beta\!=\!0.83\!<\!1$. These deviations are likely caused by finite size effects, known to be much stronger for $\phi$ than for the coordination or vibrational properties O’Hern et al. (2003), and which may thus be present for our systems of $N\!=\!484$ particles.

III.5 Vibrational properties

We computed the Hessian $H_{\mbox{\scriptsize swap}}$ and diagonalized it (see detailed in SI) to extract the density of states $D(\omega)$, as shown in Fig. S4.A for different pressures at fixed polydispersity. As expected, at low particle overlap $\Delta$ two bands appear in the spectrum. The lowest-frequency band presents a plateau above some frequency scale $\omega^{*}_{\mbox{\scriptsize swap}}$ which satisfies $\omega^{*}_{\mbox{\scriptsize swap}}\!\sim\!\sqrt{\Delta}$ as shown in the inset, as expected if the structure were marginally stable. As shown in S.I., in the absence of pre-stress the minimal eigenvalues of the Hessian increase many folds, again a signature of marginal stability Wyart et al. (2005). Further evidence appears in Fig. S4.B showing $D(\omega)$ at fixed $\Delta\!=\!10^{-4}$ for varying polydispersity. $\omega^{*}_{\mbox{\scriptsize swap}}$ essentially does not depend on $\bar{\alpha}$ as shown in the inset, as expected for marginal packings if the pressure is fixed. The cut-off frequency $\omega_{i}$ of the low-frequency plateau scales as $\omega_{i}\!\sim\!\sqrt{k_{R}}\!\sim\!1/\sqrt{\bar{\alpha}}$, as predicted above.

IV Glass transition

We now turn to the glass transition, which always takes place at a sizeable distance form the jamming transition Ikeda et al. (2012): for example for hard discs, $\phi_{g}\approx 0.78$ and $\phi_{c}\approx 0.85$. A similar difference of packing fraction occurs by compressing soft spheres at overlap $\Delta\approx 0.05$, as illustrated in Fig.S5(d). From the arguments above, we expect that if the poly-dispersity is sufficiently large, vibrational properties will be strongly affected even far away from jamming, in particular near the glass transition.

The direct consequence of this fact is that the energy landscape will be affected by swap, which will in turn affect the glass transition. At high energy configurations are unstable — they are saddles with many unstable directions — whereas below some characteristic energy minima appear. However, since stability is strictly more demanding with swap, this characteristic energy must be reduced when swap is allowed for. We prove this point in Fig.S5.(a,b), where inherent structures of energy $U_{\infty}$ are obtained after using a steepest descent for non-swap dynamics. These configurations are not stable for our generalised steepest descent that let particles deform, which leads to configurations of energy $U_{\alpha}\!<\!U_{\infty}$. This effect is stronger near jamming in relative terms as shown in Fig. S5.(c), but remains significant away from jamming if the poly-dispersity is broad enough. It corresponds for example to a reduction of energy of 25% for $\Delta=0.05$ for our $\alpha=0.06$. We show in the inset of that panel that the relative shift of energy induced by swap $(U_{\infty}-U_{\alpha})/U_{\alpha}$ is proportional to $\alpha$ and inversely proportional to $\Delta$ when $\Delta$ is large enough, in consistence with what we found for the structure in Eq.(S.10).

Thus as the temperature is lowered in these liquids, the Goldstein temperature where activation sets in will be smaller when swap is allowed for. This analysis thus predicts an entire range of temperature where the non-swap dynamics is slowed down by activation, whereas the swap dynamics can flow along unstable modes. More quantitatively, we predict the shift of glass transition temperature $\Delta T_{g}/T_{g}$ induced by swap to be proportional to $\alpha$, in consistence with the observation that very broad distributions lead to large swap effects Ninarello et al. (2017). We also predict that $\Delta T_{g}/T_{g}$ is inversely proportional to the distance to jamming $\Delta$ when this quantity is well-defined (e.g. for soft spheres, but also to some extent for Lennard-Jones potentials Wyart (2005); Xu et al. (2007)) and large enough. In real space, the unstable modes that render activation useless involve both translational degrees of freedom as well as swelling and shrinking of the particles. We show an example of such a mode in Fig. S6, corresponding to the softest mode of the generalised Hessian we obtain with parameters $\alpha=0.06$ and $\Delta\!=\!10^{-2}$. It illustrates that the particle displacements are not necessarily divergent free when swap is allowed, since the system can locally compress or expand by changing the particle sizes.

This interpretation of swap acceleration is consistent with the observation that the dynamics is less collective with swap at the temperature where the non-swap dynamics is activated, since the system can rearrange locally without jumping over barriers if there are enough unstable modes. Collective dynamics is expected only when these modes become less abundant at lower temperatures. Likewise, we expect the Debye-waller factor to be larger with swap, since the vibrational spectrum is softer. Note that these arguments are not restricted to finite range interactions. We expect them to apply as well to Lennard-Jones potentials for example, where the abundance of degrees of freedom vs. the number of strong interactions is also known to affect the vibrational spectrum Wyart (2005); Xu et al. (2007).

V Hard sphere systems

Our arguments above consider the energy landscape. For interactions potentials which are very sharp, non-linearities induced by thermal fluctuations are important, and the vibrational properties of a glassy configurations at finite temperature $T$ can differ significantly from those of its inherent structure obtained by quenching it rapidly. Here we consider the extreme case of hard spheres where the energy is always zero, and cannot be used to define vibrational modes. Instead, by averaging on vibrational time scales within a glassy configuration, a local free energy can be defined Brito and Wyart (2006, 2009); Altieri et al. (2016) where particles that are colliding within that state interact with a logarithmic potential. This description is exact near jamming and systematic deviations are expected away from it Altieri (2018). However in practice, the Hessian defined from this free energy captures well the fluctuations of particle positions and the vibrational dynamics throughout the glass phase Brito and Wyart (2009). (This procedure can be pursued to include thermal effects in soft spheres as well DeGiuli et al. (2015)).

VI Conclusion

In swap algorithms, the dynamics is governed by an effective potential ${\cal V}(\{r\},\{R\})$ that describes both the particles interaction and their ability to deform. As a result, we have shown that vibrational and elastic properties are softened when swaps are allowed for, while thermodynamic quantities are strictly preserved (when thermal equilibrium is reached). This result supports that the cross-over temperature $T_{0}$ where mechanical stability appears and dynamics becomes activated must be reduced with swap with $T_{0}^{\mbox{\scriptsize swap}}\!<\!T_{0}$, leading to a natural explanation as to why the glass transition occurs then at a lower temperature $T_{g}^{\mbox{\scriptsize swap}}\!<\!T_{g}$. Secondly, swap must strongly affect the structure of the glass phase. This is particularly striking near the jamming transition that occurs in hard and soft spheres, where we predict that well-known key properties such as isostaticity must disappear. We have confirmed these predictions numerically, and found that for rapid quenches the effective potential ${\cal V}(\{r\},\{R\})$ appears to be marginally stable throughout the glass phase.

Concerning the glass transition, our work does not specify the mechanism by which activation occurs in glasses, but it does support that swap delays the temperature where activation is required to relax, which potentially explains several previous observations of swap algorithms Berthier et al. (2016a); Gutiérrez et al. (2015); Ninarello et al. (2017); Berthier et al. (2017). Possible theories to describe the mechanism by which activation occurs in glasses include elastic Dyre (2006) and facilitation models Garrahan and Chandler (2002). We do think however that theories based on a growing thermodynamic length will be hard to reconcile with the notion that some collective modes do not see any barriers at all.

Our analysis also makes additional qualitative testable predictions. By increasing continuously the width $\alpha$ of the radii distribution $\rho(R)$, we predict that $T_{g}^{\mbox{\scriptsize swap}}(\alpha)$ will smoothly decrease, while $T_{g}(\alpha)$ should be essentially unchanged, with $(T_{g}(\alpha)-T_{g}^{\mbox{\scriptsize swap}}(\alpha))/T_{g}(\alpha)\propto\alpha$ and more specifically $\propto\alpha/\Delta$ for soft spheres. Furthermore, many studies have analyzed correlations between dynamics and vibrational modes, see e.g. Broderix et al. (2000); Grigera et al. (2002); Widmer-Cooper et al. (2008); Brito and Wyart (2009), which can be repeated to relate the swap dynamics to the spectrum of the effective potential ${\cal V}(\{r\},\{R\})$. Near $T_{g}$, we predict the latter to have more abundant modes at low or negative frequencies than the much studied Hessian of the potential energy, and its softest modes to be better predictors of further relaxation processes. Lastly, the present analysis suggests that adding additional degrees of freedom (such as changing the shape of the particles, and not only their size) will increase even further the difference between swap and non-swap dynamics.

Finally, we have shown that ultra-stable glasses can be built on the computer, simply by descending along the effective potential ${\cal V}(\{r\},\{R\})$. As illustrated in Fig.S7, these configurations must sit strictly inside the stable region of the regular dynamics (i.e. at a finite distance from the blue line). As a consequence, the usual potential energy landscape ${\cal U}(\{r\})$ around the obtained configurations does not display excess soft anomalous modes at very low frequency, even near the jamming transition: these modes are gapped. This result must hold for the ground state too (which must be stable toward swap) and by continuity also for low-temperature equilibrated states. It may explain why marginal stability (and the Gardner transition leading to it) could not be observed in protocols where a thermal quench was used from swap-generated configurations Scalliet et al. (2017). It would be very interesting to see if other well-known excitations of low-temperature glassy solids are also gapped in these configurations, including two-level-systems, reported to be almost absent in experimental ultra-stable glasses Queen et al. (2013).

References

Appendix A SUPPLEMENTAL MATERIAL

This supplemental material (SM) provides: $(i)$ descriptions of the numerical model, protocols and methods used to generate athermal packings under swap dynamics at different pressures, $(ii)$ a computation of the Hessian of the potential energy, together with explanations about how pre-stress affects the vibrational modes, and $(iii)$ a discussion about the effect of the radii distribution generated by swap dynamics on the value of the packing fraction of our athermal packings obtained while freezing the degrees of freedom associated with particles’ radii.

We employ systems of $N\!=\!484$ particles in a square box in two dimensions. The total potential energy depends upon the particles coordinates $\{r\}$ and radii $\{R\}$, as

where $k$ is a stiffness, set to unity, $r_{ij}$ is the distance between the $i^{\mbox{\tiny th}}$ and $j^{\mbox{\tiny th}}$ particles, and $\Theta(x)$ is the Heaviside step function. The chemical potential associated with the radii is

where $\bar{k}_{R}$ is the stiffness of the potential associated with the radii $\{R\}$, that serves as a parameter in our study, and is set as described below. $R_{i}^{(0)}$ denotes the intrinsic radius of the $i^{\mbox{\tiny th}}$ particle. In each configuration we randomly assigned $R_{i}^{(0)}\!=\!0.5$ for half of the particles, and $R_{i}^{(0)}\!=\!0.7$ for the other half. The mass $m$ of particles, and that associated with their fluctuating radii, are all set to unity. Vibrational frequencies should be understood as expressed in terms of $\sqrt{k/m}$, and pressures in terms of $k$.

Configurations in mechanical equilibrium at zero temperature and at a desired target pressure $p_{0}$ were generated as follows; we start by initializing systems with random particle positions at packing fraction $\phi=1.2$, and set the initial radii to be $R_{i}\!=\!R_{i}^{(0)}$. We then minimize the total potential energy ${\cal V}(\{r\},\{R\})$ at a target dimensionless pressure $\Delta_{0}\!=\!10^{-1}$ using a combination of the FIRE algorithm Bitzek et al. (2006) and the Berendsen barostat Berendsen et al. (1984), see futher discussion about the latter below. Each packing is then used as the initial conditions for sequentially generating lower pressure packings, as demonstrated in Fig. S8. Following this protocol, we generated 1000 independent packings at each target dimensionless pressure, that ranges from $\Delta_{0}\!=\!10^{-1}$ up to $\Delta_{0}\!=\!10^{-5}$. For each target pressure, we set the stiffness $\bar{k}_{R}$ of the chemical potential of the radii according to $\bar{k}_{R}\!=\!p_{0}/\bar{\alpha}$, and vary $\bar{\alpha}$ systematically between $3\times 10^{-4}$ and $1$. During minimizations we calculate a characteristic net force scale F_{\mbox{\tiny typ}}\!\equiv\!\big{(}\sum_{i}||\vec{F}_{i}||^{2}/N\big{)}^{1/2}, where $\vec{F}_{i}\!=\!-\partial{\cal V}/\partial\vec{r}_{i}$ is the net force acting on the $i^{\mbox{\tiny th}}$ particle, whose coordinates are denoted by $\vec{r}_{i}$. A packing is considered to be in mechanical equilibrium when $F_{\mbox{\tiny typ}}$ drops below $10^{-8}\Delta_{0}$.

Berendsen barostat parameter: The FIRE algorithm Bitzek et al. (2006) features equations of motion which are to be integrated as in conventional MD simulations. We exploit this feature and embed the Berendsen barostat Berendsen et al. (1984) in our Verlet integration scheme Allen and Tildesley (1989). This amounts to scaling the simulation cell volume by a factor $\chi$, calculated as

where $\delta t$ is the (dynamical) integration time step, and $\xi$ is a parameter that determines how quickly the instantaneous dimensionless pressure converges to the target dimensionless pressure Allen and Tildesley (1989). In Fig. S9 shows the $\xi$-dependence of the convergence of the instantaneous dimensionless pressure $\Delta$ to the target value $\Delta_{0}$. Below $\xi=0.01$, the behavior of $\Delta$ as a function of iteration number is similar. We therefore set $\xi=0.01$ througthout this work.

The total potential energy ${\cal V}(\{r\},\{R\})$ of our model system is spelled out in Eqs. (S.14)-(S.16). We next work out the expansion of ${\cal V}$ in term of small displacements $\delta\vec{r}_{i}$ of particle positions, and small fluctuations $\delta R_{i}$ of the radii, about a mechanical equilibrium configuration with energy ${\cal V}_{0}$, as

where $H_{ij}\!\equiv\!\partial^{2}{\cal V}/\partial\vec{r}_{i}\partial\vec{r}_{j}$, $Q\!\equiv\!\partial^{2}{\cal V}/\partial R_{i}\partial R_{j}$, and $T_{ij}\!\equiv\!\partial^{2}{\cal V}/\partial\vec{r}_{i}\partial R_{j}$. The expansion given by Eq. (S.18) can be written using bra-ket notation as

The elements of the submatrix $H_{{\scriptscriptstyle\rm Nd,Nd}}$ can be written as tensors of rank $d\!=\!2$ as

where $\vec{n}_{ij}$ is a unit vector connecting between the $i^{\mbox{\tiny th}}$ and $j^{\mbox{\tiny th}}$ particles, $\vec{n}_{ij}^{\perp}$ is a unit vector perpendicular to $\vec{n}_{ij}$, $\otimes$ is the outer product, $\delta_{\langle ij\rangle}\!=\!1$ when particles $i,j$ are in contact, $\delta_{i,j}$ is the Kronecker delta, and the sum is taken over all particles $l$ in contact with particle $i$. The elements of the submatrix Q_{\mbox{\tinyN,\!N}} are scalars given by:

The matrix T_{\mbox{\tinyN,Nd}} is not diagonal and each element can be expressed as a vector with two components given by:

The eigenvectors of $H_{\mbox{\scriptsize swap}}$ are the normal modes of the system, and the eigenvalues are the vibrational frequencies squared $\omega^{2}$. The distribution of these frequencies is known as the density of states $D(\omega)$.

Effect of the pre-stress on the vibrational modes: When a system of purely repulsive particles is at mechanical equilibrium, forces $f_{ij}$ are exerted between particles in contact. These forces give rise to a term in the expansion of the energy, of the form

often referred to as the “pre-stress term”. For plane waves, it can be shown that the energy contributed by this term is very small. However, for the soft modes present when the system is close to the marginal stability limit, it can be shown that this term reduces the energy of the modes by a quantity proportional to the pressure Wyart et al. (2005). Marginal stability corresponds to a buckling transition where the destabilising effect of pre-stress exactly compensate the stabilising effect of being over-constrained. In this scenario where two effects compensate each other, the eigenvalue of the softest (non-Goldstone) modes of the Hessian in the absence of pre-stress $\bar{\omega}^{2}$ must be much larger than $\omega^{*}{}^{2}$ computed when pre-stress is present. To demonstrate this, we have calculated the density of states for systems while including and excluding the pre-stress term. The results are shown in Fig. S10, where it is found that near jamming $\omega^{*}{}^{2}/\bar{\omega}^{2}\approx 5\%$, which is consistent with what previously found for the traditional jamming transition DeGiuli et al. (2014a) and supports that the system is very close to (but not exactly at) marginal stability.

A.3 Packing fraction

In the main text we have shown that the jamming packing fraction $\phi_{c}$ generated using the swap dynamics increases when $\rho(R)$ broadens, i.e. for smaller values of the parameter $\bar{\alpha}$ that controls the stiffness of the potential energy associated with the radii. Here we compare the dependence of the packing fraction on pressure as measured for systems in which the radii are not allowed to fluctuate. In addition, in this test we borrow the distribution of radii $\rho(R)$ from swap-packings generated at $p\!=\!10^{-4}$, and at various values of the parameter $\bar{\alpha}$, varied between 1 to $\infty$ (the latter corresponds to disallowing particle radii fluctuations). Packings were generated using the total potential energy as given by Eq. S.15 (with the radii $R_{i}$ considered to be fixed), and using the same protocol and numerical methods used to generate the swap packings. The results are shown in Fig. S11, where it can be seen that the value of $\phi_{c}$ is essentially the same for any borrowed $\rho(R)$ from the swap packings.