## Abstract

Amorphous solids, such as glasses, have complex responses to deformations, with substantial consequences in material design and applications. In this respect, two intertwined aspects are important: stability and reversibility. It is crucial to understand, on the one hand, how a glass may become unstable due to increased plasticity under shear deformations, and, on the other hand, to what extent the response is reversible, meaning how much a system is able to recover the original configuration once the perturbation is released. Here, we focus on assemblies of hard spheres as the simplest model of amorphous solids such as colloidal glasses and granular matter. We prepare glass states quenched from equilibrium supercooled liquid states, which are obtained by using the swap Monte Carlo algorithm and correspond to a wide range of structural relaxation time scales. We exhaustively map out their stability and reversibility under volume and shear strains using extensive numerical simulations. The region on the volume-shear strain phase diagram where the original glass state remains solid is bounded by the shear yielding and the shear jamming lines that meet at a yielding-jamming crossover point. This solid phase can be further divided into two subphases: the stable glass phase, where the system deforms purely elastically and is totally reversible, and the marginal glass phase, where it experiences stochastic plastic deformations at mesoscopic scales and is partially irreversible. The details of the stability-reversibility map depend strongly on the quality of annealing of the glass. This study provides a unified framework for understanding elasticity, plasticity, yielding, and jamming in amorphous solids.

## INTRODUCTION

Understanding the response of amorphous materials to deformations is a central problem in condensed matter both from fundamental and practical viewpoints. It is not only a way to probe the nature of amorphous solids and their properties but also crucial to understand a wide range of phenomena from the fracture of metallic glasses to earthquakes and landslides. Furthermore, it has important applications in material design (*1*). Although many research efforts have focused on the mechanisms leading to the formation of amorphous solids from liquids (*2*–*5*), an orthogonal approach is to study these materials deep inside their amorphous phase (*6*–*8*). In this work, we focus on this second strategy by addressing the problem of understanding the nature of the response of glasses to volume and shear strains.

To a first approximation, glasses are solids much like crystals: They deform essentially elastically for small deformations but yield under large enough shear strains and start to flow. However, glasses are fundamentally different from crystals, being out-of-equilibrium states of matter. As a consequence, the properties of glasses strongly depend on the details of the preparation protocol (*3*). As an example, the yielding of glasses prepared via a fast quench or very slow annealing is qualitatively different (*9*). Thus, in sharp contrast to ordinary states of matter such as gases, liquids, and crystals, the equations of state (EOS), or the constitutive laws, of glasses, which characterize their macroscopic properties, must depend on the preparation protocol. Understanding the mechanical properties of glasses from a unified microscopic point of view thus emerges as a challenging problem (*10*).

To this aim, a central question is to understand the degree of stability of a glass, i.e., to what extent it can resist deformations. In isotropic materials such as glasses, it is sufficient to consider two types of deformations, namely, the volume strain, which changes the volume of the system isotropically, and the shear strain, which preserves the volume but changes the shape of the container. Under volume strains, glasses melt by decompression, and in the presence of a hard-core repulsion, as in granular matter and in colloids, they exhibit jamming upon compression. The melting and the jamming transitions delimit the line where the glass remains solid. Taking a glass on that line, one can probe its stability along the other axis of deformation, i.e., shear strain. Typically, the response of a glass to shear can be either (i) purely elastic and stable (note that this does not mean that the response is purely affine, as elasticity can emerge even in the presence of a nonaffine response), (ii) partially plastic (*6*–*8*), which is accompanied by slip avalanches and might be associated to the property of marginal stability (*4*, *11*), or (iii) purely plastic and unstable, once yielding takes place (*7*, *10*, *12*). Furthermore, granular materials (*13*) and dense suspensions (*14*) may (iv) jam when they are sheared.

A question related to stability is reversibility, i.e., to what extent a glass can recover its initial configuration when the deformation is released. This question has been one of the key interests in cyclic shear experiments of colloidal suspensions (*15*). In simulations of some model glasses, it has been found that a reversible-irreversible transition accompanies the occurrence of yielding (*16*–*18*).

The purpose of this work was to study, through extensive numerical simulations, the volume and shear strain phase diagram of a model glass former, hard spheres (HSs), to unify the abovementioned phenomena, i.e., plasticity, yielding, compression and shear jamming, and the failure of reversibility. Thanks to the swap algorithm, introduced by Kranendonk and Frenkel (*19*) and recently adapted to simulate polydisperse HS systems with unprecedented efficiency (*20*), we are able to prepare initial equilibrium supercooled liquid configurations up to high densities going even beyond experimental limits (*21*). While the standard molecular dynamics (MD) algorithm mimics the real dynamics, the swap algorithm accelerates the relaxation by introducing artificial exchanges of particles at different positions. With this trick, a dense supercooled liquid state with very large relaxation time can be prepared. Given such a system, by turning off the swap moves and switching to standard MD simulations, the system is effectively confined in a glass state, because its relaxation time is much larger than the achievable MD simulation time. By perturbing this initial equilibrium state with a given rate of compression, decompression, or shear strain during MD simulations, the system is driven out of the original equilibrium supercooled liquid state. In this way, we study the out-of-equilibrium response to these external perturbations of the glass selected by the initial supercooled liquid configuration, thus realizing what in (*22*) is called adiabatic state following. Using this procedure, we completely map out the degree of stability of the HS glasses corresponding to widely different preparation protocols. We show that there is a unique mapping between different types of stability and reversibility, that the stable and the marginally stable glass phases can be well separated by sensitive measurement protocols (*23*, *24*), and that marginality is manifested by a new type of reversibility, which we denote as partial irreversibility.

The idea of establishing a phase diagram to unify the glass transition, jamming, and yielding of amorphous solids was initially proposed by Liu and Nagel (*4*, *25*) and subsequently explored by others [e.g., see (*26*, *27*)]. Here, we explicitly construct such a phase diagram for HS glasses, represented by a stability-reversibility map, which complements the conjecture in (*4*) with new ingredients, namely, the existence of the marginal glass phase and the dependence on the quality of annealing (*22*, *28*, *29*). Our phase diagram is expected to be reproducible in experiments on vibrated granular glasses (*30*, *31*) and colloids (*10*), while molecular glasses are usually described by soft potentials, for which the phase diagram needs to be modified.

The plasticity of amorphous solids has been extensively studied in both phenomenological (*12*, *32*–*34*) and first principle (*5*, *29*) theories. According to the exact mean-field (MF) solution of the HS model in infinite dimensions (*29*), the glass phase can be decomposed into stable regions where plasticity is absent, and marginally stable regions where it is expected. The two phases are separated by a line where the so-called Gardner transition takes place (*5*, *23*, *24*). Determining whether this MF Gardner transition is also present in three dimensions is an extremely hard and currently open problem. Numerical simulations in three dimensions have found consistent evidence that an HS glass changes from a stable state to a marginally stable state across a certain threshold density before reaching jamming (*23*, *24*) but are not capable to determine whether such a change corresponds to a phase transition or a crossover, because of the lack of a careful analysis of finite size effects. Here, we relate the signatures of the Gardner transition/crossover to the emergence of plastic behavior and avalanches (*11*, *35*, *36*), which can be measured in simulations via the onset of partial plasticity and the emergence of a protocol-dependent shear modulus (*24*, *37*). The Gardner threshold determined in this approach is consistent with an independent estimate based on the growth of a spin glass–like susceptibility (*23*). Because the scope of our work is not to decide on the existence of a sharp Gardner phase transition, here we keep the conventional use of the terminology “Gardner transition” but do not exclude the possibility that it may become a crossover in three dimensions. Moreover, it remains an open question whether the Gardner transition and the associated marginality are of relevance to other systems. For example, the absence of marginality has been reported in simulations of a three-dimensional soft potential model (*38*) and a system of HSs confined in a one-dimensional channel (*39*). While details may change among various systems, the approach used in this study provides an example of how to construct a stability-reversibility map for generic glasses.

## RESULTS

### Preparation of annealed glasses

We study a three-dimensional HS glass with continuous polydispersity, identical to the one in (*20*) (see Materials and Methods). Note that for HSs, the temperature is irrelevant: It only fixes the overall kinetic energy of the system, which is related to the sphere velocities, and thus to the unit of time. In our simulations, we set *k*_{B}*T* to unity. The relevant control parameters in this study are the packing fraction φ and the shear strain γ. The reduced or dimensionless pressure *p* = *P*/(*k*_{B}*T*ρ), with *P* being the pressure and ρ the number density, can be determined uniquely from the EOS for the given φ and γ. Because the jamming limit is the point where the reduced pressure of HSs diverges, it corresponds, for our system, to the infinite pressure limit for fixed temperature or the zero-temperature limit for fixed pressure.

One can consider HSs as the limit of soft repulsive particles when the interaction energy scale divided by *k*_{B}*T* goes to infinity: Then, the HS system formally corresponds to the zero-temperature limit of soft repulsive particles in the unjammed phase where particles do not overlap. The jamming limit coincides in both systems, but the overjammed phase is inaccessible by definition for HSs. As a consequence, one of the axes (the temperature axis) in the Liu-Nagel phase diagram (*4*) will be missing in our context. The HS phase diagram established here should correspond to the zero-temperature plane of the Liu-Nagel phase diagram without the overjammed part.

Our HS model is chosen in such a way that the particle swap moves (*19*) can be used in combination with standard event-driven MD to fully equilibrate the system up to very high densities, covering a very wide range of time scales for the standard MD without swap (*20*). Switching off the swap movements at volume fraction φ_{g} and leaving only MD acting on the particles, one gets effectively an HS amorphous solid, corresponding to the glass that would be formed during an annealing process that falls out of equilibrium at φ_{g}. Therefore, φ_{g} is the glass transition density. Because the system is still in equilibrium at φ_{g}, its reduced pressure *p*_{g} follows the liquid EOS (L-EOS) *p*_{g} = *p*_{liq}(φ_{g}).

The possibility to explore a wide range of glass transition densities, thanks to the swap algorithm, is crucial to our work. In the following, we choose to work on three different values of φ_{g}, representing ascending levels of annealing:

1) Weakly annealed case: φ_{g} = 0.609, corresponding to the pressure *p*_{g} = 25.9. Berthier *et al*. (*21*) fitted the data of α-relaxation time τ_{α} as a function of *p* in liquids using the standard Vogel-Fulcher-Tammann (VFT) form τ_{α} = τ_{∞} exp[*A*/(*p*_{vft} − *p*)], a generalized VFT form τ_{α} = τ_{∞} exp[*A*/(*p*_{vft} − *p*)^{2}], and the facilitation model (FM) form τ_{α} = τ_{∞} exp[*A*(*p* − *p*_{fm})^{2}] [see (*21*) for details on the fitting]. We estimate that the α-relaxation time corresponding to *p*_{g} = 25.9 is about τ_{α}/τ_{0} ~ 5 × 10^{4} for all these forms, where τ_{0} ~ 10^{4} is the α-relaxation time at the onset density φ_{0} ≈ 0.56 of glassy dynamics. Both VFT and FM forms give consistent values of τ_{α}. The time scale τ_{α}/τ_{0} ~ 5 × 10^{4} corresponds to a typical time scale measured experimentally in colloidal glasses (τ_{α}/τ_{0} ≲ 10^{5} s and τ_{0} ≈ 10^{−1} s).

2) Moderately annealed case: φ_{g} = 0.631 and *p*_{g} = 30.9. At this density, the standard VFT fitting gives an estimated time scale τ_{α}/τ_{0} ~ 3 × 10^{10}, the generalized VFT gives τ_{α}/τ_{0} ~ 10^{10}, and the FM fitting gives τ_{α}/τ_{0} ~ 10^{9}. These time scales are typically reachable in molecular glass-forming liquids (τ_{α}/τ_{0} ≲ 10^{13} s and τ_{0} ≈ 10^{−10} s).

3) Deeply annealed case: φ_{g} = 0.655 and *p*_{g} = 40.0. At this density, the relaxation time is enormously large, and both VFT and FM fittings are unreliable. Fullerton and Berthier (*40*) measured the stability ratio S (the ratio between the melting time and the equilibrium relaxation time at the melting temperature) of this system. According to the data in (*40*), the stability ratio at this density is around (the value depends on the melting pressure), which is comparable to experimental scales of vapor-deposited ultrastable glasses (*41*).

While the time scales we can access correspond to different materials, as discussed above, it is important to stress that molecular glass-forming liquids and ultrastable glasses do not display a hard-core repulsion. The repulsion between molecules in these systems is usually better described by a Lennard-Jones–like soft potential. Therefore, some of the phenomena that we will describe in the following, which are strongly related to the presence of a hard-core potential, will be absent in these materials. The most important example is jamming, which is, by definition, not present in Lennard-Jones–like soft potentials. The nature of the Gardner transition could also be markedly different in some soft materials (*38*), and the applicability of some of our results on partial irreversibility should then be checked. Yet, we believe that the HS model is a remarkable benchmark as it displays many important instability mechanisms (melting, yielding, compression and shear jamming, and the onset of marginal stability). It thus allows us to study in full detail the interplay between these instability mechanisms and their dependence on the quality of annealing.

### Stability and reversibility

Starting from the equilibrated supercooled liquid configurations at φ_{g}, we now turn off the swap moves. By doing this, the liquid relaxation time goes beyond the time scale that we can access in our numerical experiments, and the system is thus effectively trapped into a glass state. We can then follow the quasi-static evolution of the system under slow changes of the volume strain ϵ = (φ_{g} − φ)/φ and the shear strain γ (see Materials and Methods) and measure the corresponding evolution of the pressure and the shear stress. Although the system is formally out of equilibrium (from the liquid point of view), one can reach a perfectly stationary state on the time scale we explore, restricted to the glass basin (*22*). The basin can then be followed in restricted metastable “equilibrium.” We call the resulting trajectory in control parameter space metastable EOS or glass EOS (G-EOS) to distinguish it from the L-EOS. The G-EOS can be obtained by plotting the pressure and stress as functions of the volume and shear strains.

The change of volume strain ϵ can be converted to that of volume fraction φ via the relation φ = φ_{g}/(1 + ϵ). To achieve a change in volume fraction, all particle diameters *D* are uniformly changed with rate for compression and for decompression. The resulting rate of the change of volume strain is . The change of the shear strain is given at a rate . The corresponding time scales of these rates are in between the fast β- and the α-relaxation times, in such a way that the glass is followed nearly adiabatically, while the α-relaxation remains effectively frozen (*23*, *24*). The target strains (ϵ, γ) can be achieved starting from the initial point (0, 0) following various paths in the volume-shear strain plane. For example, one can apply first a shear strain followed by a volume strain or vice versa. In the following, we specify explicitly the paths that we follow and check the dependency of the final outcome on the choices of paths.

Now, let us start our analysis by considering what happens following a simple cyclic deformation: First, the system is strained normally (0, 0) → (ϵ, 0), then sheared (ϵ, 0) → (ϵ, γ), and lastly sheared back in the reversed way (ϵ, γ) → (ϵ, 0). The following three typical behaviors are found: The response of the glass can be reversible, partially irreversible, or totally irreversible, which signals stable, marginally stable, and unstable states of the glass. Typical examples of the stress-strain curves are shown in Fig. 1.

(i) Reversible regime: For small γ, the stress σ increases smoothly and monotonically with increasing γ (green lines in Fig. 1, A and D). To the first order, the stress is linear, δσ = μδγ, where μ is the shear modulus. If the strain is released with , the stress-strain curve reverses to the origin—this is a typical elastic response.

(ii) Partially irreversible regime: For larger γ above a certain threshold γ_{G} , the stress-strain curve becomes jerky, consisting of piecewise linear elastic responses followed by small but abrupt stress drops (red lines in Fig. 1, A and D). Each stress drop corresponds to a plastic event, where some particles rearrange their positions. The glass in this regime is marginally stable in the sense that a tiny δγ could make the system unstable by triggering these plastic events, but the particles immediately find another stable configuration nearby, avoiding further failure of the entire system. Although the stress-strain curve is locally irreversible for small reversed strain, globally, it eventually returns to the origin when the shear strain is released back to γ = 0 (the red lines in Fig. 1, A and D, merge with the green lines below γ_{G}). We call this behavior partial irreversibility.

(iii) Limit of existence of the solid: For even larger γ, the glass faces two kinds of consequences depending on the volume strain ϵ applied before shearing.

• Yielding: At the yielding strain γ_{Y}, a sudden and notable stress drop occurs. When this happens, the entire system breaks into two pieces that can slide with respect to each other along a fracture. As shown by the stress-strain curve (black line in Fig. 1A), yielding is irreversible—once the glass is broken, it cannot be “repaired.” In a constant volume protocol where we keep the total volume of the system unchanged, yielding can be seen only if the system is not compressed to too high packing fractions, i.e., for not too negative volume strain ϵ.

• Shear jamming: The behavior changes markedly if the system is compressed more before shearing. In this case, the system jams at the shear jamming strain γ_{J}, which is signaled by the divergence of the shear stress (black line in Fig. 1D).

To examine the reversibility more carefully, we measure the relative mean squared displacement Δ_{r} (see Materials and Methods for the definition) between the initial state at γ = 0 before the shear is applied and the final state at γ = 0 after a single cycle of shear is applied (Fig. 1, B and E). If the initial and the final configurations are identical, Δ_{r} = 0; otherwise, the more different they are, the larger Δ_{r} is. The value of Δ_{r} returns to zero in the reversible and partially irreversible cases but becomes nonzero in the irreversible case, being consistent with the above analysis based on the stress-strain curves. Note that, here, we neglect differences on the microscopic scale of vibrational cage size Δ ≲ 0.01 (see Materials and Methods), i.e., a system is called irreversible only if the difference on Δ_{r} between the initial and final configurations is larger than Δ. We have also examined that the above behaviors persist in multicycle shears (see fig. S2).

It is useful to understand our observations using a schematic picture of the free-energy landscape. Each glass state is represented by a basin of free-energy *F*(φ_{g}; ϵ, γ; *N*), which is distorted upon increasing shear strain γ (Fig. 1, C and F). The shear stress is nothing but the slope of the free-energy , with β being the inverse temperature. The associated shear modulus is obtained by taking one more derivative with respect to γ, which gives nothing but the curvature of the free-energy basin. In the stable regime, the basin is smooth; in the marginally stable regime, the basin becomes rough, consisting of many sub-basins with larger associated shear modulus, which results in the failure of pure elasticity (*33*, *36*, *37*). In this state, the system can release the stress via hopping between different sub-basins, corresponding to plastic events, which leads to emergent slow relaxation of shear stress (*24*, *37*). For very large strains, the system either yields by escaping from the glass basin (Fig. 1C) or jams by hitting the vertical wall due to the hard-core constraint (Fig. 1F).

The plastic behavior appearing in the partially irreversible regime takes place at mesoscopic scales, and it would be averaged out in a macroscopic system at large enough time scales (*37*). There is evidence which shows that the minimum strain increment δγ_{trigger}(*N*) to trigger a plastic event vanishes in the thermodynamic limit *N* → ∞ (*33*, *42*). This implies that in a macroscopic system, any small but finite increment of strain would cause a nonzero number of mesoscopic plastic events (*11*). Moreover, time-dependent aging effects associated to these plastic events were observed in stress relaxations (*24*). Therefore, in macroscopic systems at large enough time scales, the plasticity would be averaged out, and one would observe just a renormalized “elastic” response. The bare elastic response can only be seen within the piecewise linear mesoscopic response for δγ < δγ_{trigger}(*N*). This means that two different shear moduli can be defined: the bare one μ_{bare} = lim_{N → ∞} lim_{δγ → 0}δσ(φ_{g}; ϵ, γ; *N*)/δγ that takes into account the piecewise elastic behavior between two subsequent avalanches, and the macroscopic one μ_{macro} = lim_{δγ → 0}lim_{N → ∞}δσ(φ_{g}; ϵ, γ; *N*)/δγ, which represent the average behavior and is smaller than the former (*33*, *37*). Therefore, the small strain δγ → 0 limit and the thermodynamics limit do not commute in the marginal plastic phase (see text S1 for a detailed discussion).

### Stability-reversibility map and G-EOS

These three different kinds of responses of the system to simple cyclic shear, listed above as (i) to (iii), can be summarized by the stability-reversibility map in the ϵ − γ plane, as shown in Fig. 2. There we also show a typical plastic event in the marginal phase (Fig. 2B) and a yielding event (Fig. 2C), which indicate two different mechanisms that can cause a failure of stability. As long as the glass remains stable or marginally stable, its macroscopic properties can be characterized by the G-EOS for the pressure *p* = *p*_{glass}(φ_{g}; ϵ, γ) and the shear stress σ = σ_{glass}(φ_{g}; ϵ, γ), as shown in Fig. 3 (A and B). The pressure *p* and the shear stress σ are derivatives of the glass free-energy − β*F*(φ_{g}; ϵ, γ) with respect to ϵ and γ, respectively.

Along the γ = 0 line, the evolution of the system under volume strain ϵ will eventually lead the system to either jamming after sufficient compression ϵ < 0 or melting after sufficient decompression ϵ > 0. At jamming, particles form an isostatic rigid contact network such that no further compression can be applied. Decompressing the system reduces *p*, which eventually melts the system into a liquid state. The evolution of the pressure *p* follows the zero–shear strain G-EOS *p* = *p*_{glass}(φ_{g}; ϵ, 0) both upon compression and decompression. Obviously, σ = σ_{glass}(φ_{g}; ϵ, 0) = 0.

Applying a shear strain at any point on the γ = 0 G-EOSs *p* = *p*_{glass}(φ_{g}; ϵ, 0) and σ = σ_{glass}(φ_{g}; ϵ, 0) allows us to explore the volume strain versus shear strain phase diagram, and we can track both the pressure *p* = *p*_{glass}(φ_{g}; ϵ, γ) and the stress σ = σ_{glass}(φ_{g}; ϵ, γ). Under shear, the glass has two possible fates: either it yields across the shear yielding line γ = γ_{Y}(φ_{g}; ϵ) or it jams at the shear jamming line γ = γ_{J}(φ_{g}; ϵ). Yielding can be detected by analyzing the stress-strain curve, i.e., σ_{glass}(φ_{g}; ϵ, γ) versus γ, while shear jamming is signaled by a divergence of both the pressure *p*_{glass}(φ_{g}; ϵ, γ) → ∞ and the stress σ_{glass}(φ_{g}; ϵ, γ) → ∞ as γ → γ_{J}(φ_{g}; ϵ). The shear yielding and the shear jamming lines define the boundaries of the stability of the HS glass, beyond which the glass is unstable or simply forbidden. The two lines meet at a yielding-jamming crossover point [ϵ_{c}(φ_{g}), γ_{c}(φ_{g})] or [φ_{c}(φ_{g}), γ_{c}(φ_{g})].

Within the boundary of the stability-reversibility map, there are two phases: the stable (reversible) phase and the marginally stable (partially irreversible) phase. We call the line that separates the two a Gardner line. Across this line, the qualitative nature of the system’s response to deformations changes: The stress-strain curve is smooth within the stable (reversible) phase but jerky in the marginally stable (partially irreversible) phase. The stability-reversibility map shown in Fig. 2 suggests that if we choose an ϵ such that the Gardner line is not crossed along the path (ϵ, 0) → (ϵ, γ_{Y}), then no marginally stable region should be observed. Figure S1 shows such a case (with ϵ = 0.057) where we do not observe partial irreversibility all the way up to yielding. The term Gardner line is inferred from the MF glass theory (*5*, *29*), in which a continuous phase transition, the Gardner transition, occurs on this line. However, whether it is a genuine transition line or a crossover line in three dimensions is an open question, as we noted in the Introduction. In the next subsection, we will explain how we estimate this line numerically in the present system.

We made the choice in Fig. 2 to represent the stability-reversibility map in terms of strains (volume and shear). In fig. S3A, we plot it in terms of volume fraction φ and shear strain γ, which can be directly compared with the theoretical prediction in (*29*). In some experiments, the shear stress is controlled instead of the shear strain, and in that case, it is customary to represent the phase diagram in the density-stress plane. Such a figure is reported in fig. S3B, which is directly comparable to the phase diagram reported in the granular experiment of (*30*).

The stability-reversibility map and the G-EOS depend on the preparation density φ_{g} of the glass, which represents the depth of annealing. As shown in Fig. 3C, where the G-EOS and L-EOS are displayed together, the γ = 0 G-EOS and the L-EOS intersect at the point (φ, γ) = (φ_{g}, 0), which shows the intrinsic connection between G-EOS and L-EOS. The initial unperturbed glass is located at (ϵ, γ) = (0, 0) in the stability-reversibility map.

### Marginal stability and partial irreversibility

Having presented above our most important results, in the following, we show more details on how the stability-reversibility map and the G-EOS are obtained in our numerical experiments. To this end, at each φ_{g}, we prepare *N*_{s} = 100 independent equilibrium supercooled liquid configurations by the swap algorithm, which have different equilibrium positions of particles and are called samples. By switching off the swap, they become glasses. For each sample of glass, we repeat *N*_{r} ~ 50 − 200 realizations of a given protocol, which is a combination of compression (or decompression) and simple shear. Each realization starts from statistically independent initial particle velocities drawn from the Maxwell-Boltzmann distribution at φ_{g}.

The Gardner transition marks the point where the elastic behavior is replaced by a partially plastic one. Avalanches and plasticity are extremely marked in finite size systems, while they are averaged out on macroscopic length and time scales. Furthermore, in finite size systems, even though each individual stress-strain curve is jerky in the marginal phase as shown in Fig. 1, the average over different samples and realizations washes out all the sudden drops, giving rise to a smooth profile. Therefore, macroscopic G-EOSs by themselves do not allow the detection of the marginally stable phase (see text S1 and fig. S4 for a detailed discussion). To precisely locate the onset of plasticity and the marginal phase, we will examine the hysteretic response to very small shear increments.

Inspired from spin glass experiments (*43*), we compare the shear stress measured by two different protocols, the so-called zero-field compression (ZFC) and the field compression (FC) protocols (*24*). Within the FC protocol, one first compresses the system and then shears it. In the ZFC, one instead reverses the order (see Materials and Methods for more details). The FC stress σ_{FC} can be considered as the large time limit of σ_{ZFC}, as long as the yielding and the α-relaxation do not occur (*24*, *37*). For elastic solids such as crystals, the two stresses are identical. For marginally stable glasses, however, σ_{FC} is lower than σ_{ZFC}, because of the stress relaxation associated to the plastic events happening at mesoscopic scales. The origin of two responses can be attributed to the organization of free-energy landscape shown schematically in Fig. 1 (C and F). Roughly speaking, the ZFC stress σ_{ZFC} is dominated by the short time response within the small sub-basins, while the FC stress σ_{FC} reflects the renormalized, long time response within the big envelope of sub-basins. The bifurcation point between the two stresses determines the Gardner point. Note that this criterion to determine the Gardner point is the same as the one used in (*24*). Figure 4A shows the data used to obtain the Gardner points ϵ_{G}(φ_{g}; γ) for a few different values of γ. Connecting the Gardner points gives the Gardner line γ = γ_{G}(φ_{g}; ϵ) in Fig. 2. See fig. S5 for the same results obtained for other values of φ_{g}.

Alternatively, one may look at caging order parameters such as the mean squared displacement Δ and the typical separation between two replicas Δ_{AB} (*23*) (see Materials and Methods for more precise definitions). The two replicas are generated from the same initial sample in two independent realizations. They are first compressed to a target ϵ under zero shear strain and then sheared to the target shear strain γ under the fixed ϵ. When the Gardner point γ_{G}(φ_{g}; ϵ) is crossed over, Δ and Δ_{AB} should also separate. However, this is a sign of critical behavior only if the corresponding susceptibility grows (*39*). Here, 〈 … 〉 represents the average over both samples and realizations. χ_{AB} is a spin glass–like susceptibility whose growth suggests the increase of heterogeneity and cooperativity in the system, as suggested by the MF theory (*5*). The behavior of χ_{AB} can be inferred from Fig. 1B, where we plot the probability distribution *P*(Δ_{AB}) of Δ_{AB}. We see a Gaussian-like behavior below the Gardner threshold, fat tailed around it, and double peaked above it. The Gardner point inferred in this way γ_{G}(φ_{g} = 0.655; ϵ = −0.036) ≈ 0.06 is consistent with the determination from ZFC-FC protocols ϵ_{G}(φ_{g} = 0.655; γ = 0.06) ≈ −0.036.

This result provides a strong evidence that partial irreversibility and plasticity in Fig. 1 are essentially related to emerging marginal stability. We perform the following test to examine their connections more directly. Starting from a compressed glass at (ϵ, γ = 0), we first shear the glass to a target shear strain at (ϵ, γ − δγ) under constant volume and then apply an additional cycle of small shear strain δγ = 0.004, following the path (ϵ, γ − δγ) → (ϵ, γ) → (ϵ, γ − δγ). If the system is reversible, then the difference between the stresses before and after the single cyclic shear, , should be zero, otherwise not. Figure 4C confirms that begins to grow around the γ_{G} estimated from the other two approaches described before (Fig. 4, A and B). However, this kind of irreversibility is only partial because the system is reversible under a circle of shear with larger strain. Systems following the path (ϵ, δγ) → (ϵ, γ) → (ϵ, δγ), where δγ = 0.004 is fixed and γ is varying, show that the stress difference is nearly zero for any γ.

Last, it is important to stress that in our three-dimensional numerical simulations, as in previous ones (*23*, *24*), we cannot decide on whether the separation between the stable and marginally stable phase corresponds to a true phase transition. This would require, for instance, a careful study of finite size effects on χ_{AB} to extract the behavior for *N* → ∞, which is very difficult already in much simpler models such as spin glasses. The focus of our work is on relating the Gardner line, which is only a (quite sharp) crossover in our simulations, to the onset of partial irreversibility.

### Shear yielding and shear jamming

Up to now we have investigated the interior of the stability-reversibility map. Next, we turn to explore the boundaries of the stability-reversibility map by analyzing the G-EOSs both in pressure and shear stress. From now on, all data presented are averaged over different samples and realizations. Therefore, even with a finite size system, the individual plastic events will be averaged out. Furthermore, we will plot the G-EOSs on a phase diagram using 1/φ (instead of, equivalently, ϵ) and γ to better show their relations to L-EOS.

First of all, starting from an equilibrium configuration at (φ, γ) = (φ_{g}, 0) or (ϵ, γ) = (0, 0), the system melts under decompression for sufficiently large ϵ. We define the melting point as the crossover point between the γ = 0 G-EOS for the pressure and the L-EOS (see Figs. 3C and 5G). The melting point sets the upper bound of the stability-reversibility map along the γ = 0 line.

To systematically explore the stability-reversibility map, we design three specific protocols combining compression/decompression and shear, namely, CP-S, CV-S, and CS-C/D (see Materials and Methods for details). These protocols can be realized also in experiments. In principle, the EOS should not be protocol dependent, but whether it is also the case for G-EOS is not so obvious.

In the CP-S protocol, for any fixed pressure *p*, the specific volume 1/φ (or volume strain ϵ) evolves with shear strain γ, which defines a G-EOS for the pressure. Figure 5A shows the G-EOSs for a few different pressures *p* in a γ − 1/φ plot. Such a plot is essentially the projection of the three-dimensional plot of the G-EOSs for the pressure *p* = *p*(φ_{g}; ϵ, γ) in Fig. 3C onto the γ − 1/φ plane. The data show that the specific volume 1/φ expands as strain γ is increased, known as the dilatancy effect. The dilatancy is stronger for better annealed glasses, as observed previously in (*24*), and at lower pressure for a fixed quality of annealing, as shown here. Both observations are consistent with theoretical predictions [ Fig. 2 in (*22*) and Fig. 2A in (*29*), respectively]. Note that this dilatancy effect shall be distinguished from the one discussed in the context of steady flow, which is necessarily correlated to friction, as shown in (*44*). At high pressures, the isobaric lines are nearly parallel to the shear jamming line, which corresponds to the *p* = ∞ isobaric line by the definition of jamming. On the other hand, the average stress σ shown in Fig. 5B initially increases with the shear strain, but it eventually approaches a plateau after a big drop corresponding to yielding. We define the yielding point as the peak of the stress susceptibility χ_{σ} = *N*(〈σ^{2}〉 − 〈σ〉^{2}) (see Fig. 5C). The yielding point is approximately at the middle of the drop on the stress-strain curve, corresponding to the steepest decrease of stress. After yielding, the shear stress generally remains nonzero, indicating that the glass is not completely fluidized. Real-space visualization shows that the glass breaks into two pieces sliding against each other (see Fig. 2C). However, near the melting point, such a picture might change, because melting could mix with yielding, giving rise to a hybrid behavior. We will not discuss this situation in detail here. By connecting the yielding strains obtained at different *p*, we obtain the yielding line. We notice that for a certain range of pressure *p* near *p*_{g} of the initial glass, the yielding strain is nearly independent of *p*.

In Fig. 5 (D and E), we show how the inverse pressure 1/*p* and shear stress σ evolve with γ for various φ in the CV-S protocol. We find a threshold density φ_{c} (see fig. S6 for how φ_{c} is determined), which separates the shear yielding and shear jamming cases. If φ < φ_{c}, the system generally yields at large γ; otherwise, both pressure and shear stress diverge as γ is increased, indicating shear jamming. In this protocol, the yielding point can be determined again from the peak of the stress susceptibility (see Fig. 5F). In the shear jamming case, the pressure and shear stress both follow the free-volume scaling laws: *p* ~ (γ_{J} − γ)^{−1} and σ ~ (γ_{J} − γ)^{−1} (see fig. S7). The shear jamming is a natural consequence of the dilatancy effect (i.e., *p* increases with γ for fixed φ), as long as the system does not yield. Thus, φ_{c} results from the competition between the dilatancy effect and the tendency to break the system at large strains. We have checked that all the shear-jammed packings that we create satisfy the isostatic condition (*45*), i.e., the average coordination number *Z* = 6, once the ratters (particles that have less than four contacts) are excluded, and that the shear jamming transition falls in the same universality class of the usual jamming transition in the absence of shear.

Figure 5 (G and H) shows the constant-γ EOSs of the pressure and shear stress for a few different γ in the CS-C/D protocol. For small shear strains, , the system jams at a γ-dependent jamming density φ_{J} under compression. For shear strains larger than the yielding strain , however, the G-EOSs for pressure collapse onto the same curve, and consequently, the jamming density φ_{J} also does not change with γ anymore. This observation is consistent with our interpretation of yielded states: The glass just breaks into two pieces of solids at by forming a planar fracture. These planar structures should have minor effect on bulk properties like the pressure. On the other hand, the glass always melts under decompression, for any γ. We find that the melting point is independent of γ, both below and above . The stress susceptibility χ_{σ} displays a peak upon decompression, which reveals the vestige of yielding, and therefore can be used to define the yielding point in the CS-C/D protocol (Fig. 5I). For , the yielding density φ increases with γ; for , the peak does not exist anymore, and the yielding point cannot be defined as expected. In addition, we show and discuss the behavior of the pressure susceptibility χ_{p} in fig. S8.

### Dependence on protocols and system sizes

Let us discuss how the stability-reversibility map and G-EOSs depend on protocols. There are two important sources of protocol dependences. First, the stability-reversibility map and the G-EOSs depend on the glass transition point φ_{g}, and φ_{g} itself depends on the protocol parameters such as the compression rate in a standard compression annealing protocol (here, it is a function of where we stop swap moves). Figure 6A shows the stability-reversibility maps for three different φ_{g}, corresponding to three typical experimental time scales as discussed previously (see fig. S9A for the three-dimensional representations). They share common qualitative features in general. The stable regime expands with φ_{g}, as one would naturally expect that more deeply annealed glasses should be more stable. The shear jamming line becomes more vertical with decreasing φ_{g}. This trend is consistent with previous numerical observations, which show that, in the thermodynamical limit, the shear jamming line is completely vertical for infinitely rapidly quenched systems (*46*). Moreover, we point out that the Gardner transition points cannot be determined unambiguously using our approaches for the less annealed systems φ_{g} = 0.631, γ > 0.06 and φ_{g} = 0.609, γ > 0 (see fig. S5) because different activated dynamics, such as plastic rearrangements, formation of fractured structures, and α-relaxations cannot be well separated.

Second, we show in Fig. 6B how the stability-reversibility map and also the G-EOSs for the pressure depend on the exploration protocols (CP-S, CV-S, and CS-C/D); see fig. S9 (B and C) for the three-dimensional representations. We find a protocol-independent regime , where all three protocols give the same pressure. The part of the stability-reversibility map above cannot be accessed by the CS-C/D protocol. For ϵ > ϵ*, the yielding line bends down differently, depending on the protocol. The system yields most easily in the CV-S protocol, presumably because the liquid bubbles formed around melting are easier to expand in a volume-controlled protocol (*40*).

Last, we discuss briefly how the stability-reversibility maps depend on the system size *N* in fig. S10. We do not observe appreciable finite size effects on the shear jamming line. On the other hand, the yielding line exhibits strong finite size effects, but we expect it to converge at larger sizes, based on the recent results of (*9*). Using the present methods, we also do not find strong *N* dependence on the Gardner line, consistent with the data in (*24*). However, we stress that based on available numerical results, we cannot conclude on the thermodynamic behavior of the Gardner transition. Understanding whether it is a sharp transition or a crossover is an active and hot topic in the field, through numerical (*23*), experimental (*31*), and theoretical analyses (*39*, *47*–*49*). While the finite size analysis presented here shall not be considered conclusive, we leave a more detailed finite size study on yielding, shear jamming, and the Gardner transition for future work.

## DISCUSSION

Here, we investigate the stability and the reversibility of polydisperse HS glasses under volume and shear strains. We prepare equilibrium supercooled liquid states, with different degrees of stability ranging from a fast quench to an extremely slow annealing, corresponding to ultrastable configurations. Each configuration corresponds to a glass within a time scale that is shorter than the structural relaxation time. We study the stability of the glass under volume and shear strains and find that the region of stability is delimited by lines where the system can either yield or jam. We also find that within the region of stability, the system can be either a normal solid, which essentially responds elastically and reversibly to perturbations, or a marginally stable solid, which responds plastically and in a partially irreversible way. More precisely, the main outcomes of our analysis are the following:

1) Response. The response of the system to a shear strain is either purely elastic, partially plastic, or fully plastic (yielding), depending on the quality of annealing and the amount of volume and shear strains imposed to it.

2) Failure. Well-annealed glasses (large φ_{g}), when sheared at sufficiently low densities (large volume strain ϵ), behave purely elastically up to yielding, which is an abrupt process where a fracture is formed and the glass fails. At higher densities, they display a partial plastic phase before yielding is reached. At even higher densities, they display shear jamming (under constant volume shear). The shear yielding and shear jamming lines delimit the region of existence of the HS amorphous solid.

3) Marginality. Along the solid part of the stress-strain curves, the partial plastic behavior is well separated from the purely elastic one by the Gardner point. The onset of partial plasticity is accompanied by the emergence of critical behavior and marginal stability. Beyond the Gardner point, the shear modulus of the system becomes history dependent. At the same time, a growing spin glass–like susceptibility is observed.

4) Reversibility. The purely elastic phase is globally reversible: Once the shear is released, the system gets back to the original configuration. The partially plastic marginal glass phase is partially irreversible: Upon releasing the deformation by a small amount, the system is not able to get back to the previous state, while upon complete release, the system is able to get back to the original configuration. Yielding corresponds to complete irreversibility: Once broken, the system starts to flow, and it is not able to get back to the original configuration once the strain is completely released.

By collecting together the boundaries of the different regions, we obtained a complete stability-reversibility map (phase diagram), reported in Fig. 2. The stability-reversibility map obtained in the present study for three-dimensional HS glasses can be compared with the one obtained by the MF theory in the large dimensional limit (*29*). The most important features, such as the presence of the shear jamming and the shear yielding lines that delimit the stability region and the presence of the Gardner line, are qualitatively in good agreement with the predictions of the theory. There are, however, several important differences. (i) The shear yielding line in the three-dimensional system is not a spinodal line, as predicted by the MF theory (*22*). The abrupt formation of a fracture is completely missed by the MF theory, which does not describe the spatial fluctuations of stress that accumulate around the fracture. (ii) The point (ϵ_{c}, γ_{c}) where the shear yielding and the shear jamming lines meet is predicted to be a critical point in the MF theory, but it is rather a crossover point in the three-dimensional system. (iii) The marginally stable phase has larger γ_{Y} than the stable phase. This suggests that the plastic events in the marginal phase help the system to avoid total failure. In the theory, the shear yielding line bends down rather than bends up near the point (ϵ_{c}, γ_{c}) [see Fig. 3 of (*29*)].

Note that the MF predictions of (*29*) were obtained using the so-called replica symmetric (RS) ansatz. To properly consider yielding in the marginally stable phase, one should extend the computation to a full-step RS breaking (fullRSB) ansatz (*50*). This might help in solving some of the discrepancies between the analytical and the numerical results. According to the RS theory, yielding is a spinodal transition with disorder (*22*). However, it is not clear how this picture will be modified by a fullRSB theory.

Our simulation results show that a well-annealed glass (φ_{g} well above the MCT density φ_{MCT}) yields abruptly—it is brittle. However, a poorly annealed glass (φ_{g} ~ φ_{MCT}) may instead continuously yield into a plastic flow state (*10*, *12*)—it is more ductile. We expect that near the melting point, even a well-annealed glass would behave similarly to a poorly annealed one, as it would become much “softer” upon decompression. Nevertheless, the yielding point can be determined from the peak of χ_{σ} for both cases as shown here. Our approach thus provides a unified framework to study the transition between the two distinct mechanisms of yielding. The possibility of two yielding mechanisms is missed by the current MF theory. A dynamical extension of MF might account for these effects. Understanding the nature of the yielding transition (*51*–*53*) is a crucial problem that requires further analysis.

The plastic events we observe in the partially irreversible phase could correspond to two different types of soft modes: collective modes, associated to a diverging length scale, as predicted by the MF theory in the marginally stable phase (*5*, *11*), or localized modes, such as the ones that have been observed in numerical studies of low-dimensional systems (*54*–*56*). In this study, we did not investigate systematically the nature of the plastic events in our system, but the growth of the spin glass–like susceptibility in our data suggests the presence, in our HS model, of large-scale collective excitations. Note that the situation could be radically different in soft potential models (*38*, *55*, *56*). We also stress that while the existence of partial plasticity before yielding is well known (*6*–*8*), our well-annealed systems provide an example where the pure elasticity and partial plasticity regimes are well separated, allowing us to define a line (the Gardner line) that separates them in the stability-reversibility map.

Last, concerning the reversibility, here we focus on the reversibility with respect to just one cycle of simple shear (see fig. S2 for the results under a few cycles). In cyclic shear protocols, a steady state can be reached after many cycles (*17*, *57*). Very complicated dynamics should be involved in these processes. It would be interesting to systematically extend the present study to multiple cyclic shear to understand better these processes.

## MATERIALS AND METHODS

### Model

The system consists of *N* = 1000 (unless otherwise specified) HS particles with a diameter distribution *P*(*D*) ~ *D*^{−3}, for *D*_{min} ≤ *D* ≤ *D*_{min}/0.45. The continuous polydispersity was sufficient to suppress crystallization even in deep annealing and optimized the efficiency of swap algorithm. The volume fraction is , where ρ = *N*/*V* is the number density and *V* is the total volume. We define the reduced pressure *p* = *PV*/*Nk*_{B}*T* and the reduced stress σ = Σ*V*/*Nk*_{B}*T*, where *P* and Σ are the pressure and the stress of the system, respectively. For simplicity, in the rest of this paper, we refer as pressure and stress to *p* and σ instead of *P* and Σ. We set the Boltzmann constant *k*_{B}, the temperature *T*, the mean diameter , and the particle mass *m* to unity.

### Swap algorithm

At each dynamical step, the swap Monte Carlo algorithm attempts to exchange the positions of two randomly picked particles as long as they do not overlap with their new neighbors. These nonlocal Monte Carlo moves eliminate the local confinement of particles in supercooled states, which, combined with standard event-driven MD, substantially facilitates the equilibration procedure. It has been carefully examined that the swap algorithm does not introduce crystalline order in the polydisperse HS model studied here (*20*).

### Compression/decompression algorithm

We used the Lubachevsky-Stillinger algorithm (*58*) to compress and decompress the system. The particles were simulated by using event-driven MD. The sphere diameters were increased/decreased with a constant rate. The MD time was expressed in units of .

### Simple shear algorithm

At each step, we performed *N*_{collision} = 100 − 1000 collisions per particle using the event-driven MD and then instantaneously increased the shear strain by , where δ*t* is the time elapsed during the collisions. The instantaneous shear shifts all particles by *x*_{i} → *x*_{i} + δγ*z*_{i}, where *x*_{i} and *z*_{i} are the *x* − and *z* −coordinates of particle *i*. To remove the possible overlappings introduced during this shift, we switched to a harmonic interparticle potential and used the conjugated gradient (CG) method to minimize the energy. The harmonic potential was switched off after CG. The Lees-Edwards boundary conditions (*59*) were used [see (*24*) for more details].

### Protocols of ZFC and FC

In the ZFC protocol, starting from the initial equilibrium configuration at (ϵ, γ) = (0, 0), we (i) sheared the system to a target shear strain at (0, γ) while keeping the volume strain unchanged, (ii) compressed it to a target volume strain at (ϵ, γ) while keeping the shear strain unchanged, (iii) applied an additional small shear strain δγ = 0.002, and (iv) measured the stress σ_{ZFC} at the state point (ϵ, γ + δγ). In the FC protocol, the order of steps (ii) and (iii) was interchanged. The FC protocol therefore has the path (0, 0) → (0, γ) → (0, γ + δγ) → (ϵ, γ + δγ). The target shear strain was chosen such that it is below the yielding strain γ < γ_{Y}. Here, the shear strain serves as an external "field" with respect to compression, in analogy to the magnetic field in cooling experiments on spin glasses (*43*). The stress was measured on a time scale *t* = 10 ≈ 10τ_{0}, where τ_{0} is the ballistic time. This choice ensures that the ZFC protocol measures the short time response to shear, while the FC measurement corresponds to the long time response because the shear strain γ + δγ was reached before the volume strain was applied [see (*24*) for a detailed analysis on the stress relaxation dynamics]. This protocol generalizes the one used in (*24*), which corresponds to the case γ = 0.

### Protocols of CP-S, CV-S, and CS-C/D

In the CP-S protocol, the system was first compressed or decompressed (depending on whether the target *p* is higher or lower than *p*_{g}) from the equilibrium state at (p, γ) = (*p*_{g}, 0) to the state at (*p*, 0). Then, simple shear was applied under the constant-*p* condition, until the system reached the target shear strain at (*p*, γ). At each shear step, the particle diameters were adjusted to keep *p* constant. In the CV-S protocol, the system was first compressed or decompressed from φ = φ_{g} to the target density φ, and then the simple shear was applied by keeping the volume constant. In the CS-C/D protocol, the system was first sheared from (φ, γ) = (φ_{g}, 0) to a target strain at (φ_{g}, γ), and then compression or decompression was applied while keeping the shear strain γ constant.

### Caging order parameters

We considered two order parameters Δ_{r} and Δ_{AB} defined below to characterize the glass state. The relative mean squared displacement is defined as(1)where {*r*_{i}} and are the particle coordinates of the target and reference configurations. In Fig. 1, the target and reference are the configurations after and before shear, respectively. The replica mean squared displacement(2)measures the distance between two replicas of the same sample generated by two independent realizations.

One may also consider the time-dependent mean squared displacement , whose value at the ballistic time scale τ_{0} ~ 1 gives the typical vibrational cage size of particles. We found that in our systems, Δ(τ_{0}) ≲ 0.01 [see (*23*)]. The cage size is nearly unchanged under simple shear.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/12/eaat6387/DC1

Fig. S1. Singe-realization stress-strain curve for φ_{g} = 0.655 and ϵ = 0.057.

Fig. S2. Multicycle stress-strain curves.

Fig. S3. Other representations of the stability-reversibility map.

Fig. S4. Rescaled stress-strain curves.

Fig. S5. Determination of the Gardner threshold for other φ_{g}.

Fig. S6. Determination of the yielding-jamming crossover point.

Fig. S7. Free-volume scalings in shear jamming.

Fig. S8. Pressure susceptibility in the CS-C/D protocol.

Fig. S9. Dependence of the stability-reversibility map on φ_{g} and protocols.

Fig. S10. Dependence of the stability-reversibility map on the system size.

Text S1. Bare and macro shear moduli.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank L. Berthier, M. Ozawa, C. Scalliet, M. Wyart, A. Altieri, O. Dauchot, K. Miyazaki, and T. Kawasaki for the discussions.

**Funding:**This work was supported by KAKENHI (no. 25103005 “Fluctuation & Structure” and no. 50335337) from MEXT, Japan, by the CAS Pioneer Hundred Talents Program (Y.J.), by a grant from the Simons Foundation (no. 454955, Francesco Zamponi), and by “Investissements d’Avenir” LabEx PALM (ANR-10-LABX-0039-PALM) (Pierfrancesco Urbani). The computations were performed using the computing facilities in the Research Center for Computational Science, Okazaki, Japan, and in the Cybermedia Center, Osaka University.

**Author contributions:**Y.J., P.U., F.Z., and H.Y. designed the research. Y.J. wrote the code and performed the numerical simulations and the data analysis, in close collaboration with H.Y. All authors contributed to the data analysis, the theoretical interpretation of the results, and the writing of the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).