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

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 studied in great detail. Motivated by the exploration of the energy landscape that has been successfully used to investigate athermal jamming, we introduce a new method that includes the possibility of the thermally excited crossing of energy barriers. We then determine whether the ground state configurations of a soft sphere system are accessible or not and as a consequence whether the system is ergodic or effectively non-ergodic. Interestingly, we find an transition where the system becomes effectively non-ergodic if the density is increased. The transition density in the limit of small but non-zero temperatures 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, we show that the ergodicity breaking transition 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.

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 studied in great detail. Motivated by the exploration of the energy landscape that has been successfully used to investigate athermal jamming, we introduce a new method that includes the possibility of the thermally excited crossing of energy barriers. We then determine whether the ground state configurations of a soft sphere system are accessible or not and as a consequence whether the system is ergodic or effectively non-ergodic. Interestingly, we find an transition where the system becomes effectively non-ergodic if the density is increased. The transition density in the limit of small but non-zero temperatures 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, we show that the ergodicity breaking transition 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.
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 φ 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 36 . 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 37,38 .
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 39 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 φ 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 Fig. 2(e,f) a lot of unjammed configurations occur while for large packing fractions (cf. Figure 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. Figure 3 shows the number of overlapping particles as a function of the number of steps t for selected p and different φ. 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 = σ, 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 φ = 0.052 and φ = 0.056. This case is similar to the transition studied in 36,40 which has the same universality class as the random organization transition considered in 37,38,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 36,38 . 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 36,38 . 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 36 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 ≤ 10 −4 (Fig. 3(d-f)) the transition occurs at roughly the same packing fraction between φ = 0.53 and φ = 0.555. Furthermore, this transition packing fraction is much smaller than the athermal jamming transition packing fraction φ 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 φ denoting a state without remaining overlaps the smallest φ of all states that possess overlaps in the end for a given p (marked by brown error bars in Fig. 4). Similarly, for given φ 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 φ c (p) of this transition that we term the thermal jamming transition for p → 0 approaches a packing fraction p lim ( ) 0 55 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 φ 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 φ 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 ov (t → ∞) as well as the relaxation time τ for the relaxation towards the state that we find for long times as functions of φ − φ c .
In order to determine f ov (t → ∞) and τ we fit ov ov to the relaxation curves above the jamming transition. This approach assumes a critical power law f ov (t) ∝ t α at the transition density from which the curves at other densities deviate at a time τ. Such an approach already has been used in 36,42 . As shown in Supplementary Note 3, our results are in agreement with α = − 0.732 as expected for directed percolation 43 . For the fits to Eq. (1), we use a fixed α = − 0.732 and otherwise employ f ov (t → ∞), the prefactor A, and the relaxation time τ as fitting parameters. For φ < φ 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 τ below φ c . In Fig. 5(a) f ov (t → ∞) is shown as function of φ − φ c . It is zero for unjammed configurations. Above the thermal jamming transition at φ c we find that our results can be described by a power law.
ov c with a critical exponent β. As can be seen in the log-log representation shown in the inset of Fig. 5(a) the exponent β does not depend on p and therefore all curves possess the same critical behavior as directed percolation.    . Thermal jamming phase diagram. Thermal jamming phase diagram showing the transition between states that can lead to non-overlapping configurations (blue) and states where overlaps remain (red) depending on the probability p of random steps (corresponding to a temperature T) and the packing fraction φ. Note that we are especially interested in the case of small p where the blue area denotes ergodic, unjammed states while the red area corresponds to thermally jammed states that are effectively non-ergodic. The yellow area marks the packing fractions where the system would be athermally jammed in case of p = 0. The white and brown bars denote the range in which the transition occurs in case φ or p is kept constant, respectively. In addition we show the transition packing fractions that we obtain from fits of critical power laws to the steady state values f ov (t → ∞) in case of jammed states (white solid circles, cf. Fig. 5(a)) or to the relaxation times τ (yellow open squares, cf. Fig. 5(b)). Finally, the triangles indicate where we explore whether chains of touching or overlapping particles in the final configurations are unpercolated (cyan triangles), continuously percolated in space (orange triangles), or spatially directed percolated (blue triangles). These spatial percolation transitions are analyzed in Fig. 6. For comparison, the black line indicates the slope β = 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 τ as a function of φ − φ c obeys a critical power law of the form with an exponent ν. The log-log plot in the inset of Fig. 5(b) demonstrates that ν also does not depend on p and that all simulation results are in agreement with ν = − 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 0 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 → 0 (or T → 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 φ 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 φ 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 transitions might be the reason why for small p there is no significant increase of the thermal jamming transition packing fraction φ 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 → 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 → 0, the thermal jamming transition density where the system becomes effectively non-ergodic is φ G = 0.55 ± 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 φ G ≈ 0.56 is reported 7 , while another experiment with a larger polydispersity indicates φ G ≈ 0.58 44 . By fitting the power law divergence of the relaxation time predicted by mode-coupling theory 45 to data of colloidal experiments φ 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 φ ≈ 0.6 47 . Numerical studies with soft repulsive harmonic spheres predict φ 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 Figure 7. Starting configuration dependence of glass transition density. Analysis of how thermal moves affect systems that are athermally unjammed, i.e., in a local energy minimum determined with the athermal jamming protocol. In (a,b) we show how such a system relaxes if one particle crosses a barrier In (c,d) one particle is moved over a barrier and then the thermal jamming protocol is switched on with p = 10 −5 . In (a,c,d)  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 φ < φ J while the thermally jammed configurations for small p and φ > φ 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 φ = 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. Figure 7(c) shows the number overlaps as function of time after the initial barrier crossing for various realizations for p = 10 −5 and φ = 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 φ = 0.58 as plotted in Fig. 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 φ = 0.555 but still well below φ J . Note that we expect this apparent ergodicity breaking transition to depend on p even in the limit p → 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 → ∞ for p → 0. For example, in our case for T ~ 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-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 σ 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) = ε(1 − r/σ) 2 for particle distances r < σ and zero otherwise. The prefactor ε sets the energy scale. The packing fraction is given by φ = πσ 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 ε or less (denoting an unjammed state). Accordingly, we consider two particle to overlap if σ − r > 10 −7 σ. 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 − σ < 10 −7 σ) 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.