Finite-disorder critical point in the yielding transition of elasto-plastic models Saverio Rossi, 1 Giulio Biroli, 2 Misaki Ozawa, 2 Gilles Tarjus, 1 and Francesco Zamponi 2 1 LPTMC, CNRS-UMR 7600, Sorbonne Université, 4 Pl. Jussieu, F-75005 Paris, France 2 Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, F-75005 Paris, France (Dated: November 28, 2022)

Upon loading, amorphous solids can exhibit brittle yielding, with the abrupt formation of macroscopic shear bands leading to fracture, or ductile yielding, with a multitude of plastic events leading to homogeneous flow. It has been recently proposed, and subsequently questioned, that the two regimes are separated by a sharp critical point, as a function of some control parameter characterizing the intrinsic disorder strength and the degree of stability of the solid. In order to resolve this issue, we have performed extensive numerical simulations of athermally driven elasto-plastic models with long-range and anisotropic realistic interaction kernels in two and three dimensions. Our results provide clear evidence for a finite-disorder critical point separating brittle and ductile yielding, and we provide an estimate of the critical exponents in 2D and 3D.

Yielding of amorphous materials is a practically and scientifically important problem [1][START_REF] Barrat | Dynamical Heterogeneities in Glasses, Colloids, and Granular Media[END_REF][START_REF] Rodney | [END_REF][4][5]. When a material is mechanically slowly driven from an initial quiescent glassy state, two different types of yielding behavior are observed. One is brittle yielding, where the sample catastrophically breaks into pieces and displays one or several macroscopic shear bands (usually experimentally encountered in atomic and molecular glasses). The other one is ductile yielding, for which the sample deforms rather homogeneously via a series of local and mesoscopic plastic events that prevent catastrophic failure (usually experimentally encountered in soft materials like colloids and pastes). It has been established that a given material may show brittle or ductile yielding depending on the preparation history of the sample [START_REF] Rodney | [END_REF][6][7][8]. In particular, a well-annealed, hence stable, glass sample shows brittle yielding, whereas a poorly-annealed, less stable, glass sample exhibits ductile yielding. Note that in the materials science and engineering communities the yield point is traditionally defined as the end of the purely elastic branch and the onset of plastic behavior. However, because several molecular simulations and elastoplastic model (EPM) studies (see, e.g., [9][10][11]) demonstrated that a purely elastic branch does not exist in sheared amorphous solids as plasticity appears for any infinitesimal deformation in the thermodynamic limit, the statistical physics community usually adopts a different definition of yielding.

Recent theoretical studies suggest that brittle yielding corresponds to a nonequilibrium first-order transition (or spinodal [12][13][14][15][16][17]), associated with a macroscopic discontinuous stress drop at a given strain value, whereas ductile yielding corresponds to a continuous stress-strain curve, corresponding to a progressive plastic softening of the material. In athermal quasi-static (AQS) conditions [18], it was observed that these two distinct behaviors are separated by a critical value of the stability (or the disorder) [19][20][21]. It was then proposed that the brittle-to-ductile transition is a novel nonequilibrium phase transition, similar to that of an athermally driven random-field Ising model [22,23]. Further understanding the transformation from brittle to ductile yielding appears as a major challenge in many fields, from materials science to statistical physics [24][25][26][27][28][29][30][31][32][33][34][35].

The above scenario has been challenged in Refs. [36][37][38]. In particular, Ref. [36] argues that in AQS condition and provided the samples are large enough, yielding always takes place in a brittle manner. This should happen irrespectively of the stability or disorder of the samples, except for the putative infinitely disordered sample, thereby implying that the brittle-to-ductile transition does not exist in the thermodynamic limit. Largescale molecular simulations [38] seem to give some support to the statements of Refs. [36,37]. Yet, it remains hard to conclude due to the limited system sizes accessible in molecular simulations and the very small number of samples involved.

In this Letter, in order to overcome this difficulty and obtain conclusive results, we perform a thorough numerical analysis of the brittle-to-ductile transition in EPMs [1,39]. These mesoscopic models have already been successfully applied to describe several aspects of the rheology of amorphous materials, in particular the yielding transition [37,40]. Their coarse-grained lattice nature enables us to access very large system sizes and a large number of samples, allowing for a careful finitesize scaling analysis of the critical point. Our main result is a direct numerical evidence for the existence of a finite-disorder critical point separating brittle and ductile behavior.

The type of EPM we focus on provides a simple scalar description of the AQS dynamics [41][42][43] and corresponds to a cellular automaton on two-dimensional (2D) and three-dimensional (3D) cubic lattices. In particular, we consider incompressible, homogenous and isotropic materials under a simple shear deformation protocol and we focus on a single shear-stress component that we denote σ [39]. (This is an approximate treatment which ignores the other stress components.) The model describes the evolution of coarse-grained local stresses σ i in the presence of an external strain γ. Whenever one such stress goes above a stability threshold, the site yields, and the resulting stress drop is propagated through the sample via a long-range Eshelby-like propagator [44]. The initial stability of the solid, which in real systems depends on the annealing protocol, is quantified by a parameter R associated with the width of the initial stress distribution (and therefore characterizing the strength of the disorder). We vary the system size over a wide range of linear box lengths, L = 256 -4096 for 2D and L = 48 -164 (with a few samples at L = 200) for 3D. Details concerning the simulated model and the numerical simulations are presented in the SM [45]. We have checked that variations of the model corresponding to different ways of accounting for the initial stability and the force balance lead to the same results (see the SM).

We first show that the model displays the same behavior as that found in numerical simulations of particle systems [START_REF] Rodney | [END_REF]19,20], with in particular the signature of a brittle-to-ductile transition accompanied by substantial finite-size effects [37,38]. In consequence, it provides a suitable framework to address the issues discussed above. In Fig. 1, we present stress-versus-strain curves for a 3D system with L = 200 for two values of the disorder strength, R = 0.3 and R = 0.8. The former clearly shows brittle yielding characterized by a discontinuous stress drop and the appearance of a shear band (top inset), while the latter displays ductile yielding characterized by a continuous monotonic stress growth and homogeneously distributed plastic events (bottom inset). Brittle and ductile yieldings are thus qualitatively distinct, and their occurrence depends on the disorder strength R. The present results are in line with previous numerical observations in two-dimensional EPMs [40,42,46].

Our second goal is to identify the putative critical point separating brittle and ductile yielding. In previous studies [19,20,28], the maximum stress drop, ∆σ max = max γ {∆σ(γ)} where ∆σ(γ) is the stress drop due to irreversible events at the strain γ, was used as an order parameter to detect the transition. Here, • • • denotes an average over many independent realizations (or samples), and the maximum is computed for each sample and then averaged over samples. We have found that for EPMs a more efficient order parameter is obtained from the fraction of sites along a single line in 2D or plane in 3D that have yielded at least once up to the strain γ. We consider the maximum of this fraction over all horizontal and vertical lines (2D) or planes (3D) and we call this quantity n(γ). By construction it is an increasing function of γ. It shows a discontinuous jump of order O(1) when σ(γ) shows a discontinuous drop of order O(1) and it increases continuously when σ(γ) shows a continuous ductile behavior (see below). Therefore, n(γ) essentially contains the same information as σ(γ) for distinguish- ing brittle and ductile yielding behavior. We then define a new order parameter, ∆n max = max γ {∆n(γ)}, where ∆n(γ) is the jump of n(γ) that takes place in the AQS dynamics from γ to γ + ∆γ in a given sample. We have observed that ∆n max better quantifies the abrupt emergence of a system-spanning shear band and, as a result, detects the critical point in EPMs more accurately than ∆σ max (see the detailed discussion in the SM).

As seen in Figs. 2(a,b), the average order parameter ∆n max is small and essentially constant at high disorder strength R and it starts to rapidly grow below some finite value of R. Moreover, as L increases both in 2D and in 3D, the increase of ∆n max with decreasing R becomes steeper while the flat part becomes smaller, suggesting the presence of a critical point. Fig- ures 2(c,d) show the variance of ∆n max , which corresponds to the associated "disconnected susceptibility", χ dis = N Var(∆n max ), defined in analogy with an AQS driven random-field Ising model [19,21]. (It provides crisper, but similar, results than the "connected susceptibility" χ con = -∂ ∆n max /∂γ [19,20].) The disconnected susceptibility is strongly peaked and the peak becomes sharper and higher with increasing L, suggesting a divergence at some critical point. Essentially the same trend is observed for ∆σ max (see SM), in agreement with the results of molecular simulations [19].

To firmly establish the existence of the critical point, we have performed a detailed finite-size scaling analysis. We use here the same scaling ansatz as that of the AQS driven random-field Ising model, χ dis (r, L) ∼ L γ/ν Ψ(rL 1/ν ), where R c (L) locates the maximum value of χ dis , r = (R -R c (L))/R is the reduced disorder strength, Ψ(•) is a scaling function, and with γ and ν some critical exponents. According to this ansatz, the maximum over r of χ dis (r, L) should diverge as L γ/ν and its full width at half maximum should vanish as L -1/ν . The corresponding plots obtained from the data in Figs. 2(c,d) are shown in Fig. 3. We observe a good power-law behavior, and by fitting these curves we obtain γ/ν = 1.86 ± 0.02 and ν = 3.0 ± 0.3 in 2D, and γ/ν = 2.66±0.04 and ν = 2.5±0.2 in 3D, where the errors are derived from the fit. We also show in Fig. 4 the scaling collapse of the disconnected susceptibility, in which the parameters γ, ν, and R c (L) are adjusted to provide the best visual collapse of the curves for the different values of L. The displayed collapses are for γ/ν ≈ 1.82 and ν ≈ 2.9 in 2D, and γ/ν ≈ 2.61 and ν ≈ 2.2 in 3D, values that are consistent with those determined by the fitting procedure. Work is now in progress to determine whether these critical exponents are in the same universality class as an AQS driven random-field Ising model with Eshelby-like interactions [47]. Figures 2, 3 and 4 provide very strong evidence for a critical behavior around R c (L) with an estimate for the associated exponents γ and ν in 2D and 3D. However, the critical disorder R c (L) slightly shifts to larger R as L increases, as seen from Fig. 2. Understanding the fate of the critical disorder R c (L) in the thermodynamic limit is therefore a key issue. Ref.s [19,20] proposed that R c (L → ∞) stays finite in the thermodynamic limit, whereas Refs. [36][37][38]

argued that R c (L → ∞) → ∞.
Note that in this second scenario there is no ductile phase for large enough system size, i.e., all systems are brittle in the thermodynamic limit (except in the singular infinite-disorder limit). We stress that the existence of a finite-disorder brittle-to-ductile critical point in the thermodynamic limit is a separate issue from the persistence of an overshoot in the average stress-versus-strain curve for large ductile systems, which was the main concern of Ref. [37]. We show below that by disentangling these two problems one can obtain conclusive evidence in favor of the existence of the critical point in the thermodynamic limit.

We display in Fig. 5(a) the stress-versus-strain curves of typical 3D samples at fixed R for several values of L. We set R = 0.40 (> R c (L)), which belongs to the putative ductile yielding regime as determined from the above finite-size scaling analysis. The plots focus on the stress values around the overshoot. For a fixed R, the stress drop tends to become sharper with increasing L, showing the same trend as found in Refs. [36][37][38]. Instead, as shown in Fig. 5(b), for a fixed L (here, L = 128) a clear evolution between distinct yielding patterns is observed as R is decreased, from a purely monotonic increase of the stress to a continuous overshoot and then to a discontinuous drop. To characterize the asymptotic behavior when L → ∞, we locate for each given system size L the value of R at which the overshoot first appears (coming from large R) in the average stress-versus-strain curves and we denote it by R o (L) (see the SM for details). We display R o , together with the critical disorder R c , for 2D and 3D in Fig. 6. To facilitate the comparison between 2D and 3D we plot R o and R c as a function of the number of sites N = L D . We find that R o is essentially independent of N in both cases while R c increases very slowly with N . As we explain in more detail below, this is direct evidence for the existence of a ductile phase over a finite range of disorder strength in the thermodynamic limit.

The values of R o (N ) and R c (N ) define three distinct yielding regimes in the (N, R) plane, as schematically illustrated by the insets in Fig. 6. The region R > R o (N ) corresponds to a monotonic increase of the average stress, with no overshoot. The region R < R c (N ) corresponds to a discontinuous stress drop at yielding. The regime R c (N ) < R < R o (N ) corresponds to a continuous average stress curve with a mild overshoot. By construction, R c (N ) has to remain below R o (N ), which then gives an upper bound on the critical disorder. The fact that R o (N ) is essentially independent of N thus provides strong evidence that R c (N ) converges to a finite value for large N and that a finite-disorder brittle-to-ductile critical point persists in the thermodynamic limit. The fate of the overshoot as N → ∞ is instead unclear and depends on whether R c (N ) converges to R o (∞) or to R c (∞) < R o (∞) in the thermodynamic limit. In the former case the overshoot disappears at the critical point whereas a regime of ductile yielding with an overshoot exists in the latter case. We show in Fig. 6 the best fits to R c (N ) = R c (∞) -a/N b with R c (∞), a, and b free parameters. We find that R c (∞) is finite in 2D and 3D. In the critical scaling picture and assuming that 2D and 3D are below the upper critical dimension, the parameter b is related to the (correlation length) exponent ν through 1/ν = Db. The fits then yield ν ≈ 3.57 in 2D and ν ≈ 2.22 in 3D, values which, given the large uncertainties, are consistent with the previous determinations given above.

Strictly speaking, we cannot exclude an alternative scenario in which R o (N ) would start to increase with N above some size N * which is out of reach of present-day simulations and would ultimately diverge in the thermodynamic limit together with R c (N ). However, in view of the absence of any observable N -dependence of R o (N ) in the accessible range, which spans three decades in 2D, and of the lack of any sound theoretical argument supporting the existence of a critical size N * , this possibility seems extremely unlikely.

In conclusion, we have performed extensive numerical simulations of athermally driven elasto-plastic models in two and three dimensions. Thanks to the simple coarsegrained, lattice-based, nature of the modeling, we have been able to simulate substantially larger system sizes and larger number of samples than in molecular simulations, allowing us to perform a thorough finite-size scaling analysis. We have obtained clear evidence for the existence of a critical point separating brittle from ductile yielding in 2D and 3D and we have provided estimates for two associated critical exponents. Our results establish, at least for the studied elasto-plastic models, that criticality persists in the thermodynamic limit and takes place for a finite value of the disorder characterizing the samples (and corresponding to a given initial stability of the solid), as suggested in Refs. [19,20]. The alternative scenario [36,38] according to which the critical point either takes place at infinite disorder or disappears because the disorder cannot go beyond some upper bound is not plausible in view of our results from elastoplastic modeling. It is nonetheless still unclear if the overshoot in the average stress-versus-strain curve disappears in the ther-modynamic limit for the ductile regime, as advocated in Ref. [37]. If it does, the critical point would occur exactly when the overshoot associated with brittle yielding disappears and then gives way to a monotonic (albeit singular at criticality) average stress-strain curve. The other possibility is that the critical point takes place at a value of the disorder for which a smooth overshoot is still present. The latter case implies that disorder, which is not accounted for in the linear instability argument of Ref. [37], is able to pin the propagation of the instability, thereby allowing for the presence of a smooth overshoot. It is hard to go beyond the present study in terms of numerical simulations. Thus, progress is now needed on the theoretical front.

Finally, we point out that the presence of a finitedisorder critical point is not restricted to the specific rheological setting considered in this paper. Recently, the AQS cyclic shear protocol has been actively studied, in relation to other nonequilibrium phase transition phenomena such as absorbing-state phase transitions. This protocol also leads to a transition from ductile to brittlelike behavior, as a function of the disorder or stability of the initial glass samples, as shown in molecular simulations [27,28], simulated EPMs [46,48], as well as a mean-field EPM [49]. A detailed characterization of the critical point under cyclic shear and the determination of the associated exponents would be an interesting subject for future research.

SUPPLEMENTAL MATERIAL Implementation of Elasto-Plastic Models

We study two and three-dimensional lattice-based elasto-plastic models [39] with periodic boundary conditions. We consider a simple scalar description whose main quantity is the local shear stress σ i on each site i of the lattice. The macroscopic stress, σ (without index), is defined by σ = 1 N i σ i . We then assign a condition for the local stability: when a site i has |σ i | > σ th i it is considered to be unstable and it yields locally with a stress drop η i > 0; otherwise this site is stable. σ th i is the threshold stress. In this paper, we set it to one, σ th i = 1, uniformly for all sites for the model presented in the main Another choice with nonuniform σ th i 's [50] will be discussed below. The stress drop η is drawn from an exponential distribution, P (η) = 1 η e -η/η , with mean value η = 1. This exponential distribution of the local stress was observed in molecular simulations [51].

The initial condition at γ = 0, σ i (γ = 0), is independently and identically extracted from a probability distribution P 0 (σ i ), given by

P 0 (σ i ) = (1 -σ 2 i ) N e -σ 2 i /(2R 2 ) , σ i ∈ [-1, 1], (1)
where R is a parameter that characterizes the strength of the initial disorder of the sample and N is a normalization constant. Notice that when R is small this distribution has essentially the same form as a Gaussian distribution of zero mean and standard deviation R.

To ensure that each site respects the stability condition, |σ i (γ = 0)| < 1 at γ = 0, we add a correction term, 1 -σ 2 i , in front of the Gaussian distribution. When the external loading is applied, the local stress at each site is increased by the same amount ∆σ ext , which is the minimum stress increment needed to induce a local instability at a single site somewhere in the sample, triggering a plastic event [42]. When this unstable site, say i, yields, it influences all the other sites via a discrete stress propagator from site i to j, namely,

G j,i = G(r j -r i ) , (2)
leading to a distribution of new stresses on sites j = i

σ j → σ j + G j,i η i , (3)
where η i is the stress drop at site i. Note that choosing for convenience G i,i = -1 allows one to use Eq. (3) also for j = i. This stress propagation may lead to other sites reaching the condition of instability. This process goes on until all sites become stable, forming an avalanche. We next proceed to a new increase of the external stress ∆σ ext so that a single site yields, and the process is repeated. This driving mechanism decouples the timescale of the external loading and the timescale of formation of an avalanche, which amounts to the so-called athermal quasi-static driving (AQS) [18]. More precisely, the evolution of the model is given by the following algorithm:

1. Initialize the stresses σ i (γ = 0) that are i.i.d. from P 0 (σ i );

2. Find site i that is closest to its positive threshold, i.e., minimizes σ th i -σ i ;

3. Increase the stress on all sites by ∆σ ext = σ th i -σ i such that site i yields; 4. Extract η i from P (η) and evolve all sites according to σ j → σ j + G j,i η i (recall that G i,i = -1);

5. For all unstable sites with |σ i | > σ th i , choose a random stress drop δσ i = sign(σ i )(|σ i | -σ th i + η i ) with η i i.i.d. from P (η), assign δσ i = 0 to stable sites, and transform all sites in parallel according to σ j → σ j + i G j,i δσ i ; repeat until all sites are stable.

Repeat from (2).

In each step, the local stresses evolve as

σ j → σ j + ∆σ ext + i G j,i ∆σ i , (4)
where ∆σ i is the total accumulated stress change during steps 4 and 5 of the procedure. Because the macroscopic stress is the spatial average of σ i , summing Eq. ( 4) over j and noting that Ĝq=0 = j G j,i does not depend on i, it evolves as

σ → σ + ∆σ ext + Ĝq=0 1 N i ∆σ i , (5)
where q = (q x , q y ) in 2D and q = (q x , q y , q z ) in 3D denote the wave-vector. The macroscopic strain at each step of the evolution is obtained as follows. In the beginning, γ = 0. At each step, the strain γ is then increased by an amount ∆γ = ∆σ ext /µ, where µ is the local shear modulus that we set to one (µ = 1). The choice of Ĝq=0 is then related to the way the system is driven. We consider in this work the strain-controlled driving, which corresponds to Ĝq=0 = -1. With this choice, the stress evolution is

σ → σ + µ∆γ -∆σ p , (6)
where ∆σ p = 1 N i ∆σ i corresponds to the irreversible plastic stress drop at constant strain. Other driving settings including the stress-control one were for instance studied in Ref. [52].

In this work we choose the Eshelby elastic stress propagator, which in the case of a continuous material is given by

G E (r) = cos(4θ) πr 2 (2D), (7)
G E (r) = 3 4π

r 2 (x 2 + y 2 ) -10x 2 y 2 r 7 (3D), ( 8
)
where r = (x, y) in 2D and r = (x, y, z) in 3D denote the position of a site relative to the yielded site [START_REF] Mura | Micromechanics of defects in solids[END_REF]. Notice that, since our system is defined on a discrete lattice, we need to consider a carefully discretized version of this propagator. First of all, we note that the value of G E (r = 0) is not defined because the Eshelby kernel describes the long-range elastic field and has to be regularized at the origin. For the discrete version, we can thus freely choose G i,i = G(r = 0) = -1, as discussed above. Next, we compute the Eshelby propagator in the continuous Fourier space [START_REF] Mura | Micromechanics of defects in solids[END_REF], ĜE (q 2 x , q 2 y ) = -

4q 2 x q 2 y q 2 x + q 2 y 2 (2D), ĜE (q 2 x , q 2 y , q 2 z ) = - 4q 2 x q 2 y q 4 - q 2 z q 2 (3D). ( 9
)
Because we consider a finite system with a linear box length L and periodic boundary conditions, we need to discretize the Fourier space as well, following the same procedure as in Ref. [52]. First, the wave-vector components are written as q µ = (2π/L)n µ , with µ = x, y(, z) and n µ = -L/2 + 1, • • • , L/2. Second, we consider a discrete Laplacian instead of the continuous one, which amounts to replacing q 2 µ → 2 -2 cos q µ in Eq. ( 9). Moreover, Eq. ( 9) is not defined at q = 0 because the Eshelby propagator is not integrable in infinite volume, but in the discrete version we can normalize it in such a way that Ĝq=0 = -1, as discussed above. These two conditions, Ĝq=0 = -1 and G i,i = -1, can be simultaneously satisfied by normalizing the Eshelby propagator in Fourier space [40]. Our final discrete propagator is then Ĝq

= -1 , q = 0 , ĜE (q 2 µ →2-2 cos qµ) G , q = 2π L n = 0 , (10)
with the constant G determined by the condition,

G i,i = 1 N q Ĝq = -1 , (11)
where the sum is over the discretized wave-vectors. The propagator in discrete real space is finally obtained from Ĝq by a discrete inverse Fourier transform.

Comparison with molecular simulations

To perform sample averages of the various observables, we use 1000 -2000, 400 -600, 400 -600, 200 -400, and 100 -200 samples for L = N 1/2 = 256, 512, 1024, 2048, and 4096, respectively, in 2D, and 800-2000, 800-1000, 800 -1000, 400 -1000, 200 -500, and 100 -200 samples for L = N 1/3 = 48, 64, 80, 104, 128, and 164, respectively, in 3D. These numbers are larger than those used in molecular dynamics (MD) simulation studies: in previous work by some of us, we averaged over 800, 700,400,200,200,200,200, and 100 samples for N = 1000, 2000, 4000, 8000, 16000, 32000, 64000, and 128000, respectively, in 2D [20], and 800, 400, 200, 100, 100, 50, and 25 -50 samples for N = 1500, 3000, 6000, 12000, 24000, 48000, and 96000, respectively, in 3D [19]. Moreover, we study a larger number of values of the disorder (R in EPMs) than in molecular simulations (in which disorder is encoded by T ini [19,20]).

It is also instructive to compare the system sizes achieved in MD and EPM simulations. This can be done by following two distinct approaches. The first one is measuring the typical size of a shear transformation zone, which corresponds to a building block (or single site) for the EPM mesoscopic description, in an MD simulation: see, e.g., [51]. The second one is a quantitative calibration or mapping from MD to EPM, such that the statistical properties and the macroscopic responses in both studies match: see, e.g., [START_REF] Liu | [END_REF]55]. All of this indicates that a single site in an EPM corresponds to of the order of magnitude of 100 particles (atoms) in an MD simulation. Although these studies focus on a rather ductile yielding regime, we expect that this order of magnitude does not change significantly in the brittle yielding regime. This issue is further discussed in Ref. [56]. Using this conversion, 16 millions of particles in an MD simulation as in Ref. [38] correspond to 160000 sites in EPMs. Our largest EPM in 2D (N = L 2 = 16777216) is thus about 100 times larger and in 3D (N = L 3 = 4410944) about 30 times larger than the MD counterpart in Ref. [38]. Therefore, effectively, our study accesses much larger system sizes than previous MD studies.

Definition of the order parameter

In order to characterize the presence of a critical point, we need to define a proper order parameter associated with the brittle-to-ductile transition. The most direct one is based on the presence of a discontinuity in the individual stress-versus-strain curve, namely, the stress drop, ∆σ(γ). One can then define ∆σ max = max γ {∆σ(γ)}, which amounts to finding the largest stress drop in the whole dynamics [19,20]. If one considers a system with a strong disorder, i.e., in the ductile regime, this quantity goes to zero as the system size increases. However, in the brittle regime a macroscopic stress drop is present, and as a result ∆σ max remains finite. Therefore, one can study how ∆σ max changes when the initial stability (or disorder strength) R is varied in order to distinguish the two yielding regimes.

The mean and variance of ∆σ max are presented in Fig. 7 in both 2D and 3D. As in Fig. 2 in the main text, we observe the signature of a critical point: ∆σ max is of order O(1) when R is small, while it goes to nearly zero when R is large, and more so as the system size L increases. Around the onset of the growth of ∆σ max , the (disconnected) susceptibility, N ∆σ 2 max -∆σ max 2 , shows a peak that grows with the system size. However, as we decrease R further, the susceptibility seems to diverge again with the system size. This unexpected behavior can be explained by considering a simple argument in the case of very weak disorder R. Let us assume that when R → 0 + the system evolves purely elastically up to the macroscopic yielding strain γ Y , which corresponds to the value of γ at which the least stable site yields. We now define the maximum stress of the initial distribution, σ max = max i {σ i (γ = 0)}. Consequently, we have µγ Y = σ th -σ max . If we also assume that the macroscopic stress goes to a plateau value σ ss after yielding, we obtain ∆σ max = µγ Y -σ ss . Therefore, we can relate the fluctuations of ∆σ max and σ max through

∆σ 2 max -∆σ max 2 = σ 2 max -σ max 2 . ( 12
)
Since the initial distribution is nearly a Gaussian, σ max follows a Gumbel distribution, associated with extreme-value statistics [57]. We then arrive at

∆σ 2 max -∆σ max 2 = π 2 R 2 12 ln N . ( 13
)
This expression explains the growth of the susceptibility at R → 0 in Fig. 7. In real experiments and molecular simulations, such a growth would not be present because in practice the elastic branch contains multiple plastic events and the purely elastic assumption taken above does not hold in these systems [10].

In order to avoid this problem in the elasto-plastic models, we have decided to consider a different order parameter. Since the presence of a macroscopic stress drop in the brittle regime is accompanied by the abrupt formation of a shear band, one can think of characterizing yielding by looking at the number of sites that yield in each horizontal or vertical line (in 2D) or plane (in 3D). More precisely, we define a quantity n x,y (γ) in 2D that is 1 if the site at (x, y) has already had a plastic event during its evolution up to γ, and 0 otherwise (and similarly for n x,y,z (γ) in 3D). We can then define

n x (γ) = 1 L L y=1 n x,y (γ) (2D), (14)
n x (γ) = 1 L 2 L y,z=1 n x,y,z (γ) (3D), (15)
in 2D and 3D respectively. In the same way we also define n y (γ) for vertical lines in 2D, and n y (γ) and n z (γ) for vertical planes in 3D. We expect that, in the ductile regime, there is essentially no difference between n x (γ) and n y (γ) (and n z (γ)) since the plastic activity occurs rather uniformly throughout the sample and no spatial organization of the events is observed. On the other hand, in the brittle regime, there is a specific x * (or y * or z * ) such that n x * (γ) ≈ 1 at γ ≈ γ Y due to the formation of the shear band. This corresponds to a strong anisotropic localization of the plastic activity. One can then characterize whether the system is in a ductile or brittle phase by studying the evolution of these quantities. In particular we consider n(γ) = max x,y {n x (γ), n y (γ)} (and similarly in 3D). The inset of Fig. 4 in the main text shows the evolution of n(γ) for some samples. The discontinuous stress drop of σ(γ) in the brittle regime corresponds to the discontinuous jump of n(γ), whereas the mild continuous crossover of σ(γ) in the ductile regime corresponds to a continuous increase of n(γ). Therefore, n(γ) essentially contains the same information as σ(γ) as far as characterizing brittle and ductile yielding is concerned. As for ∆σ(γ), we compute the jumps, ∆n(γ), in this quantity, which allows us to define a new order parameter ∆n max , defined as ∆n max = max γ {∆n(γ)}. The results for ∆n max (γ) are presented in Fig. 2 in the main text.

Force balance

Since we use the AQS driving, imposing force balance is an important element in the simulation of the EPM. In the continuum limit, mechanical equilibrium imposes that the divergence of the stress vanishes. For example in 2D, one should have

∂σ xx (x, y) ∂x + ∂σ xy (x, y) ∂y = 0 . (16)
Integrating this relation over x, we obtain

dx ∂σ xy (x, y) ∂y = -dx ∂σ xx (x, y) ∂x . (17)
Because the right-hand side of Eq. ( 17) vanishes due to periodic boundary conditions for any y, we obtain

∂ ∂y dx σ xy (x, y) = 0 , (18)
implying that the average shear stress σ xy (x, y) must be constant over a row (or a column). A similar derivation can be obtained in 3D.

In the discrete version, averaging the local shear stress (here σ = σ xy ) evolution in Eq. ( 4) over x in 2D, we obtain

1 L x σ x,y → 1 L x σ x,y + ∆σ ext -∆σ p + 1 N qy =0
Ĝqx=0,qy ∆σ qx=0,qy e iqyy .

(

) 19
Because the Eshelby kernel in Fourier space vanishes whenever q x = 0 for all q y = 0, see Eq. ( 9), the second line vanishes and the average stress along all rows changes by the same quantity, ∆σ ext -∆σ p [START_REF] Liu | Critical dynamics at the yielding transition and creep behavior of amorphous systems: mesoscopic modeling[END_REF]. The same result is obtained for columns. This means that the difference between the average stress along a given row and along a different row or column remains constant during the dynamics. Therefore, to impose mechanical equilibrium, or force balance condition, it is sufficient to ensure that the sum of the stresses along rows and columns is the same for the initial stress distribution at γ = 0. For 3D, a similar argument holds for the average local shear stress over an x-z (or y-z) plane.

As already explained, we initialize the local stresses by a nearly Gaussian distribution with zero mean. Thus, the force balance condition is satisfied asymptotically at large L because the average stress over a single row or column goes to zero when γ = 0 and L → ∞. However, this condition does not hold exactly in finite system sizes. We have thus tested another way to initialize the local stresses, which strictly enforces the force balance condition. The initial stress at each site is drawn from a Gaussian distribution of zero mean and standard deviation R, but we follow the procedure described in Ref. [40] and implement an additional random operation to keep the sum of the stresses equal along all rows and columns. We then use this initial condition for our EPMs, and we find that the results quickly converge to the results presented in the main text with increasing L. Therefore, our conclusions on the critical point are not affected by the choice of the force balance condition at γ = 0.

Variations of the model: random thresholds

We have also considered a different way to tune the disorder of the sample. We fix R = 0.45 for the initial condition and we let instead the stress threshold σ th i be randomly distributed [40]. In this model, the local stress threshold at each site, σ th i , when γ = 0 is chosen according to a Gaussian distribution of mean M and variance 0.01. After each plastic event, σ th i is updated from another Gaussian distribution of mean 1 and variance 0.01. Notice that in this model the system shows brittle or ductile yielding depending on the average value M of the initial σ th i . A large M corresponds to a stable or weakly disordered sample, whereas a small M corresponds to a less stable or highly disordered sample. Although this model contains some ingredients that are different than those of the model we have studied in the main text, such as the threshold distribution and its softening, the precise way to enforce force balance, etc., it essentially displays the same phenomenology. This is illustrated in Fig. 8 for the mean value of the order parameter ∆n max and the associated disconnected susceptibility, with 1/M playing the role of R. These numerical observations suggest that the critical point that we have identified in this work is quite robust and universal with respect to changes in the details of the simulated elasto-plastic models.

Determination of the onset of the stress overshoot

We determine the onset value R o at which the stress overshoot appears (or disappears) in the evolution with the applied strain of the stress averaged over many independent samples. We illustrate the operational procedure in Fig. 9 for the 2D elasto-plastic model with L = 2048. Three averaged stress-versus-strain curves in the vicinity of the first appearance of the overshoot are shown. From such a plot, one can detect R o in a reliable manner. We have performed this analysis for different values of L both in 2D and 3D. The outcome is summarized in Fig. 5 in the main text.

FIG. 1 .

 1 FIG. 1. Stress versus strain curves in the 3D elasto-plastic model illustrating the brittle (R = 0.3) and the ductile (R = 0.8) cases. The linear box length is L = 200. Three independent samples are presented. Insets: Real-space configurations at γ = 0.6 in the two cases. The color bar corresponds to the number of local plastic events.

FIG. 2 .

 2 FIG. 2. Evidence for a critical point in 2D and 3D EPMs. Upper: Average value of the order parameter ∆nmax as a function of R for several system sizes in 2D (a) and 3D (b). Lower: Variance of ∆nmax multiplied by N = L D , where D is the spatial dimensions, i.e., disconnected susceptibility, in 2D (c) and 3D (d).

FIG. 3 .FIG. 4 .

 34 FIG.3. Log-log plot of the maximum (a) and the full width at half maximum (b) of the disconnected susceptibility associated with the order parameter ∆nmax as a function of N = L D for both 2D (blue) and 3D (red). The straight black lines have slopes 0.9 in (a) and -0.15 in (b).

FIG. 5 .FIG. 6 .

 56 FIG. 5. Stress-versus-strain curves for 3D samples at fixed disorder strength R = 0.40 for several system sizes (a) and at fixed system size L = 128 for different values of R (b). Insets: The corresponding n(γ) curves.

FIG. 7 .

 7 FIG.7. Mean value of the order parameter ∆σmax as a function of R for several system sizes L in 2D (a) and 3D (b). The corresponding variance of ∆σmax, which corresponds to a disconnected susceptibility, is shown in (c) for 2D and (d) for 3D.

FIG. 8 .FIG. 9 .

 89 FIG. 8. Finite-size analysis of the critical point for the model described in Ref. [40]. (a) Average value of ∆nmax. (b) Associated disconnected susceptibility.

Acknowledgments: We thank Ludovic Berthier, Chen Liu, Sylvain Patinet, and Matthieu Wyart for discussions. We also thank Chen Liu for multiple advices for the numerical implementation of the elasto-plastic models. This work was granted access to the HPC resources of the MeSU platform at Sorbonne Université and the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement n. 723955 -GlassUniversality) and from the Simons Foundation (Grant No. 454935 to G.B.; Grant No. 454955 to F.Z.).


Daniele Coslovich

Misaki Ozawa

Walter Kob

Dynamic and thermodynamic crossover scenarios in the Kob-Andersen mixture: Insights from multi-CPU and multi-GPU simulations

Dynamic and thermodynamic crossover scenarios in the Kob-Andersen mixture: Insights from multi-CPU and multi-GPU simulations

Daniele Coslovich, Misaki Ozawa, Walter Kob

The physical behavior of glass-forming liquids presents complex features of both dynamic and thermodynamic nature. Some studies indicate the presence of thermodynamic anomalies and of crossovers in the dynamic properties, but their origin and degree of universality is difficult to assess. Moreover, conventional simulations are barely able to cover the range of temperatures at which these crossovers usually occur. To address these issues, we simulate the Kob-Andersen Lennard-Jones mixture using efficient protocols based on multi-CPU and multi-GPU parallel tempering. Our setup enables us to probe the thermodynamics and dynamics of the liquid at equilibrium well below the critical temperature of mode-coupling theory, T MCT = 0.435. We find that below T = 0.4 the analysis is hampered by partial crystallization of the metastable liquid, which nucleates extended regions populated by large particles arranged in an fcc structure. By filtering out crystalline samples, we reveal that the specific heat grows in a regular manner down to T = 0.38. Possible thermodynamic anomalies suggested by previous studies can thus occur only in a region of the phase diagram where the system is highly metastable. Using the equilibrium configurations obtained from the parallel tempering simulations, we perform molecular dynamics and Monte Carlo simulations to probe the equilibrium dynamics down to T = 0.4. A temperature-derivative analysis of the relaxation time and diffusion data allows us to assess different dynamic scenarios around T MCT . Hints of a dynamic crossover come from analysis of the four-point dynamic susceptibility. Finally, we discuss possible future numerical strategies to clarify the nature of crossover phenomena in glass-forming liquids.

I. INTRODUCTION

If a liquid is cooled quickly enough, it will bypass crystallization and form an amorphous solid called glass. Such a glass transition is directly related to the rapid increase of the structural relaxation time upon cooling, whose temperature-dependence is given by the Arrhenius law for strong glass-formers and is super-Arrhenius for fragile glass-formers [START_REF] Angell | Formation of glasses from liquids and biopolymers[END_REF] . Despite being purely kinetic in nature, the glass transition is accompanied by a change in the thermodynamic properties of the system at the glass transition temperature T g . For instance, the specific heat shows a mild increase upon cooling and displays a drop at T g due to the freezing of the configurational degrees of freedom. Although no phase transition occurs at T g , it has been suggested that the slowing down of the dynamics and the resulting glass formation might be driven by a hidden phase transition [START_REF] Binder | Glassy materials and disordered solids: An introduction to their statistical mechanics[END_REF][START_REF] Berthier | Theoretical perspective on the glass transition and amorphous materials[END_REF] .

While most glass-forming liquids follow the trends outlined above, there are notable exceptions. This is the case for instance of liquids with strong, directional interactions, such as silica, silicon, and water. These systems display a thermodynamic anomaly in the form of a local maximum of the specific heat at a temperature T * > T g [START_REF] Angell | Glass-formers and viscous liquid slowdown since David turnbull: Enduring puzzles and new twists[END_REF] , a behavior which has also been observed in simulation studies of silica [START_REF] Scheidler | Frequency-dependent specific heat of viscous silica[END_REF][START_REF] Saika-Voivod | Free energy and configurational entropy of liquid silica: Fragile-to-strong crossover and polyamorphism[END_REF] , supercooled water [START_REF] Saito | Frequency dependence of specific heat in supercooled liquid water and emergence of correlated dynamics[END_REF] , and other simple models [START_REF] Moreno | Energy Landscape of a Simple Model for Strong Liquids[END_REF][START_REF] Xu | A monatomic system with a liquid-liquid critical point and two distinct glassy states[END_REF][START_REF] Gutiérrez | The static lengthscale characterizing the glass transition at lower temperatures[END_REF][START_REF] Ozawa | Tuning pairwise potential can control the fragility of glass-forming liquids: from a tetrahedral network to isotropic soft sphere models[END_REF] . Concomitantly, also the dynamic behavior of these liquids changes around T * , crossing over from a super-Arrhenius to an Arrhenius temperature dependence ("fragile-to-strong" crossover). The origin of these thermodynamic and dynamic anomalies has often been attributed to the presence, or proximity, of a liquid-liquid transition or its Widom-line, located at temperatures higher than T g 4 . Interestingly, this fragileto-strong crossover is not limited to systems with directional interactions [START_REF] Mallamace | Transport properties of glassforming liquids suggest that dynamic crossover temperature is as important as the glass transition temperature[END_REF] . Recent experimental [START_REF] Zhang | Fragile-to-strong transition in metallic glass-forming liquids[END_REF][START_REF] Zhou | Structural evolution during fragile-to-strong transition in cuzr (al) glass-forming liquids[END_REF] and simulational [START_REF] Lad | Signatures of fragileto-strong transition in a binary metallic glass-forming liquid[END_REF] studies have demonstrated that also several metallic glass-formers exhibit a notable fragile-to-strong crossover. Several interpretations have been proposed to explain the crossover in these systems, including the presence of thermodynamic anomalies such as a liquidliquid transition [START_REF] Lad | Signatures of fragileto-strong transition in a binary metallic glass-forming liquid[END_REF][START_REF] Wei | Liquid-liquid transition in a strong bulk metallic glass-forming liquid[END_REF][START_REF] Stolpe | Structural changes during a liquid-liquid transition in the deeply undercooled Zr58.5Cu15.6Ni12.8Al10.3Nb2.8 bulk metallic glass forming melt[END_REF] , anomalous crystallization [START_REF] Yang | Anomalous crystallization as a signature of the fragileto-strong transition in metallic glass-forming liquids[END_REF] , or the evolution of structural medium range order [START_REF] Zhou | Structural evolution during fragile-to-strong transition in cuzr (al) glass-forming liquids[END_REF] .

The dynamics of some molecular and polymeric liquids show even a different kind of crossover at a temperature T D , typically located about 15%-25% above T g . This dynamic crossover thus occurs when the relaxation times are around 10 -8 -10 -6 s and is subtle, but it can be revealed by temperature-derivative analysis of high-quality dynamic measurements in molecular glass-formers [START_REF] Stickel | Dynamics of glassforming liquids. I. Temperature-derivative analysis of dielectric relaxation data[END_REF][START_REF] Martinez-Garcia | The new insight into dynamic crossover in glass forming liquids from the apparent enthalpy analysis[END_REF][START_REF] Novikov | Qualitative change in structural dynamics of some glass-forming systems[END_REF] . At T D the increase of relaxation times crosses over from super-Arrhenius to a milder temperature dependence and tend to become Arrhenius at very low temperature. While some authors have considered systems with such a crossover as marginal cases [START_REF] Elmatad | Corresponding States of Structural Glass Formers[END_REF] , an alternative point of view suggests that this behavior may fairly general [START_REF] Mallamace | Transport properties of glassforming liquids suggest that dynamic crossover temperature is as important as the glass transition temperature[END_REF] , even though the location of the crossover is highly system-dependent [START_REF] Martinez-Garcia | The new insight into dynamic crossover in glass forming liquids from the apparent enthalpy analysis[END_REF][START_REF] Novikov | Qualitative change in structural dynamics of some glass-forming systems[END_REF] .

The physical origin of this dynamic crossover is not completely clear and there are diverging viewpoints on this. Some authors have pointed out the closeness of T D with the temperature at which a power law fit predicts a divergence of the relaxation time data [START_REF] Mallamace | Transport properties of glassforming liquids suggest that dynamic crossover temperature is as important as the glass transition temperature[END_REF][START_REF] Casalini | Dynamic crossover in supercooled liquids induced by high pressure[END_REF][START_REF] Casalini | Viscosity at the Dynamic Crossover in o-Terphenyl and Salol under High Pressure[END_REF] . The dynamic crossover should then be identified with the critical temperature of mode-coupling theory (MCT) [START_REF] Götze | Complex dynamics of glass-forming liquids: A modecoupling theory[END_REF] . Others have attributed this crossover to an upper limit of the arXiv:1804.04559v2 [cond-mat.stat-mech] 7 Jun 2018 activation energy [START_REF] Novikov | Qualitative change in structural dynamics of some glass-forming systems[END_REF][START_REF] Elmatad | Corresponding States of Structural Glass Formers[END_REF] , which may saturate at an arbitrarily low temperature. While these two interpretations are completely different in nature, they are difficult to disentangle in practice. We emphasize that the dynamic crossover scenario is a priori unrelated from one that accompanies liquid-liquid transitions, in that the former does not involve the thermodynamics and is purely dynamic in origin.

Also computer simulations have been used to probe the existence of anomalies in the dynamic and thermodynamic properties of glass-forming liquids. The results of these studies suggest that both the dynamic correlation length scales [START_REF] Kob | Non-monotonic temperature evolution of dynamic correlations in glass-forming liquids[END_REF] and the dynamic finite size effects 27 may show a crossover compatible with several theoretical predictions [START_REF] Biroli | The random firstorder transition theory of glasses: A critical assessment[END_REF][START_REF] Rizzo | Qualitative features at the glass crossover[END_REF][START_REF] Rizzo | Dynamical landau theory of the glass crossover[END_REF] , although these results do not seem to hold universally. On the other hand, the dynamic range accessible in conventional numerical studies is limited to only 4-5 decades in time. Thus, some of these results might be affected by insufficient equilibration.

In the present study, we focus on a simple glass-former, proposed by Kob and Andersen (KA) [START_REF] Kob | Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture I: The van Hove correlation function[END_REF] , that has so far been fairly robust against crystallization and hence has been used in many investigations of the glass-transition. Previous studies of this system have given evidence for the presence of a thermodynamic anomaly in the form of a peak in the specific heat [START_REF] Flenner | Hybrid Monte Carlo simulation of a glass-forming binary mixture[END_REF] and a possible fragile-tostrong crossover [START_REF] Ashwin | Low-temperature behaviour of the Kob-Andersen binary mixture[END_REF][START_REF] Doliwa | Energy barriers and activated dynamics in a supercooled Lennard-Jones liquid[END_REF] , although these results might have been affected by finite size and finite sampling effects. In the present work we employ an optimized simulation setup, which exploits the parallel tempering method [START_REF] Hukushima | Exchange Monte Carlo method and application to spin glass simulations[END_REF] and state-of-art molecular dynamics code running on graphics processing units (GPU), to extend the temperature range in which thermodynamic and dynamic measurements can be done at equilibrium. We find that, in contrast to previous reports [START_REF] Flenner | Hybrid Monte Carlo simulation of a glass-forming binary mixture[END_REF] , the thermodynamics of the liquid is regular. At lower temperatures, simulations are hindered by crystallization, which involves structures of fcc symmetry formed by large particle. Finally, we discuss possible ways and setups to detect numerically the presence of dynamic anomalies in model glass-formers.

II. MODEL AND METHODS

A. Model parameters

The system we consider is a binary mixture in which both species have the same mass m. The particles interact via a Lennard-Jones potential given by

u αβ (r) = 4 αβ σ αβ r 12 - σ αβ r 6 , (1)
where α, β ∈ {A, B} are species indices. The value of the parameters σ αβ and αβ are given in Ref. [START_REF] Kob | Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture I: The van Hove correlation function[END_REF] . The units of length and energy are set by the parameters σ AA = 1 and AA = 1, respectively. The potentials are cut and shifted at a distance 2.5σ αβ . We simulate systems com-posed by N particles in a cubic box of side L with periodic boundary conditions and a number density given by ρ = N/V = 1.1998. Note that even small differences in density can quantitatively affect the thermodynamic and dynamic observables. For example, in the original paper [START_REF] Kob | Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture I: The van Hove correlation function[END_REF] a different density ρ = 1.204 was used, which leads to slight differences in statics and dynamics when sufficiently high quality data are available. The system size ranges from N = 300 to 3600 for parallel tempering simulations (see below). Additional dynamic data have been obtained for a much bigger sample (N = 100000) using the LAMMPS simulation package [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF][START_REF]LAMMPS: Large-scale Atomic/Molecular Massively Parallel Simulator[END_REF] .

B. Simulation protocols

Our simulations implement the parallel tempering (PT) algorithm [START_REF] Hukushima | Exchange Monte Carlo method and application to spin glass simulations[END_REF][START_REF] Yamamoto | Replica-exchange molecular dynamics simulation for supercooled liquids[END_REF] , in which M replicas of the system of interest perform independently simulations at temperatures {T i } with the potential energies {U i }. At regular intervals, exchanges are attempted between pairs of replicas at neighboring states and the temperatures are exchanged with probability

p = min {1, exp [(U i -U j ) (β i -β j )]} , (2)
where In the first implementation, named PT-1 protocol, we perform multi-CPU parallel tempering simulations with an in-house molecular dynamics code. The MD simulations are performed in the N V T ensemble using the Nose-Hoover thermostat [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF] with a time step δt = 0.004 and a thermostat relaxation time 0.4 = 100δt. We used M = 12 -14 replicas depending on system size. The code is parallelized using MPI to handle communications between replicas, which attempt to exchange their state, i.e., the temperature of the associated thermostat, every 50000 MD steps.

β i = 1/k B T i (with k B =
The second implementation, named PT-2 protocol, relies on a multi-GPU parallel tempering code [START_REF]atooms-pt: Multi-core / multi-GPU parallel tempering[END_REF] that uses the RUMD package 41 as a simulation backend. This multi-GPU code was implemented in python building on the atooms framework [START_REF]atooms: A python framework for simulations of interacting particles[END_REF] . Multiple replicas are simulated on individual GPUs and communication between GPUs are handled at high-level via the mpi4py package [START_REF]mpi4py[END_REF] . We ran the multi-GPU simulations on a dedicated cluster of inexpensive gaming cards (GTX-980 and even GTX-750Ti). The thermostat is again of the Nose-Hoover type, the time step is δt = 0.004 and the thermostat relaxation time τ T = 0.2 = 50δt. We used the same exchange intervals and number of replicas as in PT-1. We checked that increasing the interval between exchanges did not change our results. Finally, we performed additional PT simulations (PT-3 protocol) to extend our thermodynamic measurements to even lower temperatures than the ones attained by protocols PT-1 and PT-2. The PT-3 simulations are performed using the multi-GPU code described above starting from uncorrelated configurations obtained using the PT-2 protocol at a temperature T = 0.4. These efficient multi-GPU simulations enabled us to carefully measure the waiting time dependence of the results so as to assess equilibration issues, see Sec. III A.

From these PT simulations, we evaluated the specific heat per particle c V from the fluctuations of the potential energy U ,

c V = 1 N T 2 U 2 -U 2 + 3 2 , (3)
and from the temperature derivative of the average potential energy,

c V = 1 N ∂ U ∂T + 3 2 , (4)
where (• • • ) is the thermal average. The two expressions yield identical results provided the averages are carried out over the equilibrium measure. In a simulation, the agreement between the estimates of c V obtained through the two methods is often taken as a test of equilibration. We note that, in practice, only very accurate measurements can reveal discrepancies between the two methods. The distributions P (U/N ) and the T -dependence of U obtained from well-equilibrated simulations (PT-3 protocol) are shown in Fig. 1.

Strictly speaking, parallel tempering simulations only give access to thermodynamic and static properties. However, it is possible to perform extended dynamic measurements by starting "regular" simulations from configurations sampled during the PT simulations at a given temperature. Here, again, we follow two distinct protocols to corroborate our results.

In the MC protocol we performed normal Monte Carlo simulations using an in-house code and starting from uncorrelated configurations obtained from protocol PT-1. The MC simulations are carried out in the N V T ensemble using simple displacement moves [START_REF] Berthier | The Monte Carlo dynamics of a binary Lennard-Jones glass-forming mixture[END_REF] , in which we attempt to displace a randomly selected particles over a cube of side 0.15. We used 10 -30 independent configurations depending on temperature. The length of the simulations at the lowest temperature (T = 0.4) is 10 9 Monte Carlo steps. In the following, the time unit for the MC protocol is given by one MC sweep, i.e., N attempted displacement moves.

In the MD protocol we performed molecular dynamics simulations using the RUMD package starting from uncorrelated configurations obtained from protocol PT-1 and PT-2. The MD simulations are carried out in the N V T ensemble using the Nose-Hoover thermostat, the time step is δt = 0.004 and the thermostat relaxation time τ T = 0.2 = 50δt. For the N = 1200 samples, we used 128 independent configurations down to T = 0.4 from protocol PT-2. For T = 0.39, we used only 30 configurations, namely the final configurations of the PT-2 runs. For the N = 3600 samples, we used 20 independent configurations from protocol PT-1. The duration of each of the MD simulations ranged from 4.2 × 10 6 steps (at T = 0.51) to 2.1 × 10 9 steps (at T = 0.39), thus each run was about 10 -20 times longer than the typical structural relaxation time τ α (see below for its definition). In total, for each temperature, our simulations cover over about 2500 structural relaxation times. This high-quality statistics enables us to perform a temperature-derivative analysis of the dynamic data, see Sec. III B.

For both MC and MD protocols, we checked that the initial configurations were uncorrelated from one another by measuring their mutual self overlaps 45

Q s = 1 N i Θ(a -|r α i -r β i |), (5)
where α and β denote two configurations, and their mu- tual collective overlaps

Q c = 1 N i,j Θ(a -|r α i -r β j |). (6)
A sensible choice of parameter a is a fraction of the typical interparticle distance. We chose a = 0.3. We found that both Q s and Q c are close to the values expected for uncorrelated pairs of configurations, i.e., O(1/N ) and From MD and MC simulations we extract the self part of the intermediate scattering functions

F A s (k, t) = f A s (k, t) = 1 N A j e -ik•[rj (t)-rj (0)] , (7)
where the sum runs over the particles of type A. We choose a wave-vector k = 7.25, close to the first peak of the structure factor [START_REF] Kob | Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture I: The van Hove correlation function[END_REF] . The corresponding structural relaxation time τ α is defined as usual as F A s (k, τ α ) = 1/e. In Fig. 2, we show the dynamic data obtained from the MD protocol.

C. Crystalline order detection

The study of glass-forming liquids is often hampered by crystallization and the very relation between glassy behavior and crystallization remains a matter of debate [START_REF] Kawasaki | Formation of a crystal nucleus from liquid[END_REF][START_REF] Tanaka | Bond orientational order in liquids: Towards a unified description of water-like anomalies, liquid-liquid transition, glass transition, and crystallization[END_REF] . The KA mixture, which is a simple model of a metallic glass-former, has been extensively used as a model to study the glass transition because of its stability against crystallization. Until very recently, the note added in the proofs of Ref. 48 was, to the best of our knowledge, the only report of crystallization of this model by direct simulation, achieved through runs of about 3.7×10 7 time units (7.4×10 9 steps) at T = 0.40. At this temperature, however, the nucleation time is still much larger than the structural relaxation time τ α (∼ 10 5 time units), and therefore MD / MC simulations of the metastable liquid can be carried out safely. In this work, however, we were able to equilibrate the mixture at even lower temperatures. Below T = 0.4, crystallization events become increasingly frequent, as also demonstrated by a very recent simulation study [START_REF] Ingebrigtsen | Crystallisation instability in glassforming mixtures[END_REF] . Within the studied range of system sizes, the smaller the system, the stronger the tendency to crystallization.

As a first indicator of crystallization events in our simulations, we monitored the evolution of the inherent structure (IS) energy as a function of time [START_REF] Stillinger | Hidden structure in liquids[END_REF] . However, the IS energy may also display large but reversible fluctuations, unrelated to crystal nucleation. We thus studied two additional order parameters that allow us to disentangle "amorphous" and crystalline fluctuations. The first one relies on the so-called common neighbor analysis (CNA) [START_REF] Honeycutt | Molecular dynamics study of melting and freezing of small Lennard-Jones clusters[END_REF] . In this approach, the bonds formed by neigh- boring particles are classified according to the number of shared neighbors. It has been shown that the fraction f 142 of bonds of type 142, see e.g. Ref. 52, allows one to detect crystallization in biased simulations of the KA mixture. We found that this approach allows one to detect crystallization in the bulk mixture as well. An example of a crystallization event is depicted in Fig. 3. Even though the nature of the fluctuation is not always clear-cut, we found that a threshold on f 142 is an effective criterion to filter out crystalline configurations. Note that since we run several replicas at a time, only a few of them may be affected by crystallization. When this occurs they typically remain "stuck" in the lowest portion of temperature space. As in Ref. 52, we used a threshold of 12% CNA-142 bonds. A large crystalline cluster detected in our simulations is shown in Fig. 3. The crystal nucleus is formed by fcc pockets of A particles, which implies compositional fluctuations that deplete B particles. To detect it, we introduced an even simpler order parameter, namely the concentration of cages formed by A particles only [START_REF] Hedges | Dynamic Order-Disorder in Atomistic Models of Structural Glass Formers[END_REF] . In particular, we evaluated the connected component [START_REF] Newman | Networks: An Introduction[END_REF] of clusters formed by pure-A cages. We found that in typical crystalline samples, the size of these connected clusters is about a few hundred particles. Finally, in Fig. 3(c) we show the percentage of crystalline configurations in the simulations of protocol PT-2 and PT-3 for N = 1200 particles. We see that they increase markedly below T = 0.4. To deal with data analysis in this delicate regime while retaining the maximum possible amount of statistics, we filtered our c V measurements by discarding individual configurations whose CNA-142 concentration was higher than 0.12, see Sec. III A. On the other hand, above T = 0.4 the fraction of samples above the crystalline threshold is negligible and we rarely encountered problematic runs. These runs were simply discarded all together.

III. RESULTS

A. Thermodynamics

The specific heat c V is a sensitive, although not unambiguous, indicator of thermodynamic changes in dense liquids. For instance, a sudden drop in c V as a function of the control parameters may be due to a phase transition, or more generally to a change in the topography of the underlying energy landscape, but might as well indicate incomplete equilibration, as is the case at the laboratory glass transition observed upon cooling. Previous work on the KA mixture showed the presence of a peak in c V at some temperature close to the mode-coupling temperature T MCT = 0.435 [START_REF] Kob | Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture I: The van Hove correlation function[END_REF] . In this section, we show that the thermodynamics of the KA mixture is regular, i.e., it shows no anomaly, at least down to T = 0.39 and that the peak found in Ref. 32 is due to partial equilibration. Below T = 0.4, equilibration becomes hard on current simulation time scales and analysis is further hampered by crystallization events. Finally, we discuss the possible presence of a thermodynamic anomaly in this highly metastable portion of the phase diagram.

In Fig. 4(a) we show specific heat measurements from parallel tempering simulations using protocols PT-1 and PT-2. We also include results from a previous study using PT simulations, which reported a peak in c V around T MCT [START_REF] Flenner | Hybrid Monte Carlo simulation of a glass-forming binary mixture[END_REF] . From our data, we conclude that no anomaly is observed in the specific heat, which increases monotonically down to T = 0.39. The length of our simulations is typically one order of magnitude longer than those of Ref. 32. In retrospect, our work warns that tests such as histogram reweighting or the consistency of fluctuation and derivative expressions of response functions, such as c V , are not sufficient to ensure equilibration of supercooled liquids, see Ref. 54 for a more detailed discussion of this issue. Also, our results imply that one needs at least ∼ 10 9 MD time steps to equilibrate the KA mixture below T MCT .

The empirical model of Rosenfeld-Tarazona [START_REF] Rosenfeld | Density functional theory and the asymptotic high density expansion of the free energy of classical solids and fluids[END_REF] proposes the following functional form for the potential energy U = aT 3/5 + b, which yields c V ∼ T -2/5 . In Fig. 4(b) we draw the specific heat data from protocols PT-1 and PT-2 as a function of T -2/5 , which linearizes the Rosenfeld-Tarazona law. A similar representation yields an excellent data collapse for simple liquids in both normal and moderately supercooled regime [START_REF] Ingebrigtsen | Communication: The rosenfeld-tarazona expression for liquids specific heat: A numerical investigation of eighteen systems[END_REF] . We find that this functional form provides a very good description of the data but a slight upward bending is observed around T = 0.4, suggesting the presence of additional fluctuations not accounted for by this simple liquid-state model. This is confirmed by the analysis of the third moment of the potential energy distribution. We found that the skewness remains constant at high temperature and starts to increase slightly around T = 0.4 (not shown), even when crystalline configurations are removed from the analysis. This behavior may be thus attributed to crystallization precursors or to enhanced fluctuations of the locally favored structure of the system [START_REF] Coslovich | Locally preferred structures and many-body static correlations in viscous liquids[END_REF][START_REF] Turci | Nonequilibrium phase transition in an atomistic glassformer: The connection to thermodynamics[END_REF] .

We found that the measurements using PT-1 and PT-2 start to deviate significantly from each other below T = 0.39, Moreover, in this temperature regime, the liquid has an increased tendency to crystallize, see Sec. II. In Fig. 4(c) we show measurements of c V from protocol PT-3, in which crystalline samples are removed on a per-configuration basis and equilibration is assessed directly, i.e., by measuring the waiting time dependence of c V . Specifically, we start parallel tempering simulations from previously equilibrated configurations at T = 0.4 and perform averages over restricted portions of the trajectories as follows

A tw = 1 t p -t x tw+tp tw dtA(t)W (t), ( 8
)
where t w is the waiting time and W (t) = Θ(0.12 -f 142 ) is a windowing function that removes from the averages samples identified as crystalline (see Sec. II) and t x = tw+tp tw dt(1 -W (t)). The production time t p was kept fixed to 36000 PT steps, independent of waiting time.

In Fig. 4(c) we show the results obtained by calculating the specific heat starting from time t w with increasing values of t w . We emphasize that these specific heat measurements cover a temperature regime that has never been probed before at equilibrium. For temperatures higher than 0.4, the results of the PT-3 protocol show a consistent, i.e., t w -independent, growth of c V , thus corroborating our previous analysis. Below T = 0.4, the specific heat measurements display a more marked waiting-time dependence but fluctuation and derivative formulas converge for sufficiently long waiting times and reveal a monotonic increase of c V down to T = 0.38. Only at the lowest temperature, the specific heat c V measured from energy fluctuations shows a local maximum. The dependence of this peak as a function of waiting time suggests that this feature may be actually due to lack of equilibration.

Recent simulations based on trajectory path sampling [START_REF] Turci | Nonequilibrium phase transition in an atomistic glassformer: The connection to thermodynamics[END_REF] suggest the existence of a liquid-liquid transition in the low temperature part of the phase diagram of the KA mixture. Our data narrow down the temperature range over which this hypothetical transition may occur and rule out thermodynamic anomalies for T > 0.37. These results do not exclude, however, the scenarios discussed in Ref. 58. Given the strong tendency to crystallize below T = 0.4, however, our results show that a thermodynamic anomaly, if present at all in the KA mixture, is hidden in a highly metastable portion of the phase diagram. In practice, it will be very difficult to detect it in large enough samples. Therefore, we think that future work in this context should focus on more robust models of glass-formers.

B. Dynamics

The parallel tempering algorithm allows one to accelerate the sampling of static and thermodynamic observables, but it does not provide by itself useful information on the dynamics. However, it is possible to carry out dynamic measurements by using uncorrelated configurations extracted from PT simulations as starting points of conventional molecular dynamics or Monte Carlo simulations. This approach can be parallelized in a trivial way by performing independent simulations and allows one to extend the accessible dynamic range by short-circuiting equilibration issues. In this section, we implement these ideas following using the MD and MC protocols described in Sec. II and test various dynamic crossover scenarios below the MCT temperature.

Figure 5(a) shows an Arrhenius representation of the relaxation times obtained using the MD and MC protocols for 1200 and 3600 particles. We note that equilibrium sampling is ensured by the fact that both MD and MD simulations start from previously equilibrated configurations down to T = 0.39. The time scale of Monte Carlo simulations has been scaled to match the relaxation times of MD at low temperature. We also include dynamic data from Ref. 57, which agree with the new set of simulations over the common temperature range. Very similar behavior is found also for B-type particles (not shown). Thanks to the hybrid protocols employed in the present work, the accessible dynamic range has increased by almost two orders of magnitude compared to the conventional molecular dynamics simulations of Ref. 57 and we can comfortably study the dynamics be- low T MCT . The data suggest that the temperature dependence of the relaxation times gets milder, i.e., more Arrhenius-like, at the lowest temperatures. However, it is notoriously difficult to draw firm conclusions based on analysis of the relaxation time alone, which has led to a number of controversies [START_REF] Mallamace | Transport properties of glassforming liquids suggest that dynamic crossover temperature is as important as the glass transition temperature[END_REF] . Temperature-derivative analysis of the dynamic data provides a most stringent test of analytic expressions for τ α (T ) and is particularly well-suited to reveal the presence of a dynamic crossover [START_REF] Stickel | Dynamics of glassforming liquids. I. Temperature-derivative analysis of dielectric relaxation data[END_REF][START_REF] Martinez-Garcia | The new insight into dynamic crossover in glass forming liquids from the apparent enthalpy analysis[END_REF][START_REF] Novikov | Qualitative change in structural dynamics of some glass-forming systems[END_REF] . Even though this approach requires very accurate data, it has the advantage of being parameter free and requires no data fitting. So far, it has mostly been applied to high-quality dielectric relaxation measurements. To the best of our knowledge, the only numerical study to conduct this analysis across the MCT temperature is Ref. 59 for a mixtures of harmonic spheres. Here, we push this kind of analysis even further by combining the trivial parallelism of our MD and MC protocols and the efficiency of the RUMD simulation package on small system sizes [START_REF] Bailey | RUMD: A general purpose molecular dynamics package optimized to utilize GPU hardware down to a few thousand particles[END_REF] .

The central quantity in our analysis is the apparent activation energy

E(T ) = d ln τ α d(1/T ) , (9)
which we compute by the centered difference method 60 . A graph of E(T ) versus 1/T provides a simple way to test the parabolic law proposed by Elmatald et al. [START_REF] Elmatad | Corresponding States of Structural Glass Formers[END_REF] in the context of dynamic facilitation

τ α = τ 0 exp J T 0 2 T 0 T -1 2 , (10)
where τ 0 , J, and T 0 are material parameters. The activation energy is then a linear function of 1/T given by the following expression

E p (T ) = 2J 2 T 0 (T 0 /T -1) . (11)
It should be emphasized that Eq. ( 10) is only expected to hold below the onset temperature T 0 and above an additional reference temperature T x . Outside this range of temperatures, the dynamics is expected to be Arrhenius [START_REF] Elmatad | Corresponding States of Structural Glass Formers[END_REF] .

A similar approach can be used to linearize the classic VFT equation. Following Stickel et al. [START_REF] Stickel | Dynamics of glassforming liquids. I. Temperature-derivative analysis of dielectric relaxation data[END_REF] we introduce φ = E(T ) -1/2 . By computing the derivative of the VFT expression

τ α = τ 0 exp T vft K(T -T vft ) , (12)
we obtain

φ(T ) = K T vft 1 - T vft T , (13)
where τ 0 , K, and T vft are parameters. Thus, the VFT equation appears linear in a graph where φ is shown as a function of 1/T 19 . In Fig. 5(b) we show E(T ) as a function of 1/T . Overall, the activation energy displays an approximately linear behavior at sufficiently low temperature, i.e., T < 0.5. This behavior is well captured by Eq. ( 11), although fitting the relaxation times in this range requires a value T 0 ≈ 0.72, which is lower than the usual estimates 0.8 -1.0, see e.g. Ref. 61, but consistent with a very recent analysis 62 based on Eq. (11). We also note that the value of J/T 0 = 3.7 is larger than the one used in Ref. 22. The fit parameters used in Ref. 22 appear to strike a balance between the high and low T portions of the data.

For the most accurate data set, i.e., MD with N = 1200, and at the lowest temperatures, this representation also reveals a slight inflection in E(T ). This slight saturation may suggest a dynamic crossover similar to the one observed in certain molecular glass-formers [START_REF] Novikov | Universality of the dynamic crossover in glass-forming liquids: A "magic" relaxation time[END_REF] and in the simulation of harmonic spheres [START_REF] Kob | Nonlinear dynamic response of glassforming liquids to random pinning[END_REF] . One comment on possible finite size effects is in order. In a small sample, the activation energy necessarily reaches an upper bound and a smooth crossover to Arrhenius behavior is expected at a system-size dependent temperature [START_REF] Kim | Apparent finite-size effects in the dynamics of supercooled liquids[END_REF] . Within the quality of our data the temperature dependence for different system sizes gives compatibles results. However, due to the larger scatter in the N = 3600 samples, we cannot completely rule out the possibility that the inflection of the N = 1200 is due to a finite size effect. Pushing our analysis one step further, we compute the derivative of E(T ) with respect to 1/T and show the results in Fig. 5(c). This quantity should be constant and equal to 2J 2 in the dynamic facilitation scenario, whereas other models [START_REF] Novikov | Universality of the dynamic crossover in glass-forming liquids: A "magic" relaxation time[END_REF] predict that it should peak around some crossover temperature. The inset of the figure illustrates schematically these two possible scenarios. Within the noise of the data, and in particular of our most accurate data set (MD with N = 1200), our measurements are compatible with a decrease of dE(T )/d(1/T ) below some temperature T D close to T MCT . Simulations at even lower temperature and better statistics for larger system sizes would be needed to fully confirm this behavior.

In Fig. 6(a) we test the validity of the VFT law by computing φ(T ) as a function of 1/T . We find that the VFT law holds well at intermediate temperatures, but clear deviations are visible at low temperature. This confirms the well-known observation [START_REF] Stickel | Dynamics of glassforming liquids. I. Temperature-derivative analysis of dielectric relaxation data[END_REF] that, except in rare cases, the VFT equation cannot describe the full T de- pendence of the dynamic data. In Fig. 6(b) we show a plot E(T )/E p (T ) vs 1/T . In this representation, the parabolic law would correspond to a horizontal line. We see that the data flatten out only below T = 0.52, thus the range over which this law holds in the KA mixture is actually more limited than previously thought [START_REF] Elmatad | Corresponding States of Structural Glass Formers[END_REF] .

We also computed the activation energy E D from the diffusive time τ D , defined as the time needed to reach a mean square displacement of A particles equal to 1. These diffusive times are shown in Fig. 7(a) for several protocols and system sizes. We then compute, as before, the activation energies E D associated to diffusive times τ D . Figure 7(b) indicates that the growth of E D gets milder below some crossover temperature T D 0.45. Again, our N = 1200 MD data-set suggests an inflection of E D around this temperature. This crossover temperature is remarkably close to the one reported long ago by Ashwin and Sastry in Ref. 33 for a smaller system size (N = 256). However, in small systems E D actually saturates at T D , while bigger samples only cross over to a milder temperature dependence. These results indicate that, even though the results of Ref. 33 were probably affected by finite size effects, the diffusion mechanism might actually change in a temperature range close to the MCT transition temperature. This is corroborated by the analysis of a much larger sample (N = 100000) simulated using LAMMPS, which gives a trend compatible with the ones observed for N = 1200 and N = 3600 particles. These large-scale simulations were performed in the N V T simulations using the Nose-Hoover thermostat but without parallel tempering.

Finally, we compute the dynamic fluctuations associ- Dynamic susceptibility χ4(t) from MD protcol (N V T ensemble) for N = 1200 particles. The full lines are measurements over the entire length of the runs, while the symbols are measurements restricted to the first half of the runs. The error bars are estimated from the difference between the first and second halves of the runs. The crosses correspond to χ N V E 4 (t) + T 2 dQs dT (t) 2 /cV , see Eq. ( 16), at

T = 0.4877.
ated to the time-dependent self overlap function

Q s (t) = N -1 j Θ(0.3 -|r j (t) -r j (0)|), (14)
where the sum is taken over all particles. The dynamic susceptibility is then defined as

χ 4 (t) = N {[ Q s (t) 2 ] -[ Q s (t) ] 2 }, (15)
where [. . . ] indicates an average over statistically independent initial samples and . . . a time average over a given run. The dynamic susceptibilities, shown in Fig. 8 for N = 1200, display a peak at a time τ 4 proportional to the structural relaxation time. The peak height χ * 4 is a standard proxy for the extent of dynamic heterogeneity in supercooled liquids [START_REF] Berthier | Theoretical perspective on the glass transition and amorphous materials[END_REF] .

In Figs. 9(a) and (b) we show χ * 4 as a function of 1/T and as a function of τ * 4 , respectively. The dynamic fluctuations quantified by the susceptibility in the N V T ensemble cross over around the MCT temperature to a very mild dependence on both temperature and relaxation time. Within the quality of our data, this effect cannot be attributed to the finite size of the system. We point out, however, that the MCT scaling expected in the N V T ensemble [START_REF] Berthier | Spontaneous and induced dynamic correlations in glass formers. II. Model calculations and comparison to numerical simulations[END_REF]

, χ * 4 ∼ τ γ/2 α ∼ τ γ/2
4 , where γ ≈ 2.4 is the exponent of power law fit to the relaxation time, does not hold well for this system. Therefore, the fact that the crossover in the dynamic susceptibility occurs close to the mode-coupling temperature may be coincidental. We also note that similar results may not hold for the full dynamic susceptibility 3 , which includes contributions from number and concentration fluctuations and which we have neglected here. Our results resonate with those of Flenner and Szamel 66 , who showed that some contributions to the total susceptibility, notably those associated to dynamic fluctuations in the N V E ensemble, may be more sensitive to the presence of a crossover than the full susceptibility. To address this point quantitatively, we computed the dynamic susceptibility in the N V E ensemble from the exact expression 67

χ N V E 4 (t) = χ N V T 4 (t) - T 2 c V dQ s dT (t) 2 . ( 16
)
We checked that the above expression yields results consistent with those obtained from simulations in the N V E ensemble at T = 0.4877, at which the energy drift during N V E simulations was not too severe, see Fig. 8. The temperature dependence of the peak height of χ N V E 4 (t), shown in Fig. 9, provides clear hints of a dynamic crossover around the MCT temperature and corroborates the findings of Ref. 66. Royall et al. [START_REF] Royall | Strong geometric frustration in model glassformers[END_REF] have also noted that the growth of χ * 4 observed in simulations at above T MCT is incompatible with the typical correlation lengths measured experimentally. From this point of view, our finding of a crossover to a much milder rate of growth of dynamic correlations is a welcome result, which may solve the apparent conundrum evidenced in Ref. 68.

IV. CONCLUSIONS

In this work we have performed extensive computer simulations to probe the existence of thermodynamic and dynamic anomalies in the Kob-Andersen binary Lennard-Jones mixture. Our simulations build on an efficient simulation setup, which combines algorithmic and hardware optimizations. In particular, we exploit the parallel tempering algorithm and multi-GPU acceleration to measure the specific heat at equilibrium down to unprecedented temperatures. We also appreciably extend the range of dynamic measurements by starting regular MD and MC simulations from independent configurations previously equilibrated with parallel tempering. This approach is trivially parallel and enables us to perform a temperature-derivative analysis of the dynamic data.

Thanks to these advances, we could clarify some issues related to the thermodynamic and dynamic behavior of the mixture. In particular, we found that the specific heat increases monotonically down to at least T = 0.38, in contrast with previous findings 32 , which suggested a maximum around T = 0.43 and which was probably due to out-of-equilibrium effects. Although crystallization renders the analysis difficult at low temperature, our data indicate that thermodynamic anomalies, if present at all, must occur in a highly metastable, lowtemperature portion of the phase diagram. Conversely, any possible crossover above T = 0.38 must be purely dynamic in origin. By performing a temperature-derivative analysis of the dynamic data and by analyzing the dynamic susceptibilities, we have assessed several scenarios that may involve a dynamic crossover in this temperature regime. One possible scenario suggests the presence of a crossover [START_REF] Berthier | Finite-size effects in the dynamics of glass-forming liquids[END_REF][START_REF] Rizzo | Dynamical landau theory of the glass crossover[END_REF] around the mode-coupling temperature T MCT . There are however alternative interpretations that attribute the dynamic crossover to a saturation of energy barriers, as is the case in elastic models [START_REF] Dyre | Colloquium: The glass transition and elastic models of glass-forming liquids[END_REF][START_REF] Mirigian | Elastically cooperative activated barrier hopping theory of relaxation in viscous fluids. i. general formulation and application to hard sphere fluids[END_REF] , freevolume models [START_REF] Cohen | Liquid-glass transition, a freevolume approach[END_REF] and even possibly in the dynamic facilitation scenario at sufficiently low temperature [START_REF] Elmatad | Corresponding States of Structural Glass Formers. II[END_REF] . In these latter scenarios, the crossover need not be around T MCT Although our analysis is not entirely conclusive yet, we argue that state-of-the-art simulation methods combined with temperature-derivative analysis of simulation data may hold the key to fully disentangle these scenarios in the near future. Judging from the quality of our temperature-derivative analysis and dynamic fluctuations measurements, we infer that simulations of the order of several thousands of structural relaxation times are needed to fully confirm or rule out the presence of dynamic crossovers. Our results also suggest that contributions to the dynamic susceptibility measured in ensembles where conserved quantities are not free to fluctuate may better probe the dynamic behavior than the full susceptibility [START_REF] Berthier | Spontaneous and induced dynamic correlations in glass formers. II. Model calculations and comparison to numerical simulations[END_REF] . For specific systems, specialized Monte Carlo moves, like particle swaps [START_REF] Gazzillo | Equation of state for symmetric non-additive hard-sphere fluids: An approximate analytic expression and new Monte Carlo results[END_REF][START_REF] Sindzingre | Calculation of partial enthalpies of an argon-krypton mixture by NPT molecular dynamics[END_REF][START_REF] Grigera | Fast Monte Carlo algorithm for supercooled soft spheres[END_REF][START_REF] Ninarello | Models and algorithms for the next generation of glass transition studies[END_REF][START_REF] Berthier | Equilibrium sampling of hard spheres up to the jamming density and beyond[END_REF] , may provide an additional and important efficiency improvement and facilitate thermodynamic and dynamic studies below the mode-coupling temperature.

FIG. 2 .

 2 FIG. 2. Self part of the intermediate scattering functions F A s (k, t) obtained from the MD protocol for N = 1200 particles. Errors bars are one standard deviation on the mean, calculated over about 128 runs.

4 3

 4 πa 3 ρ, respectively.

FIG. 3 .

 3 FIG. 3. Detection of crystalline configurations. (a) Crystallization event during a PT-2 simulation for N = 1200 particles. The fraction of CNA-142 bonds, f142, is shown as a function of PT steps (1 PT step=50000 MD steps). (b) The probability density calculated during the run shows a bimodal distribution with a sharp minimum around f142 ≈ 12%. (c) Percentage of crystalline configurations, detected using a 12% threshold on f142, as a function of temperature. Below T = 0.4, both PT-2 and PT-3 protocols have a large fraction of crystalline configurations. PT-2 data were not analyzed above T = 0.4. (d) Connected component of a cluster formed by particles surrounded by A-particles only during a crystallization event.

FIG. 4 .

 4 FIG. 4. (a) Specific heat from protocols PT-1 and PT-2 for N = 1200 particles in the regime T 0.4 where crystallization is negligible. Results from fluctuation and derivative expressions of the specific heat are shown as indicated in the legend. The temperature of the data from Ref. 32 (circles) is divided by the density scaling factor (1.204/1.2) γ , with γ = 5.0 to correct for the small density mismatch. (b) Rosenfeld-Tarazona scaling of the same specific heat cV data as in panel (a). (c) Specific heat from protocol PT-3 for different waiting times tw expressed in unit of PT steps. Results from fluctuation and derivative expressions of the specific heat correspond to circles and squares, respectively. Configurations with f142 > 0.12 are discarded from the calculations.

FIG. 5 .

 5 FIG. 5. Temperature-derivative analysis of relaxation times. The relaxation times obtained from MC simulations have been scaled by 5 × 10 2 MC steps. Error bars on MD data represent one standard deviation on the mean and are only shown when larger than the symbol size. (a) Relaxation times τα as a function of 1/T for various protocols and system sizes. (b) Activation energy E(T ) for the same set of data as panel (a). The dashed line indicates a linear fit to Eq. (11) in the range T ≤ 0.52 with J = 2.7 and T0 = 0.72. (c) Derivative of E(T ) with respect to 1/T . The horizontal line corresponds to 2J 2 . The inset shows a schematic representation of two possible scenarios for this derivative: the solid line depicts the behavior expected from facilitation models, while the dashed line is the qualitative behavior described in Ref.21.

FIG. 6 .

 6 FIG. 6. Critical assessment of the VFT and parabolic laws. (a) Stickel plot φ(T ) versus 1/T . (b) Plot of E(T )/Ep(T ) as a function of 1/T . The dashed lines are guides to the eye. The arrow indicates the temperature T = 0.52 from which the parabolic law begins to apply.

FIG. 7 .

 7 FIG. 7. Analysis of diffusive times τD. The diffusive times obtained from MC simulations have been scaled by 5 × 10 2 MC steps. (a) Arrhenius representation of the diffusive times τD. (b) Activation energies from the diffusive time τD. The vertical line indicates the crossover observed by Ashwin and Sastry in Ref. 33 for a system of 256 particles.

τ 4 FIG. 9 .

 49 FIG. 9. Peak height of the dynamic susceptibility χ * 4 in the N V T ensemble (full symbols) and in the N V E ensemble (empty symbols) as a function of (a) 1/T and (b) the peak time τ4. The vertical line in the two panels marks the MCT temperature. The dotted lines in panel (b) indicate the predicted MCT scaling in the N V T (exponent 2/γ) and N V E ensemble (exponent 1/γ). Squares in panel (a) are taken from from Ref. 69. Squares in panel (b) indicate results from additional, high temperature N V T simulations for N = 1200 particles.

ACKNOWLEDGMENTS

We thank E. Flenner and G. Szamel for useful discussions. We acknowledge PRACE for awarding us access to Curie at GENCI@CEA, France, which enabled the development of part of the software used in this work. Data relevant to this work can be accessed at https://doi.org/10.5281/zenodo.1227831.

The final publication is available at Springer via http://dx.doi.org/10.1140/epje/i2018-11671-2, albeit with bitmapped figures.

AUTHOR CONTRIBUTION STATEMENT

The simulations were performed as follows: protocol PT-1 and LAMMPS by WK, protocol PT-2, PT-3, and MD by DC and MC by MO. DC and MO analyzed the data. DC wrote the first draft of the article. All authors discussed the results and contributed to the final manuscript.


Misaki Ozawa

Walter Kob

Atsushi Ikeda

Kunimasa Miyazaki

Reply to Chakrabarty et al.: Particles move even in ideal glasses

In their letter, Chakrabarty et al. [START_REF] Chakrabarty | Vanishing of configurational entropy may not imply an ideal glass transition in randomly pinned liquids[END_REF], point out that their data on the relaxation dynamics are inconsistent with the thermodynamic data presented in our paper [START_REF] Ozawa | Equilibrium phase diagram of a randomly pinned glassformer[END_REF]. They argue that from their results and the predictions of the random first-order transition theory [START_REF] Cammarota | Ideal glass transitions by random pinning[END_REF] one must conclude that our configurational entropy s c is quantitatively not accurate. In the following we will show that this conclusion is not necessarily valid.

The main argument of Chakrabarty et al. [START_REF] Chakrabarty | Vanishing of configurational entropy may not imply an ideal glass transition in randomly pinned liquids[END_REF] (figure 1, Left, of Ref. [START_REF] Chakrabarty | Vanishing of configurational entropy may not imply an ideal glass transition in randomly pinned liquids[END_REF]) is that the self part of the intermediate scattering function F s (k, t) decays to zero even in the glass phase (defined by s c = 0) and that hence this phase cannot be nonergodic. Although this observation is interesting, it does not imply the mentioned conclusion. This can be seen by considering the case of a simple crystal. In a crystal the diffusion constant of a tagged particle is nonzero (because a crystal at nonzero temperatures will have a finite concentration of defects) and hence F s (k, t) decays to zero at long times. From this one can, however, not conclude that a crystal is an ergodic system because its structure is obviously independent of time (i.e., the collective correlation functions do not decay to zero).

Thus, also in the case of the ideal glass transition in the pinned system one has to study not only the time dependence of the self functions, but also the one of the collective functions. This is done in FIG. 1, where we compare the time dependence of these two functions for different values of c at the temperature at which Chakrabarty et al. [START_REF] Chakrabarty | Vanishing of configurational entropy may not imply an ideal glass transition in randomly pinned liquids[END_REF] have presented their data. The arrows indicate the static overlap function from Fig. 2C of our article [START_REF] Ozawa | Equilibrium phase diagram of a randomly pinned glassformer[END_REF], that is, the expected long time limit of the collective function, showing the consistency of statics and dynamics. These data also show that with increasing c the self and collective functions decouple from each other and that the relaxation time of the latter increases rapidly, and thus shows a behavior that is qualitatively consistent with the presence of a nonergodic transition.

Regarding the accuracy of our configurational entropy, we point out that in our paper [START_REF] Ozawa | Equilibrium phase diagram of a randomly pinned glassformer[END_REF] we obtained the entropy by using a thermodynamic integration of the potential energy from the high temperature limit to temperatures inside the ideal glass state. In a completely independent calculation, we have estimated the vibrational entropy for the glass state via a harmonic approximation that should give accurate results at lower temperatures. Because inside the glass state the two approaches give the same result, we feel confident that the configurational entropy has been obtained with good accuracy.

Finally, we mention that the concentration for which Chakrabarty et al. [START_REF] Chakrabarty | Vanishing of configurational entropy may not imply an ideal glass transition in randomly pinned liquids[END_REF] present their results is close to or beyond the endpoint of the glass transition line (Fig. 4 in our paper [START_REF] Ozawa | Equilibrium phase diagram of a randomly pinned glassformer[END_REF]). Thus, unusual dynamics can be expected, and it will be interesting to study this dynamics in more detail. [START_REF] Ozawa | Equilibrium phase diagram of a randomly pinned glassformer[END_REF] and show that the static data are consistent with the dynamic data at long times.

FIG. 1 :

 1 FIG. 1: Time dependence of the self and collective overlap correlation functions (open and filled symbols, respectively) for pinning concentrations c approaching the ideal glass transition of a pinned system at T = 0.5. The horizontal arrows indicate the value of the static overlap function from Fig. 2C of our paper[START_REF] Ozawa | Equilibrium phase diagram of a randomly pinned glassformer[END_REF] and show that the static data are consistent with the dynamic data at long times.

Acknowledgments

We thank Chandan Dasgupta for helpful discussions. MO acknowledge the financial support by Grant-in-Aid for JSPS Fellows (26.1878). WK acknowledges the Institut Universitaire de France. AI acknowledges JSPS KAKENHI No. 26887021. KM and MO acknowledge KAKENHI No. 24340098, 25103005, 25000002, and the JSPS Core-to-Core Program. The simulations have been done in Research Center for Computational Science, Okazaki, Japan, at the HPC@LR, and the CINES (grant c2014097308).