Ergodicity breaking transition in a glassy soft sphere system at small but non-zero temperatures

TL;DR

This study investigates the ergodicity breaking transition in a glassy soft sphere system at small non-zero temperatures, revealing a transition density independent of temperature and distinct from athermal jamming, with critical behavior in the directed percolation class.

Contribution

Introduces a new method to determine ergodicity in soft sphere systems, showing the transition's universality class and its difference from athermal jamming.

Findings

Ergodicity breaking occurs at a density below athermal jamming.

Transition density approaches a temperature-independent limit.

Transition belongs to the directed percolation universality class.

Abstract

While the glass transition at non-zero temperature seems to be hard to access for experimental, theoretical, or simulation studies, jamming at zero temperature has been explored in great detail. It is a widely discussed question whether this athermal jamming transition is related to the glass transition. Motivated by the exploration of the energy landscape that has been successfully used to describe athermal jamming, we introduce a new method to determine whether the configuration space of a soft sphere system can be explored within a reasonable timescale or not, i.e., whether the system is ergodic or effectively non-ergodic. While in case of athermal jamming for a given random starting configuration only the local energy minimum is determined, we allow the thermally excited crossing of energy barriers. Interestingly, we observe that a transition exists where the system becomes effectively non-ergodic if the density is increased. In the limit of small but non-zero temperatures the density where the ergodicity breaking transition occurs approaches a value that is independent of temperature and below the transition density of athermal jamming. This confirms recent computer simulation studies where athermal jamming occurs deep inside the glass phase. In addition, with our method we determined the critical behavior of the ergodicity breaking transition and show that it is in the universality class of directed percolation. Therefore, our approach not only makes the transition from an ergodic to an effectively non-ergodic systems easily accessible and helps to reveal its universality class but also shows that it is fundamentally different from athermal jamming.

Equations

$$f_{\text{ov}}(t)=f_{\text{ov}}(t\to \infty )+A{t}^{\alpha }\times \exp (-\frac{t}{\tau })$$

$$f_{\text{ov}}(t\to \infty )\sim (\varphi -{\varphi }_{\text{c}}{)}^{\beta }$$

$$\tau \sim (\varphi -{\varphi }_{\text{c}}^{\tau }{)}^{\nu }$$

Full text

Published in Scientific Reports 8, 1837 (2018), www.nature.com/articles/s41598-018-20152-3

Introduction

When increasing the density or decreasing the temperature many particulate systems reach a state where no longer any significant dynamics can be observed such that the system is in an amorphous, effectively solid state. Such a dramatic slowdown of the dynamics has been observed in many systems [1, 2, 3, 4, 5], even for particles with quite simple pair interactions as in the case of colloids [5, 6, 7, 8, 9]. In simulations, a glassy slowdown even occurs in hard spheres systems [10, 11, 12, 13]. The origin of the associated glass transition as well as its fundamental properties are still under discussion (cf. [14, 15]).

A solid amorphous state that on the first glimpse seems to be similar to a glass can be achieved by increasing the density at exactly zero temperature [16, 17]. For example, in the protocol proposed by O’Hern et al. [18, 19] one starts with a random configuration of soft spheres that interact according to a finite-ranged repulsive interaction like a Hertzian or a harmonic potential. Then the local energy minimum is determined, i.e., the energy is minimized without crossing energy barriers. Note that the energy that is minimized is given as sum of all pair interaction energies. Depending on the packing fraction of the system, either all overlaps have been removed, which corresponds to a ground state and is called an unjammed system, or the configuration at the local minimum contains overlapping particles, which is called a jammed configuration. Note that such a jammed configuration obviously is not a ground state and that as a consequence jammed systems usually are not in equilibrium, because in principle ground states might still exist but are just not accessed. In the limit of large system size, jamming occurs at a well-defined packing fraction corresponding to random closed packing [18, 19]. Note that starting with other configurations will lead to an athermal jamming transition with the same scaling behavior but a different transition packing fraction [20, 21]. In this article we explore how the athermal jamming transition changes if energy barriers are crossed due to thermal fluctuations during the quest to approach the ground state by energy minimization. We then study the transition from systems that reach a ground state to systems where the ground state is not accessible. We show that the onset of the effectively non-ergodic behavior is given by a directed percolation transition in time.

A unified jamming phase diagram [22] has been proposed where athermal jamming is the endpoint of the glass transition line, i.e., of the jamming transition at non-zero temperature. Interestingly, recent theoretical studies [23, 24, 25] and simulations [26, 27, 28, 29, 30] suggests that the athermal jamming transition might occur inside the glass phase at small but nonzero temperatures.

In a recent work [31] Morse and Corwin modified the athermal jamming protocols and force the particles to stay in contact during energy minimization. They then observe a percolation transition of clusters formed by locally rigid particles and relate it to the dynamical glass transition. To be specific, Morse and Corwin show for systems in three to six dimensions that the resulting rigidity transition occurs at packing fractions that are close to the packing fractions expected for the dynamical glass transition according to [32], e.g., in three dimensions they find a transition at the packing fraction $0.55(8)$ [31] which is far below the athermal jamming packing fraction $\phi_{J}=0.639$ [19] but close to the packing fraction of the dynamical glass transition at $0.571$ [32]. Morse and Corwin argue that particles that are not locally rigid have more degrees of freedom to explore the configuration space than particles in large clusters of locally rigid particles [31]. An open question remains, namely why should the particles stay in contact during the unjamming procedure? This question will be answered by our finding that due to rare thermal rearrangements the system can effectively be trapped in a region of the configuration space where a significant number of particles stays in contact. As a consequence, the ergodicity-breaking transition studied in this article is directly connected to the percolation transition described by Morse and Corwin.

Model: Testing the accessibility of ground states. In this work we study the competition and interplay of the slow relaxation by energy minimization within an energy basin and the rare crossing of barriers. The minimization part of our protocol is known from athermal jamming [18, 19], the hopping over barriers is often employed in models that are used in order to describe dynamical properties of glasses, e.g., dynamics heterogeneities or ageing [3, 33, 34, 35]. Note that in the later models usually an energy landscape of disconnected basins is assumed where thermalization within a basin takes place almost instantaneously. However for the case that we are interested in, i.e., for densities below athermal jamming, all energy basins are connected in case of an infinite system. The dynamics still can be slow in case of a slow relaxation within this one basin due to narrow and long pathways to the ground state.

If we explored the energy landscape by a Brownian dynamics simulation or a local Monte Carlo simulation, the trajectory in configuration space mainly would fluctuate in a valley and on average goes downhill as far as possible (see sketch in Fig. 1(a)). On rare occasions an energy barrier can be crossed. Since we want to study large systems and long timescales, we consider a simplified approach in order to determine which parts of the configuration space can be accessed within a reasonable time. We start with random configurations of systems consisting of monodisperse spheres in three dimensions and then usually employ the energy minimization protocol of athermal jamming. However, we introduce additional steps that for each particle and each minimization step can occur with a small probability $p$ and can result in the crossing of an energy barrier (see sketch in Fig. 1(b)). We have tested different implementations of the minimization as well as the barrier hopping steps. Details of the protocol that we usually use in the main text are given in the methods section and alternative protocols are discussed in the Supplementary Notes 1. All of these protocols lead to the same ergodicity-breaking transition in the limit of small $p$. Note that close to the observed transition we employ large systems of up to $N=10^{7}$ particles in order to avoid finite size effects that are explored in the Supplementary Notes 2.

The goal of the work is to find out whether the ground state can be reached or not. In analogy to the terminonoly used for athermal jamming, a system is considered to be unjammed if all overlaps can be eliminated. If for a larger packing fraction the number of particles that still possess overlaps is not decreasing, the system is called jammed. Note that we are only interested in the case of rare barrier crossing events, i.e., small $p$-values, where the dynamics is dominated by the minimization process and indeed effectively is stuck, if the system cannot reach the ground state. For large $p$ a significant number of rearrangement occurs due to the barrier crossing events such that instead of a dynamically jammed system one finds a fluid of soft, overlapping spheres. If all steps were random (for $p=1$) our protocol is similar to the one with only random displacements that we studied in [39]. The transition in this case was shown to be in the universality class of directed percolation and can be mapped onto a random organization transition, which was first observed in cyclically sheared colloidal systems [36, 37].

In this article we identify the barrier crossing probability $p$ with temperature $T$ in the sense that for $T=0$ no barriers can be crossed, i.e., $p=0$, that for $T>0$ there is a non-zero probability $p>0$ of barrier crossing, and that with increasing temperature the probability to cross barriers increases as well. We want to point out that in order to obtain an ensemble sampling at a fixed real temperature $T$ of the configurations that we find, one would have to weight the observed trajectories by appropriate factors, i.e., by using Kramers’ rate [38] for each barrier crossing. Note that in the limit of small $p$ that we are interested in only a minority of the particles cross a barrier at all. Therefore, in principle, an ensemble sampling based on our method is possible. Unfortunatally close to the glass transition such a sampling is computationally too expensive because large systems are required to avoid finite size effects. Note that Brownian dynamics or Monte Carlo simulations close to the glass transition also are too demanding if you wanted to study the critical behavior in the same system size as we do. Since our main interest in this paper is to find out, whether the ground state is accessible or not, ensemble sampling or the determination of statistical, non-zero weights for accessible configurations is not necessary, because the accessibility of a configuration would not change if it was weighted differently. Furthermore, any sampling of the power laws that we observe in the following would again lead to power laws with the same exponents. Finally, we show that the transition density that we find does not depend on $p$ in the limit of small $p$. The reason is related to a property of the configuration space, namely a spatial percolation transition that occurs at this density as we will show in the paragraph on spatial percolation and that also has been reported in [31]. If a particle crosses an energy barrier it can affect the minimization process of all particles that are part of a percolated chain of overlapping particles. As a consequence the naive expectation that p directly determines how many particles are disturbed on their way towards a non-overlapping ground state turns out to be wrong in a spatially percolated system. In such a system the whole percolated chain of particles is affected if only one particle of the chain crosses an energy barrier. Therefore, the ergodicity breaking transition then does not depend on $p$ but is given by the packing fraction of the spatial percolation transition. Note that since $p$ is given by a strongly monotonic function of $T$, the ergodicity breaking transition also cannot depend on $T$ for small $T$ no matter how the functional dependence of $p$ on $T$ actually looks like.

Results

Our new method for $p=0$, i.e., without any crossing of energy barriers, leads to the well-known athermal jamming transition at a packing fraction of $\phi_{\textnormal{J}}=0.638$ which is in agreement with the results for a monodisperse system reported in [18, 19]. As we will show in the following, for $p>0$ a different type of transition occurs which we call the thermal jamming transition.

The thermal jamming transition. To get a rough idea of the athermal and the thermal relaxation process, we sketch schematic energy landscapes in Fig. 2. The blue areas correspond to unjammed ground states. In the dilute systems shown in Figs. 2 (e,f) a lot of unjammed configurations occur while for large packing fractions (cf. Figs. 2 (a,b)) overlaps might prevail. In the energy landscapes on the left hand side of Fig. 2 the relaxation process occurs according to the athermal jamming protocol ($p=0$) where the local energy minimum is finally reached. The athermal jamming transition takes place at the packing fraction where the energy of the local minimum changes from zero to a non-zero value.

In case of thermal jamming ($p>0$) depicted on the right hand side of Fig. 2 the crossing of barriers is possible leading to a transition packing fraction that is below the one of athermal jamming.

Fig. 3 shows the number of overlapping particles as a function of the number of steps $t$ for selected $p$ and different $\phi$. For a given $p$ one finds unjammed configurations where the number of overlaps decays to zero at small packing fractions. At large packing fractions, the curves reach a steady state denoting a jammed system. Pair distribution functions $g(r)$ of jammed configurations close to the transition are analyzed in the supplementary note 3. We observe a pronounced peak of $g(r)$ close to $r=\sigma$, which is a known feature of soft sphere glasses close to the hard sphere limit [24, 26].

If one tries to cross barriers in all steps, i.e., for $p=1$ depicted in Fig. 3(a), we observe a transition between $\phi=0.052$ and $\phi=0.056$. This case is similar to the transition studied in [39, 40] which has the same universality class as the random organization transition considered in [36, 37, 41]. Note that it is known that differences in the details of the protocol lead to different transition packing fractions. It was shown that the critical behavior either corresponds to a directed percolation transition either with a conserved number of binding sites particles or to one with an unconserved number [39, 37]. Whether it is a non-conserved or a conserved directed percolation turned out to be hard to determine from an analysis of the critical exponents [39, 37]. Note that binding sites correspond to overlapping particles in our case. Since there is no reason why the number of overlapping particles should be conserved and since a detailed analysis in [39] of the random organization transition indicates that it is more likely to be a non-conserved directed percolation transition, we will compare our observed critical behavior to non-conserved directed percolation in the following.

If the probability $p$ is decreased, the observed transition packing fractions usually increases. However, in all cases with $p\leq 10^{-4}$ (Fig. 3(d-f)) the transition occurs at roughly the same packing fraction between $\phi=0.53$ and $\phi=0.555$. Furthermore, this transition packing fraction is much smaller than the athermal jamming transition packing fraction $\phi_{\textnormal{J}}=0.638$ which is obtained for $p=0$.

Thermal jamming phase diagram. In Fig. 4 we show how the transition packing fraction that separates states leading to non-overlapping configurations from states where overlaps remain depends on the probability $p$. All overlaps can be removed in the blue area while in the red area the overlaps do not vanish. For comparison we also depict where in the case $p=0$ the transition into an athermally jammed state occurs which is marked in yellow. The transition line is determined by the largest $\phi$ denoting a state without remaining overlaps the smallest $\phi$ of all states that possess overlaps in the end for a given $p$ (marked by brown error bars in Fig. 4). Similarly, for given $\phi$ the transition ranges in $p$ are determined (white error bars). In addition we obtained transition packing fractions by analyzing the critical behavior (white solid circles and yellow open square) as we will explain later. We have checked that the shown results are not affected by system size effects (see also Supplemetary Note 2.).

We find that the transition line $\phi_{\textnormal{c}}(p)$ of this transition that we term the thermal jamming transition for $p\rightarrow 0$ approaches a packing fraction $\phi_{\textnormal{G}}=\lim_{p\rightarrow 0}\phi_{\textnormal{c}}(p)=0.55\pm 0.01$ where it denotes a transition between an ergodic state and an effectively non-ergodic state as we have explained before. The packing fraction $\phi_{\textnormal{G}}$ of this ergodicity breaking in the limit of small but non-zero $p$ corresponding to small but non-zero temperatures significantly differs from the athermal jamming packing fraction $\phi_{\textnormal{J}}$. Note that if $p$ is not small the thermal rearrangements lead to a fluidization even in the case where overlaps remain. Therefore, the glass that we are interested in only occurs in the limit of small $p$.

Critical behavior. Close to the thermal jamming transition we analyze the critical behavior. We determine the fraction of overlapping particles in the long-time limit $f_{\textnormal{ov}}(t\rightarrow\infty)$ as well as the relaxation time $\tau$ for the relaxation towards the state that we find for long times as functions of $\phi-\phi_{\textnormal{c}}$.

In order to determine $f_{\textnormal{ov}}(t\rightarrow\infty)$ and $\tau$ we fit

[TABLE]

to the relaxation curves above the jamming transition. This approach assumes a critical power law $f_{\textnormal{ov}}(t)\propto t^{\alpha}$ at the transition density from which the curves at other densities deviate at a time $\tau$. Such an approach already has been used in [42, 39]. As shown in Supplementary Note 3, our results are in agreement with $\alpha=-0.732$ as expected for directed percolation [43]. For the fits to Eq. (1), we use a fixed $\alpha=-0.732$ and otherwise employ $f_{\textnormal{ov}}(t\rightarrow\infty)$, the prefactor $A$, and the relaxation time $\tau$ as fitting parameters. For $\phi<\phi_{c}$ we often observe a power law decay with exponent $-1.5$ (see Supplementary Note 2), which does not influence the transition but makes it hard to define a relaxation time $\tau$ below $\phi_{\textnormal{c}}$.

In Fig. 5(a) $f_{\textnormal{ov}}(t\rightarrow\infty)$ is shown as function of $\phi-\phi_{\textnormal{c}}$. It is zero for unjammed configurations. Above the thermal jamming transition at $\phi_{\textnormal{c}}$ we find that our results can be described by a power law

[TABLE]

with a critical exponent $\beta$. As can be seen in the log-log representation shown in the inset of Fig. 5(a) the exponent $\beta$ does not depend on $p$ and therefore all curves possess the same critical behavior as directed percolation. For comparison, the black line indicates the slope $\beta=0.81$ expected for a directed percolation transition [43]. The packing fractions that are obtained by fitting the power law in Eq. (2) to our simulation data are shown with white solid circles in Fig. 4.

In Fig. 5(b) we show that $\tau$ as a function of $\phi-\phi_{\textnormal{c}}$ obeys a critical power law of the form

[TABLE]

with an exponent $\nu$. The log-log plot in the inset of Fig. 5(b) demonstrates that $\nu$ also does not depend on $p$ and that all simulation results are in agreement with $\nu=-1.1$ (black line) as expected for directed percolation [43]. The packing fractions where we observe the divergence of fitted power laws are shown with yellow open squares in Fig. 4.

In conclusions, the critical behavior for all $p>0$ is the same as directed percolation. Physically this result can be motivated in the following way: A system is thermally jammed if there is a region of overlapping particles that does never disappear as time proceeds though it might move around. This corresponds to a path that describes the propagation of overlaps in time and that is directed percolated (directed because the time is always directed).

The athermal jamming transition is not in the universality class of directed percolation. In fact, it is a very different transition because the number of overlaps per particles jump from [math] to the value needed for isostaticity (e.g., to $4$ in two dimensions and $6$ in three dimensions) [18, 19]. Therefore, even without taking the transition packing fractions into account one can conclude from the critical behavior that the athermal jamming transition cannot be the $p\rightarrow 0$ (or $T\rightarrow 0$) limit of the thermal jamming transition. Note that also the pair correlation function $g(r)$ differs significantly in thermally jammed and athermally jammed packings as we show in Supplementary Note 3.

Spatial Percolation.

We analyze whether a spatial percolation transition occurs for the configurations obtained by the thermal jamming protocol. We consider two particles to be in contact if they touch or overlap. Starting at an arbitrary particle we determine the cluster of particles that can be connected by contacts. In Fig. 6(a,b,c) we show for $p=10^{-2}$, $10^{-3}$, and $10^{-4}$ and various $\phi$ the probability distribution $P(x)$ that a particle still is in this contact cluster if in $x$-direction (or any other given direction) it is at a distance $x$ from the starting particle. For large $\phi$ we observe that the cluster of connected particles reaches through the whole system and therefore there is a continuous percolation transition in space. If we only consider paths within the clusters that are directed in $x$-direction, i.e., if we only go from particle to particle within the cluster if this increases the $x$-coordinate, we obtain a directed percolation transition in space as is shown in Fig. 6(d,e,f).

For the three values of $p$, we mark the studied state points in Fig. 4 with triangles. Spatially unpercolated configurations are denoted by cyan triangles, continuously percolated states by orange triangles, and configurations that are also directed percolated in space by blue triangles. For large $p$ continuous and directed percolation take place deep in the thermally jammed phase. However, for small $p$ the percolation transitions both occur at a similar packing fraction as the thermal jamming transition. Therefore, our results indicates that one of these spatial percolation transition might be the reason why for small $p$ there is no significant increase of the thermal jamming transition packing fraction $\phi_{\textnormal{c}}(p)$ upon a decrease of $p$. As a consequence of the spatial percolation a single randomly moved particle can affect the whole system no matter how large it is and therefore even in the limit $p\rightarrow 0$ the relaxation into an unjammed state can be prevented.

Morse and Corwin in their work [31] find a spatial percolation of locally rigid particles at the same packing fraction up to the precision of our simulations. They claim that this transition is an echo of the dynamical glass transition [31]. As we point out, it is related to the breaking of the effective ergodicity.

Packing fraction of ergodicity breaking for different starting configurations.

We observe that in the limit $T\rightarrow 0$, the thermal jamming transition density where the system becomes effectively non-ergodic is $\phi_{\textnormal{G}}=0.55\pm 0.01$ for monodisperse spheres. To give a few examples from literature for comparison, for an experiment on colloidal suspensions with a small polydispersity $0.05$ a glass transition packing fraction $\phi_{\textnormal{G}}\approx 0.56$ is reported [7], while another experiment with a larger polydispersity indicates $\phi_{\textnormal{G}}\approx 0.58$ [44]. By fitting the power law divergence of the relaxation time predicted by mode-coupling theory [45] to data of colloidal experiments $\phi_{\textnormal{G}}$ is expected in the range $0.571$ to $0.595$ [46, 47]. In computer simulations the dynamics was studied even beyond this prediction up to $\phi\approx 0.6$ [47]. Numerical studies with soft repulsive harmonic spheres predict $\phi_{\textnormal{G}}$ at zero temperature limit at $0.637$ [27, 28] or in the zero shear stress limit at $0.59$ [29]. For emulsions, $0.589$ is reported for experiments [48] and $0.591$ for simulations [48]. Most of these differences probably are due to differences of the systems (e.g., monodisperse vs. polydisperse systems), different methods of extrapolations, or due to difficulties to determine packing fractions with high accuracy in experiments. However, as we will show in the following, different transition packing fractions might also arise due to different starting configurations.

Instead of the random initial configurations that we have used so far, for the results shown in Fig. 7 we employ athermally unjammed starting configurations that are obtained by the deterministic minimization protocol. Note that we employ the usual athermal jamming protocol and not the modification used in [31], i.e., particles are not kept in contact at the end of the minimization. Furthermore, in the following only overlapping but not touching particles are considered for the analysis such that for our athermally unjammed starting configurations there is no spatially percolation for $\phi<\phi_{\textnormal{J}}$ while the thermally jammed configurations for small $p$ and $\phi>\phi_{\textnormal{G}}$ possess spatially percolated clusters as shown in the previous subsection.

In Fig. 7(a) we show the number of overlapping particles that we observe in case of an athermally unjammed starting configuration and for the athermal protocol after a single particle is moved over a barrier for $\phi=0.58$. The different curves correspond to different realizations. Note that the system always ends up in an unjammed configuration and the number of overlapping particles always remains finite (and much smaller than the system size). In order to understand how further crossings of barriers would affect such a system (e.g., for $p>0$) we determine how many overlapping particles exist integrated over time. The value $I$ of this integral indicates the total number of overlapping events. In Fig. 7(b) the probability distribution $P(I)$ of how often a given value $I$ occurs is plotted for various packing fractions ranging from $0.55$ to $0.61$. In none of the cases an infinite $I$ was observed. However, when increasing the packing fraction the $P(I)$ are shifted to larger $I$ and they become narrower. In the inset of Fig. 7(b) the mean values of $I$ are shown as a function of the packing fraction. Interestingly there is an exponential increase.

Next, we want to find out whether we obtain a thermally jammed state if the thermal protocol with a $p>0$ is switched on. Fig. 7(c) shows the number overlaps as function of time after the initial barrier crossing for various realizations for $p=10^{-5}$ and $\phi=0.555$ where in case of a random initial configuration a thermally jammed state occurs. However, if started from an athermally unjammed configuration the system is not in a thermally jammed because for all realizations we find that the system relaxes into an unjammed configuration. Only at a larger density, e.g., for $\phi=0.58$ as plotted in Figs. 7(d), most realizations end up in a jammed state. Therefore, in case of an athermally unjammed initial configuration and for $p=10^{-5}$ a thermal jamming transition can be observed which occurs above $\phi=0.555$ but still well below $\phi_{\textnormal{J}}$. Note that we expect this apparent ergodicity breaking transition to depend on $p$ even in the limit $p\rightarrow 0$ because there is no underlying spatial percolation transition. However, the dependence on $p$ will be weak (probably logarithmic) because a small increase in packing fraction results in an exponential increase in overlapping events as shown above. Finally, we want to point out that in order to study the hard-sphere limit it might be more natural to start with the non-overlapping configurations as considered in this subsection instead of random configurations that were employed in the previous sections and that can contain large overlaps. The behavior observed for small overlaps then in principle can be mapped onto hard sphere system [49, 50].

Discussion

By adding the rare possibility to cross energy barriers to the protocol that previously was employed to study the athermal jamming transition, we obtain a powerful new method that allows the direct investigation of jamming at small but non-zero temperatures. By employing this method we determined whether the system can access the ground state or not. In the later case the system effectively is non-ergodic for small temperatures. Therefore, the observed transition is an effective ergodicity breaking transition.

We find that the ergodicity breaking transition is a directed percolation transition in time. The transition occurs at much smaller packing fractions than the athermal jamming transition and the transition density in the limit of small temperature does not depend on the temperature but is is given by the spatial percolation of particles in contact.

The specific value of the transition depends on the initial conditions. As we explained in our discussion of different starting configurations, past experiments and simulations used various protocols which might be one reason why different transition densities for the glass transition have been reported. Furthermore, as we explain in Supplementary Note 2, in too small systems the apparent ergodic to non-ergodic transition might be observed at a packing fraction that is larger than in our large system. However, for given initial conditions and large enough systems, our method can be used to directly predict the density of the ergodic to non-ergodic transition. For example, in case of a a very fast quenches from infinite to small final temperatures we find a transition at a packing fraction of 0.55 exactly as Morse and Corwin in [31].

Our protocol is constructed such that we can easily find out whether a ground state can be reached or not. Note that Brownian dynamics simulations or local Monte-Carlo simulations are superior for simulating fluctuations within a valley of the landscape or for producing an enable sampling of visited configurations. However, they are not superior in minimization and they all allow for the (maybe rare) crossing of barriers in case of non-zero temperature. As a consequence, Brownian dynamics or local Monte-Carlo simulations cannot reach a ground state that is not even accessible with our optimized approach. Note that while ensemble sampling with our method might require the use of additional weight factors, the power laws describing the critical behavior as well as the power law behavior for $g(r)$ shown in Supplementary Note 3 would not change if the resulting configurations had to be weighted differently in an ensemble sampling.

In recent years, there have been significant advances in studying and characterizing non-ergodic systems including those associated with anomalous diffusion (for reviews see, e.g., [51, 52]). Note that in case of the dynamical soft sphere glass transition that we study there might be rare random rearrangements even in the effectively non-ergodic glass phase. Due to these rearrangements the system is diffusive in the long-time limit. However, the diffusive motion occurs on a timescale that is longer than the timescale of observation and therefore the system is termed an effectively non-ergodic glass. In the language of stochastic dynamics such systems are termed weakly non-ergodic [52]. With our approach in case of a small but non-zero probability $p$ for energy barrier crossings, the rearrangements are directly associated with such rare barrier crossing events. The timescale of the rearrangements is given by $1/p$ in units of minimization steps. The relaxation time typically is a multiple of this timescale. As a consequence the Deborah number $De$ that is defined as ratio of the timescale of relaxation and the timescale of observation $T$ [53] has to be $De>1/(pT)$, i.e., $De\rightarrow\infty$ for $p\rightarrow 0$. For example, in our case for $T\sim 10^{5}$ and $p=10^{-6}$ one finds $De>10$. Since the system is ergodic and diffusive in the long time limit, asymptotically the time average and the ensemble average of the squared displacement are the same. A corresponding weakly non-ergodic system is given by a Brownian particle that moves in a random, bounded external potential. For such a system the relation between time and ensemble average at finite times has been extensively studied in simulations as well as experiments [54, 55].

We are confident that the new method can be employed to obtain more insights into the physics of glassy states, e.g., in order to study the properties of the modes in a thermally jammed but athermally unjammed state point (cf. [30]), to determine the connection to rigidity percolation [56] and contact percolation [57], to understand a possible Gardner transition [58, 59, 60], to explore the properties of the basins in the energy landscape which can be used to determine the entropy of the system [61, 62], or extend the test of Edwards’ approach to the statistics of granulates [63, 64] that recently has been tested for unjammed configurations [65]. Furthermore, it would be interesting to combine our method with alternative packing protocols [66], or to study the influence of shearing as in [41], to explore the glass transition for active particle [67, 68], and to learn more about gelation in case of particles with short ranged attraction where recently a connection to spatial directed percolation has been reported [69].

Methods

In case of athermal jamming we start with a random configuration of spheres with diameter $\sigma$ in three dimensions in a cubic box with side length $L$ and periodic boundary conditions. We instantaneously quench the system to zero temperature by minimizing the total energy. The interaction energy is given by the finite-ranged repulsive harmonic pair potential, which is $V(r)=\epsilon(1-r/\sigma)^{2}$ for particle distances $r<\sigma$ and zero otherwise. The prefactor $\epsilon$ sets the energy scale. The packing fraction is given by $\phi=\pi\sigma^{3}N/(6L^{3})$. The energy landscape depends on the position of all particles and is the sum over all pair interaction energies. For minimization we use the conjugate gradient method of the simulation package LAMMPS [70]. The number of minimization steps is denoted by $t$. The minimization is stopped either if the number of overlapping particles as a function of $t$ fluctuates around a plateau value (denoting a jammed state) or if the energy per particle is $10^{-16}\epsilon$ or less (denoting an unjammed state). Accordingly, we consider two particle to overlap if $\sigma-r>10^{-7}\sigma$. We have checked that this energy cut-off and the precision of the overlap definition do not have an influence on our results.

In order to study jamming at a non-zero temperature we employ the same system but in addition particles that are still part of the process, i.e., that are either overlapping or touching (up to our precision, $r-\sigma<10^{-7}\sigma$) are selected in each step with a probability $p$. For the selected particles a random spatial direction is chosen ant the particles are displaced in that direction until they reach a minimum or maximum of the total energy determined along the line in that direction. In the latter case particles are set slightly behind the maximum such the energy barrier is crossed. The LAMMPS software is modified to adapt it for our modified protocol. Especially, the search direction to the gradient direction during the relaxation are reset for the randomly displaced particles. After selected particles have been displaced, all particles are moved according to the minimization protocol in the same step. The stopping criterion for the modified protocol is the same as for the athermal jamming protocol.

In case of small $p$ large systems have to be considered in order to be sure that the results are not affected by system size effects. For example, for $p=10^{-6}$ close to the transition we employ systems with up to $N=10^{7}$ particles. A detailed analysis of system size effects is presented in Supplementary Note 2. We have checked that no crystallization occurs by analyzing the bond orientational order parameter $Q_{6}$ as described in [71]. Furthermore, we have tested different implementations for the minimization as well as for the random steps (see Suplementary Note 1). Such details affect how $p$ is defined, but not the transition density in the limit of small $p$ or the type of the transition anywhere.

Data availability. The data shown in this paper or the supplementary notes are available from the authors.

Supplementary Notes on “Ergodicity breaking transition in a glassy soft sphere system at small but non-zero temperatures”

**Moumita Maiti1 and Michael Schmiedeberg1,∗

1**Institut für Theoretische Physik 1, Friedrich-Alexander-Universität Erlangen-Nürnberg, Staudtstr. 7, 91058 Erlangen, Germany

∗michael.schmiedeberg@fau.de

Supplementary Note 1:

Modifications of the protocol

In this section we present different models, where we vary the protocol that is employed to access the ergodicity breaking transition in the limit of rare barrier crossing events. In Fig. 8(a) we summarize the ranges that are obtained for the transition packing fraction. The curves for the fraction of overlaps as a function of step number for cases close to this transition are displayed in Fig. 8(b,c,d,e). The employed models are

• •

Model A: The model that is used in the main text and that is explained in detail therein.

• •

Model B: In this variation of model A we only allow overlapping particles to take part in the protocol. Therefore, touching particles can no longer be randomly displaced.

• •

Model C: Variation of model A, where in a random step particles are only displaced if in the chosen direction there is a crossing of an energy barrier. Otherwise the particles is not moved. Therefore, minimization only takes place during minimization steps.

• •

Model D: The minimization is done according to the steepest descent algorithm. Otherwise the protocol of model A is used. Note that minimization using the steepest descend algorithm usually is less efficient than the conjugate gradient algorithm that is used in models A to C. As a consequence, the relaxation just below the transition is very slow and noisy as can be seen for the orange curve for $\phi=0.53$ in Fig. 8(d). As a consequence, in order to differentiate between relaxing and non relaxing systems close to the transition, we also analyzed the structure. In case of a persisting percolated cluster we are sure to be in the non-ergodic phase, while below the transition the system disassociates into distinct non-percolating clusters of overlapping particles that slowly shrink.

• •

Model E: In this model completely different protocols for minimization and random steps are used. The minimization steps corresponds to a deterministic athermal jamming step as used in [39]. In such a step two overlapping particles are deterministically displaced along the line given by the two center points of the particles until the particles are no longer overlapping but only touching. In random steps that occur with probability $p$ the touching particles in addition are rotated around the center of mass by an random angle. More details are given in [39], where we studied the individual protocols but not their mixing. Note that the minimization in this model is very ineffective such that no large systems can be studied within a reasonable computation time. For the shown results, we averaged over 25 runs with $N=25000$ particles. As a consequence, the drop after a plateau has been reached for packing fraction 0.55 in Fig. 8(e) probably is artifact due to the small system size. Therefore, in principle it is possible that the error bars shown in Fig. 8(a) for model E differ from error bars obtained in larger systems. However, up to the system size that we could consider, even in this completely different model the plateaus start to develop at the same packing fraction as for the other models.

Concerning finite size effects, we want to point out that slower minimization protocols might require to go to smaller values of $p$ in order to reach the limit where the transition packing fraction no longer depends on $p$. As a consequence, such protocols might require larger system sizes to obtain a comparable accuracy for the transition packing fraction. In general, the width of the range for the transition packing fraction is mainly limited due to finite size effects. In model A we performed an extensive finite size analysis (cf. Supplementary Note 2) and as a consequence can narrow the range significantly. For models B and C we checked two different system size to show that the given range for the transition packing fraction is not affected by the system size. In case of model E we have to point out that significant finite size effects are still present. Especially the result for $p=10^{-8}$ is probably strongly affected by the system size.

In summary, in the limit of rare barrier crossing events, all models lead to a range for the transition packing fraction that is in agreement with the transition packing fraction determined by model A in the main text.

Supplementary Note 2:

Analysis of system size effects

In the following we show how the relaxation process might be affected by finite size effects. Furthermore, at the end of this section we explain how we choose our system sizes such that the transitions that we determine do not depend on these sizes.

Fig. 9(a) demonstrates that below the transition packing fraction smaller systems relax sooner than larger systems. For the small probability $p=10^{-5}$ we point out to another interesting observation: We find that the relaxation obeys an unexpected power law with an exponent $-1.5$. Note that this power law only occurs below the transition and is not related to the critical power law decay $t^{\alpha}$ with $\alpha=-0.732$ which occurs directly at the transition.

For the relaxation above the transition, we depict different system sizes for $p=0.1$ in Fig. 9(b) and for $p=10^{-4}$ in Fig. 9(c). For small systems the fraction of overlapping particles might suddenly decay to zero instead of approaching the plateau value that should occur in the long-time limit. Therefore, in small systems it is easier to find the unjammed configuration while in large systems the zero energy states are not reached due to the larger configuration space. As a consequence of the system size analysis one can conclude, that if too small systems were considered, one would observe a ergodicity breaking only at larger densities. For example, in the limit $p\rightarrow 0$ the apparent packing fraction of the ergodic to non-ergodic transition would be much larger than $0.55$ if systems with $10^{3}$ to $10^{4}$ particles were considered as it is the case in many simulations of glassy dynamics.

In Fig. 9(d) we also show a case for $p=10^{-5}$ where an unexpected jump at intermediate system sizes occurs (red curve for $N=10^{6}$), where the relaxation first seems to decay rapidly but then fluctuates around a non-zero plateau value. Such a behavior is only observed in a few cases at small $p$ and close to the transition.

In order to determine the thermal jamming transition independent of the system size we use the following recipe: The system size for a given $p$ at first is chosen from the criteria $Np=1$ such that at each step one random move is performed, then we increase $N$ until for the curves that are closest to the transition we do not observe any significant change for two different system sizes.

The system sizes that we determined by this method and that we actually used for the curves shown in Fig. 2 are $N=50000$ for $p=1$, $N=10^{5}$ for $p=0.1$, $N=10^{6}$ for $p=10^{-3}$ and $p=10^{-4}$ as well as for $p=10^{-5}$ below the transition, $N=5\times 10^{6}$ for $p=10^{-5}$ above the transition, and $N=8\times 10^{6}$ for $p=5\times 10^{-6}$. In order to further narrow the ranges in which the transition occurs, we sometimes employ even larger systems in Fig. 3 and Fig. 4 with up to $N=10^{7}$ particles, e.g., in case of $p=10^{-6}$.

Note that the dependence on system size is most prominent close to the transition. As a consequence the system size in principle might limit how close we can get to the transition without being affected by finite size effects. For the resolution of $\phi$ that we used in our article this usually is not a problem. Only for $p>10^{-5}$ we were not able to further narrow the range of the transition packing fraction, e.g., in Fig. 4 of the main article we cannot report any plateau values $f_{\textnormal{ov}}(t\rightarrow\infty)$ close to [math] for $p>10^{-5}$ because only data that has proven to be independent of the system size is shown. Furthermore, for $p\geq 10^{-5}$ we are not able to determine values of $\tau$ that do not depend on the choice of the fitting range or how the data points of the relaxation curve are spaced. The reason is that for small $p$ we are not able to get close enough to the transition such that the range of the power law behavior that is assumed for the fitting (see also Fig. 10(a)) is too small in order to obtain a reliable fit. Therefore, there are no $\tau$-values for $p\geq 10^{-5}$ in Fig. 4(b) but only values for $f_{\textnormal{ov}}(t\rightarrow\infty)$ in Fig. 4(a) which is still independent of the way of fitting. If one wanted to obtain values of $\tau$ for smaller $p$ one would need larger systems closer to the transition such that the power law behavior would be sufficiently extended.

Supplementary Note 3:

Close to the thermal jamming transition: Relaxation curves and pair correlation functions

The fraction of overlapping particles $f_{\textnormal{ov}}(t)=N_{\textnormal{ov}}(t)/N$ as a function of step number $t$ decreases like a power law $t^{\alpha}$ close to the thermal jamming transition. Close to the transition, the relaxation curves can be fitted with the function given in Eq.(2) where we assumed $\alpha=\frac{\beta}{\nu}=-0.732$ as for directed percolation. The quality of this choice of functional behavior is demonstrated in Fig. 10(a) where we collapse our data for the relaxation curves with different $p$ and for $\phi$ just above the transition by plotting $F=\tau^{0.732}(f_{\textnormal{ov}}-f_{\textnormal{ov}}(t\rightarrow\infty))/A$ as a function of $\frac{t}{\tau}$. All curves are in good agreement with $(t/\tau)^{-0.732}\exp\left(-t/\tau\right)$ (magenta line) as expected from Eq.(2).

In Figs. 10(b) and (c) we plot the pair correlation function $g(r)$ for thermally jammed structures with $\phi$ just above the thermal jamming transition and for various $p$. In case of the non-overlapping part shown in Fig. 10(c), between the contact distance and the first minima of $g(r)$ we observe a power law decay with an exponent $-0.75$. The exponent is robust for all $p$. For $p\rightarrow 0$ the power law extends closer to the contact point. Note that for athermal jamming a power law of $g(r)$ is also observed [72]. However, in case of athermal jamming the exponent is $-0.5$ [72] which significantly differs from the exponent $-0.75$ that we observe here for all values of $p$. Therefore not only the nature of the thermal jamming transition is fundamentally different from the athermal jamming transition as demonstrated in the main text but also for the structural properties it makes a huge difference whether a system is athermally jammed or thermally jammed.

Acknowledgements

The work was supported by the German Research Foundation (DFG) by grant Schm 2657/3. We thank Lars Milz for helpful discussions.

Author contributions statement

M.M. carried out the simulations and M.S. designed the research and the model system. Both authors analyzed the results and wrote the article.

Additional information

Accession codes (where applicable).

Competing financial interests: The authors declare no competing financial interests.