Non-Equilibrium Phase Transition in an Atomistic Glassformer: the Connection to Thermodynamics

Tackling the low-temperature fate of supercooled liquids is challenging due to the immense timescales involved, which prevent equilibration and lead to the operational glass transition. Relating glassy behaviour to an underlying, thermodynamic phase transition is a long-standing open question in condensed matter physics. Like experiments, computer simulations are limited by the small time window over which a liquid can be equilibrated. Here we address the challenge of low temperature equilibration using trajectory sampling in a system undergoing a nonequilibrium phase transition. This transition occurs in trajectory space between the normal supercooled liquid and a glassy state rich in low-energy geometric motifs. Our results indicate that this transition might become accessible in equilibrium configurational space at a temperature close to the so-called Kauzmann temperature, and provide a possible route to unify dynamical and thermodynamical theories of the glass transition.


I. INTRODUCTION
Although statistical mechanics was firmly established more than a hundred years ago [1], simple liquids remain a persistent challenge when cooled to low temperatures. In particular, the dramatic super-Arrhenius increase of the relaxation time so far eludes a generally accepted explanation. A multitude of theoretical approaches have been advanced [2,3], but obtaining data that enable discrimination between these is challenging, not least because of the difficulties in handling the huge time scales required to equilibrate supercooled liquids. At some point (typically bypassing crystallization), the structural relaxation time τ α exceeds the experimentally or numerically accessible time scale, and the liquid falls out of equilibrium into a dynamically arrested state called a glass. This so-called operational glass transition is protocol dependent and distinct from equilibrium thermodynamic phase transitions.
However, the idea that at very low temperatures a genuine thermodynamic phase transition controls dynamic arrest has been around for a long time, starting with an observation by Kauzmann [4]: Extrapolating the configurational entropy of the liquid suggests that it should fall below that of the crystalline solid at a finite (Kauzmann) temperature T K . One resolution of this apparent paradox is to posit a thermodynamic transition to an "ideal glass" with very low configurational entropy. This transition is hard to access with the operational glass transition intervening because of the accompanying divergence of the structural relaxation time [5].
More recent theories of the glass transition continue to use the language of phase transitions although they disagree on even the most basic assumptions. On one hand, the emergence of slow dynamics in a rugged free-energy landscape is anticipated through freezing, not to a single crystal but to a vast number of random, aperiodic states [6]-a picture that finds justification from mean-field results obtained in higher dimensions [7,8] and recent experiments [9]. Approaches based on replica symmetry breaking also imagine a phase transition to a state where configurations are closely related to one another; i.e., they have high overlap [10]. Another example, geometric frustration, posits that the population of geometric motifs that minimize the local free energy (locally favored structures-LFS) strongly increases upon supercooling and that the arrest seen is a manifestation of an avoided phase transition. Approaching the ideal glass transition, the system may seek to minimize its entropy, which can correspond to an increase in LFS [11].
Quite in contrast, in dynamical facilitation theory [12], dynamic correlations are the fundamental mechanism for the slow-down, with configuration-based static correlations being absent or at least irrelevant. Originally developed with lattice models, evidence for the dynamical facilitation approach has been obtained in experiment [13][14][15]. Like the thermodynamic approaches noted above, the notion of a phase transition is pivotal, but now the transition is between two dynamic regimes: the supercooled liquid and a state with extremely slow dynamics [16]. The mathematical framework is that of statistical mechanics combined with large deviations [17]. The transition shares many similarities with the disorder-order transition of an Ising magnetthe crucial conceptual differences being that configurations are replaced by trajectories (i.e., time sequences of configurations) and that the order parameters consider the time-integrated quantities such as particle motion instead of magnetization. That is to say, the dynamical transition we consider (via a generalized or "dynamical" chemical potential) does not explicitly pertain to the motion of the particles but rather to a phase transition in trajectory space, based on time-integrated quantities. What is meant by dynamical is that we consider time-integrated quantities over trajectories, but those quantities themselves may be based on static information, derived from configurations along the trajectory. An intriguing prediction is the coexistence between the two dynamical regimes, which, in lattice models, is predicted to terminate at two critical points (one at high temperatures and one at low temperatures) [18]. While numerical evidence supports coexistence and a first-order nonequilibrium transition in trajectory space [16,19], the putative lower critical point remains out of reach for direct numerical investigations in atomistic models. Whereas the operational glass transition is a consequence of the exceedingly long relaxation times of deeply supercooled liquids, the nonequilibrium phase transition in trajectory space is a genuine phase transition, originating from nonanalyticities in the derivatives of the nonequilibrium equivalent of free energies (i.e., the large deviation functions of time-integrated observables [17]). Recently, such dynamical phase transitions have become accessible to particle-resolved experiments [14,15].
Returning to the idea of a thermodynamic transition, a key aspect is some kind of structural change. While such structural changes appear to be minor when looking at two-point measures like the structure factor, higher-order measures reveal a richer behavior [20][21][22][23][24][25]. In particular, the population of geometric motifs, so-called LFS, shows a strong temperature dependence [20]. Indeed, the nonequilibrium transition to a state with very low particle mobility can be driven by increasing the LFS population [26], indicating that these structural changes are not only a by-product of cooling but also play a crucial role in the dynamic arrest. This finding is corroborated by simulations showing that single-particle motion and LFS are correlated [27].
Numerical simulations are an indispensable tool to gain microscopic insights. However, accessible time scales are still 9-10 decades away in relaxation time from the operational glass transition temperature, and more indirect methods have to be devised to gain insight into the nature of the glass transition. For example, pinning (or confinement) of particles [28][29][30] shifts the putative thermodynamic transition into the time window accessible to computer simulation, as does the observation of distributions of overlaps in configurations of particles [31]. Another method to generate deeply supercooled, equilibrated configurations is that of particle swaps [32,33].
Here, we introduce a new and powerful numerical method to tackle the challenge of determining the lowtemperature equilibrium behavior of a popular atomistic glassformer, the Kob-Andersen binary mixture [34]. This model is loosely based on the metallic glassformer Ni 80 P 20 . Unlike previous approaches, our method enables a connection between nonequilibrium phase transitions (i.e., the approach of dynamic facilitation [12]) with more structurebased theoretical approaches which assume a thermodynamic phase transition. In particular, we use this method to extract configurations with exceptionally low configurational entropy and energy.
Our numerical method exploits three effects: First, the central configurations of biased trajectories contain many locally favored structures. Second, a higher population of LFS facilitates the sampling of configurations with low energies even at moderately supercooled temperatures. We develop methods to remove the simulation bias, which allows us to access not only the simulated state point but also an extended temperature range in its vicinity. The third effect is that the contribution to the total free energy from the local minima in the potential energy landscape and from the vibrational free energy decouple [35]. This allows us to access equilibrium properties at configurational temperatures T (as opposed to the conventional thermodynamic temperature) different from the sampling temperature T s .
Using the described method (along with particle swaps [33,36] adapted to the binary Kob-Andersen system with "ghost insertion" [37]), we perform extensive biased simulations of the Kob-Andersen mixture, from which we construct part of the phase diagram for the time-integrated first-order transition from a LFS-poor phase (the normal liquid) to a LFS-rich phase. We find numerical evidence that coexistence between these phases is terminated at a finite temperature, implying a lower critical point eventually accessible to the equilibrium system-were it cooled sufficiently slowly. This scenario has significant consequences because the dynamical transition at a higher temperature allows us to probe glassy configurations that are typical of much lower temperatures while circumventing the prohibitive increase of equilibration times. We are thus able to probe the low-temperature fate of this model glassformer.

A. Model
The Kob-Andersen binary mixture [34] consists of 80% large particles and 20% small particles interacting through truncated and shifted Lennard-Jones pair potentials. We employ the original potential parameters. All numerical values are reported in Lennard-Jones units with respect to the large particles, and we set Boltzmann's constant to unity. Simulations are performed for a system of N particles at number density N=V ¼ 1.2 in a periodic box with constant volume V. We employ the Andersen thermostat at sampling temperatures T s ¼ 0.74, 0.73, 0.7, 0.6, 0.55, 0.5, 048. Newton's equations of motion are solved using Verlet's velocity algorithm with time step 0.005 [38]. Trajectories are then stored as discrete sequences X ¼ fC −K=2 ; …; C 0 ; …; C K=2 g of configurations for lengths K ¼ 60 and K ¼ 100. The time between successive configurations is chosen such that trajectories have a physical duration of t obs ≈ 4.5τ α ; 7.5τ α , respectively. Upon cooling, the structural relaxation time increases faster than exponentially (super-Arrhenius behavior)-see Fig. 1(a). Whether the relaxation time diverges at a finite temperature is still debated [39,40].

B. Order parameters
In the following, we focus on two order parameters characterizing the liquid. First, the time-integrated order parameter characterizes trajectories by counting the total number of particles in locally favored structures, where nðCÞ is the fraction of particles in LFS in configuration C. For the specific model considered here, the LFS corresponds to a bicapped square antiprism formed by 11 particles as sketched in the inset of Fig. 1(a). We detect this motif by employing the topological cluster classification method [23,41]. The temperature dependence of the LFS population is fitted with the empirical Fermi function with fitted exponent α ≃ 2.5 and temperature T 1=2 ≃ 0.25 at which the population extrapolates to 1=2. The increase in LFS concentration [cf. Fig. 1(a)] has been shown to correlate well with the propensity of particle mobility on lowering the temperature [27]. Steepest descent quenches of the central configurations C 0 to the local minima of the energy landscape yield the socalled inherent states [42], which may be thought of as particle configurations with the thermal noise removed. These inherent states C is 0 constitute the local minima in the energy landscape around which particles vibrate before relaxing to another local minimum. While in a perfect crystal there would be only one inherent state, in the liquid there are many different particle arrangements, and the number of these accessible amorphous inherent states defines the configurational entropy.
Our second, static order parameter is thus the inherent state energy (ISE) per particle ϕ½X¼UðC is 0 Þ=N of the central configuration of trajectory X, where UðCÞ is the potential energy of particle configuration C and C is is the corresponding inherent state. This is distinct in nature from the timeintegrated order parameter in Eq. (1).
Practically, the ISE ϕ of central configurations C 0 is obtained using the FIRE algorithm [43] limited to 1000 iterations in order to make the generation of long sequences of trajectories computationally feasible. Ground states of equivalent energies have been obtained by employing basin hopping techniques [44], confirming that the reported values for ϕ should be understood as upper (although tight) bounds to the true inherent state energies. C. Biased simulations We go beyond the temperature regime that is accessible in conventional molecular dynamics simulations. To this end, biased simulations with N ¼216 and N ¼ 400 particles are run at a moderately supercooled sampling temperature T s to faster explore phase space while at the same time sampling configurations with very low potential energy and a high population of LFS representative of low temperatures. To improve the sampling, we employ replica exchange [19]. We simultaneously extract time-integrated and static information (central configuration) by harvesting trajectories of length t obs , which we choose to be a few structural relaxation times t obs ≈ 4.5τ α and t obs ≈ 7.5τ α for runs with K ¼ 60 and K ¼ 100 configurations, respectively.
We employ the multistate Bennet acceptance ratio estimator [45,46] in order to calculate expectations and distributions from the biased numerical data. For an arbitrary observable A½X, this amounts to evaluating the expression where h·i indicates the average over the sampled trajectories at sampling temperature T s . Here, T denotes the target configurational temperature and μ the dynamical chemical potential.
In practice, our numerical approach allows us to improve the sampling of trajectories (and configurations) that are rare for a given sampling temperature T s , explicitly favoring trajectories with exceptionally large (or small) overall concentrations of LFS n by performing importance sampling in trajectory space. Transition path sampling [47] is performed according to the structural bias in n; however, any observable (including both ϕ and n) can be simultaneously tracked, and its correct probability distribution, and, in particular, its mean value, can be recovered by reassigning the correct weight to each trajectory, as illustrated by Eq. (3).

D. Temperature reweighting
The harvested trajectories also yield equilibrium expectation values beyond the sampling temperature for observables A½X ¼ AðC is Þ that depend on inherent state configurations C is only, provided that (i) configurations are sampled according to the Boltzmann weight ∝ e −UðCÞ=T s (at the sampling temperature T s ) and (ii) vibrations are independent of the inherent state. Previous work [35] has demonstrated via thermodynamic integration that for sufficiently low temperatures (T < 0.8), the contribution of the inherent state energies to the free energy decouples from the vibrational contribution. This satisfies condition (ii), while condition (i) is ensured by our sampling technique.
The key aspects of our technique are then as follows. To show how it is possible to extract the equilibrium statistics of inherent state energies from the trajectory sampling at temperatures distinct to T s , we split the total potential energy into two contributions: where δUðCjC is Þ indicates the extra potential energy of a configuration C compatible with the inherent state configuration C is . Using such a decomposition, we can evaluate the thermal average h·i of the following expression: with partition sum ZðTÞ¼ P C e −UðCÞ=T and ϑ¼1=T −1=T s . We then define the restricted partition sum sampling the fluctuations out of the inherent state C is . If the restricted partition sum is approximately independent of the inherent state,ZðTjC is Þ ≈ZðTÞ, we can rewrite Eq. (5) as Finally, this expression can be used to compute the reweighted average h·i T;μ evaluated at equilibrium μ ¼ 0 in Eq. (3), and we can obtain corresponding to the equilibrium expectation of A at temperatures T ≠ T s , different from the sampling temperature.

Ghost-particle Monte Carlo technique
In order to compare our trajectory sampling results with an alternative equilibrium technique, we perform Monte Carlo simulations where we perturb the Hamiltonian of the system in order to favor local relaxation of caged particle arrangements.
The method is inspired by ghost-particle insertion [37] and swap Monte Carlo algorithms [33,36]. In our implementation, a subset of the particles in the system is associated with a stochastic variable that determines the maximum repulsive forces acting on the particles of the subset. This facilitates local cage breaks and relaxation. The particles that interact via this capped potential are reselected randomly after a certain number of iterations and are called "ghost" particles, hence the name of the method. The random walk of the stochastic variable satisfies detailed balance and includes states for which the Hamiltonian is that of the unperturbed system. We record those states and compute the average inherent state energy hϕi from an equilibrium distribution of configurations [37]. For more details on the method and the algorithm, see Appendix A.

III. RESULTS
Our results are twofold. First, we describe our reweighting procedure to obtain configurations representative of exceedingly low temperatures, far beyond the regime accessible to conventional techniques. This shows that the configurational entropy and inherent state energy are both compatible with a transition to a state of very low entropy at a finite temperature T K . We then proceed to consider the behavior of the system in the space of the nonequilibrium phase transition ðT; μÞ. We show that the nonequilibrium phase transition appears (upon extrapolation) to be bound from below, with a lower critical point T c . This lower critical point lies close to T K (or perhaps even at T K ). We discuss how this proximity of the lower critical point and T K may allow some union of the different approaches to the glass transition.

A. Configurational entropy
We first focus on the behavior of the central configurations. Our simulations sample configurations with a wide range of inherent state energies, which we compile into distributions pðϕÞ for the different sampling temperatures.
In agreement with previous studies [35,48], we find that these distributions are well described by a Gaussian [see Fig. 1(b)]. Moreover, it has been demonstrated that at low enough temperatures (including the sampling temperature 0.48 ≤ T s ≤ 0.7 used here), the vibrational free energy is independent of ϕ to a very good approximation [35].
Let ΩðϕÞ ∝ Ω ∞ e NσðϕÞ δϕ be the number of amorphous inherent states with energy per particle ϕ within an interval ϕ AE δϕ=2. Here, σðϕÞ is the enumeration function and Ω ∞ ≃ e Ns ∞ is the maximal available volume in configuration space, i.e., the high-temperature case where s ∞ is the configurational entropy. The (configurational) temperature T corresponds to the inverse slope, dσ=dϕ ¼ 1=T. In the limit of large N, the number Ω is either extensive or it becomes exponentially small. Hence, in the thermodynamic limit, the extensive configurational entropy becomes ln ΩðϕÞ ≃ NsðϕÞ, with sðϕÞ ¼ s ∞ þ σðϕÞ for ϕ ≥ ϕ K and sðϕÞ ¼ 0 for extreme inherent state energies ϕ < ϕ K . Here, ϕ K ≃ −7.82 is the inherent state energy of the system at the Kauzmann temperature T K .
The Gaussian shape of ΩðϕÞ implies that the enumeration function we extract, σðϕÞ, should be quadratic, σðϕÞ ¼ −ðϕ − ϕ ∞ Þ 2 =J 2 , from the measured distribution with fitted maximal energy ϕ ∞ ≃ −7.38 and energy scale J ≃ 0.502. The resulting configurational entropy sðϕÞ ¼ s ∞ þ σðϕÞ is plotted in Fig. 1(c) using the value s ∞ ≃ 0.74 reported previously [35]. Our numerical scheme is thus able to cover a wide range of inherent states, and the excellent agreement demonstrates that the configurational entropy of the liquid is indeed very well described by a quadratic function.
For the thermal average of the inherent state energy, one finds as a function of temperature T. While we make the distinction between configurational and conventional thermodynamic temperature, for quantities based on inherent states, such as the inherent state energy itself, the configurational temperature corresponds to the conventional temperature. The scenario that hϕi ¼ ϕ K is reached at a finite temperature T K is shown in Fig. 1(c) by extrapolating simulation data. Figuratively speaking, at this temperature, the liquid would "run out" of amorphous configurations and would undergo a thermodynamic phase transition to an ideal glass with constant inherent state energy ϕ K and very low entropy. Although the existence of such an ideal glass state is debated on several grounds, e.g., in Refs. [49,50], we keep ϕ K as a convenient estimate for the lowest energy accessible to amorphous configurations. Extrapolating the quadratic form for sðϕÞ to lower energies, one derives 3 for the present system, in agreement with previous estimates in the range 0.297 < T K < 0.325 [26,35,51]. This value is compatible with an extrapolation of molecular dynamics data for the average ISE hϕi, which we show in the inset of Fig. 1(c). Comparing with the increase of τ α [ Fig. 1(a)], it is clear that T K is far below temperatures for which the liquid can be equilibrated in computer simulations. In Fig. 2, the average inherent state energy hϕi T;μ¼0 is shown as a function of (configurational) temperature T employing our scheme, showing good agreement with data obtained from molecular dynamics simulations and with data obtained from the ghost-particle Monte Carlo technique that promotes relaxation via the usage of staged, capped potentials (see Appendix A). We observe that the results from trajectory sampling for the lowest sampling temperatures are in good agreement with Eq. (10) down to T ¼ 0.37, where the relative deviations start to exceed 0.1%. At T s ¼ 0.48, 0.50, the average inherent state energy departs from Eq. (10) for T < 0.37, and glassy, amorphous, very-low-energy states are obtained. These are inaccessible to conventional techniques since the relaxation time is expected to be, via extrapolation, τ α ðT ¼ 0.37Þ > 10 4 τ α ðT MC ¼ 0.435Þ. Assuming that the Kauzmann energy ϕ K represents the "bottom" of the energy landscape of amorphous states, then we expect that for T < T K , there should be no further reduction in hϕi, which is indeed supported by our data for successively lower sampling temperatures T s . B. Trajectory-space phase transition from LFS-poor to LFS-rich states We now consider the full trajectories rather than the central configuration as above. Formally, one can treat the particles in LFS and non-LFS (free) particles as two chemical species with the same chemical potential, which interconvert freely. The extension to ensembles of trajectories is, at least formally, straightforward, with μ the dynamical analog of the chemical potential difference. Unbiased equilibrium dynamics corresponds to μ ¼ 0. In line with physical intuition, positive μ > 0 gives trajectories with larger populations a higher weight, which is demonstrated in Fig. 3(a). Here, we show the average population hni T s ;μ of LFS as a function of μ for the three sampling temperatures. At some μ Ã > 0, there is a sharp increase of the population, indicating a transition from the normal supercooled liquid (LFS poor) to a state composed of trajectories with a large number of LFS (LFS rich). Indeed, it has been shown that the transition becomes sharper for larger systems [26], where "larger" can be both larger t obs (longer trajectories) and larger N (more particles while holding the density constant). This implies a first-order transition in trajectory space, where the jump of the order parameter (the population of LFS) is rounded by finite-size effects [52]. Practically, the value of μ at which the transition occurs, μ Ã , is determined from the peak of the susceptibility, maximizing the fluctuations of the order parameter n as shown in Fig. 3(a). The susceptibility measures the sensitivity of the population of LFS along the trajectory to a change in the field strength μ. This becomes very large close  to the phase transition and is expected to diverge in the limit of infinite system size. Note that when decreasing the sampling temperature, the peak position moves towards μ ¼ 0, which corresponds to the unbiased, equilibrium system, i.e., the one encountered in experiment. Moreover, upon decreasing temperature, the peak height χ Ã ¼ χðμ Ã Þ grows. We return to the intriguing question of the lowtemperature behavior of μ Ã below.
Having identified a LFS-poor and a LFS-rich phase, in Fig. 3(a) we show that the transition is accompanied by a complementary drop of the inherent state energies of central configurations [53]. Indeed, looking at the joint distribution of n and ϕ in Fig. 3(b), two distinct populations of trajectories can be identified. In particular, trajectories with a large number of LFS feature central configurations that typically have low inherent state energy.
From the marginal distribution of ISE [inset of Fig. 3(c)], we extract the two typical values ϕ poor and ϕ rich for LFSpoor and LFS-rich states, respectively, as a function of configurational temperature T, using Eq. (3). We plot these in Fig. 3(c). The inherent state energies for the LFS-poor liquid follow the prediction of Eq. (10) for the average inherent state energy and are therefore indistinguishable, up to our numerical precision, from the equilibrium liquid. In contrast, the much lower ϕ rich shows a more complex behavior: At high temperatures, the different T s converge at the same energy, which is a slow monotonically increasing function of the configurational temperature T. For the lowest temperatures T, there is a systematic dependence on the sampling temperature T s , with the lower T s allowing us to better sample lower energy states and thus leading to lower average ISE. Because of the much weaker dependence on T of ϕ rich , the energy gap ϕ poor − ϕ rich is reduced with decreasing T, which suggests the existence of a finite temperature at which the two energies reach the same value and the gap vanishes. In other words, we expect a low but finite temperature at which the equilibrium liquid becomes indistinguishable from the low-energy, LFS-rich phase.

C. Nonequilibrium phase diagram
We draw together our results in a nonequilibrium phase diagram where both the equilibrium liquid and the nonequilibrium coexistence are represented [see Figs. 4(a) and 4(b)]. To construct the nonequilibrium phase diagram, we determine the populations n poor and n rich of the coexisting LFS-poor and LFS-rich states, respectively. To this end, we produce histograms of the LFS population at μ ¼ μ Ã and determine the positions of the two peaks corresponding to each phase [see Fig. 4(b), inset]. These populations delimit the coexistence region, as shown in Fig. 4(b), for several sampling temperatures. They agree within errors for different trajectory lengths and system sizes (N ¼ 216 and N ¼ 400 particles). The unbiased equilibrium dynamics with μ ¼ 0 corresponds to a line in the ðT; nÞ plane as shown in red in Fig. 4(b); any point away from this line corresponds to a nonequilibrium state with μ ≠ 0. For the temperatures that we can sample directly, the LFS population in the equilibrium liquid lies close to the spread of the data for the LFS-poor region.

IV. DISCUSSION
A. Lower critical point of the nonequilibrium phase transition?
The existence of a thermodynamic phase transition for glass-forming fluids is widely debated [2,3,12,49]. In the case of the Kob-Andersen mixture, there are several pieces of numerical evidence pointing towards a transition at temperatures substantially below the mode-coupling transition temperature T MC ≃ 0.435: (i) Thermodynamic integration [35] provides, through an extrapolation, an estimate for a finite Kauzmann temperature T K at which the configurational entropy vanishes. (ii) As previously mentioned, relaxation times can be thought to diverge in the vicinity of T K , providing a dynamical description of the transition [see Fig. 1(a)]. (iii) Recent equilibrium simulations [29] in the presence of a variable concentration of pinned particles c show that, in the limit of vanishing c → 0, a transition between a more and a less "glassy" thermodynamic state (as indicated by the so-called overlap order parameter-see Ref. [29] and references therein) can be inferred at a nonzero temperature close to T K . (iv) Finally, the same overlap order parameter [54] has been used to show that static fluctuations exist for T ≲ 0.5, consistent with the existence of a random-field-like critical point predicted by effective field theories [55,56].
However, it is well known that the Vogel-Fulcher-Tamman fit is not the only-nor even the best-fit to the structural relaxation time [39]. Others, which do not imply a thermodynamic transition, give as good of a description [11,57]. Furthermore, the dynamical aspects, in the sense of particle mobility of the nonequilibrium transition that we study, cannot be overstated, whether they are generated dynamically [16] or by structural averaging along trajectories, as in our case. These nonequilibrium phase transitions lead to an inactive state whose time correlations do not decay on the simulation time scale. In other words, the inactive phase is a glass, consistent with the dynamic facilitation scenario, which posits that there is no thermodynamic phase transition.
In the scenario we consider here, the thermodynamic behavior of the system follows from the more general description of its nonequilibrium phase diagram, where the equilibrium limit corresponds to μ → 0. We find compelling numerical evidence for a coexistence region between LFSrich and LFS-poor trajectories, at least in the temperature interval 0.48 ≤ T s ≤ 0.74, directly probed by the biased simulations. The temperatures at which the distinction between LFS-rich and LFS-poor trajectories ceases would identify an upper and a lower critical point to the nonequilibrium phase transition, the locations of which are inferred by extrapolation. In particular, the precise location of the lower critical point T c with respect to the line corresponding to the equilibrium liquid allows for a connection between the putative thermodynamic transition at T K and the dynamical phase transition determined in the present work. In the following, we discuss our evidence for a thermodynamic transition and also consider alternative scenarios.

B. Numerical evidence for a thermodynamic transition and critical point
To interpret our numerical results, we first consider the extrapolation of the dynamical chemical potential μ Ã at coexistence. In Fig. 4(a), we illustrate the phase diagram in the T − μ plane, with the coexistence region being represented by a line separating the LFS-poor from the LFS-rich trajectories. Here, we have plotted the different values of chemical potential μ Ã at coexistence for different sampling temperatures T s , scaled with the trajectory length t obs , and we observe that they collapse onto a common master curve.
Down to the coldest sampling temperature T s ¼ 0.48, μ Ã ðTÞ follows a linear behavior, whose extrapolation to μ Ã ¼ 0 provides a second estimate for a crossing temperature T × ðμ Ã → 0Þ ≃ 0.37, 23% higher than the estimate for the Kauzmann temperature T K ≃ 0.3. We observe [inset of Fig. 4(a)] that the position of the maximum of the susceptibility χðμÞ moves from very large values of μ Ã at high temperatures (i.e., very rare fluctuations) towards μ ¼ 0 (i.e., typical equilibrium fluctuations) for lower temperatures. This is accompanied by a nonmonotonic behavior of the peak amplitude, suggesting critical-like fluctuations at some finite high temperature (which, for physical reasons, we relate to the onset of glassy dynamics) and a low temperature T × . If we only fit the coexistence points in the T − μ plane inside the regime for which inherent states are well defined and the susceptibility peak χ is monotonically increasing as we decrease the temperature [T ≤ 0.60, see inset in Fig. 4(a)], the linear fit provides a lower estimate for T × ¼ 0.33 closer to the estimate of the Kauzmann temperature.
Our second piece of evidence for a thermodynamic transition comes from the population of LFS at coexistence between the LFS-rich and LFS-poor phases, n Ã ðTÞ ¼ hni T;μ Ã . To do this, we cast the phase diagram in the T − n plane. As shown in Fig. 4(b) (green dashed line), the dependence of n Ã upon temperature is linear to a very good degree. Extrapolating towards lower temperatures, it crosses the equilibrium (μ ¼ 0) line at a "crossing" temperature T × with n Ã ðT × Þ ¼ hni T × ;μ¼0 , which implies n poor ¼ n rich . The crossing of both lines occurs at T × ðn poor ¼ n rich Þ ≃ 0.31, close to the estimated value for the Kauzmann temperature T K ≃ 0.30.
Our third piece of evidence stems from the inherent state energies. We note that the inherent state energy of the LFSrich configurations in Fig. 3(c) is approximately linear, with fitted ϕ 0 ≃ −7.82 (note that ϕ 0 ≈ ϕ K ) and γ ≃ 0.18 for the lowest sampling temperature. Assuming that the vibrational free energy and that of the inherent states still decouple, for any transition, the configurational temperatures (i.e., the slopes of the configurational entropies of the LFS-rich and LFS-poor phases) have to agree. For a putative continuous transition, we equate the inherent state energy [using the equilibrium liquid, Eq. (10) for the LFS-poor and Eq. (12) for the LFS-rich phase] and find ϕðT × Þ ≃ −7.76, with corresponding temperature T × ðϕ rich ¼ ϕ poor Þ ≃ 0.33, which also provides an estimate for the temperature at which the two phases are indistinguishable, the nonequilibrium critical temperature T c . Finally, the information contained in the narrowing of the gap between the inherent state energies of the LFS-rich and LFS-poor phases can be translated into the T − n phase diagram through additional numerical nonequilibrium simulations. The estimate of the inherent state energy of the LFS-rich phase at low temperatures is used (see Appendix B) to select a representative population of configurations. We then use these as initial configurations for molecular dynamics simulations at temperature T, which explore a metastable basin and eventually relax to the equilibrium liquid state. We assume that the LFS-poor liquid nucleates within the LFS-rich liquid; therefore, we determine the characteristic time for nucleation τ melting and the LFS populations before and after nucleation. These populations are then employed as estimates at temperatures for the coexisting LFS-rich and LFS-poor densities, respectively [see Fig. 4(b)]. We notice that our estimate for the gap n rich − n poor narrows from n rich − n poor ¼ 0.15ð1Þ at T ¼ 0.70 to n rich − n poor ¼ 0.09ð1Þ at T ¼ 0.40.
Although relying on extrapolations, these four pieces of numerical evidence indicate that the dynamical transition from LFS poor to LFS rich would become accessible for the equilibrium supercooled liquid, provided that it is cooled sufficiently slowly and that crystallization does not interfere. From the numerical data, we can bound the crossing temperature T × at which the equilibrium (LFS-poor) liquid becomes undistinguishable from the low-energy, LFS-rich amorphous phase, obtaining 0.31 < T × < 0.33. This range is intriguingly close to previous estimates of the Kauzmann temperature T K ≃ 0.30 at which a thermodynamic phase transition has been predicted [29,35]. Moreover, from the narrowing of the gap between the LFS-poor and LFS-rich coexisting phases, we infer that a lower critical point exists at a temperature T c , terminating the nonequilibrium coexistence region between LFS-poor and LFS-rich trajectories.

C. Low-temperature fate of the supercooled liquid
Our numerical precision does not enable us to distinguish the equilibrium liquid and the LFS-poor phase, that is to say, T × ≈ T c . However, considering the precise location of the nonequilibrium critical point T c with respect to T × , we obtain three possible scenarios: First, the critical temperature T c > T × is higher, in which case no equilibrium transition occurs but is closely passed (an avoided transition) [Figs. 4(d) and 4(g)]. At T × , the equilibrium system then crosses the line emanating from the nonequilibrium critical point where the fluctuations are maximal (sometimes called the "Widom line" [58,59]) [see Fig. 4(d)], analogously to some scenarios in softened kinetically constrained models [60]. Second, in the case where T c ¼ T × , the equilibrium unbiased dynamics crosses through the critical point, and a continuous transition occurs, as depicted in Fig. 4(b). The final possibility is that the nonequilibrium critical point lies below the crossing temperature, T c < T × . In this case, the equilibrium transition, corresponding to the crossing of the coexistence region with a small but finite increase of the LFS concentration, would be (weakly) first order, as sketched in Figs. 4(e) and 4(h).
Before concluding, we emphasize that our discussion concerns the linear extrapolation of the dynamical chemical potential μ Ã and LFS population n Ã at coexistence, a supposition that the inherent state energy continues as indicated in Fig. 2 and that the relaxation time follows the VFT prediction. Other alternatives are possible, such as a more gradual approach to a state of very low configurational entropy, which could be met only at T ¼ 0. This could be related to "defects" in the amorphous order [49]. Other scenarios have been suggested in which μ Ã ðTÞ approaches equilibrium [μ Ã ðTÞ ¼ 0] more slowly at lower temperatures, only reaching zero at T ¼ 0 [61], as depicted in Figs. 4(c) and 4(f). These would imply deviations in the relaxation time from the VFT behavior. It is even possible that these effects lead to behavior similar to the one we discuss but at a lower temperature than our extrapolations suggest.

V. CONCLUSIONS
By means of numerical simulations employing trajectory sampling, we have explored the connection between a dynamical phase transition and the low-temperature thermodynamics of a model atomistic glassformer. We have used locally favored structures as an order parameter to bias the simulations in trajectory space at a given sampling temperature. Thus, we have shown that this is a powerful and reliable method to obtain equilibrium and thermodynamic information (such as the inherent state energies and the configurational entropy) down to exceptionally low temperatures. This is because trajectory sampling facilitates the exploration of states that are exponentially unlikely to be accessed in a conventional molecular dynamics simulation.
Employing the probability distributions obtained, we show that the dynamical phase transition in trajectory space implies a bimodal distribution of inherent state energies: For the equilibrium dynamics, most of the configurations are compatible with a high-energy, locally favored structurepoor state: the normal liquid. However, we also reveal that an amorphous, LFS-rich, low-energy state exists and is accessible when the dynamical chemical potential μ is tuned away from zero [26]. We also determine the temperature dependence of the energy gap between the two states, which provides numerical evidence for a lower critical point at a nonzero temperature T c , where the high-energy (LFS-poor) phase would be indistinguishable from the low-energy (LFS-rich) phase. This is further corroborated by the temperature at which the first-order dynamical phase transition appears to cross the unbiased dynamics (μ ¼ 0). Numerically, we find that the temperature T c is close to the Kauzmann temperature estimated from our configurational entropy measurements. We confirm our findings with particle swap methods that we have tailored for this system.
Our results support the idea that the low-temperature fate of supercooled atomistic glassformers is determined by the fluctuations of the population of local structures, which couples with the reduction of configurational entropy.
There is now ample numerical and experimental evidence that particle dynamics, at least for a class of model glassformers [20,23,27,62], is correlated with the presence of LFS. We observe a nonequilibrium phase transition between two phases in trajectory space, LFS-poor and LFSrich phases. The former is, given our numerical precision, indistinguishable from the equilibrium liquid. The latter exhibits very slow dynamics [26]. Glassy phenomenology, such as dynamical heterogeneities and the emergence of slow dynamics, may be related to this nonequilibrium phase transition in trajectory space. Our numerical data provide evidence that the transition may terminate at a lower critical point close to the Kauzmann temperature T K .
The connections we identify between dynamical coexistence, the distinct basins in the energy landscape, and the location of a low-temperature nonequilibrium critical point close to the equilibrium (μ → 0) path in the nonequilibrium phase diagram provide the ingredients for a novel understanding of the low-temperature fate of supercooled liquids. In particular, our work opens a potential avenue to unify the competing theories that have been developed in order to understand the microscopic origin of slow glassy dynamics: On one hand, the thermodynamic picture of a decreasing configurational entropy, vanishing at a finite temperature around T K , is compatible with the emergence of the low-energy, low-entropy, LFS-rich state that we have identified. On the other hand, our results suggest that down to T K , this state may be in coexistence in trajectory space with the normal equilibrium liquid and that the glassy phenomenology may be related to the vicinity of the equilibrium line to this nonequilibrium phase coexistence. Within this approach, structural and dynamical aspects are combined, and they play a complementary role. We emphasize that our conclusions rely on extrapolation. Other possibilities, including a continuous decrease of the configurational entropy such that it approaches zero only as temperature also approaches zero, are possible. That is to say, our results do not exclude the possibility that there is no vanishing of configurational entropy or indeed any transition at T K [61]. Further alternatives include a suppression of any transition around T K by "defects" [49]. Finally, despite the Kob-Andersen mixture being a good glassformer, at low temperatures T < 0.40, it is prone to crystallization, with a microphase separation that promotes the formation of a fcc crystal of large particles coexisting with the fluid [63]. In our calculations, this is anticorrelated with the growth of the LFS population, so it identifies a different regime; however, the mechanisms that lead to crystallization may play a role in the metastability of the LFS-rich phase.
Applying the method we have developed to other model glassformers with different fragilities and different structural signatures [27], and comparing it with existing theories for critical behavior in glasses [55,64] will be the challenging subjects of future investigations. Particleresolved experiments with colloids and granular materials [65] provide the natural means by which our predictions may be verified experimentally, as these experiments provide the same data as the simulations we have used here. Progress has been made in this direction, with identification of the dynamical phase transition in trajectory space to a LFS-rich phase [15]. It is thus possible to directly investigate the phenomena we have identified here.
More generally, provided a LFS can be identified with a given system, we expect it to be possible to observe the kind of behavior we see here. In particular, we expect that other Lennard-Jones and hard sphere systems should exhibit similar behavior, as LFS are known for these systems [66]. In metallic glasses, LFS have also been found [67][68][69]. Other glassformers in which local structure has been investigated include oxides and other inorganic materials [70]. A challenging direction of research would be to apply trajectory sampling to simple models of silica (such as the so-called BKS model [71] or the NTW model [72]), where a crossover from fragile to strong behavior is observed [73], in order to identify its relationship with the dynamical phase transition discussed in the present work. Beyond that, it is important to observe that organic molecules, and even polymers at the present time, form an open challenge as LFS have yet to be identified. For some coarse-grained models (such as the Lewis-Wahnström model for ortho-terphenyl [74]), it would be possible to identify local structure and its relation to the dynamics following an analysis similar to the present work. However, the greatly increased number of degrees of freedom of molecular systems means that convincingly determining a LFS still constitutes a considerable task, and other forms of amorphous order may prove more amenable to identify the kind of phenomena we present here [24,[75][76][77]. We stress, though, that we expect the behavior we see to be generic to a wide range of materials.

ACKNOWLEDGMENTS
The authors are grateful to Ludovic Berthier, Juan P. Garrahan, Robert L. Jack, Walter Kob, and John Russo for helpful comments. Stephen Williams is gratefully acknowledged for helpful discussions and preliminary results. C. P. R. acknowledges the Royal Society for funding and Kyoto University SPIRITS fund. F. T. and C. P. R. acknowledge the European Research Council (ERC Consolidator Grant NANOPRS, Project No. 617266). This work was carried out using the computational facilities of the Advanced Computing Research Centre, University of Bristol.

APPENDIX A: GHOST-PARTICLE MONTE CARLO TECHNIQUE
We compare the inherent state energies obtained from trajectory sampling to equilibrium values obtained via a modified Monte Carlo technique inspired by the ghostparticle insertion method [37] and the employment of swap Monte Carlo techniques in order to equilibrate supercooled liquids closer to the glass transition.
The Kob-Andersen binary mixture is characterized by high density and relatively hard interactions which make traditional particle-swap techniques highly inefficient. Particle swaps appear to be much more effective for ternary mixtures [32], soft interactions [36], or very polydisperse hard spheres [33].
Here, we choose to improve the sampling efficiency in the Kob-Andersen mixture by introducing a third species in our system. At any time step, the system is split into N − 1 particles interacting via the standard nonadditive Kob-Andersen interactions (see Sec. II A) and a so-called "ghost" particle g. This can be both a large particle A or a small particle B, and it interacts with any other particle i via a capped potential: where γ g and γ i correspond to the A or B types of particles g and i, while s is a discrete variable associated with the state of particle g. We define a discrete space of N stages states corresponding to different values of the maximum allowed energy V s max . During the simulation, we perform a Markov chain in the space of states, with transitions between state s and s 0 ¼ s AE 1 governed by a Metropolis acceptance probability: We then implement the following algorithm: (1) We perform N attempts of local Monte Carlo moves following the usual Metropolis prescription. (2) We perform a Monte Carlo step according to Eq. (A2) and modify the discrete state variable s and the corresponding potential cap V s max associated with the ghost particle g.
(3) When the discrete variable s assumes the value N stages , the ghost particle becomes a normal particle. We therefore store the corresponding configuration, and the ghost status s ¼ 0 is assigned to a new, randomly selected particle. (4) We repeat steps 1-3 in order to obtain equilibrated configurations, whose inherent state energies are evaluated using the FIRE algorithm (see Sec. II A above), while trajectory sampled configurations appear to sample energies that are significantly lower. Our choice for the discrete set of states is the following: The idea behind the algorithm is to facilitate the particle motion by lowering the energy barriers for cage breaks locally while preserving detailed balance. As shown in Fig. 2, the technique produces reliable estimates of the inherent state energies down to T ¼ 0.40. At T ¼ 0.375, the system still appears to be trapped in a glassy state, within the considered computation time (256 CPU hours, 4 × 10 7 Monte Carlo sweeps), and trajectory sampling provides significantly lower inherent state energies.

APPENDIX B: NONEQUILIBRIUM RELAXATION FROM LFS-RICH CONFIGURATIONS
In order to provide an estimate for the coexisting LFSpoor and LFS-rich populations at temperatures lower than the lowest sampling temperature T s ¼ 0.48 employed in the trajectory sampling technique, we perform conventional molecular dynamics simulations starting from a selected sample of initial configurations.
We first determine the expected average inherent state energy for the LFS-rich phase ϕ rich ðTÞ at temperature T via reweighting: See Sec. II D and Fig. 2. From our population of sampled trajectories, we then extract a number L ¼ 70 of central configurations of N ¼ 216 particles whose energy is within the interval ½ϕ rich − σ; ϕ rich þ σ, where σ ¼ 0.006 is chosen as the typical width at half-height of the LFSrich peak in the reweighted probability distribution pðϕ; μ ¼ μÃÞ [see inset of Fig. 3(c)]. We measure the relaxation from the LFS-rich phase to the LFS phase with a constant temperature protocol, where the system is coupled to a Nose-Hoover thermostat at temperature T whose characteristic damping time is set to a tenth of the structural relaxation time.
Following this protocol, we produce trajectories of duration 4τ α temperatures T ¼ 0.44, 0.42, 0.40 and measure the population of LFS n along the trajectories. For every trajectory, we model the time evolution of the LFS population as a two-state function nðtÞ¼AþBtanhðt−τ melting Þ= C and fit the model to the data to obtain per-trajectory estimates of the time τ melting needed to melt the LFS-rich state. For every trajectory, we then estimate the LFS-rich coexisting population with the time average of the LFS population in the interval ½0; τ melting , and also the LFS-poor population from the remaining part of the nonequilibrium trajectory.
Averaging over all the initial conditions, we obtain mean and standard errors as represented in Fig. 4(b) by pink (NVT protocol) symbols. The estimates for the LFS-poor phase (despite not resulting from equilibrium runs) are, at most, 3% higher than the prediction for the equilibrium liquid, Eq. (2).