Creating equilibrium glassy states via random particle bonding

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.


I. INTRODUCTION
Since glasses are out-of-equilibrium materials, their properties depend on their 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 the deep minima inside this landscape is a crucial challenge for experiments as well as 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][5][6][7][8][9][10].On the atomistic simulation side, various sampling techniques have been developed, such as replica exchange methods [11][12][13], Monte-Carlo simulations with smart updates [14][15][16][17], random pinning methods [18,19], as well as machine learning-assisted sampling techniques [20][21][22][23][24][25][26].In particular, the swap-Monte Carlo [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 vary in a systematic manner 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, accelerating significantly the relaxation dynamics.While this method 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, and hence the system composed of 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][30][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, it has been demonstrated that the dynamics is significantly altered from the bulk state, possibly due to the violation of translational invariance, leading to a strong decoupling between self and collective behavior [30,37,38], a decrease of fragility as the concentration of pinned particles is increased [39,40], and a suppression of dynamical heterogeneity when the glass transition is approached [39,[41][42][43], in stark contrast to the behavior of standard bulk materials that are approaching their glass transition.
Recently, we have proposed the 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 stereolitography [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, see e.g.[46,[54][55][56][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 Ref. [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 while the system is in equilibrium if the bonded particles are chosen completely randomly (in this case, the bond lengths are thus arbitrary), the bonding from pairs from nearest neighbor particles (which corresponds to a more realistic situation) does not ensure strict equilibration.However, it turn 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 the ones from the dynamics of randomly pinned systems.We find that self and collective correlation functions for the translational degrees of freedom, as well as the rotational correlation function, are strongly coupled.Besides, we observe that the kinetic fragility does not change by increasing the concentration of bonds.Finally, we find that dynamical heterogeneity keeps growing with approaching 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.

II. STATISTICAL MECHANICS OF BONDED SYSTEMS
In this section, we discuss the statistical mechanics of randomly-bonded glass formers composed of monomers and dimers.We randomly choose pairs of neighboring particles from an equilibrium configuration and bond then together by introducing a rigid-body constraint [44].The main goal of this section is to demonstrate that the creation of these bonds does not perturb significantly the thermal equilibrium of the system, and we do this by using similar ideas as in the random pinning protocol that fixes the positions of particles [58,59].However, since the case of random bonding is quite subtle, we begin with a review of the random pinning approach and we then discuss how to extend the idea to random bonding.
A. Quiet freezing of variables Consider a system whose degrees of freedom are arbitrarily split into two distinct vectors x and y with Hamiltonian H(x, y).The equilibrium probability distribution of the total system is ρ(x, y) = exp[−βH(x, y)]/Z, where Z is the partition function and β is the inverse temperature.Let us assume that an equilibrium configuration of the system can be generated by sampling from ρ(x, y).Consider now the degrees of freedom 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) Next, consider a setting in which one first generates a configuration of the y degrees of freedom from their marginal probability distribution ρ f (y), and then considers a system in which y are "frozen" and whose dynamical variables are the x degrees of freedom, with equilibrium probability Here, the 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 x.
Using the chain rule of probabilities, we have Hence, a pair {x, y} generated from the joint distribution ρ(x, y) can be considered either as an equilibrium configuration of the full system, or as a realization y of the quenched disorder of the frozen system obtained from ρ f (y) together with a typical equilibrium realization x of that system obtained from ρ(x|y).We conclude that one can generate a pair {x, y} from ρ(x, y), then freeze (or pin) y, and automatically obtain an equilibrium configuration x (but only a single one) of the distribution ρ(x|y) that describes the frozen system with quenched disorder y.
We now define an "annealed" average ⟨• • • ⟩, a "thermal" average ⟨• • • ⟩ y , and a "disorder" average respectively.Furthermore we define the "quenched" average as ⟨• • • ⟩ y and because of the chain rule in Eq. ( 3), the annealed and quenched averages coincide: Conversely, the validity of Eq. ( 5) for an arbitrary observable implies Eq. (3).Equations ( 3) and ( 5) provide the core idea behind the random pinning and swap Monte-Carlo simulations, as we will discuss in detail below.
Since the freezing of degrees of freedom eliminates possible relaxation channels [17,28,[60][61][62], it can be expected that sampling x from ρ(x|y) at fixed y is harder than sampling x and y together from ρ(x, y).(In some cases, however, the conditional sampling gives rise to a faster sampling, see, e.g., Ref. [63].)In other words, a local dynamics for x that satisfies detailed balance with respect to ρ(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 ρ(x, y).Generating additional independent configurations from ρ(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 ρ 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 ρ(x|y) as in Eq. ( 2), but the distribution of the quenched disorder y, i.e., ρ q (y), is chosen independently, and in general ρ q (y) ̸ = Z f (y)/Z, where Z f (y) is defined in Eq. ( 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 ρ q (y) fixed while changing the temperature associated to x.The choice ρ q (y) = ρ f (y) = Z f (y)/Z, in which the quenched disorder depends on temperature, guarantees the equality of annealed and quenched averages in Eq. (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 Ref. [64] for a pedagogical discussion.
• The thermal average ⟨O(x)⟩ y of an extensive observable that is the sum of local terms of the form 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, Furthermore, note that the observable O(x, y) can also depend on the frozen degrees of freedom y (e.g. if O = H), in which case its thermal average also depends explicitly on the y through O.This does not affect the proof of equivalence of the annealed and quenched averages in Eq. ( 5), and it also does not affect the self-averaging property.Note that self-averaging only holds for extensive variables, i.e., 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 the random pinning procedure, the pair correlation function measured in a large configuration will be independent of the choice of the pinned particles, and is self averaging.On the other hand, the density in a specific point of space will depend on this choice, and is not self averaging.
• A local dynamics on x that samples from ρ(x|y) could be so slow that it becomes non-ergodic at an ideal glass transition.This has been shown to happen in randomly pinned mean field spin glasses [19,65] and finite-dimensional glasses [29][30][31].In this case, the system remains stuck forever around the initial equilibrium configuration x, and generating independent additional configurations is impossible using standard simulation approaches such as simple Monte Carlo simulations or molecular dynamics.
It is crucial to stress that the above construction, which we will call quiet freezing of variables, supposes that the division of the degrees of freedom into x and y is done before the thermalized configuration is constructed.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 Eq. (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.

B. Random pinning and wall pinning
The construction of Sec.II A can be applied in a straightforward manner to the random pinning procedure [18,58,59].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 additional degrees of freedom and sum up all the possible permutations in the partition function [27,[67][68][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.
1. Random pinning: Choosing at random the particles that will be pinned The Hamiltonian H and the partition function Z are given by where r i , p i , and m i are the position, momentum, and mass of the i-th particle, respectively.U is the potential energy, and β = 1/T is the inverse temperature.Particles are supposed 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, e.g., r N = (r 1 , r 2 , ..., 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 Eq. ( 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 y = {r N f , p N f }, and a set of N −N f unpinned (or mobile) particles with the remaining coordinates, x.Because the labeling of particles is arbitrary, we can consider that the first N f particles are the pinned ones.Note that the splitting is here performed before any thermalization, as we emphasized in Sec.II A. One can then thermalize the whole system of N particles, hence the joint set {x, y}, in the liquid phase.Because in such phase particles can freely diffuse, we expect the pinned particles to be uniformly distributed in the volume V .At this point, the positions of the N f particles, 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 particles, setting the momenta of pinned particles to zero does not affect the distribution of x (in other words, only the configurational part of the Hamiltonian matters).Hence, the particles x can be considered to be in equilibrium with the pinned particles 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 chooses at random the N f particles to be pinned.But, 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 choosing 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. Wall pinning: Choosing the particles to be pinned according to a geometrical criterion In Ref. [59], 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 in order to probe the equilibrium properties of liquids [70][71][72][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 Fig. 1.The shape of these regions is completely arbitrary.We can then generate an equilibrium 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 Sec.II A 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, the average of an arbitrary observable O(r N ) that is invariant under permutations of the particle labels can be written as where 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(r N ) and U (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 We conclude, as in Sec.II A, that we can define an "annealed" average ⟨• • • ⟩, a "thermal" average ⟨• • • ⟩ W conditioned to the wall, and a "disorder" average • • • over the realizations of the wall, from and because of the chain rule in Eq. ( 10), the annealed and quenched averages coincide, i.e., ⟨• 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 Ω 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 Ω W region.Yet this is not the case for dilute systems, and hence one has to impose an additional constraint on the dynamics of fluid particles such as a hard wall condition.
Note that the only requirement on the observable O(r N ) is that it is invariant under permutations of particles.We also mention that one can construct an observable O(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 FIG. 1. Illustration of the wall pinning procedure, of which we give in the main text two different proofs.The first, as in Ref. [59], 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 Eq. ( 10).The second corresponds to ordering the particles according to their z component, and freezing the k particles with the smallest z, as in Eq. ( 13).The number k has to be chosen such that the wall region has the desired size.
where I[E] is the indicator function of event E, which is one if E is realized and zero otherwise.Then, any observable of the form in Eq. ( 12) can be used to characterize the fluid system only, while keeping the quenched and annealed averages coincident.

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 r i ≺ r j if r i • e < r j • e, where e is an arbitrary unit vector.Similarly, r i ⪯ r j if r i • e ≤ r j • 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, {r 1 , • • • , r N }.We now define the vector for the k first ordered particles and the vector for the rest of the ordered particles as r ⪯k = (r 1 , • • • , r k ) and r ≻k = (r k+1 , • • • , r N ), respectively.Considering again only permutation-symmetric observables O(r N ), and using the permutation symmetry, we can order the N particles and write where We thus obtain a relation similar to Eq. ( 10).
The core idea both in the "wall pinning" proof (Sec.II B 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 r ⪯k and r ≻k in Eq. ( 13) can be released.This results in a similar expression to the wall pinning case, but using a separation between the k-th and (k + 1)-th ordered particles: dr ≻k e −βU (r ⪯k ,r ≻k ) . ( This relation is identical to Eq. ( 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 while the first proof corresponds to a "fixed-boundary" ensemble.As usual, the two ensembles become equivalent in the thermodynamic limit if the wall and fluid regions are both macroscopic, i.e., if k ∼ N .An illustration is given in Fig. 1.
While the derivations presented here and in Sec.II B 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 r N .At that instant, we can perform the sorting of the particles (e.g., according to their z component, setting e = 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 ↔ 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 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 , • • • , z k }.This hard constraint can be implemented in a Monte Carlo 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 ⟨• • • ⟩ W over the fluid particles.One should then either perform the disorder average over the wall particles, i.e., • • •, 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 one single run is representative of the average.
Note that, as in the wall pinning case, we can construct an observable O(r N ) = O(r ≻k ) that depends only on the last N − k ordered vector r ≻k = (r k+1 , • • • , 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.

Variable transformation to create virtual dimers
We now show that the ideas presented in Sec.II B 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 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 respectively, with dimer i being composed by particles i and i+1.We then denote r 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 Ω d for the dimers and Ω m for the monomers.Because of permutation invariance, there are distinct equivalent ways of constructing the dimers 1 , and the sum over all these equivalent possibilities reconstructs 1 There are N  the whole integration volume of the original monomer system2 .Taking advantage of the permutation symmetry of the problem, we can then write 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 Fig. 2(a).Once again, we emphasize that this is a virtual operation, and the actual system is not altered at all, i.e., 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 (Fig. 2(b)), using the center of mass and relative position for the dimer i ∈ D = {1, 3, 5, • • • , 2N d − 1} composed of the i-th and (i + 1)-th particles.These are given by R i = miri+mi+1ri+1 mi+mi+1 and ri = r i − r i+1 , respectively.The corresponding momenta are denoted by P i and pi , respectively.Also, the dimer total mass and the reduced mass are given by M i = m i + m i+1 and µ i = mimi+1 mi+mi+1 , respectively.This variable transformation is formally written as (r Making this variable transformation, we can rewrite the Hamiltonian H and the partition function Z as Note that the integration spaces Ω d and Ω m in Eq. ( 21) are not exactly the same as the ones in Eq. ( 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 (r . We note that because of Liouville's theorem the Jacobian is unity for this transformation, namely, dr i dp i = dr i dθ i dφ i dp ri dp θi dp φi .Also, the kinetic part of the Hamiltonian can be written as where I i = µ i r2 i is the moment of inertia, p ri = µ i ṙi , p θi = I i θi , and p φi = I i (sin 2 θ i ) φi are the new momenta and the dot denotes a time derivative.We note that as one can see in Eq. ( 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: We can then define the total (or annealed) average as 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 Sec.II B 2, one should imagine the following procedure: First, an equilibrium configuration of 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 into their respective domains3 , r 2N d ∈ Ω d and r Nm ∈ Ω m , one can keep running the dynamics.If at some point one variable goes out of its domain, then one should run again the procedure 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 the annealed average, as in Sec.II B 2.

Random bonding
Using the notations introduced in Sec.II C 1, it is now straightforward to perform random bonding, i.e., freeze the dimer bond lengths to define a quenched average.We consider the dimer lengths rN d as being frozen, i.e., rN d plays the same role as the positions of the pinned particles in the random pinning approach, and of the frozen degrees of freedom y in the general discussion of Sec.II A. 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 rN d , that play the role of the thermal degrees of freedom x in Sec.II A. Note that the kinetic term associated to the frozen bond lengths is decoupled from all the other degrees of freedom.Hence, like in the particle pinning case, we can set the momenta p ri 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 rN d , is then given by and the thermal average of an observable for a particular realization rN d is defined by Comparing this with Eq. ( 25), we can write From this expression, we can deduce the probability distribution of the frozen variables rN d , which is given by and define a "disorder" average over the realization of rN d as Finally, as it can be expected from the discussion of Sec.II A, we can obtain the identity that is at the basis of the random pinning procedure [58,59]: 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 to the frozen ones, corresponds exactly to the annealed average of the bulk system without bonding.By using Eqs.( 26) -( 30), we get 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 II C 1 to identify the N d dimers.The bond lengths are finally frozen, i.e., 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 II A are fulfilled, and we can rigorously state that an equilibrium configuration of a system containing N d dimers and N m monomers has 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 correlation function of the initial system of monomers, and arbitrarily long dimers (up to the system size) will be present 4 .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 (see section II A) will therefore not be strictly respected, and the consequences of this choice have to be assessed.In the following we discuss several possible choices, and study two of them numerically.

'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 |r ij | = |r i − r j | and we relabel them as r 1 and 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 among particles, which we relabel as r 3 and 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 will show 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 for creating 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. Ω d is defined by the constraint (as above, the tilde notation indicates a distance between a pair of particles, rather than an absolute position): i.e., 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 Ω d on the distances between particles belonging to distinct dimers.For example, in order to guarantee that |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 Ω 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 are 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 Ω m , i.e., the smallest monomer-monomer and monomer-dimer distance must stay larger than the largest dimer bond length R max .
In order to perform the quenched average dynamically, the above constraints must be implemented either by rejecting moves that violate them in a Monte Carlo 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 Ω 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, since, e.g., the interactions between the particles are non-zero even for arbitrarily large distances.

'Random bonding' within a cutoff
We now consider a different procedure, close to the rigorous "random bonding" described at the end of section II C 2, but now imposing that the length of the formed dimers does not exceed a maximum value R b (see Fig. 3(a)).This is the method used in Ref. [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 Ω d is defined by the only constraint that the dimer bond lengths ri ≤ R b , which is trivially satisfied if the bond lengths are frozen.For the rest, dimers are free to translate and rotate.
2. The monomer space Ω 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 a property 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 where g inter is the contribution of pairs of particles that are not bonded together by a frozen bond, and g intra (r) is the contribution of the pairs of particles that are connected.By construction, g intra (r) will have a discontinuity at R b , above which it will jump to zero.On the other hand, g(r) is the pair correlation of an equilibrium system of monomers, and is continuous.We conclude that g 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.

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 to the orientation of the frozen dimers, which are restricted to an angular sector θ b around the z direction (see Fig. 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, in comparison to the isotropic case.

A. Model
We employ the Kob-Andersen binary mixture [66], in which particles interact through the Lennard-Jones pair potential, where α, β ∈ {A, B} are species indexes.Both species have the same mass, which is set to m = 1.The value of the parameters σ αβ and ϵ αβ are given in Ref. [66].The units of length and energy are set by the parameters σ = σ AA = 1 and ϵ = ϵ AA = 1, respectively, and we put the Boltzmann constant k B = 1.The potentials are cut and shifted at a distance 2.5σ αβ .We simulate systems composed of N particles in a cubic box of side L with periodic boundary conditions at a number density ρ = 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].

B. 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 composed of 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 Fig. 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, ri = |r i − r i+1 |, which means that the particles i and i + 1 now form a dimer by a rigid body constraint.We repeat , such that c = 0 corresponds to the system with only monomers (hence the original bulk 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 percent 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 tend to align along the unit vector of the z-direction, e z = (0, 0, 1).We define a unit vector of the orientation of each dumbbell molecule composed of particles i and i + 1 by n i = ri /|r i |, where ri = (x i − x i+1 , y i − y i+1 , z i − z i+1 ).In the 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 neighbor region by taking into account the orientation of the molecule.To this end, we define a magnitude of alignment (or polarization) for each molecule, P i , given by We then introduce another cut-off threshold P b such that particles having P i > P b = cos θ b can be considered as candidate neighbors for bonding.This new neighbor region is schematically shown in Fig. 3(b).As in procedure 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 Fig. 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).

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

A. Radial distribution function
We first characterize the static structure by the radial distribution function g(r) for all particles, which is defined by  where ∆r is the bin size for computation.We use ∆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 ∈ [0, 10]) starting right after bonding (t = 0) and a longer time window (t ∈ [10000, 20000]) starting after a waiting time t w = 10000.The former time scale corresponds to a vibrational one, whereas the latter timescale corresponds to the time scale for escaping from the cage (see Fig. 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 ∈ [0, 10] and t ∈ [10000, 20000] superimpose very well, i.e., the system ends up in a new stationary state very quickly.
In order to understand better 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 intra (r)) and the rest (inter molecule contribution, g inter (r)), following the discussion in Sec.II C 4. By construction, g(r) = g intra (r) + g inter (r).In Fig. 5(b) we present g(r), g intra (r), and g inter (r) for c = 0.95 in the time window t ∈ [0, 10].This graph shows that g 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 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 inter (r), at strictly t = 0, has a small dip at r = R b (not shown), compensating the sharp edge in g intra (r), such that the total g(r) is smooth as in the original system.At t > 0, only g inter (r) evolves with time, which removes the dip around r = R b and ends up in a smooth curve as shown in Fig. 5(b).Consequently, the total contribution, g(r) = g intra (r) + g 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.

B. 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 where R i is the center of mass of molecule composed of the i-th and (i + 1)-th particles.In Fig. 5(c), we present g 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.

C. 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 Fig. 3(b) (see also the discussion in Sec.II C 5).In Fig. 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 ∈ [0, 10]) starting right after the bonding (t = 0) and a longer time window (t ∈ [10000, 20000]) after a waiting time of t w = 10000.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 seen in g ori (r) presented in Fig. 6(b).At the shorter time scale, g ori (r) is overall larger, due to 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) ≈ 0.5 at r ≫ 1 and is overall similar to that of Fig. 5(c) for a spherically-bonded system.This aging process can be directly quantified by the total molecule polarization, given by Figure 6(c) shows the time evolution of P tot .The system with P b = 0 has random molecular orientations, producing P tot ≈ 0.5 during the entire simulation.On the contrary, the system with P b = 0.9 takes a higher value, P tot ≈ 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 ≈ 0.5.

V. DYNAMICAL PROPERTIES
In the previous section we have confirmed that the random bonding protocol with a spherical cut-off produces a nearly 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.

A. 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 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, the location of the main peak in the static structure factor [66].  F Di s (q, t) displays a higher value of the plateau than the one found in F Mono s (q, t) or F All s (q, t), which is reasonable since 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 is evidence that the monomers and dimers have a very similar relaxation dynamics in terms of translational motions.We define the relaxation time τ α as the time at which the intermediate scattering function decays to 1/e and present the τ α versus 1/T plot for c = 0.5 in Fig. 7(b).The relaxation time for the center of mass of dimers exceeds τ α for the monomers by a factor around 2, independent of T , which shows 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 in the following on F All s (q, t) and drop the superscript unless otherwise specified.In Fig. 7(a), we also include data for a different waiting time t w , i.e., 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 in 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 by simulations right after bonding, akin to random pinning.
The temperature and c−dependence of F s (q, t) is presented in Fig. 8.We find that the influence of bonding is very small if 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 allows to access an extremely slow glassy dynamics in almost equilibrium.By using the data for τ α (see Fig. 11(a) below), we estimate that the bonded system at T = 0.424 and c = 0.95 (that are prepared by making bonds from the original T = 0.424 configurations) has an equilibrium relaxation time τ α ≈ 10 12 , thus about a factor of 10 7 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.
In order to see the influence of T and c together we present in Fig. 9 the iso-τ α 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 τ α , hinting at a simple functional relation between τ α , 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 intermediate scattering functions.Previous studies have reported that in randomly pinned fluids the collective part shows apparent freezing while the self part did not, rendering the analysis of the system dynamics difficult [30,37,38].In Fig. 10(a), we show 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 Fig. 10(b).This suggests that self and collective correlators do not decouple and therefore we can conclude that the self part gives reliable dynamical information about the system.

B. Dynamical scaling
In the following we use a scaling analysis to examine how τ α depends on T and c.In Fig. 11(a), we present an Arrhenius plot for the structural relaxation time τ α (c, T ) for different c.For the c = 0 system, which we will call the "original" system, τ α follows the well-known non-Arrhenius temperature dependence [66]: τ α (c = 0, T ) = τ 0 exp E(c=0,T ) T , where E(c = 0, T ) is a T -dependent activation energy accounting for the non-Arrhenius behavior.With increasing c, τ α (c, T ) increases as expected, which implies that the activation energy E(c, T ), defined by τ α (c, T ) = τ 0 exp E(c,T ) T , grows due to the addition of bonds.Interestingly, we find that the temperature dependence of τ α (c, T ) can be rescaled by an unknown function m(c), namely, τ α (c, T ) = τ α (c = 0, T /m(c)), as shown in Fig. 11(b).The inset shows m(c), which has been determined manually for each c such that τ α (c, T ) superimposes with τ α (c = 0, T ).One sees that to a first approximation m(c) is linear, but a slight upward bending can be noticed.The existence of a master curve is so far empirical, and it holds at least for the T -and c-range probed by our simulations.Yet the scaling 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 m(c).This implies that the fragility of randomly-bonded glass-forming liquids does not change inherently across different c, in contrast to randomly pinned systems [39,40].This scaling for the randomly bonded systems also implies that E(c, T ) = m(c)E(c = 0, T /m(c)).On the other hand, in the random pinning case, it was argued that E(c, T ) = q(c)E(c = 0, T ), where q(c) is an increasing function of c with q(c = 0) = 1 [38], which explains the fact that fragility decreases with increasing c.Moreover, Ref. [38] argued that if E(c = 0, T ) has a singularity at a finite T K , e.g., E(c = 0, T )/T ∼ 1/(T − T K ), pinned systems at c > 0 inherit the same singular temperature T K irrespective of c. Comparing these results with the (T, c)-dependence of our bonded system, one can conclude that the dynamics of randomly bonded glass-forming liquids is qualitatively quite different from the one of randomly pinned systems also in terms of the temperature dependence of structural relaxation.
The dynamical scaling shown in Fig. 11(b) reminds us of simple systems where structure and dynamics are invariant to a good approximation along isomorphs in the phase diagram [79,80].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 in 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 /m(c).Although the relaxation time τ α 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.

VI. DYNAMICAL HETEROGENEITY
The dynamics of glassy liquids is accompanied by strong dynamical heterogeneities, the intensity of which grows with decreasing temperature [81][82][83][84].Since our bonding procedure allows to generate configurations in the deeply supercooled regime, we 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.

A. Positional degrees of freedom
First, we compute the standard four-point correlation function χ Q 4 (t) associated with positional degrees of freedom, which is given by where Q(t) = 1 N N i=1 θ(a − |r i (t) − r i (0)|) is an overlap function taking into account all particles and θ(x) is the Heaviside step function [85].We set the distance a to the often used value a = 0.3.We note that χ Q 4 (t) defined in Eq. (43) does not contain contributions from sample-to-sample fluctuations associated with different realizations of bonds [42].Figure 13 shows χ Q 4 (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 glass transition, while the relaxation time increases significantly [39,[41][42][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, while in the present system the heterogeneities can grow unhindered.liquids demonstrate growing dynamical heterogeneities approaching the glass transition in terms of both positional and rotational degrees of freedom.We note that dynamical heterogeneities in rotational motions have so far not be studied widely in computer simulations [87,88], while these are relevant for most molecular experiments [89].

VII. CONCLUSION AND DISCUSSION
We have studied randomly bonded glass-forming liquids where pairs of neighbor particles chosen from an equilibrium configuration are bonded permanently.We confirmed 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 Ref. [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 by increasing the concentration of bonds, and 3) dynamical heterogeneity keeps growing with 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 liquids experiments characterize orientational relaxation processes probed by dielectric measurements, making a conceptual gap between the 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, e.g., test the validity of the Stoke-Einstein-Debye relation in the deeply supercooled state [90,91].On the more theoretical side, it would be interesting to compute a phase diagram of randomly bonded glass formers based on the framework developed in Ref. [92] since the obtained results will be useful to connect the dynamics of molecular systems to the ones 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][50][51][52][53][54][55][56][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 in contrast to methods that have been proposed earlier.This freedom will thus permit in the future to study the (nearly) equilibrium properties of glassforming systems at thermodynamic state points which have so far been inaccessible to simulations.

FIG. 2 .
FIG. 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.

FIG. 3 .
FIG. 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 θ b associated with directional alignment.

FIG. 5 .
FIG. 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 ∈ [0, 10] (red curve) and t ∈ [10000, 20000] (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 intra molecule contribution gintra(r) and inter molecule contribution ginter(r) for c = 0.95 at T = 0.6 in the time window t ∈ [0, 10].(c): Orientational correlation function gori(r) for the bonded systems presented in (a), measured in two different time windows.

FIG. 6 .
FIG. 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 ∈ [0, 10] (red curve) and t ∈ [10000, 20000] (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 gori(r) for the bonded systems presented in (a), measured in two different time windows.(c): Time evolution of the total polarization Ptot 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.

FIG. 7 .
FIG.7.Dynamical correlations for dimers and monomers.(a): Self-intermediate scattering function Fs(q, t) for T = 0.6 and c = 0.95 for monomers, dimers, and all particles.Dashed and solid curves indicate Fs(q, t) computed from trajectories with the waiting time tw = 0 and tw = 10 5 , respectively.(b): Relaxation time τα versus the inverse of temperature 1/T for c = 0.5.

Figure 7 (
Figure7(a) shows these intermediate scattering functions at T = 0.6 and c = 0.95.F Di s (q, t) displays a higher value of the plateau than the one found in F Mono

4 FIG. 9 . 6 FIG. 10 .
FIG. 9.Relaxation time as a function of temperature T and dimer concentration c.We report iso-τα curves in the T − c plane, obtained from the dynamical correlations shown in Fig.8.

FIG. 11 .
FIG. 11.Dynamical scaling of the relaxation time with dimer concentration.(a): Relaxation time τα(c, T ) obtained from the self-intermediate scattering function for all particles.(b): The same data using a normalized abscissa, m(c)/T .The inset shows m(c) versus c.

FIG. 14 .
FIG. 14. Dynamical correlations of the orientational degrees of freedom.(a, c): Mean rotational angle Φ(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 Φ(t) = 0.2, the threshold used to define the overlap function R(t).(b, d): The corresponding rotational overlap function R(t) for (b) constant T and (d) constant c.The inset in (b) compares the relaxation times τR measured by R(τR) = 0.3 and τα measured by Fs(q, τα) = 1/e (computed from all particles) for T = 0.60.

FIG. 15 .
FIG. 15.Dynamical heterogeneity of the orientational degrees of freedom.(a, b): Four-point correlation function χ R 4 (t) associated with the rotational degrees of freedom computed from the overlap function R at a constant temperature (T = 0.65) while varying c, panel (a), and a constant concentration (c = 0.95), varying T , panel (b).