# Creating equilibrium glassy states via random particle bonding

Misaki Ozawa, Jean-Louis Barrat, Walter Kob, Francesco Zamponi

► To cite this version:

Misaki Ozawa, Jean-Louis Barrat, Walter Kob, Francesco Zamponi. Creating equilibrium glassy states via random particle bonding. Journal of Statistical Mechanics: Theory and Experiment, 2024, 2024 (1), pp.013303. 10.1088/1742-5468/ad17b6. hal-04754917

HAL Id: hal-04754917

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

Submitted on 21 Mar 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>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

Creating equilibrium glassy states via random
particle bonding

To cite this article: Misaki Ozawa et al J. Stat. Mech. (2024) 013303

View the [article online](#) for updates and enhancements.

You may also like

- Ionic liquid glasses: properties and applications
Mikhail Yu. Ivanov, Nikolay V. Surovtsev and Matvey V. Fedin

- Solving the spherical p-spin model with the cavity method: equivalence with the replica results
Giacomo Gradenigo, Maria Chiara Angelini, Luca Leuzzi et al.

- Exact theory of dense amorphous hard spheres in high dimension. III. The full replica symmetry breaking solution
Patrick Charbonneau, Jorge Kurchan, Giorgio Parisi et al.
PAPER: Disordered systems, classical and quantum

# Creating equilibrium glassy states via random particle bonding

Misaki Ozawa¹, Jean-Louis Barrat¹, Walter Kob² and Francesco Zamponi³,*

¹ Univ. Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France

² Department of Physics, University of Montpellier and CNRS, F-34095 Montpellier, France

³ Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 2, 00185 Rome, Italy

E-mail: francesco.zamponi@uniroma1.it

Received 17 November 2023

Accepted for publication 30 November 2023

Published 24 January 2024

<figure><img src="image_2.png" /></figure>

CrossMark

Online at stacks.iop.org/JSTAT/2024/013303

https://doi.org/10.1088/1742-5468/ad17b6

**Abstract.** Creating amorphous solid states by randomly bonding an ensemble of dense liquid monomers is a common procedure that is used to create a variety of materials, such as epoxy resins, colloidal gels, and vitrimers. However, the properties of the resulting solid do a *priori* strongly depend on the preparation history. This can lead to substantial aging of the material; for example, properties such as mechanical moduli and transport coefficients rely on the time elapsed since solidification, which can lead to a slow degradation of the material in technological applications. It is therefore important to understand under which conditions random monomer bonding can lead to stable solid states, that is, long-lived metastable states whose properties do not change over time. This work presents a theoretical and computational analysis of this problem and introduces a random bonding procedure that ensures the proper equilibration of the resulting amorphous states. Our procedure also provides a new route to investigate the fundamental properties of glassy energy landscapes by producing translationally invariant ultrastable glassy states in simple particle models.

<figure><img src="image_3.png" /></figure>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
Contents

<table><tr><td>1. Introduction</td><td>3</td></tr><tr><td>2. Statistical mechanics of bonded systems</td><td>4</td></tr><tr><td>2.1. Quiet freezing of variables</td><td>5</td></tr><tr><td>2.2. Random pinning and wall pinning</td><td>7</td></tr><tr><td>2.2.1. Random pinning: choosing at random the particles that will be pinned</td><td>8</td></tr><tr><td>2.2.2. Wall pinning: choosing the particles to be pinned according to a geometrical criterion</td><td>8</td></tr><tr><td>2.2.3. Sorted pinning: choosing the particles to be pinned according to an ordering</td><td>11</td></tr><tr><td>2.3. Random bonding procedures</td><td>12</td></tr><tr><td>2.3.1. Variable transformation to create virtual dimers</td><td>12</td></tr><tr><td>2.3.2. Random bonding</td><td>15</td></tr><tr><td>2.3.3. 'Sorted bonding' from the shortest to the largest bond</td><td>17</td></tr><tr><td>2.3.4. 'Random bonding' within a cutoff</td><td>18</td></tr><tr><td>2.3.5. Random bonding with directional alignment</td><td>20</td></tr><tr><td>3. Simulation methods</td><td>20</td></tr><tr><td>3.1. Model</td><td>20</td></tr><tr><td>3.2. Making randomly bonded systems</td><td>20</td></tr><tr><td>4. Structural properties</td><td>22</td></tr><tr><td>4.1. Radial distribution function</td><td>22</td></tr><tr><td>4.2. Orientational correlation function</td><td>23</td></tr><tr><td>4.3. Random bonding with directional alignment</td><td>23</td></tr><tr><td>5. Dynamical properties</td><td>24</td></tr><tr><td>5.1. Intermediate scattering functions</td><td>25</td></tr><tr><td>5.2. Scaling of the dynamics and fragility</td><td>28</td></tr><tr><td>6. Dynamical heterogeneity</td><td>29</td></tr><tr><td>6.1. Positional degrees of freedom</td><td>30</td></tr><tr><td>6.2. Rotational degrees of freedom</td><td>31</td></tr><tr><td>7. Conclusion and discussion</td><td>33</td></tr><tr><td>Acknowledgments</td><td>34</td></tr><tr><td>References</td><td>34</td></tr></table>
# 1. Introduction

Since glasses are out-of-equilibrium materials, their properties depend on the history of formation [1]. A glass quenched from its liquid state with a slow cooling rate will have higher kinetic and mechanical stability with lower energy (or enthalpy) than if cooled quickly, which reflects the fact that the system resides in a lower region of its rugged energy landscape [2, 3]. Producing stable glasses and hence accessing deep minima inside this landscape is a crucial challenge for experiments and computer simulations, as it allows one to devise better materials and/or gain insight into the nature of the glassy state. In experiments, vapor deposition techniques with controlled substrate temperature can produce glass samples with extraordinary kinetic stability, a major breakthrough in the last decades [4–10]. On the atomistic simulation side, various sampling techniques, such as replica exchange methods [11–13], Monte-Carlo (MC) simulations with smart updates [14–17], random pinning methods [18, 19], and machine learning-assisted sampling techniques [20–26], have been developed. In particular, the swap MC [17], its generalizations [27, 28], and the random pinning approach [19, 29, 30] allow to generate equilibrium configurations deep inside the glassy landscape. The key idea behind these approaches is to systematically vary certain degrees of freedom while maintaining the thermal equilibrium properties of the system. In the swap MC algorithm [17], the diameters of the particles are the new dynamical variables that are allowed to fluctuate, significantly accelerating the relaxation dynamics. Although this method has revolutionized computational glass physics, it has so far remained in the realm of computer simulations, and extending it to real experiments has turned out to be challenging. In the random pinning approach [19], the positions of a fraction of the particles are permanently frozen; hence, the system comprising the remaining mobile particles enters a very glassy state with strong confinements due to the presence of the pinned particles. Higher concentrations of pinned particles lead to ideal glasses [29–31], whose thermodynamic behavior has been found to be consistent with the random first-order transition theory [32]. A big advantage of this method is that it can be applied to experimental systems, such as colloids [33, 34] and molecular liquids [35, 36]. However, the dynamics is significantly altered from the bulk state, possibly because of the violation of translational invariance, leading to a strong decoupling between self and collective behavior [30, 37, 38], a decrease in fragility with the increased concentration of pinned particles [39, 40], and a suppression of dynamical heterogeneity when the glass transition is approached [39, 41–43], in stark contrast to the behavior of standard bulk materials that are approaching their glass transition.

Recently, we have proposed a random bonding method, in which one creates a bond between a pair of randomly chosen nearest neighbor particles [44]. This idea is inspired by random pinning [19], but random bonding has the big advantage that it preserves the translational invariance of the system and that it can also be realized without much difficulty in real experiments. In fact, random bonding is routinely used to prepare amorphous solids such as epoxy resins (via curing agents [45] or stereolithography [46]), and vitrimers [47, 48]. It has also been used to prepare colloidal or emulsion clusters of various shapes via programmable bond activation [49, 50] using temperature control [51], salt addition [52], or UV light [53]. The problem is that, to the best of our knowledge,
it is not clear whether and under which conditions these techniques can produce stable glasses [46, 54–57]. In this work, using theoretical analysis and computer simulations to probe the kinetic and mechanical properties of the system, we demonstrate that the random bonding method does indeed create ultrastable glasses in the bulk [44]. We confirm and theoretically support the preliminary molecular dynamics simulations of [44], which suggested that the relaxation dynamics does not show aging within numerical accuracy. This implies that, right after bonding, the resulting configuration is indeed close to equilibrium [44], akin to random pinning [58].

More specifically, in this work, we demonstrate that although the system is in equi-
librium if the bonded particles are selected completely randomly (in this case, the bond
lengths are thus arbitrary), the bonding between pairs of nearest neighbor particles
(which corresponds to a more realistic situation) does not ensure strict equilibration.
However, it turns out that in practice the deviation from equilibrium is very small, and
it can be negligible for most practical purposes, which will be demonstrated by detailed
molecular simulations.

Subsequently, we present results on the (almost) equilibrium dynamics of randomly bonded glass-forming liquids and contrast our findings with those from the dynamics of randomly pinned systems. We find that the self and collective correlation functions for the translational degrees of freedom, as well as the rotational correlation functions, are strongly coupled. We also observe that the kinetic fragility does not change with increasing concentration of bonds. Finally, we find that the dynamical heterogeneity keeps growing as we approach the glass transition. These trends are thus opposite to the behavior found in the dynamics of randomly pinned systems, highlighting the importance of the nature of the quenched disorder.

Our paper is organized as follows. Section 2 presents a general statistical mechanics treatment of the process of freezing some degrees of freedom. We apply this treatment to the random bonding process and show under which conditions it produces equilibrated configurations. Section 3 discusses the simulation methods and procedures to generate randomly bonded configurations. Section 4 describes the structural properties of our generated configurations and shows that they are compatible with them being essentially equilibrated. Section 5 explores the dynamical properties of the generated glasses and shows that no aging is detected. Section 6 demonstrates that the randomly bonded systems display dynamical heterogeneities similar to standard glasses. Finally, section 7 presents our conclusions.

## 2. Statistical mechanics of bonded systems

In this section, we present a general treatment of the statistical mechanics of systems with frozen degrees of freedom [17, 19, 27–30, 59, 60]. In particular, we show that there are subtle conditions that need to be fulfilled in order to ensure that once some degrees of freedom are frozen, the remaining ones are still in equilibrium conditioned to the frozen ones [61]. We then apply this framework to randomly bonded glass formers comprising monomers and dimers. We randomly choose the pairs of neighboring particles from an equilibrium configuration and bond them together by introducing a rigid-body constraint [44]. The main goal of this section is to demonstrate that while the creation
of these bonds does not strictly enforce equilibrium, it does not perturb the thermal equilibrium of the system in a significant manner.

Because the nature of the problem is quite subtle, we begin by a pedagogical presentation of the general idea (section 2.1) and its application to the random pinning protocol (section 2.2) that fixes the positions of particles [58, 62]. Subsequently, we discuss how to extend similar ideas to random bonding (section 2.3).

## 2.1. Quiet freezing of variables

Consider a system whose degrees of freedom are arbitrarily split into two distinct vectors $\mathbf{x}$ and $\mathbf{y}$ with Hamiltonian $H(\mathbf{x}, \mathbf{y})$. The equilibrium probability distribution of the total system is $\rho(\mathbf{x}, \mathbf{y}) = \exp[-\beta H(\mathbf{x}, \mathbf{y})]/Z$, where $Z$ is the partition function and $\beta$ is the inverse temperature. Let us assume that an equilibrium configuration of the system can be generated by sampling from $\rho(\mathbf{x}, \mathbf{y})$. Consider now the degrees of freedom $\mathbf{y}$. Their statistics is described by the marginal probability distribution (we add a subscript ‘f’ because these are the degrees of freedom that will be frozen in the following):

$$ \rho_f(\mathbf{y}) = \int d\mathbf{x} \rho(\mathbf{x}, \mathbf{y}) = \frac{Z_f(\mathbf{y})}{Z} \quad \text{with} \quad Z_f(\mathbf{y}) = \int d\mathbf{x} e^{-\beta H(\mathbf{x}, \mathbf{y})}. \quad (1) $$

Next, consider a setting in which one first generates a configuration of the $\mathbf{y}$ degrees of freedom from their marginal probability distribution $\rho_f(\mathbf{y})$, and then considers a system wherein $\mathbf{y}$ are ‘frozen’ and whose dynamical variables are the $\mathbf{x}$ degrees of freedom, with equilibrium probability:

$$ \rho(\mathbf{x}|\mathbf{y}) = \frac{e^{-\beta H(\mathbf{x},\mathbf{y})}}{Z_f(\mathbf{y})}. \quad (2) $$

Here, the $\mathbf{y}$ variables play the role of a frozen or pinned quenched disorder, and one is interested in the thermal properties of the system described by $\mathbf{x}$.

Using the chain rule of probabilities, we have:

$$ \rho(\mathbf{x}, \mathbf{y}) = \rho(\mathbf{x}|\mathbf{y}) \rho_f(\mathbf{y}). \quad (3) $$

Hence, a pair $\{\mathbf{x}, \mathbf{y}\}$ generated from the joint distribution $\rho(\mathbf{x}, \mathbf{y})$ can be considered either as an equilibrium configuration of the full system or as a realization $\mathbf{y}$ of the quenched disorder of the frozen system obtained from $\rho_f(\mathbf{y})$ together with a typical equilibrium realization $\mathbf{x}$ of that system obtained from $\rho(\mathbf{x}|\mathbf{y})$. We conclude that one can generate a pair $\{\mathbf{x}, \mathbf{y}\}$ from $\rho(\mathbf{x}, \mathbf{y})$, then freeze (or pin) $\mathbf{y}$ and automatically obtain an equilibrium configuration $\mathbf{x}$ (but only a single one) of the distribution $\rho(\mathbf{x}|\mathbf{y})$ that describes the frozen system with quenched disorder $\mathbf{y}$.

We now define an ‘annealed’ average $\langle \cdots \rangle$, a ‘thermal’ average $\langle \cdots \rangle_y$, and a ‘disorder’ average $\overline{\cdots}$ by

$$ \langle \cdots \rangle = \int d\mathbf{x} d\mathbf{y} \rho(\mathbf{x}, \mathbf{y}) (\cdots), \quad \langle \cdots \rangle_y = \int d\mathbf{x} \rho(\mathbf{x}|\mathbf{y}) (\cdots), \quad \overline{\cdots} = \int d\mathbf{y} \rho_f(\mathbf{y}) (\cdots), \quad (4) $$
respectively. Furthermore, we define the ‘quenched’ average as $\langle \cdots \rangle_y$ and because of the chain rule in equation (3), the annealed and quenched averages coincide:

$$ \langle \cdots \rangle = \overline{\langle \cdots \rangle}_y . \qquad (5) $$

Conversely, the validity of equation (5) for an arbitrary observable implies equation (3). Equations (3) and (5) provide the core idea behind the random pinning and swap MC simulations, as we will discuss in detail below.

Since the freezing of degrees of freedom eliminates possible relaxation channels [17, 28, 59, 60, 63], it can be expected that sampling $x$ from $\rho(x|y)$ at fixed $y$ is harder than sampling $x$ and $y$ together from $\rho(x,y)$. (In some cases, however, the conditional sampling gives rise to faster sampling, e.g. [64]). In other words, a local dynamics for $x$ that satisfies a detailed balance with respect to $\rho(x|y)$ will in general have a much larger decorrelation time than a similar local dynamics that acts on $x$ and $y$ and satisfies detailed balance with respect to $\rho(x,y)$. Generating additional independent configurations from $\rho(x|y)$ for the same $y$ might thus be a hard task. Still, having a single equilibrium configuration $x$ of the pinned system allows one to run local dynamics starting from that configuration and thus obtain equilibrium dynamical properties without having to worry about the process of equilibration itself.

A few remarks on this construction are in order at this point:

* The distribution $\rho_f(y)$ depends on temperature and hence a different ensemble of pinned systems is obtained at each preparation temperature. The configuration $x$ is in equilibrium at the same temperature. Once $y$ is frozen, one is still allowed to change the temperature of $x$, but in that case equilibrium is not guaranteed anymore.

* In systems with quenched disorder, the thermal degrees of freedom $x$ are usually described by $\rho(x|y)$ as in equation (2), but the distribution of the quenched disorder $y$, that is, $\rho_q(y)$, is selected independently, and in general $\rho_q(y) \neq Z_f(y)/Z$, where $Z_f(y)$ is defined in equation (1). Physically, this corresponds to the fact that the quenched disorder $y$ represents impurities in the Hamiltonian $H(x,y)$ (e.g., the location of the magnetic atoms in a magnetic alloy that forms a spin glass) whose dynamics is extremely slow, thus preventing them to equilibrate with the other degrees of freedom $x$ and as a consequence the annealed and quenched averages do not coincide. Note that usually one keeps $\rho_q(y)$ fixed while changing the temperature associated with $x$. The choice $\rho_q(y) = \rho_f(y) = Z_f(y)/Z$, in which the quenched disorder depends on temperature, guarantees the equality of annealed and quenched averages in equation (5). It is a very special choice and is called *Nishimori condition* in the physics literature, *quiet planting* in optimization, and *Bayes optimal condition* in statistical inference. See [61] for a pedagogical discussion.

* The thermal average $\langle O(x) \rangle_y$ of an extensive observable, which is the sum of local terms of the form

$$ O(x) = \sum_{i} o_{1}(x_{i}) + \sum_{i,j} o_{2}(x_{i}, x_{j}) + \cdots \qquad (6) $$

is a random variable that depends on the realization of the disorder $y$ (the ‘sample’). However, in the thermodynamic limit, disordered systems usually display the so-called
self-averaging property: the disorder-induced fluctuations of thermal averages vanish and as a consequence the thermal average for a single typical sample coincides with the average over samples, namely, $\langle \cdots \rangle_y = \langle \cdots \rangle_y$. Furthermore, note that the observable $O(\mathbf{x}, \mathbf{y})$ can also depend on the frozen degrees of freedom $\mathbf{y}$ (e.g., if $O=H$), in which case its thermal average also depends explicitly on $\mathbf{y}$ through $O$. This does not affect the proof of equivalence of the annealed and quenched averages in equation (5), and it also does not affect the self-averaging property. Note that self-averaging only holds for extensive variables, that is, the averages of local terms over the whole system or at least a finite fraction of it in the thermodynamic limit. As an example, in a system that has been subjected to a random pinning procedure, the pair correlation function measured in a large configuration is independent of the choice of the pinned particles, and is self averaging. In contrast, the density at a specific point in space will depend on this choice, and is not self averaging.

* A local dynamics on $\mathbf{x}$ that samples from $\rho(\mathbf{x}|\mathbf{y})$ could be so slow that it becomes non-ergodic at the ideal glass transition. This has been shown to happen in randomly pinned mean field spin glasses [19, 65] and finite-dimensional glasses [29–31]. In this case, the system remains stuck forever around the initial equilibrium configuration $\mathbf{x}$, and generating independent additional configurations is impossible using standard simulation approaches such as simple MC simulations or molecular dynamics.

It is crucial to stress that the above construction, which we will call the *quiet freezing* of variables, supposes that the division of the degrees of freedom into $\mathbf{x}$ and $\mathbf{y}$ is done *before* the construction of a thermalized configuration. In other words, the proof described above does not hold if we first thermalize a system and, then, choose which degrees of freedom to freeze by using properties of the thermalized configuration because this would introduce a bias that is not described by equation (1). Keeping this in mind, we now discuss how this construction can be applied to randomly pinned particle systems, and how it can be generalized when pinning is not random.

## 2.2. Random pinning and wall pinning

The construction in section 2.1 can be applied in a straightforward manner to the random pinning procedure [18, 58, 62]. We will consider a system of $N$ point particles in $d$ dimensions, such as the Kob-Andersen model [66], and we will make use of the particle permutation symmetry as a crucial ingredient. Note that in most cases, the systems of interest are polydisperse, including binary or ternary mixtures, and the particle permutation symmetry does not exist because particles have distinct sizes. Yet, one can always recover it if one considers the permutation of the particle species as an additional degree of freedom and sum up all the possible permutations in the partition function [27, 67–69]. This treatment does not change the thermodynamic properties that we discuss in this paper. A formally equivalent way of reintroducing the particle permutation symmetry is to assign to each particle an additional degree of freedom describing its size. Thus, for simplicity, we focus only on positions (and momenta) as relevant degrees of freedom.
2.2.1. Random pinning: choosing at random the particles that will be pinned. The Hamiltonian $H$ and the partition function $Z$ are given by:

$$H(\mathbf{r}^N, \mathbf{p}^N) = \sum_{i=1}^{N} \frac{\mathbf{p}_i^2}{2m_i} + U(\mathbf{r}^N), \quad Z = \int \left( \prod_{i=1}^{N} \mathrm{d}\mathbf{r}_i \mathrm{d}\mathbf{p}_i \right) \exp[-\beta H(\mathbf{r}^N, \mathbf{p}^N)], \quad (7)$$

where $\mathbf{r}_i$, $\mathbf{p}_i$, and $m_i$ are the position, momentum, and mass of the $i$th particle, respectively. $U$ is the potential energy, and $\beta = 1/T$ is the inverse temperature. Particles are assumed to be confined in a volume $V$ with some boundary conditions that do not need to be specified at this stage. In this paper, we use a shorthand notation for a vector of $N$ variables, for example, $\mathbf{r}^N = (\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_N)$. Note that we omit the combinatorial factors such as $N!$ and the Planck's constant $h$ in the definition of the partition function in equation (7) because we are not concerned with the absolute value of the free energy or entropy.

One can then split the $N$ particles into a set of $N_f$ pinned particles, with coordinates $\mathbf{y} = \{\mathbf{r}^{N_f}, \mathbf{p}^{N_f}\}$, and a set of $N - N_f$ unpinned (or mobile) particles with the remaining coordinates, $\mathbf{x}$. Because the labeling of particles is arbitrary, we consider the first $N_f$ particles as pinned ones. Note that the splitting is here performed before any thermalization, as we emphasized in section 2.1. One can then thermalize the whole system of $N$ particles; hence, the joint set $\{\mathbf{x}, \mathbf{y}\}$ is in the liquid phase. Because particles can freely diffuse in such phase, we expect the pinned particles to be uniformly distributed in the volume $V$. At this point, the positions of the $N_f$ particles, $\mathbf{r}^{N_f}$, are frozen and their velocities are set to zero. Because the distribution of momenta is a product of independent distributions for each particle, setting the momenta of pinned particles to zero does not affect the distribution of $\mathbf{x}$ (in other words, only the configurational part of the Hamiltonian matters). Hence, the particles $\mathbf{x}$ can be considered to be in equilibrium with the pinned particles $\mathbf{y}$, and the procedure allows one to generate an equilibrium configuration of the pinned system even in the case in which the density of pinned particles is so high that the pinned system is glassy. Note that in usual implementations of the random pinning procedure, one first generates a configuration of the full system of $N$ particles and then randomly chooses the $N_f$ particles to be pinned. However, since from the point of view of the full system, this is just a labeling that does not depend on the equilibrated configuration, the two operations (randomly selecting the particles that are going to be pinned and equilibrating the full system) can be safely exchanged. We thus conclude that random pinning belongs to the class of quiet freezing procedures, as it is well known and numerically validated.

2.2.2. Wall pinning: choosing the particles to be pinned according to a geometrical criterion. In [62], the random pinning construction has been extended to introduce the possibility of choosing which particles have to be frozen after the equilibrated configuration of the unpinned system is prepared. This ‘wall pinning’ construction has been widely applied in the context of glass physics to probe the equilibrium properties of liquids [70–73]. We briefly describe the proof for completeness. Suppose that the total volume $V$ is split into a ‘wall’ region $W$, in which particles will eventually be frozen, and a ‘fluid’ region $F$, such that $V = W + F$, as schematically shown in figure 1. The shape of these regions is completely arbitrary. We can then generate an equilibrium
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Illustration of the wall pinning procedure, of which we give two different proofs in the main text. The first, as in [62], is based on the idea of separating the volume in the wall (bottom part) and fluid (top part) regions, and freezing particles in the former region; the corresponding averages are given in equation (10). The second corresponds to ordering the particles according to their z component, and freezing the k particles with the smallest z, as in equation (13). The number k has to be selected such that the wall region has the desired size.</figcaption></figure>

configuration of the full N-particle system and freeze the particles that fall into the W region. Because the choice of degrees of freedom that are going to be frozen depends on the equilibrated configuration, the proof of section 2.1 does not apply directly and must be generalized. We note that the momenta are always irrelevant (their distribution is a product over particles) and thus can be ignored. Given the total configurational partition function

$$Z_c = \int \left( \prod_{i=1}^{N} \mathrm{d}\mathbf{r}_i \right) e^{-\beta U(\mathbf{r}^N)}, \qquad (8)$$

the average of an arbitrary observable $O(\mathbf{r}^N)$ that is invariant under permutations of the particle labels can be written as:

$$ \begin{aligned} \langle O \rangle &= \frac{1}{Z_c} \int_{W+F} \mathrm{d}\mathbf{r}_1 \cdots \int_{W+F} \mathrm{d}\mathbf{r}_N O(\mathbf{r}^N) e^{-\beta U(\mathbf{r}^N)} \\ &= \frac{1}{Z_c} \left[ \int_W \mathrm{d}\mathbf{r}_1 + \int_F \mathrm{d}\mathbf{r}_1 \right] \cdots \left[ \int_W \mathrm{d}\mathbf{r}_N + \int_F \mathrm{d}\mathbf{r}_N \right] O(\mathbf{r}^N) e^{-\beta U(\mathbf{r}^N)} \\ &= \frac{1}{Z_c} \sum_{k=0}^{N} \binom{N}{k} \int_W \mathrm{d}\mathbf{r}_1 \cdots \int_W \mathrm{d}\mathbf{r}_k \int_F \mathrm{d}\mathbf{r}_{k+1} \cdots \int_F \mathrm{d}\mathbf{r}_N O(\mathbf{r}^N) e^{-\beta U(\mathbf{r}^N)}, \quad (9) \end{aligned} $$
where $\binom{N}{k}$ is the binomial coefficient. In the last step, we defined $k$ as the number of particles in the $W$ region, and we used the permutation symmetry of both $O(\mathbf{r}^N)$ and $U(\mathbf{r}^N)$, hence of the whole integrand, to relabel the first $k$ particles as being those in the $W$ region and the remaining $N-k$ as being those in the $F$ region.

Using the shorthand notation $\mathbf{r}^W = (\mathbf{r}_1, \dots, \mathbf{r}_k) \in W^k = \Omega_W$ and $\mathbf{r}^F = (\mathbf{r}_{k+1}, \dots, \mathbf{r}_N) \in F^{N-k} = \Omega_F$, we can write

$$ \begin{aligned} \langle O \rangle &= \frac{1}{Z_c} \sum_{k=0}^{N} \binom{N}{k} \int_{\Omega_W} d\mathbf{r}^W \int_{\Omega_F} d\mathbf{r}^F O(\mathbf{r}^W, \mathbf{r}^F) e^{-\beta U(\mathbf{r}^W, \mathbf{r}^F)} \\ &= \sum_{k=0}^{N} \binom{N}{k} \int_{\Omega_W} d\mathbf{r}^W \frac{Z_W(\mathbf{r}^W)}{Z_c} \int_{\Omega_F} d\mathbf{r}^F O(\mathbf{r}^W, \mathbf{r}^F) \frac{e^{-\beta U(\mathbf{r}^W, \mathbf{r}^F)}}{Z_W(\mathbf{r}^W)} \quad \text{with} \\ Z_W(\mathbf{r}^W) &= \int_{\Omega_F} d\mathbf{r}^F e^{-\beta U(\mathbf{r}^W, \mathbf{r}^F)}. \end{aligned} \quad (10) $$

We conclude, as in section 2.1 that we can define an ‘annealed’ average $\langle \cdots \rangle$, a ‘thermal’ average $\langle \cdots \rangle_W$ conditioned to the wall, and a ‘disorder’ average $\overline{\cdots}$ over the realizations of the wall, from

$$ \begin{aligned} \langle \cdots \rangle &= \int d\mathbf{r}^N \frac{e^{-\beta U(\mathbf{r}^N)}}{Z_c}(\cdots) , \quad \langle \cdots \rangle_W = \int_{\Omega_F} d\mathbf{r}^F \frac{e^{-\beta U(\mathbf{r}^W, \mathbf{r}^F)}}{Z_W(\mathbf{r}^W)}(\cdots) , \\ \overline{\cdots} &= \sum_{k=0}^{N} \binom{N}{k} \int_{\Omega_W} d\mathbf{r}^W \frac{Z_W(\mathbf{r}^W)}{Z_c}(\cdots) , \end{aligned} \quad (11) $$

and because of the chain rule in equation (10), the annealed and quenched averages coincide, that is $\langle \cdots \rangle = \langle \overline{\cdots} \rangle_W$. Clearly, this proof can be generalized to whatever situation in which (i) there is permutation symmetry over $N$ degrees of freedom (recall that this also holds for polydisperse systems if one considers the particle species as additional degrees of freedom) and (ii) the integration space of each individual degree of freedom can be split a priori (i.e. independently of the configuration of the system) into two distinct regions $W$ and $F$. One can then generate a full equilibrium configuration of the $N$ degrees of freedom and freeze those falling into region $W$, which produces an equilibrium configuration of the remaining degrees of freedom conditioned to the frozen ones. We note that the unfrozen (‘fluid’) degrees of freedom should be constrained inside the $\Omega_F$ region during the thermal average. This constraint is nearly satisfied in dense particle systems where an excluded volume effect prevents the fluid particles from entering the $\Omega_W$ region. Yet, this is not the case for dilute systems; hence, one has to impose additional constraints on the dynamics of fluid particles such as the hard wall condition.

Note that the only requirement on the observable $O(\mathbf{r}^N)$ is that it is invariant under permutations of particles. We also mention that one can construct an observable $O(\mathbf{r}^N)$
that only depends on the particles that are in the F region while keeping the global permutation symmetry of the observable. This can be done, for example, by writing:

$$O(\mathbf{r}^N) = \sum_{i=1}^{N} o_1(\mathbf{r}_i) \mathbb{I}[\mathbf{r}_i \in F] + \sum_{i=1}^{N} \sum_{j=1}^{N} o_2(\mathbf{r}_i, \mathbf{r}_j) \mathbb{I}[\mathbf{r}_i \in F] \mathbb{I}[\mathbf{r}_j \in F] + \cdots, \quad (12)$$

where $\mathbb{I}[\mathcal{E}]$ is the indicator function of event $\mathcal{E}$, which is one if $\mathcal{E}$ is realized and zero otherwise. Then, any observable of the form in equation (12) can be used to characterize the fluid system only while keeping the quenched and annealed averages coincident.

### 2.2.3. Sorted pinning: choosing the particles to be pinned according to an ordering.
We now consider an alternative proof for quiet freezing in a similar spirit by introducing ‘sorted pinning’ described as follows. Suppose that an ordering relation between the N particle positions can be defined. One example is to use the binary ordering relation $\mathbf{r}_i \prec \mathbf{r}_j$ if $\mathbf{r}_i \preceq \mathbf{r}_j$ if $\mathbf{r}_i \cdot \mathbf{e} < \mathbf{r}_j \cdot \mathbf{e}$, where $\mathbf{e}$ is an arbitrary unit vector. Similarly, $\mathbf{r}_i \preceq \mathbf{r}_j$ if $\mathbf{r}_i \cdot \mathbf{e} \le \mathbf{r}_j \cdot \mathbf{e}$. We will see that the proof given below can be generalized to any ordering operation that is able to uniquely sort a set of positions, $\{\mathbf{r}_1, \cdots, \mathbf{r}_N\}$. We now define the vector for the k first ordered particles and the vector for the rest of the ordered particles as $\mathbf{r}^{\preceq k} = (\mathbf{r}_1, \cdots, \mathbf{r}_k)$ and $\mathbf{r}^{\succ k} = (\mathbf{r}_{k+1}, \cdots, \mathbf{r}_N)$, respectively. Considering again only permutation-symmetric observables $O(\mathbf{r}^N)$, and using the permutation symmetry, we can order N particles and write:

$$ \begin{aligned} \langle O \rangle &= \frac{N!}{Z_c} \int_V d\mathbf{r}_1 \int_{\mathbf{r}_2 \preceq \mathbf{r}_1} d\mathbf{r}_2 \cdots \int_{\mathbf{r}_N \preceq \mathbf{r}_{N-1}} d\mathbf{r}_N O(\mathbf{r}^N) e^{-\beta U(\mathbf{r}^N)} \\ &= \int_{\Omega^{\preceq k}} d\mathbf{r}^{\preceq k} \frac{N! Z_k(\mathbf{r}^{\preceq k})}{Z_c} \int_{\Omega^{\succ k}} d\mathbf{r}^{\succ k} \frac{e^{-\beta U(\mathbf{r}^{\preceq k}, \mathbf{r}^{\succ k})}}{Z_k(\mathbf{r}^{\preceq k})} O(\mathbf{r}^{\preceq k}, \mathbf{r}^{\succ k}) \quad \text{with} \\ Z_k(\mathbf{r}^{\preceq k}) &= \int_{\Omega^{\succ k}} d\mathbf{r}^{\succ k} e^{-\beta U(\mathbf{r}^{\preceq k}, \mathbf{r}^{\succ k})}, \end{aligned} \quad (13) $$

where

$$ \int_{\Omega^{\preceq k}} d\mathbf{r}^{\preceq k} = \int_V d\mathbf{r}_1 \int_{\mathbf{r}_2 \preceq \mathbf{r}_1} d\mathbf{r}_2 \cdots \int_{\mathbf{r}_k \preceq \mathbf{r}_{k-1}} d\mathbf{r}_k, \quad (14) $$

$$ \int_{\Omega^{\succ k}} d\mathbf{r}^{\succ k} = \int_{\mathbf{r}_{k+1} \preceq \mathbf{r}_k} d\mathbf{r}_{k+1} \int_{\mathbf{r}_{k+2} \preceq \mathbf{r}_{k+1}} d\mathbf{r}_{k+2} \cdots \int_{\mathbf{r}_N \preceq \mathbf{r}_{N-1}} d\mathbf{r}_N. \quad (15) $$

We thus obtain a relation similar to equation (10).

The core idea both in the ‘wall pinning’ proof (section 2.2.2) and in the ‘sorted pinning’ proof (this section) is that the integral regions for the pinned and unpinned (fluid) particles have a separation, either by a wall that we specify or by ranking via some ordering operation. Thus, as long as a clear-cut separation is enforced, the internal ordering constraint on the two vectors $\mathbf{r}^{\preceq k}$ and $\mathbf{r}^{\succ k}$ in equation (13) can be released. This results in a similar expression for the wall pinning case, but using a separation between the kth and $(k+1)$th ordered particles:
Creating equilibrium glassy states via random particle bonding

$$ \langle O \rangle = \binom{N}{k} \int_{V^k} dr^{\le k} \frac{Z_k(\mathbf{r}^{\le k})}{Z_c} \int_{\substack{\mathbf{r}_i \ge \mathbf{r}_{\max} \\ (\forall i > k)}} dr^{> k} \frac{e^{-\beta U(\mathbf{r}^{\le k}, \mathbf{r}^{> k})}}{Z_k(\mathbf{r}^{\le k})} O(\mathbf{r}^{\le k}, \mathbf{r}^{> k}) \quad \text{with} $$

$$ \mathbf{r}_{\max} = \max\{\mathbf{r}_1, \dots, \mathbf{r}_k\} \quad \text{and} \quad Z_k(\mathbf{r}^{\le k}) = \int_{\substack{\mathbf{r}_i \ge \mathbf{r}_{\max} \\ (\forall i > k)}} dr^{> k} e^{-\beta U(\mathbf{r}^{\le k}, \mathbf{r}^{> k})}. \quad (16) $$

This relationship is identical to equation (10) with the only difference that $k$ is now fixed and the boundary between the wall and fluid regions is fluctuating and determined by the largest of the first $k$ vectors. Hence, this second proof corresponds to a ‘fixed-$k$’ ensemble, whereas the first proof corresponds to a ‘fixed-boundary’ ensemble. Typically, the two ensembles become equivalent in the thermodynamic limit if the wall and fluid regions are both macroscopic, that is, if $k \sim N$. An illustration is given in figure 1.

While the derivations presented here and in section 2.2.2 concern the thermodynamics of the system, it is useful to also discuss their dynamical meaning. Suppose that we are performing some dynamics of the $N$-particle system that results, at a given time $t=0$, in an equilibrium configuration $\mathbf{r}^N$. At that instant, we can perform the sorting of the particles (e.g., according to their $z$ component and setting $\mathbf{e} = \mathbf{e}_z$) and identify the first $k$ ‘wall’ particles and the last $N-k$ ‘fluid’ ones. Now, if we let the system evolve in an unconstrained way, the $z$ coordinate of the $k$th particle might cross the $(k+1)$th one. But if this happens, we would just relabel the particles by exchanging $k \leftrightarrow k+1$, such that at any time, the first $k$ particles would have the smallest $z$ values. This dynamics would result in the ‘annealed’ thermodynamic average. Alternatively, at time $t=0$, we can freeze the positions of the first $k$ particles, and only let the remaining $N-k$ particles evolve. We know that these $N-k$ particles start in equilibrium with the wall, but then, we have the additional hard constraint that the $z$ coordinates of the evolving particles must be at any time $t>0$ larger than $z_{\max} = \max\{z_1, \dots, z_k\}$. This hard constraint can be implemented in a MC simulation by rejecting moves that would bring a particle at $z < z_{\max}$, or in Molecular Dynamics by adding a reflecting wall at $z = z_{\max}$. In both cases, the resulting dynamics will lead to the ‘thermal’ average $\langle \cdots \rangle_W$ over the fluid particles. One should then either perform the disorder average over the wall particles, that is $\cdots$, by repeating the freezing procedure many times, or invoke the self-averaging properties for large $N$ (assuming $k$ to be of order $N$) to claim that a single run is representative of the average.

Note that, as in the wall pinning case, we can construct an observable $O(\mathbf{r}^N) = O(\mathbf{r}^{> k})$ that depends only on the last $N-k$ ordered vector $\mathbf{r}^{> k} = (\mathbf{r}_{k+1}, \dots, \mathbf{r}_N)$. Such an observable is still invariant under permutations of the $N$ particles, and it can describe the fluid region without an explicit dependence on the frozen particles.

## 2.3. Random bonding procedures

### 2.3.1. Variable transformation to create virtual dimers.
We now show that the ideas presented in section 2.2 can be implemented within the random bonding approach. Suppose that we want to create $N_d$ dimers and $N_m$ monomers from our system with $N = N_m + 2N_d$ particles. To this aim, given a configuration $\mathbf{r}^N$ of the $N$ monomers, what we need is a procedure that creates $N_d$ dimers in a permutation-invariant way. More precisely, the criteria that are used to decide which particles are going to be in the
dimers should not depend on the particle labeling itself. The remaining $N_m$ particles are left unbonded. We can then exploit the permutation symmetry to indicate the indices of dimers and monomers as belonging to sets $i \in \mathcal{D} = \{1, 3, 5, \dots, 2N_d - 1\}$ and $i \in \mathcal{M} = \{2N_d + 1, \dots, N\}$, respectively, with dimer $i$ having particles $i$ and $i+1$. We then denote $\mathbf{r}^{2N_d} = (\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_{2N_d})$ and $\mathbf{r}^{N_m} = (\mathbf{r}_{2N_d+1}, \dots, \mathbf{r}_N)$ the dimer and monomer coordinates, respectively. We stress that this is just a sorting operation of the particle labels, and no physical constraint is imposed on the system at this stage.

The splitting of all particles into dimers and monomers defines an integration space for both sets of variables, which we denote as $\Omega_d$ for the dimers and $\Omega_m$ for the monomers. Due to permutation invariance, there are

$$ \Pi_{m,d} = \frac{N!}{2^{N_d} N_d! N_m!} \qquad (17) $$

distinct equivalent ways of constructing the dimers<fnref n="4" />, and the sum over all these equivalent possibilities reconstructs the whole integration volume of the original monomer system<fnref n="5" />. Taking advantage of the permutation symmetry of the problem, we can then write:

$$ \langle O \rangle = \frac{\Pi_{m,d}}{Z} \int \left( \prod_{i=1}^{N} d\mathbf{p}_i \right) \int_{\Omega_d} d\mathbf{r}^{2N_d} \int_{\Omega_m} d\mathbf{r}^{N_m} O(\mathbf{r}^N) e^{-\beta H(\mathbf{r}^N, \mathbf{p}^N)}. \qquad (19) $$

We next consider a variable transformation from Cartesian coordinates to Jacobi coordinates for the dimers. We assign virtual bonds between particles in each dimer, as schematically shown in figure 2(a). Once again, we emphasize that this is a virtual operation, and the actual system is not altered at all, that is, this is merely a variable transformation. For the monomers, we continue to use the Cartesian coordinates. For the dimers, instead, we use the Jacobi coordinates for the two-body problem (figure 2(b)), using the center of mass and relative position for the dimer $i \in \mathcal{D} = \{1, 3, 5, \dots, 2N_d - 1\}$ comprising the $i$th and $(i+1)$th particles. These are given by $\mathbf{R}_i = \frac{m_i \mathbf{r}_i + m_{i+1} \mathbf{r}_{i+1}}{m_i + m_{i+1}}$ and $\tilde{\mathbf{r}}_i = \mathbf{r}_i - \mathbf{r}_{i+1}$, respectively. The corresponding momenta are denoted by $\mathbf{P}_i$ and $\tilde{\mathbf{p}}_i$, respectively. Also, the dimer total mass and the reduced mass are given by $M_i = m_i + m_{i+1}$ and $\mu_i = \frac{m_i m_{i+1}}{m_i + m_{i+1}}$, respectively. This variable transformation is formally written as $(\mathbf{r}^N, \mathbf{p}^N) \rightarrow (\mathbf{R}^{N_d}, \mathbf{P}^{N_d}, \tilde{\mathbf{r}}^{N_d}, \tilde{\mathbf{p}}^{N_d}, \mathbf{r}^{N_m}, \mathbf{p}^{N_m})$.

By making this variable transformation, we can rewrite the Hamiltonian $H$ and the partition function $Z$ as:

$$ V^N = (W + F)^N = \sum_{k=0}^{N} \binom{N}{k} W^k \times F^{N-k}. \qquad (18) $$

<footnote n="4">⁴ There are $\binom{N}{2N_d} = \frac{N!}{(2N_d)! N_m!}$ ways to choose particles associated with dimers. Then there are $(2N_d - 1)!! = \frac{(2N_d)!}{2^{N_d} N_d!}$ ways to construct pairs among $2N_d$ particles. Thus $\Pi_{m,d} = \binom{N}{2N_d}(2N_d - 1)!!$, giving equation (17).</footnote>

<footnote n="5">⁵ To fix ideas by an example, consider the ‘wall pinning’ setting of section 2.2.2. The criterion to decide whether a particle belongs to the wall or to the fluid regions is $\mathbf{r}_i \in W$ and $\mathbf{r}_i \in F$, respectively. So the integration space for wall particles is $\mathbf{r}^W \in \Omega_W = W^k$ and that for fluid particles is $\mathbf{r}^F \in \Omega_F = F^{N-k}$. Of course, the union of all the spaces obtained by permuting the particle identities reconstructs the original space $\mathbf{r}^N \in V^N$, because of the trivial identity</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Illustration of the creation of virtual bonds and the associated change of coordinates. (a) Schematic plot for making virtual bonds. (b) The center of mass and relative position describe a dimer connecting particles *i* and *i* + 1.</figcaption></figure>

$$H(\mathbf{R}^{N_d}, \mathbf{P}^{N_d}, \tilde{\mathbf{r}}^{N_d}, \tilde{\mathbf{p}}^{N_d}, \mathbf{r}^{N_m}, \mathbf{p}^{N_m}) = \sum_{i \in \mathcal{D}} \left( \frac{\mathbf{P}_i^2}{2M_i} + \frac{\tilde{\mathbf{p}}_i^2}{2\mu_i} \right) + \sum_{i \in \mathcal{M}} \frac{\mathbf{P}_i^2}{2m_i} + U(\mathbf{R}^{N_d}, \tilde{\mathbf{r}}^{N_d}, \mathbf{r}^{N_m}), \quad (20)$$

$$Z = \Pi_{m,d} \int \left( \prod_{i \in \mathcal{D}} d\mathbf{R}_i d\mathbf{P}_i \right) \int_{\Omega_d} \left( \prod_{i \in \mathcal{D}} d\tilde{\mathbf{r}}_i d\tilde{\mathbf{p}}_i \right) \int_{\Omega_m} \left( \prod_{i \in \mathcal{M}} d\mathbf{r}_i d\mathbf{p}_i \right) \times \exp \left[ -\beta H(\mathbf{R}^{N_d}, \mathbf{P}^{N_d}, \tilde{\mathbf{r}}^{N_d}, \tilde{\mathbf{p}}^{N_d}, \mathbf{r}^{N_m}, \mathbf{p}^{N_m}) \right]. \quad (21)$$

Note that the integration spaces $\Omega_d$ and $\Omega_m$ in equation (21) are not exactly the same as the ones in equation (19), they are the images of those spaces under the change of variable. Yet, with a little abuse of notation, we keep the same notation for both.

We further proceed with the variable transformation for the relative movements of dimers by using spherical coordinates in $d=3$, which is formally written as $(\tilde{\mathbf{r}}^{N_d}, \tilde{\mathbf{p}}^{N_d}) \rightarrow (\tilde{\mathbf{r}}^{N_d}, \theta^{N_d}, \varphi^{N_d}, p_{\tilde{\mathbf{r}}}^{N_d}, p_{\theta}^{N_d}, p_{\varphi}^{N_d})$. We note that because of Liouville's theorem, the Jacobian is unity for this transformation, namely, $d\tilde{\mathbf{r}}_i d\tilde{\mathbf{p}}_i = d\tilde{\mathbf{r}}_i d\theta_i d\varphi_i dp_{\tilde{\mathbf{r}}_i} dp_{\theta_i} dp_{\varphi_i}$. Also, the kinetic part of the Hamiltonian can be written as:

$$\frac{\tilde{\mathbf{P}}_i^2}{2\mu_i} = \frac{p_{\tilde{\mathbf{r}}_i}^2}{2\mu_i} + \frac{p_{\theta_i}^2}{2I_i} + \frac{p_{\varphi_i}^2}{2I_i \sin^2 \theta_i}, \quad (22)$$

where $I_i = \mu_i \tilde{\mathbf{r}}_i^2$ is the moment of inertia, $p_{\tilde{\mathbf{r}}_i} = \mu_i \dot{\tilde{\mathbf{r}}}_i$, $p_{\theta_i} = I_i \dot{\theta}_i$, and $p_{\varphi_i} = I_i (\sin^2 \theta_i) \dot{\varphi}_i$ are the new momenta and the dot denotes a time derivative. We note that as one can see in equation (22), the kinetic term also depends on the coordinates. Thus, we cannot treat the momenta and positions separately, unlike in the random pinning setting.

We finally arrive at the expressions for the Hamiltonian $H$ and the partition function $Z$ that are suitable for our study:

$$H(\mathbf{R}^{N_d}, \mathbf{P}^{N_d}, \tilde{\mathbf{r}}^{N_d}, \theta^{N_d}, \varphi^{N_d}, p_{\tilde{\mathbf{r}}}^{N_d}, p_{\theta}^{N_d}, p_{\varphi}^{N_d}, \mathbf{r}^{N_m}, \mathbf{p}^{N_m}) = \sum_{i \in \mathcal{D}} \left( \frac{\mathbf{P}_i^2}{2M_i} + \frac{p_{\tilde{\mathbf{r}}_i}^2}{2\mu_i} + \frac{p_{\theta_i}^2}{2I_i} + \frac{p_{\varphi_i}^2}{2I_i \sin^2 \theta_i} \right) + \sum_{i \in \mathcal{M}} \frac{\mathbf{P}_i^2}{2m_i} + U(\mathbf{R}^{N_d}, \tilde{\mathbf{r}}^{N_d}, \theta^{N_d}, \varphi^{N_d}, \mathbf{r}^{N_m}), \quad (23)$$
$$Z = \Pi_{\text{m,d}} \int \left( \prod_{i \in \mathcal{D}} \mathrm{d}\mathbf{R}_i \mathrm{d}\mathbf{P}_i \right) \int_{\Omega_{\mathrm{d}}} \left( \prod_{i \in \mathcal{D}} \mathrm{d}\tilde{r}_i \mathrm{d}\theta_i \mathrm{d}\varphi_i \mathrm{d}p_{\tilde{r}_i} \mathrm{d}p_{\theta_i} \mathrm{d}p_{\varphi_i} \right) \int_{\Omega_{\mathrm{m}}} \left( \prod_{i \in \mathcal{M}} \mathrm{d}\mathbf{r}_i \mathrm{d}\mathbf{p}_i \right) \times \exp \left[ -\beta H \left( \mathbf{R}^{N_{\mathrm{d}}}, \mathbf{P}^{N_{\mathrm{d}}}, \tilde{r}^{N_{\mathrm{d}}}, \theta^{N_{\mathrm{d}}}, \varphi^{N_{\mathrm{d}}}, p_{\tilde{r}}^{N_{\mathrm{d}}}, p_{\theta}^{N_{\mathrm{d}}}, p_{\varphi}^{N_{\mathrm{d}}}, \mathbf{r}^{N_{\mathrm{m}}}, \mathbf{p}^{N_{\mathrm{m}}} \right) \right]. \quad (24)$$

We can then define the total (or annealed) average as:

$$\langle O \rangle = \frac{\Pi_{\text{m,d}}}{Z} \int \left( \prod_{i \in \mathcal{D}} \mathrm{d}\mathbf{R}_i \mathrm{d}\mathbf{P}_i \right) \int_{\Omega_{\mathrm{d}}} \left( \prod_{i \in \mathcal{D}} \mathrm{d}\tilde{r}_i \mathrm{d}\theta_i \mathrm{d}\varphi_i \mathrm{d}p_{\tilde{r}_i} \mathrm{d}p_{\theta_i} \mathrm{d}p_{\varphi_i} \right) \int_{\Omega_{\mathrm{m}}} \left( \prod_{i \in \mathcal{M}} \mathrm{d}\mathbf{r}_i \mathrm{d}\mathbf{p}_i \right) \times O(\mathbf{r}^N) \exp \left[ -\beta H \left( \mathbf{R}^{N_{\mathrm{d}}}, \mathbf{P}^{N_{\mathrm{d}}}, \tilde{r}^{N_{\mathrm{d}}}, \theta^{N_{\mathrm{d}}}, \varphi^{N_{\mathrm{d}}}, p_{\tilde{r}}^{N_{\mathrm{d}}}, p_{\theta}^{N_{\mathrm{d}}}, p_{\varphi}^{N_{\mathrm{d}}}, \mathbf{r}^{N_{\mathrm{m}}}, \mathbf{p}^{N_{\mathrm{m}}} \right) \right]. \quad (25)$$

We note again that up to this point, we have just discussed an exact variable transformation in the statistical mechanics expectation values. We did not modify the system itself at all.

From a dynamical point of view, similarly to what has been discussed at the end of section 2.2.2, one should imagine the following procedure: first, an equilibrium configuration of $\mathbf{r}^N$ is generated at time $t=0$. Second, the $N_d$ dimers are defined using the permutationally invariant procedure described in the first paragraph of this section. Third, some dynamics is run and the configuration of the system evolves in time. As long as the dimer and monomer variables remain in their respective domains$^6$, $\mathbf{r}^{2N_d} \in \Omega_d$ and $\mathbf{r}^{N_m} \in \Omega_m$, one can keep running the dynamics. If at some point one variable goes out of its domain, then one should run the procedure again to relabel the monomers and dimers according to the new particle positions. Once again, this is just a relabeling procedure that does not alter the system in any way. Doing this results in an annealed average as shown in section 2.2.2.

### 2.3.2. Random bonding.
Using the notations introduced in section 2.3.1, it is now straightforward to perform random bonding, that is, freeze the dimer bond lengths to define a quenched average. We consider the dimer lengths $\tilde{r}^{N_d}$ as being frozen, that is, $\tilde{r}^{N_d}$ plays the same role as the positions of the pinned particles in the random pinning approach, and of the frozen degrees of freedom $\mathbf{y}$ in the general discussion of section 2.1. We can thus consider the statistical mechanics of the remaining degrees of freedom, namely a composite of monomers and dimers for a particular realization of $\tilde{r}^{N_d}$, that play the role of the thermal degrees of freedom $\mathbf{x}$ in section 2.1. Note that the kinetic term associated to the frozen bond lengths is decoupled from all the other degrees of freedom. Hence, as in the particle pinning case, we can set the momenta $p_{\tilde{r}_i}$ to zero at the instant at which the bonds are frozen. Nevertheless, we include in the partition function the kinetic energy associated with the non-frozen variables that, for a given realization $\tilde{r}^{N_d}$, is then given by

<footnote>$^6$ The reader should keep in mind that while the dimer constraints are always the same, they are formally represented in different ways depending on the choice of coordinates, hence leading to formally different spaces $\Omega_d$ and $\Omega_m$.</footnote>
$$Z_{\tilde{r}}(\tilde{r}^{N_d}) = \int \left( \prod_{i \in \mathcal{D}} d\mathbf{R}_i d\mathbf{P}_i \right) \int \left( \prod_{i \in \mathcal{D}} d\theta_i d\varphi_i dp_{\tilde{r}_i} dp_{\theta_i} dp_{\varphi_i} \right) \int_{\Omega_m} \left( \prod_{i \in \mathcal{M}} dr_i d\mathbf{p}_i \right) \times \exp \left[ -\beta H \left( \mathbf{R}^{N_d}, \mathbf{P}^{N_d}, \tilde{r}^{N_d}, \theta^{N_d}, \varphi^{N_d}, p_{\tilde{r}}^{N_d}, p_{\theta}^{N_d}, p_{\varphi}^{N_d}, \mathbf{r}^{N_m}, \mathbf{p}^{N_m} \right) \right], \quad (26)$$

and the thermal average of an observable for a particular realization $\tilde{r}^{N_d}$ is defined by:

$$\langle \cdots \rangle_{\tilde{r}^{N_d}} = \frac{1}{Z_{\tilde{r}}(\tilde{r}^{N_d})} \int \left( \prod_{i \in \mathcal{D}} d\mathbf{R}_i d\mathbf{P}_i \right) \int \left( \prod_{i \in \mathcal{D}} d\theta_i d\varphi_i dp_{\tilde{r}_i} dp_{\theta_i} dp_{\varphi_i} \right) \int_{\Omega_m} \left( \prod_{i \in \mathcal{M}} dr_i d\mathbf{p}_i \right) \times \exp \left[ -\beta H \left( \mathbf{R}^{N_d}, \mathbf{P}^{N_d}, \tilde{r}^{N_d}, \theta^{N_d}, \varphi^{N_d}, p_{\tilde{r}}^{N_d}, p_{\theta}^{N_d}, p_{\varphi}^{N_d}, \mathbf{r}^{N_m}, \mathbf{p}^{N_m} \right) \right] (\cdots) . \quad (27)$$

Comparing this with equation (25), we can write:

$$\langle O \rangle = \frac{\Pi_{\text{m,d}}}{Z} \int_{\Omega_d} \left( \prod_{i \in \mathcal{D}} d\tilde{r}_i \right) \int \left( \prod_{i \in \mathcal{D}} d\mathbf{R}_i d\mathbf{P}_i \right) \int \left( \prod_{i \in \mathcal{D}} d\theta_i d\varphi_i dp_{\tilde{r}_i} dp_{\theta_i} dp_{\varphi_i} \right) \int_{\Omega_m} \left( \prod_{i \in \mathcal{M}} dr_i d\mathbf{p}_i \right) \times O(\mathbf{r}^N) \exp \left[ -\beta H \left( \mathbf{R}^{N_d}, \mathbf{P}^{N_d}, \tilde{r}^{N_d}, \theta^{N_d}, \varphi^{N_d}, p_{\tilde{r}}^{N_d}, p_{\theta}^{N_d}, p_{\varphi}^{N_d}, \mathbf{r}^{N_m}, \mathbf{p}^{N_m} \right) \right] \\ = \frac{\Pi_{\text{m,d}}}{Z} \int_{\Omega_d} d\tilde{r}^{N_d} Z_{\tilde{r}}(\tilde{r}^{N_d}) \langle O \rangle_{\tilde{r}^{N_d}} . \quad (28)$$

From this expression, we can deduce the probability distribution of the frozen variables $\tilde{r}^{N_d}$, which is given by:

$$\rho_f(\tilde{r}^{N_d}) = \Pi_{\text{m,d}} \frac{Z_{\tilde{r}}(\tilde{r}^{N_d})}{Z}, \quad (29)$$

and define a 'disorder' average over the realization of $\tilde{r}^{N_d}$ as:

$$\overline{\cdots} = \int_{\Omega_d} d\tilde{r}^{N_d} \rho_f(\tilde{r}^{N_d}) (\cdots). \quad (30)$$

Finally, as it can be expected from the discussion of section 2.1, we can obtain the identity that is at the basis of the random pinning procedure [58, 62]: The quenched average, i.e. the disorder average over the frozen degrees of freedom, taken after the thermal average over all remaining degrees of freedom conditioned on the frozen ones, corresponds exactly to the annealed average of the bulk system without bonding. By using equations (26)–(30), we get:

$$\langle O \rangle = \int_{\Omega_d} d\tilde{r}^{N_d} \rho_f(\tilde{r}^{N_d}) \langle O \rangle_{\tilde{r}^{N_d}} = \overline{\langle O \rangle_{\tilde{r}^{N_d}}}. \quad (31)$$

The random bonding procedure then goes as follows. First, one generates an equilibrium configuration of $N$ particles. Then, one performs the permutationally-invariant procedure described at the beginning of section 2.3.1 to identify the $N_d$ dimers. The bond lengths are finally frozen, that is taken as a quenched disorder, and one considers the statistical mechanics of the remaining degrees of freedom as thermal. With this
construction, all the hypotheses discussed in section 2.1 are fulfilled, and we can rig-
ously state that an equilibrium configuration of a system containing $N_d$ dimers and
$N_m$ monomers have been created. However, it should also be noted that the resulting
system is a rather unphysical one: due to the unbiased manner in which the choice of
the dimers has been made (based purely on the labels and not on a property of the
initial configuration), their length distribution will be proportional to the pair correl-
ation function of the initial system of monomers, and arbitrarily long dimers (up to
the system size) will be present$^7$. It is clear that in order to produce an equilibrium, or
nearly equilibrium, configuration of a more realistic system, a bias must be introduced
in the choice of the frozen dimers. The Nishimori (quiet freezing) condition (section 2.1)
will therefore not be strictly respected, and the consequences of this choice have to be
assessed. In the following section, we discuss several possible choices, and study two of
them numerically.

2.3.3. ‘Sorted bonding’ from the shortest to the largest bond. We start by considering a seemingly ‘natural’ procedure, which will, however, turn out to give unsatisfactory results. In this procedure, the dimers are created between those pairs of particles that have the smallest inter-particle distance. Specifically, we find the two particles with the smallest distance $|\tilde{\mathbf{r}}_{ij}| = |\mathbf{r}_i - \mathbf{r}_j|$ and we relabel them as $\mathbf{r}_1$ and $\mathbf{r}_2$ (which particle is 1 and which one is 2 is irrelevant). Then, we look at the remaining particles (excluding 1 and 2) and find once again the pair with the smallest inter-particle distance between particles, which we relabel as $\mathbf{r}_3$ and $\mathbf{r}_4$. We continue iterating this procedure until $N_d$ pairs have been sorted. This is a permutation-invariant procedure. We then make $N_d$ dimers by freezing the bond lengths.

We note that the configuration that has been created is *not* an equilibrium configuration of the systems of dimers and monomers interacting with the original Hamiltonian. To see this, one simply has to imagine the evolution of the system starting from the initial configuration just after freezing. In this initial configuration, all monomer-monomer distances or monomer-dimer distances are, by construction, larger than the largest dimer size. The dynamics will, obviously, not preserve this constraint and the ‘hole’ created by the freezing procedure in the (say) monomer-monomer correlation function, which is an intensive, self averaging quantity, will disappear progressively. The system shows the aging of a static observable, a hallmark of nonequilibrium dynamics.

A short reflection shows that the above procedure for creating dimers can in fact be
used to create an equilibrium monomer-dimer mixture if one allows a modification of
the interaction energies between the particles. This energy function should preserve the
following properties:

1. $\Omega_d$ is defined by the constraint (as above, the tilde notation indicates a distance between a pair of particles, rather than an absolute position):

$$|\tilde{\mathbf{r}}_{12}| < |\tilde{\mathbf{r}}_{34}| < |\tilde{\mathbf{r}}_{56}| < \cdots < |\tilde{\mathbf{r}}_{2N_d-1,2N_d}| = R_{\max}, \quad (32)$$

<footnote>$^7$ The case in which such long range interactions are soft has been considered in [74, 75] and it was found that such a system does show a strong slowing down of the dynamics when the number of links is increased.</footnote>
that is, the integration over the $N_d$ dimers (first $2N_d$ particles) is constrained in such a way that each dimer has a larger inter-particle distance than the previous one. Here, $R_{max}$ is the maximal bond length of the dimers.

2. Additional constraints have to be imposed in $\Omega_d$ on the distances between particles belonging to distinct dimers. For example,

$$|\tilde{\mathbf{r}}_{12}| < |\tilde{\mathbf{r}}_{13}|, \quad |\tilde{\mathbf{r}}_{12}| < |\tilde{\mathbf{r}}_{14}|, \qquad (33)$$

in order to guarantee that $|\tilde{\mathbf{r}}_{12}|$ remains the shortest distance. Similar constraints are needed for particles with higher labels.

3. The integration over the remaining $N_m$ monomers is over the space $\Omega_m$ such that all of the $N_m(N_m - 1)/2$ monomer-monomer distances and each of the $N_m N_d$ monomer-dimer distances is larger than $R_{max}$. This must be true because, otherwise, in the sorting procedure, we would have chosen one of the monomers to belong to a dimer. The positions of the monomers are then constrained to the space $\Omega_m$, that is, the smallest monomer-monomer and monomer-dimer distances must stay larger than the largest dimer bond length $R_{max}$.

To perform the quenched average dynamically, the above constraints must be implemented either by rejecting moves that violate them in a MC simulation or by adding a hard wall term that would reflect the relative velocity of two monomers or a monomer-dimer pair if they reach the minimal distance in a Molecular Dynamics simulation. The time-average under this constraint would then result in a proper quenched average over the $\Omega_m$ space.

The positive aspect of this construction is that we can construct permutationally invariant observables $O(r^N)$ that only depend on the unfrozen degrees of freedom. However, if the number of bonds is large, such that $R_{max}$ reaches the first peak of the radial distribution of monomer distances $g(r)$, the resulting constraints on the monomers will be quite strong, and the corresponding dynamics becomes unphysical, because the interactions between the particles are non-zero even for arbitrarily large distances.

2.3.4. 'Random bonding' within a cutoff. We now consider a different procedure, close to the rigorous 'random bonding' described at the end of section 2.3.2, but now imposing that the length of the formed dimers does not exceed a maximum value $R_b$ (figure 3(a)). This method was used in [44], where it was assumed to produce an equilibrium configuration of the dimer-monomer mixture.

The procedure is as follows: first, we select at random the first particle in each of the $N_d$ dimer. Second, for each of these, we choose a particle at random among those at distance smaller than $R_b$ from it to be its partner in the dimer. If there is no viable partner, we discard that candidate and choose at random a new one. This procedure is also permutation invariant with the following properties:

1. The space $\Omega_d$ is defined by the only constraint that the dimer bond lengths $\tilde{r}_i \le R_b$, which is trivially satisfied if the bond lengths are frozen. For the rest, dimers are free to translate and rotate.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Schematic illustration of the construction of a bond connecting particles $i$ and $i+1$. (a) Random bonding with spherical cut-off. The sphere at the center of the particle $i$ with the radius $R_b$ defines its neighborhood. (b) Random bonding with directional alignment. An additional constraint is imposed by a cut-off $P_b = \cos\theta_b$ associated with directional alignment.</figcaption></figure>

2. The monomer space $\Omega_m$ is, seemingly, unconstrained because in the dimer construction procedure, nothing is implied about the monomers.

However, the choice of the $N_d$ dimers is now not only based on the labels of the particles but depends on the properties of the equilibrated configuration of the system of monomers so that the Nishimori condition is, again, not strictly respected. The distribution of the frozen variables (bond lengths) is not strictly the one in an equilibrium, unconstrained system, but will converge to it as $R_b$ is increased. As a result, the initial configuration will also display subtle correlations between non-bonded particles, which will decay with time if the system is propagated with the original interaction potential. To see this, let us consider the total pair correlation function of the system just after the freezing of the bonds, $g(r)$. This correlation function can be written as:

$$g(r) = g_{\text{inter}}(r) + g_{\text{intra}}(r), \qquad (34)$$

where $g_{\text{inter}}(r)$ is the contribution of pairs of particles is not bonded together by a frozen bond, and $g_{\text{intra}}(r)$ is the contribution of the pairs of particles that are connected. By construction, $g_{\text{intra}}(r)$ will have a discontinuity at $R_b$, above which it will jump to zero. In contrast, $g(r)$ is the pair correlation of an equilibrium system of monomers and is continuous. We conclude that $g_{\text{inter}}(r)$ will initially have a discontinuity at $R_b$, which will not be preserved by the dynamics. The situation is similar to the one described above in the 'sorted bonding' case, but more subtle. In the following numerical study, we will see that the nonequilibrium character of the initial configuration is extremely small and disappears very rapidly. The choice of this bond distribution, for $R_b$ of the order of the interparticle distance, leads to a physically reasonable system which is nearly at equilibrium and can be considered to be close to optimal.
2.3.5. Random bonding with directional alignment. As a final example, we will consider a situation in which it is intuitively more obvious that the initial configuration is out of equilibrium and will display aging over a measurable time scale. The bonding scheme is similar to the random bonding within a cutoff $R_b$ discussed in the previous section; however, an additional constraint is imposed on the orientation of the frozen dimers, which are restricted to an angular sector $\theta_b$ around the z-direction figure 3(b)). The initial configuration therefore displays nematic order of the dimers, which obviously will not be preserved by the dynamics. The numerical study will allow us to estimate the persistence time of this nonequilibrium feature compared to the isotropic case.

# 3. Simulation methods

## 3.1. Model

We employ the Kob-Andersen binary mixture [66], in which particles interact through the Lennard-Jones pair potential:

$$u_{\alpha\beta}(r) = 4\epsilon_{\alpha\beta}\left[\left(\frac{\sigma_{\alpha\beta}}{r}\right)^{12} - \left(\frac{\sigma_{\alpha\beta}}{r}\right)^{6}\right], \quad (35)$$

where $\alpha, \beta \in \{A, B\}$ are species indexes. Both species have the same mass, which is set to $m=1$. The values of the parameters $\sigma_{\alpha\beta}$ and $\epsilon_{\alpha\beta}$ are given in [66]. The units of length and energy are set by the parameters $\sigma = \sigma_{AA} = 1$ and $\epsilon = \epsilon_{AA} = 1$, respectively, and we use the Boltzmann constant $k_B = 1$. The potentials are cut and shifted at a distance $2.5\sigma_{\alpha\beta}$. We simulate systems comprising $N$ particles in a cubic box of side $L$ with periodic boundary conditions at a number density $\rho = N/V = 1.2$. We use the system size $N = 1200$. We perform constrained molecular dynamics simulations via the RATTLE algorithm [76] with a simple Nosé-Hoover thermostat [44].

## 3.2. Making randomly bonded systems

Starting from an equilibrium configuration of the original (bulk) KA model with $N$ particles described above, we generate a randomly bonded system comprising monomers and dimers. We consider the following two protocols to do so.

(1) Random bonding with spherical cut-off: First, we choose a particle randomly, say particle $i$. We then randomly pick another particle $j$ among the neighboring particles of particle $i$, located inside a sphere with a cut-off radius $R_b$ and which is not yet bonded. We then relabel the particle $j$ as $i+1$. This process is schematically shown in figure 3(a). We set $R_b = 1.5$, which is near the first minimum of the radial distribution function, thus corresponding roughly to the boundary of the first coordination shell. We then permanently freeze the distance between the two particles, $\tilde{r}_i = |\mathbf{r}_i - \mathbf{r}_{i+1}|$, which means that the particles $i$ and $i+1$ now form a dimer with a rigid body constraint. We repeat the above process for the remaining monomer particles until the number of dimers, $N_d$, reaches the target value. By construction, we have $N = N_m + 2N_d$, where $N_m$ is the number of monomers. We introduce the control parameter $c = \frac{2N_d}{N} = \frac{N - N_m}{N}$, such that $c=0$ corresponds to a system with only monomers (hence the original bulk
Creating equilibrium glassy states via random particle bonding

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Examples of initial bonded configurations, here for $T = 0.6$, $c = 0.5$, and $R_b = 1.5$. Left: $P_b = 0$ (random bonding with spherical cut-off). Right: $P_b = 0.9$ (random bonding with directional alignment).</figcaption></figure>

model), whereas $c = 1$ corresponds to a system having only dimers. Using the algorithm explained above, it is difficult in practice to reach $c = 1$, because at some point one runs out of neighboring pairs, leaving a few percentage of monomer particles. Thus, in the present work, we use $c = 0.95$ as the maximum value.

(2) Random bonding with directional alignment: In this protocol, we aim to prepare an initial bonded configuration such that the orientation of the dumbbell molecules tends to align along the unit vector of the z-direction, $e_z = (0,0,1)$. We define a unit vector for the orientation of each dumbbell molecule composed of particles $i$ and $i+1$ by $\mathbf{n}_i = \tilde{\mathbf{r}}_i/|\tilde{\mathbf{r}}_i|$, where $\tilde{\mathbf{r}}_i = (x_i - x_{i+1}, y_i - y_{i+1}, z_i - z_{i+1})$. In protocol (1) we make a bond between the $i$th and $(i+1)$th particles randomly among the neighbor particles inside a spherical shell with the radius $R_b$ centered at the position of the $i$th particle. We now impose a further limitation on the neighboring region by considering the orientation of the molecule. To this end, we define a magnitude of alignment (or polarization) for each molecule, $P_i$, given by

$$P_i = |e_z \cdot \mathbf{n}_i| = |z_i - z_{i+1}|/|\tilde{\mathbf{r}}_i| = |\cos \theta_i| \quad (36)$$

We then introduce another cut-off threshold $P_b$ such that particles having $P_i > P_b = \cos \theta_b$ can be considered as candidate neighbors for bonding. This new neighborhood region is schematically shown in figure 3(b). As in protocol (1), we repeat the above process for the remaining monomer particles until either the number of dimers reaches the target value or we run out of candidate pairs for bonding. In practice, we set $c = 0.75$ as the maximum value for this protocol with $P_b = 0.9$.

In figure 4, we show initial bonded configurations for $T = 0.6$, $c = 0.5$, and $R_b = 1.5$. The left panel is $P_b = 0.0$ (random bonding with spherical cut-off), and the right panel is $P_b = 0.9$ (random bonding with directional alignment).
Creating equilibrium glassy states via random particle bonding

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Static properties in the random bonding with spherical cut-off. (a) Radial distribution function $g(r)$ for the original system with $c=0$ at $T=0.6$ (black curve), and bonded systems generated by the spherical cut-off protocol ($P_b=0$) with $c=0.95$ at $T=0.6$ measured in time windows $t \in [0, 10]$ (red curve) and $t \in [10 000, 20 000]$ (blue curve). The arrow indicates the location of the cut-off, $R_b=1.5$. The inset shows a zoomed plot near $R_b$. (b) Decomposition of $g(r)$ into the intramolecule contribution $g_{\text{intra}}(r)$ and intermolecule contribution $g_{\text{inter}}(r)$ for $c=0.95$ at $T=0.6$ in the time window $t \in [0, 10]$. (c) Orientational correlation function $g_{\text{ori}}(r)$ for the bonded systems presented in (a), measured in two different time windows.</figcaption></figure>

# 4. Structural properties

In this section, we study how the bonding process affects the equilibrium properties of the system in terms of the static structure.

## 4.1. Radial distribution function

We first characterize the static structure by the radial distribution function $g(r)$ for all particles, which is defined by:

$$g(r) = \frac{L^3}{4\pi r^2 \Delta r N (N-1)} \overline{\left\langle \sum_{\substack{i,j \\ (i \neq j)}} \delta(r - |\mathbf{r}_i - \mathbf{r}_j|) \right\rangle}_{\tilde{r}^{N_d}}, \quad (37)$$

where $\Delta r$ is the bin size for computation. We use $\Delta r = 0.0125$. Figure 5(a) shows $g(r)$ for the original system ($c=0$) in equilibrium and a bonded system ($c=0.95$) with $P_b=0$ measured in a short time window ($t \in [0, 10]$) starting right after bonding ($t=0$) and a longer time window ($t \in [10 000, 20 000]$) starting after a waiting time $t_w=10 000$. The former time scale corresponds to a vibrational one, whereas the latter timescale corresponds to the time scale for escaping from the cage (figure 7(a)). We note that the first and second peaks correspond to the nearest neighbor contacts for A-B and A-A pairs, respectively. Overall, the three $g(r)$'s superimposed very well. However, around the cut-off distance $R_b=1.5$, there are tiny but distinct differences between the original and bonded system, as emphasized in the inset. This means that the bonding process disturbs the initial equilibrium state, and the system enters an out-of-equilibrium state and thus will age. Yet this aging process is very fast, within the
vibrational timescale, as can be recognized from the fact that in the inset the curves for $t \in [0, 10]$ and $t \in [10 000, 20 000]$ superimpose very well, i.e. the system ends up in a new stationary state very quickly.

In order to better understand the differences between the bulk and bonded systems around $r = R_b$, we decompose $g(r)$ into the contributions from the bonded pair (intra molecule contribution, $g_{\text{intra}}(r)$) and the rest (inter molecule contribution, $g_{\text{inter}}(r)$), following the discussion in section 2.3.4. By construction, $g(r) = g_{\text{intra}}(r) + g_{\text{inter}}(r)$. In figure 5(b), we present $g(r)$, $g_{\text{intra}}(r)$, and $g_{\text{inter}}(r)$ for $c = 0.95$ in the time window $t \in [0, 10]$. This graph shows that $g_{\text{intra}}(r)$ is only a small contribution to the total radial distribution and, as expected, vanishes beyond $R_b$. Since the bonding we employ in this work is a rigid body constraint, $g_{\text{intra}}(r)$ is completely frozen at time $t = 0$, and it never changes during the simulation for any $t > 0$. We also find that $g_{\text{inter}}(r)$, at strictly $t = 0$, has a small dip at $r = R_b$ (not shown), compensating the sharp edge in $g_{\text{intra}}(r)$, such that the total $g(r)$ is smooth as in the original system. At $t > 0$, only $g_{\text{inter}}(r)$ evolves with time, which removes the dip around $r = R_b$ and ends up in a smooth curve as shown in figure 5(b). Consequently, the total contribution, $g(r) = g_{\text{intra}}(r) + g_{\text{inter}}(r)$, has a bump around $r = R_b$ at $t > 0$. Indeed, this out-of-equilibrium effect originates from the bonding process with a cut-off $R_b$, but we show that it is a very minor perturbation of the whole system even for a very high value of the dimer concentration $c = 0.95$.

## 4.2. Orientational correlation function

The radial distribution function $g(r)$ considers only the particle positions irrespective of molecular orientations. To characterize the configuration of molecules in more detail, we compute an orientational correlation function given by:

$$g_{\text{ori}}(r) = \frac{\left\langle \sum_{\substack{i,j \in D \\ (i \neq j)}} |\mathbf{n}_i \cdot \mathbf{n}_j| \delta(r - |\mathbf{R}_i - \mathbf{R}_j|) \right\rangle_{\tilde{r}^{N_d}}}{\left\langle \sum_{\substack{i,j \in D \\ (i \neq j)}} \delta(r - |\mathbf{R}_i - \mathbf{R}_j|) \right\rangle_{\tilde{r}^{N_d}}}, \quad (38)$$

where $\mathbf{R}_i$ is the center of mass of molecule composed of the $i$th and $(i+1)$th particles. In figure 5(c), we present $g_{\text{ori}}(r)$ for shorter and longer time windows, respectively. We find that the two curves are superimposed very well. No aging is thus detected right after a quick relaxation within a vibrational timescale in terms of molecular orientational correlations.

## 4.3. Random bonding with directional alignment

The above simulation results suggest that although it is not strictly in equilibrium, the bonding protocol with a spherical cut-off produces a nearly equilibrium state, and its aging process is very fast. We now test whether this is also the case for an arbitrary bonding protocol. To this end, we study the bonding protocol with directional alignment illustrated in figure 3(b) (see also the discussion in section 2.3.5). In figure 6(a), we show $g(r)$ for the original system ($c = 0$) in equilibrium and a bonded system ($c = 0.75$) with $P_b = 0.9$ measured in a shorter time window ($t \in [0, 10]$) starting right after the bonding
Creating equilibrium glassy states via random particle bonding

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Random bonding with directional alignment. (a) Radial distribution function $g(r)$ for the original system with $c=0$ at $T=0.6$ (black curve), bonded systems generated by the spherical and directional cut-off protocol ($P_b=0.9$) with $c=0.75$ measured in time windows $t \in [0, 10]$ (red curve) and $t \in [10 000, 20 000]$ (blue curve). The arrow indicates the location of the spherical cut-off, $R_b=1.5$. The inset shows a zoomed plot near $R_b$. (b) Orientational correlation function $g_{ori}(r)$ for the bonded systems presented in (a), measured in two different time windows. (c) Time evolution of the total polarization $P_{tot}$ for the spherical cut-off ($c=0.95$ and $T=0.6$ with $P_b=0$) and directional alignment ($c=0.75$ and $T=0.6$ with $P_b=0.9$) protocols. Three representative trajectories are shown for each protocol.</figcaption></figure>

$(t=0)$ and a longer time window $(t \in [10 000, 20 000])$ after a waiting time of $t_w=10 000$. Similar to the spherical cut-off case, all $g(r)$'s superimpose well and the tiny discrepancy can be recognized only in the zoomed-in plot near $R_b$. However, a strong aging effect can be observed in $g_{ori}(r)$ presented in figure 6(b). At shorter time scale, $g_{ori}(r)$ is overall larger because of the initially aligned molecular configuration. It then decays with time, as expected from the fact that randomly oriented configurations are entropically more favored. At longer timescale, $g_{ori}(r)$ converges to $g_{ori}(r) \approx 0.5$ at $r \gg 1$ and is overall similar to that of figure 5(c) for a spherically bonded system.

This aging process can be directly quantified by the total molecule polarization, given by:

$$P_{tot}(t) = \frac{1}{N_d} \sum_{i \in D} P_i(t) \qquad (39)$$

Figure 6(c) shows the time evolution of $P_{tot}$. The system with $P_b=0$ has random molecular orientations, producing $P_{tot} \approx 0.5$ during the entire simulation. In contrast, the system with $P_b=0.9$ takes a higher value, $P_{tot} \approx 0.95$ at $t=0$, and it decays slowly with time, demonstrating a strong aging effect. At a longer time, it then reaches $P_{tot} \approx 0.5$.

## 5. Dynamical properties

In the previous section, we have confirmed that the random bonding protocol with a spherical cut-off produces a near equilibrium configuration right after bonding, and the subsequent aging process lasts only a short period of time. Thus, in practice, this procedure allows us to study equilibrium dynamical properties on a longer time scale. In the present section, we hence study the (equilibrium) dynamical properties of randomly bonded glass-forming liquids generated by the spherical cut-off protocol.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Dynamical correlations for dimers and monomers. (a) Self-intermediate scattering function $F_s(q,t)$ for $T=0.6$ and $c=0.95$ for monomers, dimers, and all particles. Dashed and solid curves indicate $F_s(q,t)$ computed from trajectories with the waiting time $t_w=0$ and $t_w=10^5$, respectively. (b) Relaxation time $\tau_\alpha$ versus the inverse of temperature $1/T$ for $c=0.5$.</figcaption></figure>

## 5.1. Intermediate scattering functions

To characterize the dynamic properties of the system, we compute the self part of the intermediate scattering functions for monomers, the center of mass of dimers, and all particles, given by:

$$F_s^{\text{Mono}}(q,t) = \left\langle \frac{1}{N_m} \sum_{j \in \mathcal{M}} e^{-i \mathbf{q} \cdot (\mathbf{r}_j(t) - \mathbf{r}_j(0))} \right\rangle_{\tilde{r}^{N_d}}, \quad (40)$$

$$F_s^{\text{Di}}(q,t) = \left\langle \frac{1}{N_d} \sum_{j \in \mathcal{D}} e^{-i \mathbf{q} \cdot (\mathbf{R}_j(t) - \mathbf{R}_j(0))} \right\rangle_{\tilde{r}^{N_d}}, \quad (41)$$

$$F_s^{\text{All}}(q,t) = \left\langle \frac{1}{N} \sum_{j=1}^{N} e^{-i \mathbf{q} \cdot (\mathbf{r}_j(t) - \mathbf{r}_j(0))} \right\rangle_{\tilde{r}^{N_d}}, \quad (42)$$

respectively. We have averaged over 5–20 different realizations to calculate these time correlation functions. The wave-vector $q$ is chosen to be $q=7.25$, and the location of the main peak in the static structure factor [66].

Figure 7(a) shows these intermediate scattering functions at $T=0.6$ and $c=0.95$. $F_s^{\text{Di}}(q,t)$ displays a higher value of the plateau than the one found in $F_s^{\text{Mono}}(q,t)$ or $F_s^{\text{All}}(q,t)$, which is reasonable because the effective cage size of the dimers is smaller than that of the monomers. More important is the observation that these functions relax on essentially the same timescale, which proves that the monomers and dimers have a very similar relaxation dynamics in terms of translational motions. We define the relaxation time $\tau_\alpha$ as the time at which the intermediate scattering function decays to $1/e$ and present the $\tau_\alpha$ versus $1/T$ plot for $c=0.5$ in figure 7(b). The relaxation time for the center of mass of dimers exceeds $\tau_\alpha$ for the monomers by a factor around
Creating equilibrium glassy states via random particle bonding

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Intermediate scattering function for various temperatures and dimer concentrations. $F_s(q, t)$ for all particles for several values of $c$ for $T=2.0$ (a), $T=0.8$ (b), $T=0.6$ (c), and $T=0.424$ (d).</figcaption></figure>

2, independent of $T$, demonstrating that the dynamics of the two types of particles stays coupled in the whole accessed $T$-range, and we have checked that this is also the case for the other values of $c$. Hence one can conclude that the three definitions of the intermediate scattering functions provide, essentially, the same information in terms of structural relaxation. Therefore, we focus on the following on $F_s^{\text{All}}(q, t)$ and drop the superscript unless otherwise specified.

In figure 7(a), we also include data for a different waiting time $t_w$, that is, the time interval between the initial configuration ($t=0$), and the time when we start measuring the correlation ($t=t_w$). The graph demonstrates that there is no detectable waiting-time dependence on the time scale of structural relaxation, as expected from the structural analysis presented in the previous section. Thus, the random bonding (using the spherical cut-off protocol) allows us to probe equilibrium dynamics relevant to structural relaxation in simulations right after bonding, akin to random pinning.

The temperature and $c$-dependence of $F_s(q, t)$ is presented in figure 8. We find that the influence of bonding is very small if the temperature is high, $T=2.0$, panel (a). However, once $T$ is decreased, panels (b)–(d), the bonding affects the dynamics very strongly, akin to the behavior of randomly pinned systems [38, 41]. For example, at the lowest temperature $T=0.424$, one can see that for $c > 0.5$ the dynamics is completely frozen on the timescale of our simulation, demonstrating that the random bonding
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Relaxation time as a function of temperature *T* and dimer concentration *c*. We report iso-τ<sub>α</sub> curves in the *T* − *c* plane, obtained from the dynamical correlations shown in figure 8.</figcaption></figure>

allows to access an extremely slow glassy dynamics in almost equilibrium. By using the data for τα (figure 11(a)), we estimate that the bonded system at T = 0.424 and c = 0.95 (that are prepared by forming bonds from the original T = 0.424 configurations) has an equilibrium relaxation time τα ≈ 10¹², about a factor of 10⁷ larger than the largest τα accessed in our simulations [44]. This demonstrates that the random bonding protocol indeed provides us with a huge gain in terms of computational time for the preparation of the initial equilibrium state.

To observe the influence of *T* and *c* together we present in figure 9 the iso-τ<sub>α</sub> curves
in the *T* − *c* plane. These curves increase with increasing *c*, which is again similar to the
results found for randomly pinned systems [29, 30]. The shape of the curves depends
only mildly on τ<sub>α</sub>, hinting at a simple functional relation among τ<sub>α</sub>, *T*, and *c*. This
point will be investigated in more detail below.

Finally, we present a comparison between the self and collective parts of the inter-
mediate scattering functions. Previous studies have reported that in randomly pinned
fluids the collective part shows apparent freezing while the self part did not render the
analysis of the system dynamics difficult [30, 37, 38]. In figure 10(a), we show that the
self- and collective parts for the randomly-bonded glass formers at T = 0.6 for different
values of c. One recognizes that, in contrast to the pinned systems, the collective part
also relaxes to zero and that the relaxation time is slightly larger than the one for the self
part, at least for the wave vector considered. Interestingly, however, the ratio between
the two timescales is about a factor two irrespectively of c, as shown in figure 10(b).
This suggests that the self and collective correlators do not decouple and therefore we
can conclude that the self part gives reliable dynamical information about the system.
Creating equilibrium glassy states via random particle bonding

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Comparison of the self and collective dynamics. (a) Self (dashed curves) and collective (solid curves) intermediate scattering functions at $T = 0.6$ for different $c$. (b) Relaxation time $\tau_{\alpha}$ versus $c$, computed from the data in (a).</figcaption></figure>

## 5.2. Scaling of the dynamics and fragility

In the following, we study how the relaxation time $\tau_{\alpha}$ depends on $T$ and $c$ and in particular we examine the scaling of the dynamics and kinetic fragility of the system.

Figure 11(a) presents an Arrhenius plot of the structural relaxation time $\tau_{\alpha}(c, T)$ for different $c$ and one recognizes that for all investigated values of $c$ there is a pronounced non-Arrhenius dependence. Interestingly, we find in figure 11(b) that the different curves can collapse onto a master curve by introducing the scaled temperature $T/\mu(c)$, i.e. $\tau_{\alpha}(c, T) = \tau_{\alpha}(c=0, T/\mu(c))$, where $\mu(c)$ is an unknown function. The inset of figure 11(b) shows $\mu(c)$, which for each $c$ has been determined manually such that $\tau_{\alpha}(c, T)$ superimposes with $\tau_{\alpha}(c=0, T)$. One sees that to a first approximation $\mu(c)$ is linear, but a slight upward bending can be noticed. The existence of a master curve is so far just an empirical observation, and it holds at least for the $T$- and $c$-range probed by our simulations. The presence of this scaling behavior suggests that all relevant temperature scales, such as the mode-coupling crossover $T_{mct}$ [77] and the Kauzmann transition temperature $T_K$ [78] (if it exists) for the original system ($c=0$), are just scaled by $\mu(c) \ge 1$. This is in contrast to the random pinning case, where the value of $T_K$ as determined from the dynamics is essentially unchanged (or slightly decreases) with increasing the pinning density [40]. (Note however, that instead the value of $T_K$ as estimated from the thermodynamics increases [30, 31].)

We now discuss the consequences of the observed empirical scaling in terms of kinetic fragility [79, 80]. The kinetic fragility $m_F(c)$ (for a given $c$) can be defined by

$$m_F(c) = \frac{\partial \log_{10} \tau_{\alpha}(c, T)}{\partial (T_g(c)/T)} \Big|_{T=T_g(c)}, \quad (43)$$

where $T_g(c)$ is a kinetic glass transition temperature. Conventionally, this temperature is defined by $\tau_{\alpha}(c, T_g(c)) = \tau_g$, where $\tau_g$ is a typical experimental timescale, e.g. 100 s. The dynamical scaling implies that $\tau_g = \tau_{\alpha}(c, T_g(c)) = \tau_{\alpha}(c=0, T_g(c)/\mu(c)) = \tau_{\alpha}(c=0, T_g(c=0))$, and hence $T_g(c) = \mu(c)T_g(c=0)$. These relations and equation (43) give:
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Dynamical scaling of the relaxation time with dimer concentration. (a) Relaxation time $\tau_{\alpha}(c, T)$ obtained from the self-intermediate scattering function for all particles. (b) The same data using a normalized abscissa, $\mu(c)/T$. The inset shows $\mu(c)$ versus $c$.</figcaption></figure>

$$m_F(c) = \left. \frac{\partial \log_{10} \tau_{\alpha}(c=0, T/\mu(c))}{\partial(\mu(c)T_g(c=0)/T)} \right|_{T=\mu(c)T_g(c=0)} = \left. \frac{\partial \log_{10} \tau_{\alpha}(c=0, T')}{\partial(T_g(c=0)/T')} \right|_{T'=T_g(c=0)} = m_F(c=0), \quad (44)$$

irrespective of $c$. This means that the fragility of randomly-bonded glass-forming liquids is independent of $c$. This is in stark contrast to randomly pinned systems [38–40], where the fragility decreases with increasing the pinning density.

From these observations, one concludes that the dynamics of randomly bonded glass-forming liquids is qualitatively quite different from that of the randomly pinned systems also in terms of the temperature dependence of structural relaxation.

The dynamical scaling shown in figure 11(b) reminds us of simple systems where structure and dynamics are invariant to a good approximation along isomorphs in the phase diagram [81, 82]. In these systems, the dynamics at different state points can be rescaled by a uniform scaling of space and time. In contrast, the randomly bonded systems introduce strong constraints on the system (in the form of bonds) that alter dynamical relaxation processes significantly as $c$ is increased. The origin of the observed empirical scaling must then be different from the isomorph invariance. Figure 12 shows $F_s(q, t)$ having similar $T/\mu(c)$. Although the relaxation time $\tau_{\alpha}$ is similar, there is a trend that the correlators with higher $c$ show a lower plateau compared to those with smaller $c$. This suggests that the observed scaling collapse cannot be understood as a simple uniform space-time rescaling. Further investigations are needed to understand the origin of the empirical dynamical scaling.

## 6. Dynamical heterogeneity

The dynamics of glassy liquids is accompanied by strong dynamical heterogeneities, the intensity of which grows with decreasing temperature [83–86]. Since our bonding procedure allows to generate configurations in the deeply supercooled regime, we
Creating equilibrium glassy states via random particle bonding

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Comparison of dynamical correlations at corresponding state points under dynamical scaling. $F_s(q, t)$ for different state points with similar $T/\mu(c)$. (a) $T/\mu(c) \approx 0.55$. (b) $T/\mu(c) \approx 0.42$.</figcaption></figure>

can thus access these heterogeneities in the randomly bonded glass-forming liquids in (nearly) thermal equilibrium. Since bonded dimer systems involve rotational motion as an additional relaxation channel, we consider dynamical heterogeneities not only for the positional degrees of freedom but also the rotational ones.

## 6.1. Positional degrees of freedom

First, we compute the standard four-point correlation function $\chi_4^Q(t)$ associated with the positional degrees of freedom, which is given by:

$$ \chi_4^Q(t) = N \left( \langle \hat{Q}^2(t) \rangle_{\hat{r}^{N_d}} - \langle \hat{Q}(t) \rangle_{\hat{r}^{N_d}}^2 \right), \quad (45) $$

where $\hat{Q}(t) = \frac{1}{N} \sum_{i=1}^N \theta(a - |\mathbf{r}_i(t) - \mathbf{r}_i(0)|)$ is an overlap function taking into account all particles and $\theta(x)$ is the Heaviside step function [87]. We set the distance $a$ to the often used value $a=0.3$. We note that $\chi_4^Q(t)$ defined in equation (45) does not contain contributions from sample-to-sample fluctuations associated with different realizations of bonds [42]. Figure 13 shows $\chi_4^Q(t)$ for different values of $c$ at a constant temperature, panel (a) and with decreasing $T$ at a constant $c$, panel (b). We find systematic growth with increasing glassiness in both cases, which is in contrast to randomly pinned particle systems. It has been reported that the four-point correlation function of randomly pinned particle systems does not grow systematically or decreases approaching the glass transition, whereas the relaxation time increases significantly [39, 41–43]. This difference in the dynamical behavior is directly related to the fact that in pinned systems the size of the dynamical heterogeneities is hindered by the presence of the pinned particles, whereas, in the present system, the heterogeneities can grow unhindered.
Creating equilibrium glassy states via random particle bonding

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Dynamical heterogeneity of positional degrees of freedom. (a) Four-point correlation function $\chi_4^Q(t)$ associated with positional degrees of freedom computed from the overlap function $\hat{Q}$ at a constant temperature ($T = 0.65$), panel (a), and at different $T$ for constant dimer concentration $c = 0.950$, panel (b).</figcaption></figure>

## 6.2. Rotational degrees of freedom

We next consider the rotational degrees of freedom for the dimers. We first compute the average dynamics with the mean rotational angle [88], which is given by

$$ \Phi(t) = \left\langle \frac{1}{N_d} \sum_{i \in \mathcal{D}} \varphi_i(t) \right\rangle_{\hat{r}^{N_d}}, \qquad (46) $$

where $\varphi_i(t) = \arccos(\mathbf{n}_i(t) \cdot \mathbf{n}_i(0))$. Figure 14(a) shows $\Phi(t)$ for different $c$ at constant temperature. One finds that the correlator has a two-step relaxation with a plateau on an intermediate timescale, akin to the mean-squared displacement. At sufficiently long times, the correlator approaches the asymptotic value $\pi/2$, which is expected when $\mathbf{n}_i(t = 0)$ and $\mathbf{n}_i(t \to \infty)$ are uncorrelated. In short, $\Phi(t)$ can separate vibrational motion and structural rearrangement in terms of rotational relaxation. In addition, the separation becomes more distinct when $c$ is increased, which is also in qualitative agreement with the mean-squared displacement.

We can now define a mean overlap function associated with the rotational motion via:

$$ R(t) = \left\langle \frac{1}{N_d} \sum_{i \in \mathcal{D}} \theta(b - \varphi_i(t)) \right\rangle_{\hat{r}^{N_d}}, \qquad (47) $$

where $b$ is a threshold separating vibrational motion and structural relaxation. In practice, we have chosen the value $b = 0.2$, see figure 14(a). Figure 14(b) shows that $R(t)$ presents a two-step relaxation, similar to $F_s(q, t)$. This correlation function allows to compute a characteristic timescale for rotational relaxation, $\tau_R$, defined by $R(\tau_R) = 0.3$, which is presented in the inset, together with $\tau_\alpha$ obtained by $F_s(q, \tau_\alpha) = 1/e$. Both $\tau_R$ and $\tau_\alpha$ track each other very well, particularly at larger $c$, suggesting that positional and
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Dynamical correlations of the orientational degrees of freedom. (a) and (c): mean rotational angle $\Phi(t)$ for (a) a constant temperature ($T=0.65$) while varying $c$ and (c) a constant concentration ($c=0.95$) varying $T$. The horizontal dashed line is at $\Phi(t) = 0.2$, the threshold used to define the overlap function $R(t)$. (b) and (d): the corresponding rotational overlap function $R(t)$ for (b) constant $T$ and (d) constant $c$. The inset in (b) compares the relaxation times $\tau_R$ measured by $R(\tau_R) = 0.3$ and $\tau_\alpha$ measured by $F_s(q, \tau_\alpha) = 1/e$ (computed from all particles) for $T=0.60$.</figcaption></figure>

rotational relaxations are strongly coupled in the deep glassy regime. If instead of varying the concentration of the dimers one changes the temperature, one finds qualitatively the same glassy slowing down phenomenology, see figures 14(c) and (d).

Finally, we define the corresponding four-point correlation function associated with the rotational degrees of freedom by:

$$ \chi_4^R(t) = N_d \left( \langle \hat{R}^2(t) \rangle_{\hat{r}^{N_d}} - \langle \hat{R}(t) \rangle_{\hat{r}^{N_d}}^2 \right), \quad (48) $$

where $\hat{R}(t) = \frac{1}{N_d} \sum_{i \in D} \theta(b - \varphi_i(t))$. Figure 15(a) shows the time evolution of $\chi_4^R(t)$ varying $c$ at fixed $T$. We find that $\chi_4^R(t)$ grows systematically with increasing glassiness, that is here the concentration of dimers. We find the same trend when $T$ is decreased while $c$ is fixed, see figure 15(b). We thus conclude that randomly bonded glass-forming liquids demonstrate growing dynamical heterogeneities approaching the glass transition
Creating equilibrium glassy states via random particle bonding

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Dynamical heterogeneity of the orientational degrees of freedom. (a) and (b): four-point correlation function $\chi_4^R(t)$ associated with the rotational degrees of freedom computed from the overlap function $\hat{R}$ at a constant temperature ($T = 0.65$) while varying $c$, panel (a), and a constant concentration ($c = 0.95$), varying $T$, panel (b).</figcaption></figure>

in terms of both positional and rotational degrees of freedom. We note that dynamical heterogeneities in rotational motions have so far not been studied widely in computer simulations [89, 90], whereas these are relevant for most molecular experiments [91].

## 7. Conclusion and discussion

We have studied randomly bonded glass-forming liquids where pairs of neighbor particles selected from an equilibrium configuration are bonded permanently. We confirm theoretically and numerically that random bonding with a neighbor cut-off is not in strict equilibrium right after bonding. However, if one generates the bonds using a spherical cut-off as in [44], the deviation from equilibrium is very small, and the aging process stops soon after the timescale of vibrations. Therefore, this random bonding method can be used to probe the (almost) equilibrium dynamics of stable bonded glass-forming liquids deep inside the energy landscape.

Our detailed computer simulations demonstrated that (1) there is no decoupling between self and collective correlation functions, (2) fragility does not change with increasing concentration of bonds, and (3) dynamical heterogeneity keeps growing as approaching the glass transition. All these features are thus in contrast to the behavior found in the dynamics of randomly pinned systems. These discrepancies are (partly) related to the preservation of the translational invariance in the random bonding process, emphasizing the importance of the details on how the quenched disorder is generated.

Most previous studies on low-temperature glassy dynamics have been performed in simple spherical particle systems, whereas most real molecular liquid experiments characterize orientational relaxation processes probed by dielectric measurements, making a conceptual gap between simulation and experiment. Instead, our randomly bonded system with rotational relaxation allows us to study phenomena that are relevant for real
experiments. As a future investigation, this approach permits thus to measure various rotational observables and test the validity of the Stoke–Einstein–Debye relationship in the deeply supercooled state [92, 93]. On the more theoretical side, it would be interesting to compute the phase diagram of randomly bonded glass formers based on the framework developed in [94] because the obtained results will be useful to connect the dynamics of molecular systems to those of gels. On the more applied side, it would be extremely interesting to revisit a series of random bonding protocols that are routinely used to prepare amorphous solids such as epoxy resins [45, 46], vitrimers [47, 48], colloidal or emulsion clusters [46, 49–57], in order to understand to what extend these protocols give rise to equilibrated glass samples whose properties remain stable over time.

In conclusion, we emphasize that the bonding approach presented in this work is not limited to the creation of dimers, since it can easily be extended to trimers, oligomers, etc, and this is in contrast to methods that have been proposed earlier. This freedom will thus permit in the future to study the (nearly) equilibrium properties of glass-forming systems at thermodynamic state points, which have so far been inaccessible to simulations.

## Acknowledgments

We thank Y Iwashita for many useful discussions and for previous collaboration on this subject [44]. We also thank J Brujic, S Corezzi, D Fioretto, F Sciortino, and E Zaccarelli for useful discussions and support. W K is a senior member of the Institut Universitaire de France.

## References

[1] Ediger M D, Angell C A and Nagel S R 1996 *J. Phys. Chem.* **100** 13200

[2] Debenedetti P G and Stillinger F H 2001 *Nature* **410** 259

[3] Rodney D, Tanguy A and Vandembroucq D 2011 *Model. Simul. Mater. Sci. Eng.* **19** 083001

[4] Swallen S F, Kearns K L, Mapes M K, Kim Y S, McMahon R J, Ediger M D, Wu T, Yu L and Satija S 2007 *Science* **315** 353

[5] Queen D, Liu X, Karel J, Metcalf T and Hellman F 2013 *Phys. Rev. Lett.* **110** 135901

[6] Yu H-B, Luo Y and Samwer K 2013 *Adv. Mater.* **25** 5904

[7] Yoon H and McKenna G B 2018 *Sci. Adv.* **4** eaau5423

[8] Raegen A N, Yin J, Zhou Q and Forrest J A 2020 *Nat. Mater.* **19** 1110

[9] Ediger M D 2017 *J. Chem. Phys.* **147** 210901

[10] Rodriguez-Tinoco C, Gonzalez-Silveira M, Ramos M A and Rodriguez-Viejo J 2022 *Riv. Nuovo Cimento* **45** 325

[11] Marinari E and Parisi G 1992 *Europhys. Lett.* **19** 451

[12] Hukushima K and Nemoto K 1996 *J. Phys. Soc. Japan* **65** 1604

[13] Yamamoto R and Kob W 2000 *Phys. Rev. E* **61** 5473

[14] Santen L and Krauth W 2000 *Nature* **405** 550

[15] Grigera T S and Parisi G 2001 *Phys. Rev. E* **63** 045102

[16] Gutiérrez R, Karmakar S, Pollack Y G and Procaccia I 2015 *Europhys. Lett.* **111** 56009

[17] Ninarello A, Berthier L and Coslovich D 2017 *Phys. Rev. X* **7** 021039

[18] Kim K 2003 *Europhys. Lett.* **61** 790

[19] Cammarota C and Biroli G 2012 *PNAS* **109** 8850

[20] Noé F, Olsson S, Köhler J and Wu H 2019 *Science* **365** eaaw1147

[21] Wu D, Wang L and Zhang P 2019 *Phys. Rev. Lett.* **122** 080602
Creating equilibrium glassy states via random particle bonding

[22] McNaughton B, Milošević M, Perali A and Pilati S 2020 *Phys. Rev. E* **101** 053312

[23] Wu D, Rossi R and Carleo G 2021 *Phys. Rev. Res.* **3** L042024

[24] Hibat-Allah M, Inack E M, Wiersema R, Melko R G and Carrasquilla J 2021 *Nat. Mach. Intell.* **3** 952

[25] Gabrié M, Rotskoff G M and Vanden-Eijnden E 2022 *PNAS* **119** e2109420119

[26] Ciarella S, Trinquier J, Weigt M and Zamponi F 2023 *Mach. Learn.: Sci. Technol.* **4** 010501

[27] Brito C, Lerner E and Wyart M 2018 *Phys. Rev. X* **8** 031050

[28] Hagh V F, Nagel S R, Liu A J, Manning M L and Corwin E I 2022 *PNAS* **119** e2117622119

[29] Kob W and Berthier L 2013 *Phys. Rev. Lett.* **110** 245702

[30] Ozawa M, Kob W, Ikeda A and Miyazaki K 2015 *PNAS* **112** 6914

[31] Ozawa M, Ikeda A, Miyazaki K and Kob W 2018 *Phys. Rev. Lett.* **121** 205501

[32] Biroli G and Bouchaud J-P 2023 *C. R. Phys.* **24** 1

[33] Gokhale S, Nagamanasa K H, Ganapathy R and Sood A 2014 *Nat. Commun.* **5** 1

[34] Williams I, Turci F, Hallett J E, Crowther P, Cammarota C, Biroli G and Royall C P 2018 *J. Condens. Matter Phys.* **30** 094003

[35] Kikumoto G, Torii N, Fukao K, Royall C P, Yao H, Saruyama Y and Tatsumi S 2020 arXiv:2003.06089

[36] Das R, Bhowmik B P, Puthirath A B, Narayanan T N and Karmakar S 2023 *PNAS Nexus* **2** gad277

[37] Charbonneau P and Tarjus G 2013 *Phys. Rev. E* **87** 042305

[38] Chakrabarty S, Das R, Karmakar S and Dasgupta C 2016 *J. Chem. Phys.* **145** 034507

[39] Kim K, Miyazaki K and Saito S 2011 *J. Condens. Matter Phys.* **23** 234123

[40] Chakrabarty S, Karmakar S and Dasgupta C 2015 *Sci. Rep.* **5** 1

[41] Jack R L and Fullerton C J 2013 *Phys. Rev. E* **88** 042304

[42] Kob W and Coslovich D 2014 *Phys. Rev. E* **90** 052305

[43] Li Y-W, Zhu Y-L and Sun Z-Y 2015 *J. Chem. Phys.* **142** 124507

[44] Ozawa M, Iwashita Y, Kob W and Zamponi F 2023 *Nat. Commun.* **14** 113

[45] Johari G P 1994 Dynamics of irreversibly forming macromolecules *Disorder Effects on Relaxational Processes* ed R Richert and A Blumen (Springer) pp 627–57

[46] Corcione C E, Greco A and Maffezzoli A 2006 *Polym. Eng. Sci.* **46** 493

[47] Kloxin C J and Bowman C N 2013 *Chem. Soc. Rev.* **42** 7161

[48] Denissen W, Winne J M and Du Prez F E 2016 *Chem. Sci.* **7** 30

[49] Duguet E, Désert A, Perro A and Ravaine S 2011 *Chem. Soc. Rev.* **40** 941

[50] Peng B, Smallenburg F, Imhof A, Dijkstra M and van Blaaderen A 2013 *Angew. Chem.* **125** 6841

[51] McMullen A, Muñoz Basagoiti M, Zeravcic Z and Brujic J 2022 *Nature* **610** 502

[52] McMullen A, Holmes-Cerfon M, Sciortino F, Grosberg A Y and Brujic J 2018 *Phys. Rev. Lett.* **121** 138002

[53] Yuan Q, Gu J, Zhao Y-n, Yao L, Guan Y and Zhang Y 2016 *ACS Macro Lett.* **5** 565

[54] Goldbart P M, Castillo H E and Zippelius A 1996 *Adv. Phys.* **45** 393

[55] Carbas R, Marques E, Da Silva L and Lopes A 2014 *J. Adhes.* **90** 104

[56] Corezzi S, Fioretto D and Rolla P 2002 *Nature* **420** 653

[57] Mereu I, Liotta A, Comez L and Corezzi S 2015 *J. Chem. Phys.* **142** 154905

[58] Krakoviack V 2010 *Phys. Rev. E* **82** 061501

[59] Szamel G 2019 *J. Stat. Mech.* 104016

[60] Ikeda H, Zamponi F and Ikeda A 2017 *J. Chem. Phys.* **147** 234506

[61] Zdeborová L and Krzakala F 2016 *Adv. Phys.* **65** 453

[62] Scheidler P, Kob W and Binder K 2004 *J. Phys. Chem.* **108** 6673

[63] Kapteijns G, Ji W, Brito C, Wyart M and Lerner E 2019 *Phys. Rev. E* **99** 012106

[64] Marchand T, Ozawa M, Biroli G and Mallat S 2022 arXiv:2207.04941

[65] Cammarota C and Biroli G 2013 *J. Chem. Phys.* **138** 12A547

[66] Kob W and Andersen H C 1995 *Phys. Rev. E* **51** 4626

[67] Hiroike K 1960 *Prog. Theor. Phys.* **24** 317

[68] Morita T and Hiroike K 1961 *Prog. Theor. Phys.* **25** 537

[69] Ozawa M, Parisi G and Berthier L 2018 *J. Chem. Phys.* **149** 154501

[70] Biroli G, Bouchaud J-P, Cavagna A, Grigera T S and Verrocchio P 2008 *Nat. Phys.* **4** 771

[71] Hocky G M, Markland T E and Reichman D R 2012 *Phys. Rev. Lett.* **108** 225506

[72] Hocky G M, Berthier L, Kob W and Reichman D R 2014 *Phys. Rev. E* **89** 052311

[73] Yaida S, Berthier L, Charbonneau P and Tarjus G 2016 *Phys. Rev. E* **94** 032605

[74] Nandi U K, Kob W and Bhattacharyya S M 2021 *J. Chem. Phys.* **154** 094506

[75] Nandi U K, Patel P, Moid M, Nandi M K, Sengupta S, Karmakar S, Maiti P K, Dasgupta C and Bhattacharyya S M 2022 *J. Chem. Phys.* **156** 014503
[76] Andersen H C 1983 *J. Comput. Phys.* **52** 24

[77] Götze W 2008 *Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory (International Series of Monographs on Physics)* (Oxford University Press) (https://doi.org/10.1093/acprof:oso/9780199235346.001.0001)

[78] Cammarota C, Ozawa M and Tarjus G 2023 The Kauzmann transition to an ideal glass phase *Spin Glass Theory and Far Beyond: Replica Symmetry Breaking After 40 Years* ed P Charbonneau, E Marinari, M Mézard, G Parisi, F Ricci-Tersenghi, G Sicuro and F Zamponi (World Scientific) pp 203–18

[79] Angell C A 1995 *Science* **267** 1924

[80] Alba-Simionesco C and Tarjus G 2022 *J. Non-Cryst. Solids* **X** **14** 100100

[81] Gnan N, Schrøder T B, Pedersen U R, Bailey N P and Dyre J C 2009 *J. Chem. Phys.* **131** 234504

[82] Schrøder T B and Dyre J C 2014 *J. Chem. Phys.* **141** 204502

[83] Hurley M and Harrowell P 1995 *Phys. Rev. E* **52** 1694

[84] Kob W, Donati C, Plimpton S J, Poole P H and Glotzer S C 1997 *Phys. Rev. Lett.* **79** 2827

[85] Yamamoto R and Onuki A 1998 *Phys. Rev. E* **58** 3515

[86] Karmakar S, Dasgupta C and Sastry S 2009 *PNAS* **106** 3675

[87] Donati C, Franz S, Glotzer S C and Parisi G 2002 *J. Non-Cryst.* **307** 215

[88] Shiraishi K, Mizuno H and Ikeda A 2023 *PNAS* **120** e2215153120

[89] Kawasaki T and Kim K 2019 *J. Stat. Mech.* 084004

[90] Kou B, Cao Y, Li J, Xia C, Li Z, Dong H, Zhang A, Zhang J, Kob W and Wang Y 2018 *Phys. Rev. Lett.* **121** 018002

[91] Cicerone M T and Ediger M D 1995 *J. Chem. Phys.* **103** 5684

[92] Tarjus G and Kivelson D 1995 *J. Chem. Phys.* **103** 3071

[93] Kawasaki T and Kim K 2019 *Sci. Rep.* **9** 8118

[94] Hall R W and Wolynes P G 2003 *Phys. Rev. Lett.* **90** 085505