Theory for Swap Acceleration near the Glass and Jamming Transitions for Continuously Polydisperse Particles

SWAP algorithms can shift the glass transition to lower temperatures, a recent unexplained observation constraining the nature of this phenomenon. Here we show that SWAP dynamics is governed by an effective potential describing both particle interactions as well as their ability to change size. Requiring its stability is more demanding than for the potential energy alone. This result implies that stable configurations appear at lower energies with SWAP dynamics, and thus at lower temperatures when the liquid is cooled. The magnitude of this effect is predicted to be proportional to the width of the radii distribution, and to decrease with compression for finite-range purely repulsive interaction potentials. We test these predictions numerically and discuss the implications of our findings for the glass transition. These results are extended to the case of hard spheres where SWAP is argued to destroy metastable states of the free energy coarse grained on vibrational timescales. Our analysis unravels the soft elastic modes responsible for the speed-up induced by SWAP , and allows us to predict the structure and the vibrational properties of glass configurations reachable with SWAP . In particular, for continuously polydisperse systems we predict the jamming transition to be dramatically altered, as we confirm numerically. A surprising practical outcome of our analysis is a new algorithm that generates ultrastable glasses by a simple descent in an appropriate effective potential.


I. INTRODUCTION
Understanding the mechanisms underlying the slowing down of the dynamics near the glass transition is a long-standing challenge in condensed-matter [1,2]. Unexpectedly, swap algorithms [3,4] (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 [5][6][7][8]. For judicious choice of poly-dispersity, one finds that: (i) the glass transition is shifted to lower temperatures: with swaps the α-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 swap 0 . (ii) The spatial extent of dynamical correlations, which are significant near T g , are greatly reduced with swap and only occur at T swap 0 . (iii) The mean square displacement of the particles on vibrational time scales is increased significantly in this temperature range [7]. 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 [9]. A theory of the glass transition should explain both swap and non-swap dynamics. Goldstein [10] 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 [11][12][13][14] 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 [9], as proven and confirmed numerically in a mean-field model of glasses [15]. 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 [16,17], to show that they are equivalent to a system of identical particles that can individually deform according to a chemical potential µ(R), where R is the particle radius. In the (practically important) case where poly-dispersity is continuous, µ(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: hyperstaticity is found with an excess number of contacts δz with respect of the Maxwell bound δz ∼ α 1/2 > 0, where α characterises the width of the radii distribution ρ(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 [18][19][20]. These results are verified numerically by introducing a novel algorithm performing a steepest descent in the generalised potential energy that includes µ(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 swap 0 < T < T g . We predict this shift to be proportional to α 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 [19,[21][22][23]. We use this approach to provide a simple phase diagram where the Goldstein transition and the emergence of marginality [18] (referred to as a Gardner transition in infinite dimension [23]) can be related to structure for both swap and non-swap dynamics.

A. Effective potential
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 [16,17]. We consider a system of N particles with continuous polydispersity ρ(R), of width α = ( R 2 − R 2 ) 1/2 / R . Here {R} indicates the set of particle radii and {r} their positions. In what follows we denote R ≡ R 0 , and 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 P({R}) of the particle radii. In the thermodynamic limit, a grandcanonical formulation is equivalent, in which particles of different radii correspond to different species. The associated partition function writes: (S.2) where µ(R) is the chemical potential at radius R. It is chosen such that in the thermodynamic limit, the distri-bution 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:

B. Mechanical stability under swap
Let us consider first an athermal system. In the thermodynamic limit, mechanical stability under swap dynamics requires 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 where p is the pressure and d the spatial dimension. To achieve a distribution of radii of width α, the stiffness k R acting on each particle radius must thus be of order: where the average is made on all particles i.

C. Generalised vibrational modes
Stability also requires the Hessian H swap (the matrix of second derivatives of V) to be positive definite. Since they are now N (d + 1) degrees of freedom, H swap is a N (d + 1) × N (d + 1) symmetric matrix, of eigenvalues ω 2 swap . It contains a block of size N d×N d which is the regular Hessian H ij = ∂ 2 U/∂r i ∂r j . We denote by ω 2 its eigenvalues. Because hybridisation with additional degrees of freedom can only lower the minimal eigenvalues of the Hessian, H 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 ≫ k, where k is the characteristic stiffness of the interaction potential U. In general, the eigenvalues of H are functions of the set of stiffnesses {k ij }, but also of the interaction forces {f ij } [24]. 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 ∼ ω 2 x 2 and change forces by a characteristic amount δf satisfying δf 2 /k ∼ x 2 ω 2 . Because of such change, Eq. (S.5) is not satisfied anymore. Thus the potential V can be reduced further by an amount of order δf 2 /k R ∼ ω 2 x 2 k/k R if the radii are allowed to adapt. This reduced energy can be approximatively written as x 2 ω 2 swap , where ω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 Θ(x) is the Heaviside step function.

A. Jamming transition for soft spheres under swap
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 [25,26], 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 δz with respect to the Maxwell threshold where the numbers of degrees of freedom and constraints match. Effective medium [27] or a variational argument [28] 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 α << ∆, where we introduced the dimensionless particle overlap ∆ ≡ pR d−2 0 /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 ω * swap 2 ∼ δz 2 k(1 − C 0 α/∆) where C 0 is a numerical constant. Pre-stress can be shown to shift eigenvalues of the Hessian by some amount ≈ −C 1 k∆ [18,29], leading to eigenvalues satisfying ω 0 swap 2 ∼ δz 2 k(1 − C 0 α/∆) − C 1 k∆. 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 α/∆. 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 δz which is above (but very close to in the limit α << ∆) the bound for non-swap dynamics of [18], recovered by setting α = 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 α >> ∆, 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 = ∞), whose number is N z/2. The dimension of the vector space satisfying such hard constraints is N (d+1)− N z/2 = N(1 −δ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 δz is small the number of degrees of freedom N (1 − δ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 prestress 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 ω i ∼ √ k R . When pre-stress is accounted for, eigenvalues of the Hessian are again shifted by ∼ −k∆. Mechanical stability then implies ∆/αδz 2 > C 2 ∆ and: In this regime, marginal stability (the saturation of the stability bound of Eq. S.12) corresponds to a coordination independent of pressure, with δz ∼ √ α and ω * swap 2 ∼ k∆. We thus predict that swap dynamics destroys isostacity, and significantly affect structure and vibrations. For sufficiently large α, 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 [30]) 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 [31], 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 → ∞, above which our predictions should apply.

B. 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: wherek R is a characteristic stiffness. We considered a potential diverging as R i → 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 This choice leads to a bimodal distribution of size ρ(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 [7], which would be similar to having a very large α 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 ∆ = 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 thatk R = p/ᾱ, wherē α is a parameter that controls the width α as shown in the inset of Fig. S1. For this bimodal distribution, α is defined as α = (α 1 + α 2 )/2, where α 1 , α 2 are the relative width of each peak in ρ(R). In the limit where the non-swap dynamics is recovered -which happens when α → 0 -α andᾱ are proportional.

C. 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 δz v.s. p for various values ofc. Isostaticity is indeed lost and the coordination converges to a plateau as ∆ decreases. Strikingly, we find for the plateau value δz ∼ √ α, 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 ∆ ∼ α. Overall, these results supports that the numerical curves δz(∆) in Fig. S2.A correspond to the marginal stability lines under swap dynamics, shown for different polydispersity (see more on that below).

D. Packing fraction
For traditional dynamics, polydispersity tends to have very limited effects on the value of jamming packing fraction φ J . We have confirmed this result in the S.I., by showing that although our model can generate very different distributions ρ(R), the values we obtain for φ 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 φ c ≡ lim ∆→0 φ(∆). The inset of Fig. S3.A confirms that φ c increases significantly as ρ(R) broadens. To quantify this effect we consider φ(∆,ᾱ), 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 ∆ ≫ᾱ is some scaling function and β = 1. Since the coordination does not change for ∆ ≪ᾱ, we expect that it is true for the structure overall and for φ, implying that f (x) ∼ x 0 as x → 0. These predictions are essentially confirmed in Fig. S3.B. Note however that the best scaling collapse is found for β = 0.83 < 1. These deviations are likely caused by finite size effects, known to be much stronger for φ than for the coordination or vibrational properties [32], and which may thus be present for our systems of N = 484 particles.

E. Vibrational properties
We computed the Hessian H swap and diagonalized it (see detailed in SI) to extract the density of states D(ω), as shown in Fig. S4.A for different pressures at fixed polydispersity. As expected, at low particle overlap ∆ two bands appear in the spectrum. The lowest-frequency band presents a plateau above some frequency scale ω * swap which satisfies ω * swap ∼ √ ∆ 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 [18]. Further evidence appears in Fig. S4.B showing D(ω) at fixed ∆ = 10 −4 for varying polydispersity. ω * swap essentially does not depend onᾱ as shown in the inset, as expected for marginal packings if the pressure is fixed. The cut-off frequency ω i of the low-frequency plateau scales as ω i ∼ √ k R ∼ 1/ √ᾱ , 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 [30]: for example for hard discs, φ g ≈ 0.78 and φ c ≈ 0.85. A similar difference of packing fraction occurs by compressing soft spheres at overlap ∆ ≈ 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 Uα (open symbols) as a function of the dimensionless pressure ∆ imposed during the quenches as a function ofᾱ as shown in the legend. The ratio U∞/Uα is shown in (c), it is stronger near jamming but remains significant even far for jamming is the system is sufficiently poly-disperse. As shown in the inset, in relative terms the shift of the energy of inherent structures (U∞ − Uα)/Uα is proportional to α and inversely proportional to ∆ when ∆ is large enough. (d) The packing fraction φc obtained after swap is turned on is larger than φ∞ obtained for nonswap dynamics, an effect that is stronger near jamming. reduced when swap is allowed for. We prove this point in Fig.S5.(a,b), where inherent structures of energy U ∞ 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 α < U ∞ . 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 ∆ = 0.05 for our α = 0.06. We show in the inset of that panel that the relative shift of energy induced by swap (U ∞ − U α )/U α is proportional to α and inversely proportional to ∆ when ∆ 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 ∆T g /T g induced by swap to be proportional to α, in consistence with the observation that very broad distri- butions lead to large swap effects [7]. We also predict that ∆T g /T g is inversely proportional to the distance to jamming ∆ when this quantity is well-defined (e.g. for soft spheres, but also to some extent for Lennard-Jones potentials [26,33]) 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 α = 0.06 and ∆ = 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 [26,33].

V. HARD SPHERE SYSTEMS
Our arguments above consider the energy landscape. For interactions potentials which are very sharp, nonlinearities 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 [21][22][23] 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 [34]. 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 [22]. (This procedure can be pursued to include thermal effects in soft spheres as well [35]).
Stability and vibrational properties can be computed in terms of this Hessian for non-swap dynamics [21,22]. Salient results are shown in the simplified diagram of Fig  S7. Once again, two key determinant of stability are the typical gap between interacting particles∆ ≡ T /(pR d 0 ) relative to the particle radius, and the excess coordination δz, where the coordination is defined from the network of particles that are colliding within a glassy state. A marginal stability line separates stable and unstable configurations, as illustrated in Fig.S7, whose asymptotic behaviour follows δz ∼∆ (2+2θ)/(6+2θ) where θ ≈ 0.41 [35]. (Strictly speaking, this line will depend slightly on the system preparation, but this dependence is expected to be modest, and is irrelevant for the present discussion).
Under a slow compression the system follows a line (in red) in the (∆, δz) plane. Mechanical stability is reached only for some∆ <∆ 0 , a characteristic onset gap where the dynamic crosses-over to an activated regime where vibrational modes become stable, in consistence with Goldstein's proposal. In these materials, deeper in the glass phase the system eventually returns to the stability line, and undergoes a sequence of buckling events that leave it marginally stable [18,19,22]. Marginal stability implies the presence of soft elastic modes (that differs from Goldstone modes) up to nearly zero frequency, and fixes the scaling properties of both structure and vibrations as jamming is approached [18,19]. These results, valid in finite dimensions, have been quantitatively confirmed in infinite dimension calculations [20,23,36]. In that case, the point where buckling sets in was argued to be a sharp transition, coined Gardner, where the free energy landscape fractures in a hierarchical way [36], as supported by numerical studies [37]. For very rapid quenches, it was argued that the entire glass phase should be marginal [22,36].
How is this picture affected by swap? Our arguments for the generalised Hessian of the soft sphere system essentially go through unchanged for the generalized Hessian of the free energy in the hard sphere system. Once again, stability becomes more demanding with swap, and the marginal stability line is shifted to higher coordination in the (∆, δz) plane as represented in the right panel of Fig. S7. Thus the glass transition is shifted FIG. S7: Log-log representation of the stability diagram in the coordination δz = z − zc and∆ plane for continuously poly-disperse thermal hard spheres with non-swap (Left) and swap (Right) dynamics. Note that for hard spheres∆ is independent of temperature and vanishes at jamming. In the Left panel, the blue line separates mechanically stable and unstable configurations. The red line indicates the trajectory of a system under a slow compression. When∆ decreases toward the onset gap∆0, metastable states appear and the dynamic becomes activated and spatially correlated. In the glass phase, the red trajectory will depend on the compression rate, but will eventually reach the blue line at some (ratedependent)∆G. When this occurs, a buckling or Gardner transition takes place where the material becomes marginally stable, leading to a power-law relation between∆ and δz. Right: for swap dynamics, stability is more demanding and is achieved only on the green line, which differs strictly from the blue one. Thus the onset gap decreases to some valuẽ ∆ s 0 <∆0: the dynamics become activated and correlated at larger densities, shifting the position of the glass transition. Marginality is still expected beyond some pressure∆ s G , but leads to plateau value for the coordination, indicating that isostaticity is lost.
toward higher packing fractions. At smaller gap∆ (corresponding to the approach of jamming), stability implies δz ≥ α (2+2θ)/(6+2θ) , again implying that isostaticity is lost with swap. We conjecture that, just as for soft particles, marginal stability is reached in the glass phase, which would correspond to a Gardner transition in infinite dimensions.

VI. CONCLUSION
In swap algorithms, the dynamics is governed by an effective potential 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 swap 0 < T 0 , leading to a natural explanation as to why the glass transition occurs then at a lower temperature T swap g < 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 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 [5][6][7][8]. Possible theories to describe the mechanism by which activation occurs in glasses include elastic [38] and facilitation models [39]. 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 α of the radii distribution ρ(R), we predict that T swap g (α) will smoothly decrease, while T g (α) should be essentially unchanged, with (T g (α) − T swap g (α))/T g (α) ∝ α and more specifically ∝ α/∆ for soft spheres. Furthermore, many studies have analyzed correlations between dynamics and vibrational modes, see e.g. [12,13,22,40], which can be repeated to relate the swap dynamics to the spectrum of the effective potential 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 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 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 [41]. 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 [42]. The pairwise potential term reads (S.2) where k is a stiffness, set to unity, r ij is the distance between the i th and j th particles, and Θ(x) is the Heaviside step function. The chemical potential associated with the radii is wherek 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 (0) i denotes the intrinsic radius of the i th particle. In each configuration we randomly assigned R 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 φ = 1.2, and set the initial radii to be R i = R (0) i . We then minimize the total potential energy V({r}, {R}) at a target dimensionless pressure ∆ 0 = 10 −1 using a combination of the FIRE algorithm [43] and the Berendsen barostat [44], 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 ∆ 0 = 10 −1 up to ∆ 0 = 10 −5 . For each target pressure, we set the stiffnessk R of the chemical potential of the radii according tok R = p 0 /ᾱ, and varȳ α systematically between 3 × 10 −4 and 1. During minimizations we calculate a characteristic net force scale F typ ≡ i || F i || 2 /N 1/2 , where F i = −∂V/∂ r i is the net force acting on the i th particle, whose coordinates are denoted by r i . A packing is considered to be in mechanical equilibrium when F typ drops below 10 −8 ∆ 0 . FIG. S8: Dimensionless pressure ∆ as a function iteration number for a packing-generating simulation starting from one particular initial condition. Each step of the staircase shape of the signal corresponds to the production of a packing at some desired target pressure. The criterion for convergence to mechanical equilibrium at each pressure is explained in the text. We produced packings ranging from ∆ = 10 −1 to ∆ = 10 −5 . Berendsen barostat parameter: The FIRE algorithm [43] features equations of motion which are to be integrated as in conventional MD simulations. We exploit this feature and embed the Berendsen barostat [44] in our Verlet integration scheme [45]. This amounts to scaling the simulation cell volume by a factor χ, calculated as where δt is the (dynamical) integration time step, and ξ is a parameter that determines how quickly the instantaneous dimensionless pressure converges to the tar-get dimensionless pressure [45]. In Fig. S9 shows the ξ-dependence of the convergence of the instantaneous dimensionless pressure ∆ to the target value ∆ 0 . Below ξ = 0.01, the behavior of ∆ as a function of iteration number is similar. We therefore set ξ = 0.01 througthout this work. Note that for these curves ω * 2 /ω 2 ≈ 5%, indicating that the system is very close to marginal stability.

Computation of Hswap
The total potential energy V({r}, {R}) of our model system is spelled out in Eqs. (S.1)-(S.3). We next work out the expansion of V in term of small displacements δ r i of particle positions, and small fluctuations δR i of the radii, about a mechanical equilibrium configuration with energy V 0 , as The expansion given by Eq. (S.5) can be written using bra-ket notation as where |δℓ is a (d+1)N -dimensional vector which concatenates the spatial displacements δ r i and the fluctuations of the radii δR i : (δ r 1 , δ r 2 , .., δ r N , δR 1 , ... δR N ). The operator H swap can be written as: The elements of the submatrix H Nd,Nd can be written as tensors of rank d = 2 as where n ij is a unit vector connecting between the i th and j th particles, n ⊥ ij is a unit vector perpendicular to n ij , ⊗ is the outer product, δ ij = 1 when particles i, j are in contact, δ 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 N,N are scalars given by: (S.7) The matrix T N, N d is not diagonal and each element can be expressed as a vector with two components given by: The eigenvectors of H swap are the normal modes of the system, and the eigenvalues are the vibrational frequencies squared ω 2 . The distribution of these frequencies is known as the density of states D(ω). 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 FIG. S11: Packing fraction φ as function of dimensionless pressure ∆ measured for packings in which radii are not allowed to fluctuate, and whose distribution of radii is borrowed from swap packings generated at ∆ = 10 −4 and at different values ofᾱ, as indicated by the legend. The inset shows a zoom into very small dimensionless pressures, demonstrating that the value of φc depends very weakly on the borrowed distribution of radii of the swap packings, as determined by the parameterᾱ.
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 [18].
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ω 2 must be much larger than ω * 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 ω * 2 /ω 2 ≈ 5%, which is consistent with what previously found for the traditional jamming transition [19] and supports that the system is very close to (but not exactly at) marginal stability.

Packing fraction
In the main text we have shown that the jamming packing fraction φ c generated using the swap dynamics increases when ρ(R) broadens, i.e. for smaller values of the parameterᾱ 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