# Irreversible swap algorithms for soft sphere glasses

Yoshihiko Nishikawa, Federico Ghimenti, Ludovic Berthier, Frédéric van Wijland

► To cite this version:

Yoshihiko Nishikawa, Federico Ghimenti, Ludovic Berthier, Frédéric van Wijland. Irreversible swap algorithms for soft sphere glasses. Physical Review E, 2025, 111 (4), pp.045416. 10.1103/Phys-RevE.111.045416. hal-05201182

HAL Id: hal-05201182

https://hal.science/hal-05201182v1

Submitted on 6 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Irreversible swap algorithms for soft sphere glasses

Yoshihiko Nishikawa,¹,² Federico Ghimenti,³,⁴ Ludovic Berthier,⁵ and Frédéric van Wijland³,⁶

¹Department of Physics, Kitasato University, Sagamihara, Kanagawa, 252-0373, Japan

²Graduate School of Information Sciences, Tohoku University, Sendai, Miyagi 980-8579, Japan

³Laboratoire Matière et Systèmes Complexes (MSC),
Université Paris Cité & CNRS (UMR 7057), 75013 Paris, France

⁴Department of Applied Physics, Stanford University, 348 Via Pueblo, Stanford, CA 94305, USA

⁵Gulliver, UMR CNRS 7083, ESPCI Paris, PSL Research University, 75005 Paris, France

⁶Yukawa Institute for Theoretical Physics, Kyoto University,
Kitashirakawaoiwake-cho, Sakyo-ku, Kyoto 606-8502, Japan

(Dated: January 20, 2025)

We extend to soft repulsive interaction potentials a recently proposed irreversible swap algorithm originally designed for polydisperse hard spheres. The original algorithm performs rejection-free, irreversible, collective swap moves. We show that event-driven cluster updates of particle diameters can also be performed in continuous potentials by introducing a factorised Metropolis probability. However, the Metropolis factorisation needed to deal with continuous potentials decreases the efficiency of the algorithm and mitigates the benefits of breaking detailed balance. This leads us to propose another irreversible swap algorithm using the standard Metropolis probability that accelerates the relaxation of soft sphere glasses at low temperatures, compared to the original swap algorithm. We apply these efficient swap algorithms to produce very stable inherent structures with vibrational density of states lacking the quasi-localised excitations observed in conventional glasses.

# I. INTRODUCTION

Supercooled liquids experience a dramatic slowing down of their dynamics upon a mild reduction in temperature or increase in density. A fundamental understanding of this phenomenon is the central goal of theories of the glass transition [1]. The slow dynamics of glass-formers poses serious practical challenges for their numerical investigation [2]. When approaching the glass transition, numerical simulations performed by integrating Newton's equations of motion experience the glassy slowing down at low temperature, and thermal equilibrium cannot be achieved within practically accessible timescales. Similarly to Newtonian dynamics, Monte Carlo simulations implementing local translational moves reproduce the characteristic dynamics of supercooled liquids on long timescales and experience slow dynamics as well [3]. However, Monte Carlo algorithms that sample a given target distribution need not be restricted to local translations, with irreducibility and aperiodicity being the only additional constraints the Markov chain must abide by [4–6]. Thus any ‘unphysical’ move, as far as it satisfies these two conditions without affecting the target long-time distribution, can be incorporated into a Monte Carlo algorithm. Devising efficient dynamical rules that reduce the dynamical slowing down extends the range of temperature and densities that can be accessed by numerical simulations, leading to a better understanding of the system properties. Beyond glasses, extended ensemble and cluster algorithms [7–14] are successful examples of how efficient Monte Carlo algorithms can bypass slow physical pathways, and can bring new insights into the thermodynamics of many complex systems.

In the physics of glasses and supercooled liquids, clever Monte Carlo algorithms have also been used to com-

pute the nontrivial static properties, including the Franz–Parisi potential and the point-to-set length [15–20]. The swap Monte Carlo algorithm [21], simply referred to as ‘Swap’ in the following, which allows for the exchange of particle diameters in polydisperse mixtures, has been particularly impactful on the numerical study of glasses in the last decade. When applied to size polydisperse systems [22–25], this relatively simple algorithm can accelerate the equilibration of certain models by several orders of magnitude, allowing the simulations to explore equilibrium states even below the experimental glass transition temperature. Starting from stable configurations obtained with the Swap algorithm, simulations employing the physical dynamics have revealed important correlations between the dynamic and static structures of ultrastable glasses [26, 27]. The Swap algorithm was numerically applied to other systems, such as crystals with complex structures and lattice models [19, 28–31]. However, the Swap itself eventually becomes slow at very low temperatures [32], and the quest for even more efficient Monte Carlo algorithms needs to be continued.

The Swap algorithm samples the Boltzmann distribution through reversible moves, i.e. that obey detailed balance: it randomly selects two particles and attempts an exchange of their diameters with the standard Metropolis probability. A line of research dating back to the nineties [33–37] points instead to the fact that giving up on the reversibility of a Markov chain while preserving its stationary distribution, can be rewarded with faster convergence. A most prominent family of algorithms exploiting this idea are the so-called lifted Markov chains [38, 39]. In a nutshell, an additional degree of freedom, the lifting variable, is introduced which drives the system away from equilibrium and breaks detailed balance. The dynamics is neverthe-
<figure><img src="image_2.png" /><figcaption>FIG. 1. Schematic description of various swap algorithms. Each panel shows a configuration (top) and the one-dimensional array, $s$, listing particle labels in order of increasing diameters. (a) In Swap, two particles are randomly selected along $s$ and a reversible exchange of their diameters is performed using a Metropolis acceptance rule. (b) In cSwap, a directed chain of swaps is built starting from an active particle (in bright red). Consecutive particles in $s$ are successively swapped resulting in a collective move. The activity label (red arrow) is updated either to the left (backward cSwap) or to the right (forward cSwap). (c) In kSwap, the consecutive swaps have a fixed size $k \ge 1$ (here $k = 2$), and the activity label is updated according to the forward prescription. Typically kSwap clusters contain less particles but span a larger interval in diameter space.</figcaption></figure>

less carefully designed to sample the Boltzmann distribution in its steady state, upon marginalising over the lifting degrees of freedom. This approach was successfully applied in a variety of physical settings [40–44]. Despite the recent progress in the mathematical understanding of lifted Markov chains [45–47], the design of an appropriate scheme for a given physical system remains a nontrivial task. It has recently been shown [48, 49] how the lifting approach can be employed to improve the efficiency of Swap for hard particles in two and three dimensions. The introduction of lifted degrees of freedom allows to perform driven, collective updates of the particle diameters. The resulting collective Swap (named ‘cSwap’) algorithm for hard sphere glasses achieves a significant speedup compared to Swap in the dense glassy regime. However, cSwap is tailored for the hard sphere potential and it cannot be applied for continuous potentials.

Here, we extend the algorithm to soft repulsive potentials, and study its efficiency. We design and compare algorithms that perform driven particle swaps in models of soft repulsive spheres in two and three dimensions. The algorithms studied in this work are illustrated in Fig. 1. In addition to the standard Swap in Fig. 1(a), we first discuss a variant of the cSwap algorithm for hard spheres, which differs from the original version in Ref. [48] in the manner the lifting degrees of freedom are updated. We then generalise these two cSwap algorithms for soft sphere continuous potentials, which exploit the factorised Metropolis probability [43] to perform directed collective

moves, see Fig. 1(b). Our simulations show that, contrary to hard spheres, the cSwap algorithms become less efficient than the reversible Swap algorithm for polydisperse soft sphere models. We prove that the probability factorisation is the cause of this somewhat disappointing result. This leads us to design another irreversible swap algorithm, which we name ‘kSwap’. This algorithm employs the full Metropolis probability and performs consecutive directed swap of a fixed jump size $k$, see Fig. 1(c). We find that kSwap is more efficient than Swap and cSwap for these continuous potentials. As an application, we use these new algorithms to produce inherent glass states with original physical properties.

The manuscript is organised as follows. We present the Swap and cSwap algorithms, and their extensions to continuous potentials in Sec. II. The models used to benchmark these algorithms are introduced in Sec. III. The efficiency of these algorithms is studied in Sec. IV. The new kSwap algorithm is introduced and studied in Sec. V. Inherent glassy states are produced and analysed in Sec. VI, while Sec. VII concludes the manuscript.

## II. COLLECTIVE SWAP ALGORITHMS

### A. Revisiting the hard sphere system

The collective swap (cSwap) algorithm [48, 49] performs driven cluster updates in the space of particle di-
ameters. This is done by extending the configuration space of the system through an extra degree of freedom. The system then evolves through a driven Markovian dynamics in the extended phase space. The dynamical rules are designed in such a way that the nonequilibrium stationary state is the Boltzmann distribution. This general construction can be rewarded with a faster convergence toward the target stationary distribution than equilibrium dynamics, and it is known in the literature as a lifted Markov chain [34, 35, 43–45, 50, 51].

In the collective swap algorithm the particle labels are sorted in a one-dimensional array by increasing diameter. An activity label, $v$, is initialised at random among the $N$ possible particle labels. This activity label $v$ constitutes the additional lifting degree of freedom. The active particle serves to initiate the collective swap move. An exchange of diameters between particle $v$ and the particle $i$ immediately to its right in $s$ is proposed, where $s$ is the one-dimensional array that lists particle labels in order of increasing diameters. This exchange is accepted with a probability given by the Metropolis rule. For hard spheres, the move is accepted with probability 1 if the exchange of diameters generates no overlaps, and rejected otherwise. If the move is accepted, the diameters of $v$ and $i$ are exchanged and the array is updated. A new attempt at updating the diameter of $v$ with the nearest right neighbour is again performed, following the same procedure. The accumulation of a number $n_c$ of accepted moves of this type is equivalent to an overall permutation where particle $v$ is maximally inflated, while all other particles involved in the move are simultaneously deflated and replace their left neighbour in the array $s$.

The inflation of the active particle $v$ stops when an exchange in the chain produces an overlap. When this happens, the lifting degree of freedom $v$ is updated instead of simply rejecting the exchange, thus terminating a single collective Swap move. In the original version of the algorithm [48, 49], the particle immediately to the left of the cluster becomes active for the next step, see Fig. 1(b), thus moving backward along the diameter array. We refer to this choice as the backward collective Swap. It is also possible to choose the particle immediately to the right of the updated cluster, see Fig. 1(b). In this case, the activity label moves in the same forward direction as the inflating particle along the sorted array. We christen this algorithm the forward collective Swap.

In a simulation, sets of cSwap moves alternate with sets of translational Metropolis Monte Carlo moves. A set of $N_{cSwap}$ cSwap moves is performed with probability $p_{cSwap}$, and a sweep made by $N$ translational moves is performed otherwise. Here $N_{cSwap}$ is a number comparable with the number of particles in the system, $N_{cSwap} \sim N$. Finally we observe that in both forward and backward versions of cSwap, the activity label $v$ is updated deterministically. One way to warrant ergodicity is to intersperse this deterministic evolution with a resetting process: After a cluster update, the activity label $v$ is resampled uniformly along the array $s$ with a

small probability $p_r$. The precise value of $p_r$ does not impact the efficiency of the algorithm, as long as $p_r \sim N^{-1}$. Here we take $p_r = \frac{1}{N_{cSwap}}$.

The computational cost of the cluster update can be noticeably reduced by resorting to an event-driven approach. Instead of attempting the pairwise swaps of the active particle $v$ one by one, the largest particle $i^*$ with which $v$ can be exchanged without generating overlaps is easily identified through a search operation, e.g. binary search [52], along $s$. The identification of particle $i^*$ defines the entire cluster to be updated in the collective Swap move, namely particles between $v$ and $i^*$. The cluster move is then performed by updating the diameter of $v$ to the one of $i^*$, and by updating the diameter of all the other particles in the cluster, $i^*$ included, to the diameter of the particle on their left along the sorted array. The resulting cascade of deflation can be safely performed, as the shrinking of a hard sphere never generates overlaps with its surrounding neighbours. This event-driven approach makes the computational cost of a collective update comparable with the one of a single equilibrium swap move in Swap, thus allowing cSwap to maintain an edge over equilibrium Swap in practical implementations.

## B. Generalisation to continuous potentials

The cSwap algorithm described above needs reconsideration when continuous potentials are used, as finite energy differences now need to be handled correctly. In soft spheres, Monte Carlo moves which produce particle overlaps involving finite positive increments in energy have a finite probability to occur. Given an assignment of the particle diameters, the probability to accept an exchange of the diameters of two particles $i$ and $j$ in a system at inverse temperature $\beta$, is given by the Metropolis probability

$$P_{\text{Met}} = \min[1, e^{-\beta \Delta E(i,j|s)}] = e^{-\beta \max[0, \Delta E(i,j|s)]}, \quad (1)$$

where $\Delta E(i, j|s)$ is the change of energy of the system via the proposed swap move. Cluster moves can of course be constructed out of this Metropolis acceptance probability, but the evaluation of Eq. (1) would become numerically very expensive when the size of the cluster becomes large as $\Delta E$ would now involve a large number of energy evaluations. As observed in Ref. [44], this problem can be alleviated by replacing the acceptance probability in Eq. (1) with a factorised Metropolis probability $P_{fMet}$, defined as

$$P_{\text{fMet}}(i, j|s) = e^{-\beta \max[0, \Delta E_{ij}(s)]} \times e^{-\beta \max[0, \Delta E_{ji}(s)]}, \quad (2)$$

where $\Delta E_{ij}(s)$ is the energy difference of particle $i$, and of particle $i$ only, before and after the exchange of its diameter with particle $j$ on the sorted array $s$. The total energy change in the swap move is thus $\Delta E(i, j|s) = \Delta E_{ij}(s) + \Delta E_{ji}(s)$. The key observation is that $P_{fMet}$ also obeys detailed balance and can be used to sample
the Boltzmann distribution. For repulsive soft interactions, which we use below, the shrinking particle has a negative energy change while the inflating one has a positive change. If the diameter of particle $i$ inflates, the factorised Metropolis probability in Eq. (2) reduces to

$$P_{\text{fMet}}(i, j|s) = e^{-\beta \Delta E_{ij}(s)}. \quad (3)$$

Using the factorised Metropolis probability, we propose a generalisation to soft repulsive potential of the forward cSwap algorithm that is amenable to an event-driven formulation. As in the hard sphere case, the activity label is initialised at random along the sorted array $s$. Particle $v$ attempts to exchanging its diameter with the particle on its right along $s$. We denote this particle by $v_1$. The exchange is accepted according to the factorised Metropolis probability $P_{\text{fMet}}(v, v_1|s)$. If accepted, the diameters and the array $s$ are updated, after which $v$ attempts a new exchange. This procedure is iterated until the first rejection. The activity label is then updated according to the forward cSwap rule. As for hard spheres, a uniform resampling of the activity label is implemented to ensure aperiodicity. This procedure defines the forward cSwap algorithm for a soft repulsive potential.

The steady state of the forward cSwap algorithm is given by the Boltzmann distribution, since the latter verifies the stationarity, or global balance, condition

$$\pi_{ss}(s|v) = \sum_{s', v'} \phi((s', v') \rightarrow (s, v)). \quad (4)$$

Here $\pi_{ss}(s|v)$ is the stationary probability distribution for an instance of ordered diameters $s$ and of the lifted label $v$, defined as

$$\pi_{ss}(s|v) = \frac{1}{N}\pi_{B}(s), \quad (5)$$

with $\pi_B(s) \propto e^{-\beta E(s)}$ the Boltzmann distribution. The quantity $\phi((s', v') \rightarrow (s, v))$ is the steady state probability flux from a configuration $(s', v')$ to the configuration $(s, v)$. The forward cSwap rules imply that

$$\begin{aligned} \sum_{s', v'} \phi((s', v') \rightarrow (s, v)) &= \pi_{ss}(s^*|v)e^{-\beta \Delta E_{v, v_1}(s^*)} \\ &\quad + \pi_{ss}(s|v_{-1})(1 - e^{-\beta \Delta E_{v_{-1}, v}(s)}) \\ &= \frac{1}{N}\pi_B(s) = \pi_{ss}(s|v). \end{aligned} \quad (6)$$

The configuration $(s^*, v)$ is the configuration that flows into $(s, v)$ upon a successful exchange of particle $v$ with its right neighbour along $s$. The configuration $(s, v_{-1})$ flows into $(s, v)$ after a rejection of the exchange of the particle diameters $v_{-1}$ with $v$. Here, and in what follows, $v_n$ denotes the particle $n$ steps away from $v$ along the ordered array $s$ (if $n$ is positive $v_n$ lies to the right of $v$, if $n$ is negative $v_n$ lies to the left of $v$). The second

equality follows from the definition of the steady state and Boltzmann distributions in Eq. (5) and the fact that

$$\pi_B(s^*)e^{-\beta \Delta E_{v, v_1}(s^*)} = \pi_B(s)e^{-\beta \Delta E_{v_{-1}, v}(s)}. \quad (7)$$

By plugging Eq. (6) into Eq. (4) we confirm that the stationarity condition is indeed satisfied.

The use of the factorised Metropolis filter allows us to to easily calculate the probability to accept $n$ successive swaps along the array $s$ as

$$\prod_{i=v}^{v_n} P_{\text{fMet}}(i, i_1|s) = e^{-\beta \Delta E_{v, v_n}(s)}. \quad (8)$$

The simplification of the Boltzmann factors in this product is the main rationale behind the introduction of the factorised Metropolis filter. As for hard particles, we can find $i^* = v_n$ efficiently in an event-driven way. To this end, we draw a uniform random number $\Upsilon \in (0, 1)$, and find $i^*$ such that

$$e^{-\beta \Delta E_{v, i^*} (s)} < \Upsilon < e^{-\beta \Delta E_{v, i^*} (s)}. \quad (9)$$

Finding $i^*$ can be efficiently done with computational complexity $O(\log N)$ through algorithms such as binary search. A cluster move of all the particles between $v$ and $i^*$ on the array $s$ can then be performed.

An event-driven backward cSwap algorithm for soft potentials is constructed analogously, by changing the update rule for the lifting degree of freedom. A proof of stationarity for the event-driven backward cSwap is obtained as follows. Given a configuration $(s|v)$, we denote by $m$ the integer such that $v_m$ is the rightmost label along the ordered array $s$. In the backward cSwap algorithm, there are $m$ probability flows into state $(s|v)$. The right hand side of the stationarity condition, Eq. (4), becomes

$$\pi_{ss}(s|v) + \sum_{n=1}^{m-1} [\phi_{\text{in}}(n+1) - \phi_{\text{out}}(n)], \quad (10)$$

where

$$\phi_{\text{out}}(n) = \pi_{ss}(s_n|v_n)e^{-\beta \Delta E_{v_1, v_{n+1}}(s_n)}, \quad (11)$$

$$\phi_{\text{in}}(n) = \pi_{ss}(s_n|v_n)e^{-\beta \Delta E_{v_1, v_n}(s_n)}, \quad (12)$$

and $s_n$ is the configuration that flows into $s$ after performing an event-driven backward Swap move from the active particle $v_1$. Since

$$E(s_n) + \Delta E_{v_1, v_{n+1}}(s_n) = E(s_{n+1}) + \Delta E_{v_1, v_{n+1}}(s_n), \quad (13)$$

we have $\phi_{\text{in}}(n+1) - \phi_{\text{out}}(n) = 0$ for any $n$, and the total probability flow is $\pi(s|v_1) = \pi(s|v)$, which completes the proof.

Note that the backward algorithm balances the probability flows of multiple cluster moves, not of single exchanges. Thus, interrupting a cluster update violates the balance condition, whereas it does not cause any problem in the forward algorithm, as the latter satisfies the balance condition through every single diameter exchange during the construction of the cluster.
### III. GLASS-FORMING MODELS USED FOR BENCHMARKING

To test these new cSwap algorithms, we explore their efficiency in the context of two models of glass-formers. We use polydisperse soft spheres in two [53] and three dimensions [23] with number density $\phi = 1$, interacting through the following potential:

$$V(r_{ij})/v_0 = \left(\frac{d_{ij}}{r_{ij}}\right)^{12} + c_0 + c_2 \left(\frac{r_{ij}}{d_{ij}}\right)^2 + c_4 \left(\frac{r_{ij}}{d_{ij}}\right)^4, \quad (14)$$

if $r_{ij} < r_{ij}^{\text{cutoff}}$ with $r_{ij}^{\text{cutoff}} = 1.25d_{ij}$ and $v_0 > 0$ an energy scale, and $V(r_{ij}) = 0$ otherwise. Here, $r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|$ with $\mathbf{r}_i$ the position of particle $i$. The particle diameter is non-additive, $d_{ij} = \frac{d_i+d_j}{2} (1 - \epsilon |d_i - d_j|)$ with $\epsilon = 0.2$, which prevents crystallisation and demixing. We use the following power-law distribution for the diameters:

$$P(d) = \begin{cases} Ad^{-3} & : d_{\min} < d < d_{\max}, \\ 0 & : \text{otherwise,} \end{cases} \quad (15)$$

where $d_{\min}/d_{\max} = 0.45$. The constant factor $A$ is a normalisation constant. This distribution leads to an overall polydispersity $\delta = (d^2 - \bar{d}^2)^{1/2}/\bar{d} = 0.23$. The units of length and temperature are $\bar{d}$ and $v_0$.

The particle positions are updated with the standard Metropolis translational dynamics, where we first displace a particle by a random vector $\mathbf{\delta} \in (-\delta_{\max}, \delta_{\max})^d$ with $d$ the spatial dimension, and then accept it with the Metropolis probability. We set $\delta_{\max} = 0.175\bar{d}$ and $\delta_{\max} = 0.1\bar{d}$ for the two- and three-dimensional models, respectively. For the swap algorithms, we perform a set of diameter swaps with probability $p_{\text{swap}} = 0.2$ and otherwise perform a set of $N$ translational moves. The value of $p_{\text{swap}}$ is approximately optimised to minimise the relaxation time at low temperatures. When using the event-driven schemes for the collective algorithms, $N_{\text{cSwap}} = 512$ clusters are updated per unit time, while for the swap algorithms using the full Metropolis probability, one unit time consists of $N$ swap attempts. We verified that this choice of time units faithfully mirrors, up to an overall rescaling, the CPU time required by each algorithm to run. The system size is $N = 1024$ throughout the paper, unless explicitly mentioned otherwise.

For the two-dimensional model, we measure the time correlation of the local orientational order [54] at time $t$,

$$C_6(t) = \frac{\left\langle \sum_j \psi_j(0) \psi_j^*(t) \right\rangle}{\left\langle \sum_j |\psi_j|^2 \right\rangle}, \quad (16)$$

$$\psi_j(t) = \frac{1}{|\partial j(t)|} \sum_{k \in \partial j(t)} e^{-6i\theta_{jk}(t)}, \quad (17)$$

where $\theta_{jk}(t)$ is the angle of $\mathbf{r}_k(t) - \mathbf{r}_j(t)$ with respect to the $x$ axis, and $\langle \cdot \rangle$ represents an average over initial equilibrium configurations. The set of neighbouring particles

<figure><img src="image_3.png" /><figcaption>FIG. 2. Time correlations of the two dimensional (a) and three-dimensional (b) models at $T = 0.08$ comparing Metropolis translations, with Swap, forward and backward cSwap algorithms. In both cases, Swap remains marginally faster than cSwap algorithms.</figcaption></figure>

$\partial j(t)$ for particle $j$ is found using the Voronoi tessellation for each configuration at time $t$. On the other hand, the relaxation of the three-dimensional model is studied using the overlap function

$$Q(t) = \left\langle \frac{1}{N} \sum_i w(|\mathbf{r}_i(0) - \mathbf{r}_i(t)|) \right\rangle, \quad (18)$$

where $w(r)$ is 0 when $r > 0.2\bar{d}$ and 1 otherwise.

In all cases, we define the relaxation time $\tau_\alpha$ as the time where the time correlation functions decrease to the value $e^{-1}$. While this is a measurement of a specific correlation time, its temperature dependence faithfully reflects the evolution of the equilibration timescale that characterises either the approach to stationarity or the ergodic exploration of the configuration space at thermal equilibrium.

### IV. EQUILIBRATION TIMES OF COLLECTIVE SWAP ALGORITHMS

We test the two cSwap algorithms on the models at $T = 0.08$, much lower than the onset temperatures $T_o \approx 0.3$ in two dimensions and $T_o \approx 0.2$ in three dimensions, below which the heterogeneous glassy dynamics sets in [55, 56]. The data in Fig. 2 show the time correlations for the simple Metropolis, Swap and both cSwap algorithms. At this temperature, the Metropolis translational dynamics is completely frozen and the corresponding time correlation shows a plateau covering the longest simulated timescale, showing that relaxation is impossible within our simulation window using simple Monte Carlo simulations. On the other hand, all three
<figure><img src="image_4.png" /><figcaption>FIG. 3. Swap probability $P(d_1, d_2)$ between particles with diameters $d_1$ and $d_2$ at equilibrium in the three-dimensional model at $T = 0.08$, using (a) the full Metropolis expression in Eq. (1) and (b) the factorised Metropolis probability, Eq. (2), respectively. Off-diagonal elements corresponding to large diameter exchanges are strongly suppressed by the factorised probability in (b).</figcaption></figure>

swap algorithms can decorrelate the particle positions or relative orientations, and all time correlations decay to a small value after $t \approx 10^5$. Strikingly, however, the two cSwap algorithms are both slightly slower than the simple Swap despite the large number of clusters updated per unit time. This result is in stark contrast with the hard sphere situation [48, 49], where cSwap is roughly 10 times faster than Swap in both spatial dimensions.

The main difference between the two families of algorithms is the use of the factorised Metropolis probability for the continuous potentials, in order to satisfy the stationarity condition. We may therefore expect the loss of efficiency of cSwap for soft spheres to be caused by the factorisation. For hard spheres, both expressions for the Metropolis probabilities become identical, whereas for soft spheres we generally have $P_{\text{Met}} < P_{\text{Met}}$.

We test this hypothesis in Fig. 3, which compares full and factorised swap acceptance probabilities. It can be seen that the factorisation significantly reduces the acceptance of swap moves at low temperatures. The full Metropolis probability can be close to 1 even when the two diameters are dissimilar, as the negative and positive energy changes can cancel each other. By contrast, the factorised probability only uses the positive energy change of the inflating particle and rapidly decreases with the diameter difference.

Since the cSwap algorithms gradually increase the diameters of particles by successive moves between neighbouring particles in the $s$ array, the factorised probability results in clusters with a relatively small number of particles ($\approx 20$ at the lowest temperature), leading to slow decorrelation in diameter space. In fact, the factorised Metropolis probability also slows down the event-chain algorithm for translations for continuous potentials compared with hard particles [43, 50]. We have explicitly checked the influence of the factorised probability in our algorithm by comparing the equilibrium Swap against a version of forward cSwap that uses the full Metropolis

probability, with the same number of cluster updates per unit time as in the factorised case. This unfactorised algorithm is faster than Swap, thus demonstrating that the factorised probability is indeed responsible for the efficiency loss. However, a forward cSwap with a large number of cluster updates using the full Metropolis acceptance probability is not practically useful, because of the large overhead needed to estimate $P_{\text{Met}}$.

V. A NEW IRREVERSIBLE SWAP ALGORITHM WITH LARGE JUMPS

The generalisation of cSwap to continuous potentials involves a factorisation of the Metropolis acceptance rule. When applied to large clusters, the factorisation decreases the acceptance so much that the algorithm becomes inefficient. Yet the efficiency of cSwap over Swap relies on the rapid inflation of the active particle to the maximum value allowed by the Boltzmann distribution. Large changes of the particle are indeed the source of the speedup offered by swap moves [23].

These arguments and the search for an efficient algorithm suggest to introduce collective and directed swap moves that use the Metropolis probability and allow for large moves in diameter space. To this end, we generalise the forward cSwap so that it now incorporates larger jumps, of size $k$ along the $s$ array, see Fig. 1(c). We christen this algorithm ‘kSwap’. kSwap is a modification of cSwap where an active particle $v$ is swapped with another particle $v_k$, $k \ge 1$, further away along the one-dimensional array $s$ with the full Metropolis probability (recall that cSwap swaps adjacent particles $v$ and $v_1$ with the factorised Metropolis probability). The jump size $k$ is sampled from a uniform distribution between 1 and $k_{\text{max}}$ every time we start the swap update, while it is kept fixed until the update stops. The parameter $k_{\text{max}}$ thus controls the typical jump size. We incorporate these large jumps into the forward algorithm. Since it uses the full Metropolis probability, the unit time of kSwap consists of $N$ swap attempts. Its computational complexity is thus equivalent to that of Swap.

Large jumps in diameter space typically yield a smaller acceptance probability, and kSwap is not as collective as the cSwap algorithms, with the typical cluster size $\langle n_c \rangle \approx 3$ when $k_{\text{max}} = 100$, while $\langle n_c \rangle \approx 20$ for the cSwap algorithms. Nevertheless, kSwap is faster than cSwap and even faster than the original Swap. In Fig. 4, we show the relaxation time $\tau_\alpha$ of the two models as a function of inverse temperature $1/T$, for all algorithms. All algorithms incorporating swap moves are much faster than the original Metropolis Monte Carlo. For soft spheres, kSwap is faster than Swap, which is itself slightly faster than both forward and backward cSwap. The speedup of kSwap with respect to Swap weakly depends on temperature and increases towards low temperatures.

The relaxation time of kSwap depends on $k_{\text{max}}$. For $k_{\text{max}} = 1$, the algorithm reduces to a forward cSwap with
<figure><img src="image_5.png" /><figcaption>FIG. 4. Relaxation time $\tau_{\alpha}$ as a function of inverse temperature for (a) the two- and (b) three-dimensional systems. The backward and forward cSwap algorithms are slower than Swap, whereas kSwap with $k_{max} = 100$ is faster.</figcaption></figure>

the full Metropolis probability, and is slower than Swap. We find, for $N = 1024$ and $k_{max} \in [50, 150]$, that kSwap is faster than Swap, that the value $k_{max} = 100$ is roughly optimal and that it yields the shortest relaxation time at low temperatures. For continuously polydisperse models, the typical interval between two neighbouring diameters in the array decays as $N^{-1}$, thus we expect $k_{max}$ to linearly scale with $N$ for optimal performance. We have checked $k_{max} = 50$ indeed yields almost the same relaxation time for kSwap when $N = 512$, supporting the above scaling of $k_{max}$ with $N$.

## VI. USING SWAP ALGORITHMS TO PRODUCE NON-CONVENTIONAL STABLE GLASSES

Algorithms such as gradient descent and FIRE (for Fast Inertial Relaxation Engine) [57] starting from an equilibrium configuration at temperature $T$ produce an energy minimum, also called inherent structure (IS). Inherent structures and the related concept of a potential energy landscape are frequently used to discuss the properties of supercooled liquids and glasses [58–68]. Two quantities are of particular interest. First, the value of the energy density in the inherent state reached from a given $T$ reveals the stability of the corresponding amorphous structures. Typically, an IS has a lower energy when it is reached from a lower $T$, indicating more stable glassy structures lying deeper in the potential energy landscape. The second quantity is the vibrational density of states (DOS), which describes the range of eigen-

<figure><img src="image_6.png" /><figcaption>FIG. 5. Inherent structure energy $e_{IS}/N$ as a function of initial temperature $T_{init}$ for (a) the two- and (b) the three-dimensional systems, respectively. The dashed line in each panel represents $e_{IS}/N$ at the lowest $T_{init}$ for the Metropolis quench.</figcaption></figure>

frequencies characterising eigenmodes of the vibrational matrix at the energy minimum. Roughly speaking, an abundance of low frequency modes would correspond to an energy minimum with many soft vibrational modes. Remarkably, the vibrational density of states $D(\omega)$ of inherent glassy states displays a power law scaling at small frequencies, $D(\omega) \sim \omega^4$ [69–72], in addition to the Debye law. The stability of the glass is usually encoded into the prefactor of this power law, which decreases sharply with increasing stability [72]. The $\omega^4$ scaling of the DOS is thus a universal property of amorphous solids, which suggests an abundance of low-frequency soft modes resulting from the disordered structure of the IS.

We generate inherent structures from equilibrium configurations at temperature $T_{init}$ using zero-temperature Monte Carlo algorithms, in which only updates lowering the total energy are accepted. For the Metropolis translational dynamics, we reduce the jump size to $\delta_{max} = 0.0005\bar{d}$. We refer to this as the Metropolis quench. We also use hybrid Monte Carlo quenches using both translational and swap moves, using either Swap or kSwap. Note that cSwap algorithms do not work for quenches, as the factorised probability is always zero for collective swaps at zero temperature. The number of Monte Carlo sweeps per particle in our simulations ranges from $10^6$ to $2 \times 10^6$. The number of negative eigenvalues of the Hessian matrix for a generated configuration after the quench is typically zero, indicating that the Monte Carlo approach indeed produces energy minima.

In Fig. 5, we show the IS energy density $e_{IS}/N$ as a function of initial temperature $T_{init}$ for the Metropolis,
Swap, and $k_{max}$ = 100 kSwap algorithms. The Metropolis quench reproduces the behaviour observed previously in glass-formers, with a very weak dependence of the IS energy at high $T_{init}$, and a sharper dependence below a crossover temperature, corresponding to the onset temperature [55]. Below the onset, IS energies depend more sensitively on the initial temperature, lower $T_{init}$ ending deeper in the potential energy landscape. The crossover correspond roughly to the onset of slow dynamics, near $T_{init} \approx 0.3$ in two dimensions and $T_{init} \approx 0.2$ in three dimensions [26].

A clear observation in Fig. 5 is that quenching from the same configuration using a zero-temperature algorithm that performs both translational and swap moves yields much lower energies than Metropolis quenches. The performance of the two swap algorithms are nearly equivalent. The temperature dependence of the IS energies is less pronounced than for conventional quenches, especially at high temperatures. and it is essentially the same for the two swap algorithm. There is no obvious qualitative difference between two and three dimensions. In both cases, we observe that instantaneous quenches starting from very high temperatures yield inherent states that are extremely deep in the potential energy landscape, and in fact much deeper than what could be obtained without swap algorithms.

These findings for continuous potentials confirm the results obtained for fast compressions of hard particles with the swap Monte Carlo algorithm [22, 48, 49], and rationalise the efficiency of approaches where gradient descent involving not just translational degrees of freedom are used to produce stable glassy structures as developed in [73–76]. Qualitatively, using swap moves during annealing seems much more efficient than very slowly cooling the system without swap moves. We note in Fig. 5(a) a qualitative difference between hard and soft potentials. For two dimensional hard disks, fast compressions with swap moves seemingly all terminate at the same density, independently of the initial state [48]. This observation led Bolton-Lum et al. to conclude that this end-point represents an ‘ideal’ disk packing [76]. For continuous potentials, we do not find evidence for such an ‘ideal’ amorphous structure (or ground state) that would be systematically found in the potential energy landscape.

Using swap moves to prepare inherent structures not only produces much lower energies, but also affects the physical properties of the energy minima that are found. This is revealed by the vibrational DOS. In Fig. 6, we show the DOS for two and three dimensional glasses and compare in each case inherent structures that have the exact same energy (dashed lines in Fig. 5) but are produced either with or without swap moves. Clearly, the DOSs are very different for the two quench protocols, showing that the IS energy does not uniquely characterise amorphous structures. For the two-dimensional system, Fig. 6(a), the Metropolis quench has several isolated modes extending towards low frequencies with an envelope compatible with the $\omega^4$ power law reported be-

<figure><img src="image_7.png" /><figcaption>FIG. 6. Vibrational density of states of inherent structures for (a) two and (b) three dimensions. Temperatures are chosen so that IS energies are the same in the two algorithms, see dashed lines in Fig. 5. The swap moves deplete the soft modes at low frequencies.</figcaption></figure>

fore. These few low-frequency modes are nearly removed by the swap quenches. This effect becomes more obvious in three dimensions, Fig. 6(b), where the smooth tail found for Metropolis quenches is replaced by a gapped distribution. This again shows that swap moves during the approach to the energy minimum remove the soft vibrational modes and change the nature of the energy minima that are reached. The observation of gapped vibrational spectra was also reported using augmented gradient descent techniques [74, 75].

VII. CONCLUSION

This paper extends to soft continuous potentials our recent efforts [48, 49] to develop irreversible versions of the swap Monte Carlo algorithm to speedup even further the equilibration of glass-forming models. The broader rationale is the belief that any algorithm developed for hard spheres should have an efficient counterpart for continuous potentials, as illustrated before by the event-chain and geometric cluster algorithms, and even the event-driven molecular dynamics [6, 11, 43, 50, 77–81].

Here we successfully generalised the irreversible cSwap algorithms developed for hard spheres to soft potentials, following the path used for event-chain translational moves. However, the performance of the resulting algorithms applied to soft sphere models are not as good as expected. We have identified the factorisation of the Metropolis probability as the root cause of the degradation of the algorithmic performance, which led us to
devise an irreversible swap algorithm that uses the full Metropolis probability and borrows from both cSwap and Swap. This new algorithm accelerates equilibrium relaxation at low temperature compared to Swap and cSwap, despite the fact that it is less collective than cSwap. Our results provide instructive insights into the nature of collective moves needed in Monte Carlo algorithms for dense glassy systems, and the edge that can be achieved by tuning the range of the irreversible moves. Future research directions should aim at improving further Monte Carlo algorithms, for instance by exploiting, among others, reinforcement learning [82] and normalising flows strategies [83].

## ACKNOWLEDGMENTS

We thank Manon Michel and Ernest van Wijland for useful discussions. Y.N. acknowledges support from JSPS KAKENHI (Grant No. 22K13968). F.G., L.B. and F.v.W. acknowledge the financial support of the ANR THEMA AAPG2020 grant.

[1] L. Berthier and G. Biroli, Theoretical perspective on the glass transition and amorphous materials, Rev. Mod. Phys. 83, 587 (2011).

[2] L. Berthier and D. R. Reichman, Modern computational studies of the glass transition, Nat. Rev. Phys. 5, 102 (2023).

[3] L. Berthier and W. Kob, The Monte Carlo dynamics of a binary Lennard-Jones glass-forming mixture, J. Phys. Condens. Matter 19, 205130 (2007).

[4] D. Levin and Y. Peres, *Markov Chains and Mixing Times* (American Mathematical Society, Providence, Rhode Island, 2017).

[5] M. E. Newman and G. T. Barkema, *Monte Carlo methods in statistical physics* (Clarendon Press, 1999).

[6] W. Krauth, *Statistical mechanics: algorithms and computations*, Vol. 13 (OUP Oxford, 2006).

[7] R. H. Swendsen and J.-S. Wang, Nonuniversal Critical Dynamics in Monte Carlo Simulations, Phys. Rev. Lett. 58, 86 (1987).

[8] U. Wolff, Collective Monte Carlo Updating for Spin Systems, Phys. Rev. Lett. 62, 361 (1989).

[9] B. A. Berg and T. Neuhaus, Multicanonical ensemble: A new approach to simulate first-order phase transitions, Phys. Rev. Lett. 68, 9 (1992).

[10] E. Marinari and G. Parisi, Simulated Tempering: A New Monte Carlo Scheme, EPL 19, 451 (1992).

[11] C. Dress and W. Krauth, Cluster algorithm for hard spheres and related systems, J. Phys. A: Math. Gen. 28, L597 (1995).

[12] K. Hukushima and K. Nemoto, Exchange Monte Carlo Method and Application to Spin Glass Simulations, J. Phys. Soc. Jpn. 65, 1604 (1996).

[13] L. Santen and W. Krauth, Absence of thermodynamic phase transition in a model glass former, Nature 405, 550 (2000).

[14] F. Wang and D. P. Landau, Efficient, Multiple-Range Random Walk Algorithm to Calculate the Density of States, Phys. Rev. Lett. 86, 2050 (2001).

[15] Q. Yan, T. Jain, and J. J. de Pablo, Density-of-States Monte Carlo Simulation of a Binary Glass, Phys. Rev. Lett. 92, 235701 (2004).

[16] L. Berthier, Overlap fluctuations in glass-forming liquids, Phys. Rev. E 88, 022313 (2013).

[17] R. L. Jack and J. P. Garrahan, Phase Transition for Quenched Coupled Replicas in a Plaquette Spin Model of Glasses, Phys. Rev. Lett. 116, 055702 (2016).

[18] L. Berthier, P. Charbonneau, and S. Yaida, Efficient measurement of point-to-set correlations and overlap fluctuations in glass-forming liquids, J. Chem. Phys. 144, 024501 (2016).

[19] Y. Nishikawa and K. Hukushima, Lattice Glass Model in Three Spatial Dimensions, Phys. Rev. Lett. 125, 065501 (2020).

[20] B. Guiselin, L. Berthier, and G. Tarjus, Random-field Ising model criticality in a glass-forming liquid, Phys. Rev. E 102, 042129 (2020).

[21] T. S. Grigera and G. Parisi, Fast Monte Carlo algorithm for supercooled soft spheres, Phys. Rev. E 63, 045102 (2001).

[22] L. Berthier, D. Coslovich, A. Ninarello, and M. Ozawa, Equilibrium Sampling of Hard Spheres up to the Jamming Density and Beyond, Phys. Rev. Lett. 116, 238002 (2016).

[23] A. Ninarello, L. Berthier, and D. Coslovich, Models and Algorithms for the Next Generation of Glass Transition Studies, Phys. Rev. X 7, 021039 (2017).

[24] L. Berthier, E. Flenner, C. J. Fullerton, C. Scalliet, and M. Singh, Efficient swap algorithms for molecular dynamics simulations of equilibrium supercooled liquids, J. Stat. Mech. 2019, 064004 (2019).

[25] A. D. S. Parmar, M. Ozawa, and L. Berthier, Ultrastable Metallic Glasses In Silico, Phys. Rev. Lett. 125, 085505 (2020).

[26] B. Guiselin, C. Scalliet, and L. Berthier, Microscopic origin of excess wings in relaxation spectra of supercooled liquids, Nature Physics 18, 468 (2022).

[27] C. Scalliet, B. Guiselin, and L. Berthier, Thirty Milliseconds in the Life of a Supercooled Liquid, Phys. Rev. X 12, 041028 (2022).

[28] P. K. Bommineni, N. R. Varela-Rosales, M. Klement, and M. Engel, Complex Crystals from Size-Disperse Spheres, Phys. Rev. Lett. 122, 128005 (2019).

[29] P. K. Bommineni, M. Klement, and M. Engel, Spontaneous Crystallization in Systems of Binary Hard Sphere Colloids, Phys. Rev. Lett. 124, 218003 (2020).

[30] G. Alfaro Miranda, L. F. Cugliandolo, and M. Tarzia, SWAP algorithm for lattice spin models, Phys. Rev. E 110, L043301 (2024).

[31] Y. Nishikawa and L. Berthier, Collective Relaxation Dynamics in a Three-Dimensional Lattice Glass Model, Phys. Rev. Lett. 132, 067101 (2024).
[32] K. Shiraishi and L. Berthier, Characterizing the Slow Dynamics of the Swap Monte Carlo Algorithm, J. Phys. Chem. B **128**, 12279 (2024).

[33] C.-R. Hwang, S.-Y. Hwang-Ma, and S.-J. Sheu, Accelerating Gaussian Diffusions, Ann. Appl. Probab. **3** (1993).

[34] F. Chen, L. Lovász, and I. Pak, Lifting Markov chains to speed up mixing, in *Proceedings of the thirty-first annual ACM symposium on Theory of Computing* (ACM, New York, NY, USA, 1999) pp. 275-281.

[35] P. Diaconis, S. Holmes, and R. M. Neal, Analysis of a nonreversible Markov chain sampler, Ann. Appl. Probab. **10**, 726 (2000).

[36] C.-R. Hwang, S.-Y. Hwang-Ma, and S.-J. Sheu, Accelerating diffusions, Ann. Appl. Probab. **15**, 1433 (2005).

[37] F. Ghimenti, L. Berthier, G. Szamel, and F. van Wijland, Sampling Efficiency of Transverse Forces in Dense Liquids, Phys. Rev. Lett. **131**, 257101 (2023).

[38] K. S. Turitsyn, M. Chertkov, and M. Vucelja, Irreversible Monte Carlo algorithms for efficient sampling, Physica D **240**, 410 (2011).

[39] Lifting—A nonreversible Markov chain Monte Carlo algorithm, Am. J. Phys. **84**, 958 (2016).

[40] E. P. Bernard, W. Krauth, and D. B. Wilson, Event-chain Monte Carlo algorithms for hard-sphere systems, Phys. Rev. E **80**, 056704 (2009).

[41] E. P. Bernard and W. Krauth, Two-step melting in two dimensions: First-order liquid-hexatic transition, Phys. Rev. Lett. **107**, 155704 (2011).

[42] Y. Sakai and K. Hukushima, Dynamics of one-dimensional ising model without detailed balance condition, J. Phys. Soc. Jpn. **82**, 064003 (2013).

[43] M. Michel, S. C. Kapfer, and W. Krauth, Generalized event-chain Monte Carlo: Constructing rejection-free global-balance algorithms from infinitesimal steps, J. Chem. Phys. **140**, 054116 (2014).

[44] M. Michel, J. Mayer, and W. Krauth, Event-chain Monte Carlo for classical continuous spin models, EPL **112**, 20003 (2015).

[45] W. Krauth, Event-Chain Monte Carlo: Foundations, Applications, and Prospects, Front. Phys. **9**, 663457 (2021).

[46] T. Guyon, A. Guillou, and M. Michel, Necessary and sufficient symmetries in Event-Chain Monte Carlo with generalized flows and application to hard dimers, J. Chem. Phys. **160**, 024117 (2024).

[47] A. Monemvassitis, A. Guillou, and M. Michel, PDMP characterisation of event-chain Monte Carlo algorithms for particle systems, J. Stat. Phys. **190**, 66 (2023).

[48] F. Ghimenti, L. Berthier, and F. van Wijland, Irreversible Monte Carlo Algorithms for Hard Disk Glasses: From Event-Chain to Collective Swaps, Phys. Rev. Lett. **133**, 028202 (2024).

[49] L. Berthier, F. Ghimenti, and F. van Wijland, Monte Carlo simulations of glass-forming liquids beyond Metropolis, J. Chem. Phys. **161** (2024).

[50] E. P. Bernard, W. Krauth, and D. Wilson, Event-chain Monte Carlo algorithms for hard-sphere systems, Phys. Rev. E **80**, 056704 (2009).

[51] Y. Nishikawa, M. Michel, W. Krauth, and K. Hukushima, Event-chain algorithm for the Heisenberg model: Evidence for $z \simeq 1$ dynamic scaling, Phys. Rev. E **92**, 063306 (2015).

[52] D. E. Knuth, *The art of computer programming, volume 3: (2nd ed.) sorting and searching* (Addison Wesley Longman Publishing Co., Inc., USA, 1998).

[53] L. Berthier, P. Charbonneau, A. Ninarello, M. Ozawa, and S. Yaida, Zero-temperature glass transition in two dimensions, Nat. Commun. **10**, 1508 (2019).

[54] E. Flenner and G. Szamel, Fundamental differences between glassy dynamics in two and three dimensions, Nat. Commun. **6**, 7392 (2015).

[55] S. Sastry, P. G. Debenedetti, and F. H. Stillinger, Signatures of distinct dynamical regimes in the energy landscape of a glass-forming liquid, Nature **393**, 554 (1998).

[56] A. S. Keys, L. O. Hedges, J. P. Garrahan, S. C. Glotzer, and D. Chandler, Excitations are localized and relaxation is hierarchical in glass-forming liquids, Phys. Rev. X **1**, 021013 (2011).

[57] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, Structural Relaxation Made Simple, Phys. Rev. Lett. **97**, 170201 (2006).

[58] M. Goldstein, Viscous liquids and the glass transition: a potential energy barrier picture, The Journal of Chemical Physics **51**, 3728 (1969).

[59] F. H. Stillinger and T. A. Weber, Hidden structure in liquids, Phys. Rev. A **25**, 978 (1982).

[60] F. H. Stillinger, Supercooled liquids, glass transitions, and the Kauzmann paradox, J. Chem. Phys. **88**, 7818 (1988).

[61] F. H. Stillinger, A Topographic View of Supercooled Liquids and Glass Formation, Science **267**, 1935 (1995).

[62] J. Kurchan and L. Laloux, Phase space geometry and slow dynamics, J. Phys. A: Math. Gen. **29**, 1929 (1996).

[63] S. Büchner and A. Heuer, Potential energy landscape of a model glass former: Thermodynamics, anharmonicities, and finite size effects, Phys. Rev. E **60**, 6507 (1999).

[64] F. Sciortino, W. Kob, and P. Tartaglia, Inherent Structure Entropy of Supercooled Liquids, Phys. Rev. Lett. **83**, 3214 (1999).

[65] F. Sciortino, Potential energy landscape description of supercooled liquids and glasses, J. Stat. Mech. **2005**, P05015 (2005).

[66] A. Heuer, Exploring the potential energy landscape of glass-forming systems: from inherent structures via metabasins to macroscopic transport, J. Phys. Condens. Matter **20**, 373101 (2008).

[67] D. Khomenko, C. Scalliet, L. Berthier, D. R. Reichman, and F. Zamponi, Depletion of two-level systems in ultra-stable computer-generated glasses, Phys. Rev. Lett. **124**, 225901 (2020).

[68] Y. Nishikawa, M. Ozawa, A. Ikeda, P. Chaudhuri, and L. Berthier, Relaxation Dynamics in the Energy Landscape of Glass-Forming Liquids, Phys. Rev. X **12**, 021001 (2022).

[69] E. Lerner, G. Düring, and E. Bouchbinder, Statistics and Properties of Low-Frequency Vibrational Modes in Structural Glasses, Phys. Rev. Lett. **117**, 035501 (2016).

[70] H. Mizuno, H. Shiba, and A. Ikeda, Continuum limit of the vibrational properties of amorphous solids, Proc. Natl. Acad. Sci. U.S.A. **114**, E9767 (2017).

[71] G. Kapteijns, E. Bouchbinder, and E. Lerner, Universal Nonphononic Density of States in 2D, 3D, and 4D Glasses, Phys. Rev. Lett. **121**, 055501 (2018).

[72] L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel, and E. Flenner, Low-frequency vibrational modes of stable glasses, Nat. Commun. **10**, 26 (2019).

[73] C. Brito, E. Lerner, and M. Wyart, Theory for Swap Acceleration near the Glass and Jamming Transitions for Continuously Polydisperse Particles, Phys. Rev. X **8**,
031050 (2018).

[74] G. Kapteijns, W. Ji, C. Brito, M. Wyart, and E. Lerner, Fast generation of ultrastable computer glasses by minimization of an augmented potential energy, Phys. Rev. E **99**, 012106 (2019).

[75] V. F. Hagh, S. R. Nagel, A. J. Liu, M. L. Manning, and E. I. Corwin, Transient learning degrees of freedom for introducing function in materials, Proc. Natl. Acad. Sci. U.S.A. **119**, e2117622119 (2022).

[76] V. Bolton-Lum, R. C. Dennis, P. Morse, and E. Corwin, The Ideal Glass and the Ideal Disk Packing in Two Dimensions, arXiv preprint arXiv:2404.07492 (2024).

[77] B. J. Alder and T. E. Wainwright, Phase Transition for a Hard Sphere System, J. Chem. Phys. **27**, 1208 (1957).

[78] J. Liu and E. Luijten, Rejection-Free Geometric Cluster Algorithm for Complex Fluids, Phys. Rev. Lett. **92**,

035504 (2004).

[79] E. A. J. F. Peters and G. de With, Rejection-free Monte Carlo sampling for general potentials, Phys. Rev. E **85**, 026703 (2012).

[80] A. Bouchard-Côté, S. J. Vollmer, and A. Doucet, The Bouncy Particle Sampler: A Nonreversible Rejection-Free Markov Chain Monte Carlo Method, J. Am. Stat. Assoc. **113**, 855 (2018).

[81] D. C. Rapaport, *The art of molecular dynamics simulation* (Cambridge university press, 2004).

[82] L. Galliano, R. Rende, and D. Coslovich, Policy-guided Monte Carlo on general state spaces: Application to glass-forming mixtures, J. Chem. Phys. **161**, 064503 (2024).

[83] G. Jung, G. Biroli, and L. Berthier, Normalizing flows as an enhanced sampling method for atomistic supercooled liquids, Mach. Learn.: Sci. Technol. **5**, 035053 (2024).