Yoshihiko Nishikawa

Koji Hukushima

Lattice glass model in three spatial dimensions

The understanding of thermodynamic glass transition has been hindered by the lack of proper models beyond mean-field theories. Here, we propose a three-dimensional lattice glass model on a simple cubic lattice that exhibits the typical dynamics observed in fragile supercooled liquids such as two-step relaxation, super-Arrhenius growth in the relaxation time, and dynamical heterogeneity. Using advanced Monte Carlo methods, we compute the thermodynamic properties deep inside the glassy temperature regime, well below the onset temperature of the slow dynamics. The specific heat has a finite jump towards the thermodynamic limit with critical exponents close to those expected from the hyperscaling and the random first-order transition theory for the glass transition. We also study an effective free energy of glasses, the Franz-Parisi potential, as a function of the overlap between equilibrium and quenched configurations. The effective free energy indicates the existence of a first-order phase transition, consistent with the random first-order transition theory. These findings strongly suggest that the glassy dynamics of the model has its origin in thermodynamics.

A thermodynamic (or ideal) glass transition in finite low dimensions has been actively discussed both theoretically and experimentally for decades since the seminal works by Kauzmann [1], and Adam and Gibbs [2]. A mean-field theory, the random first-order transition (RFOT) theory of structural glasses, proposed and developed in Refs. [3][4][5][6] indeed shows that a thermodynamic glass transition at finite temperature exists with vanishing configuration entropy, or complexity [7]. In finite dimensions, numerical simulation of models of fragile supercooled liquids is a promising way to theoretically explore a possibility of the thermodynamic glass transition. However, the notoriously long relaxation time prevents us to access directly low-temperature thermodynamics of glass-forming supercooled liquids. Although recent progress on particle models simulated using the swap Monte Carlo method [8,9] allows us to get much more stable glass configurations at lower temperature or higher density, the thermodynamic glass transition is still inaccessible.

In an analogy to phase transitions into long-range ferromagnetic and crystal states, the lower critical dimension for the thermodynamic glass transition in models with discrete symmetry may be lower than that in models with continuous symmetry. Thus, it would be crucial to explore a possibility of a thermodynamic glass transition in finite-dimensional lattice models with discrete symmetry. Although several simple lattice models with the mean-field thermodynamic glass transition have already been proposed [10][11][12], they are not fully suitable to study the finite-dimensional glass transition: For the models in Refs. [10,12], their mean-field glass transitions are turned into a crossover in finite dimensions or their low-temperature glassy states are unstable due to crystallization. The other model [11], while having the typical glassy dynamics, is a monodisperse model, where an efficient algorithm such as the swap method [8,9] is not known. Another lattice glass model was proposed in Ref. [13], which is shown to have irregular configurations as an ordered state. However, its autocorrelation function decays without any plateau even at low temperature. Thus the model does not have the essential features of the structural glasses.

To find finite-dimensional models that allow us to access equilibrium low-temperature states without crystallization is still actively discussed. In this letter, we propose a simple lattice glass model to study glassy behaviors in three dimensions. Our model, in contrast to several lattice models [10,12,13], shows typical two-step relaxation dynamics as observed in fragile supercooled liquids. We study a binary mixture of the model, where the non-local swap dynamics explained below provide the benefits for equilibration. By large-scale Monte Carlo simulations, we equilibrate the system at temperature well below that where two-step relaxation emerges. We study the effect of a coupling field conjugate to the overlap between the system and a quenched configuration using the Wang-Landau algorithm [14,15], and compute the quenched version of the Franz-Parisi potential [16][17][18], an effective free energy of the glass transition. Our results show that the system indeed has a thermodynamic behavior consistent with the RFOT theory.

A lattice glass model we study in this letter is a binary mixture of particles defined by the Hamiltonian with a positive constant A

H (n, σ) = A L 3 i=1 n i j∈∂i n j -σi 2 , (1)
where the first summation runs over all the lattice sites. where a crystal structure is mainly composed of particles with 1 = 3 and 2 = 5, respectively. The gas, liquid, crystalline phases are separated by first-order phase transitions. In the glassy region, the system shows an aging effect and its autocorrelation function relaxes in two steps. We study the system with ρ2/ρ = 0.6 in this work.

The lattice is the simple cubic lattice with linear dimension L. We denote the occupancy of site i with n i ∈ {0, 1}, and the type of a particle at site i with σ i ∈ {1, 2}. The parameter σi determines the most favorable number of neighboring particles of type σ i particle. Here, we consider a binary mixture of 1 = 3 and 2 = 5 particles. The boundary conditions in all the directions are periodic. We study this model at finite temperature T by Monte Carlo simulation that preserves the number of each type of particles [19]: Two randomly chosen particles are swapped or a randomly chosen particle is moved to a vacant site of the lattice chosen also randomly with the Metropolis probability. Depending on the density of each type of particle, ρ 1 and ρ 2 , with fixed total density ρ = ρ 1 + ρ 2 = 0.75, our model has crystal phases at low temperature, see Fig. 1. When ρ 1 = 0.3 and ρ 2 = 0.45, we confirmed the absence of a drop in the energy and a large peak in the specific heat in the glassy region, indicating no crystallization.

Our model is somewhat related to a softened model of a lattice glass proposed by Biroli and Mézard (BM) [10,20]. The Hamiltonian of their model H BM = A i=1 n i θ j∈∂i n jσi j∈∂i n jσi with the Heaviside step function θ(x). The crucial difference between the two models is the number of neighboring particles that each particle energetically favors: In the soft BM model, any number of neighboring particles lower than σ achieves the lowest energy while only σ does in our model. This slight change makes the entropy of our model at low temperature smaller than that of the BM model and the dynamics more similar to fragile supercooled liquids.

The dynamics of our model is studied by a local Monte Carlo algorithm [19], which usually gives dynamics qualitatively similar to molecular and Brownian dynamics. In our simulation, particles can move to only neighboring sites. While any physical quantity of our model relaxes rapidly at high temperature, the relaxation gets extremely slow with decreasing temperature. To quantify the slow dynamics of the model, we study the autocorrelation function

C (t; t w ) = 1 N i δ σi(tw),σi(tw+t) -C 0 1 -C 0 , ( 2
)
where t w is the waiting time,

C 0 = ρ(ρ 2 1 + ρ 2 2
), N = ρL 3 the total number of particles, and the summation is taken over sites with n i (t w ) = 1. We also measure the dynamical susceptibility χ 4 (t; t w ) characterizing the dynamical heterogeneity observed in supercooled liquids. In the limit t w → ∞, the system is in equilibrium and we denote C (t) = C (t; t w → ∞) and χ 4 (t) = χ 4 (t; t w → ∞). We set t w to a sufficiently large value to study the equilibrium dynamics of the model so that C (t; t w ) and C (t; t w /10) agree with each other. Typical values of t w range from 10 2 to 10 8 Monte Carlo sweeps depending on temperature. Whereas the autocorrelation function of the system decays rapidly at high temperature, a two-step relaxation emerges in C(t) at temperature lower than T /A 0.6, see Fig. 2 (see also [START_REF]See Supplemental Material for a further discussion on the dynamics, the finite-size dependence of the relaxation time, and the determination of the dynamical transition temperature T d[END_REF] for a physical interpretation of the plateau in C(t)). The dynamical susceptibility χ 4 (t) shows a peak at finite time, indicating the emergence of heterogeneous dynamics of the system [7,[START_REF] Kob | [END_REF][23][24][25]. The peak value grows with decreasing temperature (see inset of Fig. 2), and it suggests that the dynamics gets more heterogeneous at lower temperature.

The RFOT theory of the glass transition predicts slow dynamics frozen at a dynamical transition temperature without any thermodynamic anomaly. In the vicinity of the dynamical transition temperature, the relaxation time diverges algebraically. However, an exponentially growing relaxation time well described by the Vogel-Fulcher-Tammann (VFT) law has been observed in experiments of fragile super-cooled liquids. Activated process in finite dimensions is supposed to wipe out the mean-field dynamical transition. Here, we study temperature dependence of the relaxation time τ α measured as a time when C (t) decays to e -1 . Around temperature where the two-step relaxation emerges, the relaxation time shows a super-Arrhenius growth with decreasing temperature (see Fig. 3). We find that the VFT law fits our data at lower temperatures very well with T VFT /A = 0.177(6) (see also [START_REF]See Supplemental Material for a further discussion on the dynamics, the finite-size dependence of the relaxation time, and the determination of the dynamical transition temperature T d[END_REF] for estimation of the dynamical transition temperature). Note that, as the relaxation time at low temperature increases with the system size, T VFT would shift to higher temperature in larger systems [START_REF]See Supplemental Material for a further discussion on the dynamics, the finite-size dependence of the relaxation time, and the determination of the dynamical transition temperature T d[END_REF].

In the RFOT theory, the overlap between independent replicas has been studied as an order parameter for the thermodynamic glass transition. The overlap successfully detects the glass transition in systems without spatial symmetry [26][27][28][29]. The overlap can detect even the slow dynamics at low temperature through its effective free energy [16][17][18]. In our model, however, the overlap and its distribution function can have no anomaly at any temperature due to its spatial symmetry [START_REF] Mézard | Glasses and replicas[END_REF]. The temperature dependence of the long-time limiting value of the dynamical susceptibility χ 4 (t → ∞), equivalent to the spin glass susceptibility χ SG , is indeed almost independent of temperature (see inset of Fig. 2). This is in contrast to a Potts spin glass model [START_REF] Takahashi | [END_REF] where χ SG in-creases with approaching its transition temperature [32].

To study rare fluctuations and the effective free energy of the overlap appropriately, we introduce a field ε in our model: A system at temperature T is coupled by the field to a quenched configuration randomly sampled at the same temperature in equilibrium [16][17][18]. The field explicitly breaks the invariance associated with the spatial symmetry, and thus the overlap can have nontrivial distributions. In the RFOT theory, the field induces a first-order transition at low temperature and the transition terminates at a critical point belonging to the random-field Ising model universality class [33]. The Franz-Parisi potential reveals the existence of metastable states at low-temperature and allows us to compute the configurational entropy as a free-energy difference between two minima [34,35]. Here, we use the Wang-Landau algorithm [14,15] with the multi-overlap ensemble [36] to compute the density of states Ω T (Q) as a function of the overlap Q at a given temperature T with a fixed reference configuration. The Franz-Parisi potential V q (Q) is computed from Ω T (Q) by averaging over quenched reference configurations. We prepared reference configurations by equilibrating the system at each temperature for 10 8 Monte Carlo sweeps with non-local swaps of particles. The non-local swap dynamics is faster with a factor of ∼ 10 2 at low temperature than the local swap dynamics while the factor slightly depends on temperature and the system size. The number of reference configurations is 912 for T /A = 0.31 and 0.33, and 48 for higher temperatures. At high temperature, the potential V q (Q) is convex and the overlap [ Q ] as a function of the field ε grows gradually, see Fig. 4 and Fig. 5 (c). Here, the brackets • and [•] represent the thermal and the reference-configuration averages, respectively. At T /A = 0.33 and 0.31, the Franz-Parisi potential is slightly non-convex (see Fig. 4), and the probability distribution of the overlap [P q (Q)] with finite ε has two

V q (Q) Q FIG. 4.
The Franz-Parisi potential Vq(Q) of the system with N = 1296 (L = 12). At high temperature the potential is convex while it is slightly non-convex at low temperature, T /A = 0.31. Inset shows an enlarged view. The broken line is a guide to eyes to show the non-convexity. separated peaks whereas that at temperature higher than T /A = 0.35 shows a clear single peak at any ε, see Fig. 5 (a) and (b), respectively [37]. We thus conclude that the coupling field induces a first-order transition into our model that terminates at a critical point at finite temperature as in mean-field and particle models in finite dimensions [17,18,33,38,39] (see Fig. 5

(d) for T -ε phase diagram).

When we assume the existence of the thermodynamic glass transition at T K > 0 for a system with spatial symmetry, we find the system-size dependence of the effective transition point ε c (L) and the configurational entropy density s conf (L) as follows. In the thermodynamic limit, an infinitesimal coupling field should make them have a finite overlap at T < T K whereas the overlap is strictly zero when ε = 0 due to spatial symmetry [START_REF] Mézard | Glasses and replicas[END_REF]. In finite systems, the first-order transition at finite field thus never goes to the expected T K , even in the limit ε → +0, while it may approach zero temperature in the limit. At temperature lower than T K , the effective transition point ε c should scale as ε c = O(L -a ) while at higher temperature ε c -ε

(∞) c = O(L -a ) where ε (∞) c
= lim L→∞ ε c (L) > 0 and a > 0, which depends on systems [40,41]. Regarding that the configurational entropy density s conf measured as the free energy difference in the Franz-Parisi potential is almost equivalent to ε c , we expect s conf of finite sys- tems to be finite even below expected T K , but decreases with ∼ L -a , as ε c (L) does. Studying the finite-size dependence of the first-order transition line and the configurational entropy would be a decisive test for the RFOT theory.

To measure the thermodynamic properties at low temperature, we use the non-local swap dynamics of randomly chosen pairs of particles and the exchange Monte Carlo (or parallel tempering) method [42] to further enhance equilibration. We also utilize the multipletemperature reweighting technique [43][44][45], which significantly improves the accuracy of Monte Carlo results. Typical number of Monte Carlo sweeps for equilibration range from 10 8 to 10 11 per site depending on the system size.

At temperature T /A 0.6 where the two-step relaxation emerges, the specific heat is rather smooth, and it starts to drop without a divergent behavior at lower temperature (see Fig. 6 (a)). To study further the systemsize dependence of the specific heat, we compute the temperature-derivative of the specific heat ∂c/∂T . We find that ∂c/∂T grows with increasing the system size, indicating that the drop gets steeper with the system size. Similar size-dependence has been observed in a meanfield model with RFOT [46] and a three-dimensional Potts glass model that has a RFOT-like spin-glass transition [START_REF] Takahashi | [END_REF].

Assuming the hyperscaling relation of the critical exponents dν = 2 -α with d the spatial dimension, we find that the peak value of ∂c/∂T follows the finite-size scaling relation

∂c ∂T * ∝ L θ (3)
with θ = (α + 1)/ν. We thus can evaluate the two conventional critical exponents from θ as ν = 3/ (θ + d) and α = 2 -dν = 2 -3d/ (θ + d). If the specific heat does not diverge but has a finite jump at a thermodynamic glass transition as in the RFOT theory, the critical exponent α = 0, suggesting ν = 2/d. The peak value (∂c/∂T ) * thus grows algebraically with exponent θ = 1/ν = d/2.

This provides a useful way for studying thermodynamic anomaly with a finite-size scaling analysis even when an order parameter is unknown. Here, in our model, the exponent θ = 1.44(2) (see inset of Fig. 6 (b)), implying ν = 0.68(1) and α = -0.03(1), marginally compatible with the RFOT theory.

In summary, we proposed a lattice glass model that is very stable against crystallization and shows the typical dynamics observed in fragile supercooled liquids at low temperature. Our numerical results show that a glass transition is detected by a singularity of the temperaturederivative of the specific heat with exponents compatible with the RFOT theory. We also numerically computed the quenched version of the Franz-Parisi potential using the Wang-Landau algorithm. The potential shows very large overlap fluctuations at low temperature, and a firstorder transition was found in the coupled system that terminates at a critical point. The large overlap fluctuations strongly suggest that low-temperature glassy dynamics emerge from thermodynamics of the system, similar to the mean-field models with the RFOT. We thus conclude that our model is useful to study and further examine the mean-field RFOT predictions. Thanks to the lattice nature of our model, it is easy to study its mean-field solution using the cavity method and compare it with our results. To determine precisely the phase diagram both in mean-field and finite dimensions will constitute a crucial test of the validity of the theory.

Supplemental Item 1: Dynamics at low temperature In off-lattice particle models for supercooled liquids and glasses, the dynamics at low temperature shows two-step relaxation as well as our model. The physical interpretation of the plateau regime in the particle models is well known that particles vibrate locally inside cages effectively formed by surrounding particles [7]. However, since particles of our model defined on a lattice have discrete degrees of freedom, the appropriate physical interpretation analogous to the local vibration in the plateau regime is unclear. In this supplemental item, we show evidence of local vibrations of particles during the plateau regime.

Thanks to the discrete nature of our model, it is easy to find lattice sites that have no contribution to the autocorrelation function Eq. ( 2), where δ σi(tw),σi(tw+t) = 0, i.e. the lattice sites that have different types of particles at t = 0 and t > 0 or that are vacant at t > 0 but occupied at t = 0. Roughly speaking, a cluster represents a movable region during time interval t. We identify connected clusters of those lattice sites at time t, and compute the normalized distribution of the cluster size n c . Fig. S1 shows the distribution of the connected-cluster size depending on the time t .

At temperature T /A = 0. 10, and the distribution is almost independent of time, see Fig. S1. Although the size distribution does not depend on time very much, the positions of the clusters change with time during the plateau, and some of them disappear with time (see Fig. S2 for real-space visualization of the changed lattice sites). We thus conclude, at short time scale t 10 3 , the system has local vibrations that involve only O(1) number of particles. These local vibrations produce the plateau in the autocorrelation function, even in the lattice model.

As time passes, the distribution of n c has a broad tail towards larger n c , and eventually becomes double-peaked at large time t τ α 10 7 . The emergence of two peaks in the distribution with small clusters of n c = O(1) and giant clusters of n c = O(N ) indicates the heterogeneous dynamics. At time t = 10 7 , the largest cluster spreads over the system while there are still some regions where only small clusters occupy, see Fig. S2. In the RFOT theory, the relaxation time diverges algebraically at temperature T d > T K rather than exponentially as observed in experiments. Whereas the dynamical transition turns into a crossover in finite dimensions, a well-defined localization transition of the potential energy landscape was shown to control the dynamical cross-over [48]. However, their method to identify the localization transition is not available in our model. We thus estimate the dynamical transition temperature T d by a more traditional method, that is, power-law fitting of the relaxation time as done in Ref. [10].

We show in Fig. S3 the best fit of the inverse relaxation time to a power-law behavior ∼ (T -T d ) γ , with T d = 0.348(4) and γ = 4.18 (12). In mean-field models, the dynamical transition temperature is close to the critical temperature of the ε-coupled system [17,18,33]. Our estimation of T d is compatible with the temperature-coupling phase diagram Fig. 6. The relaxation time of glass forming liquids has a decreasing behavior at moderately low temperature [47,49]. The magnitude of the decrease is especially remarkable just below the dynamical transition temperature (or the modecoupling temperature) [49]. In our model, with the local dynamics, the relaxation time does not decrease with the system size even around the dynamical transition temperature T d = 0.348(3), see Fig. S4. At lower temperature, the relaxation time rather increases with the system size. The finite-size effect is more prominent in the relaxation time of the non-local dynamics we used for equilibrium computation, see Fig. S4. Intuitively, the finite-size effect is clearer because of the following reason. The relaxation process at moderately low temperature is mixed with many modes, with each contribution. The slowest mode of the relaxation process, which diverges towards the transition temperature and has strong finite size effects, is buried in many other modes due to its small contribution in the dynamics of inefficient algorithms, including simple local dynamics. The relaxation time shows a finite-size effect only when the slowest mode contribution becomes very large, which is very low temperature. Efficient algorithms, on the other hand, can reduce contributions of modes faster than the slowest mode but still very slow, and thus a finite-size effect that the slowest mode has is seen more clearly. This is actually quite general in other models, even in simple ferromagnetic spin models, by comparing the dynamics of the simple Metropolis algorithm and more efficient algorithms such as the over-relaxation, for instance.

FIG. 1 .

 1 FIG. 1. Phase diagram of the model in the plane of temperature T and type-2 particle density ρ2 with the total density ρ = 0.75. Crystal3 and Crystal5 stand crystalline phases where a crystal structure is mainly composed of particles with

  FIG. 2.(a) The autocorrelation function of the system. Inset shows an enlarged view around the plateau region. (b) The dynamical susceptibility χ4 (t) of the system. Inset shows temperature dependence of the peak and the plateau values of χ4 (t) that are denoted as χ * 4 and χ p 4 , respectively. The system size N = 6000 (L = 20).

  FIG.3. The Arrhenius plot of the relaxation time τα of the system with N = 6000 (L = 20). The broken and dotted lines are the Arrhenius law with ∆E = 4 and the VFT law with TVFT/A = 0.177 respectively for guides to eyes.

FIG. 5 .

 5 FIG. 5. (a, b) The probability distribution of the overlap averaged over quenched reference configurations [Pq(Q)] with finite coupling field ε at T /A = 0.33 (a) and 0.45 (b), respectively. (c) The overlap [ Q ] as a function of the field ε. With increasing ε, [ Q ] grows smoothly at high temperature while it shows a rapid increase at low temperature. The system size N = 1296 (L = 12). (d) Temperature-coupling phase diagram of the system. The empty points and their errors are estimated by the peak of the spin glass susceptibility χSG from the data for L = 12. The filled circle point represents an expected location of the critical point that is not determined precisely here. The triangular point TVFT is from the VFT-law fit of the relaxation time for L = 20, see Fig. 3. The diamond point indicates the dynamical transition temperature Tc 0.348(4) estimated by fitting [21].

FIG. 6 .

 6 FIG.6. Temperature dependence of (a) the specific heat c and (b) its temperature derivative ∂c/∂T of the system. Inset of (b) shows the peak value of ∂c/∂T as a function of L.

7 FIG

 7 FIG. S1. The normalized distribution function of the size of connected clusters at T /A = 0.35 for (a) 10 1 ≤ t ≤ 5 × 10 5 and (b) 3 × 10 6 ≤ t ≤ 8 × 10 7 . The system size L = 20 (N = 6000).

  FIG. S2.Real-space visualization of typical configurations at time t = 10 1 , 10 2 , 5 × 10 2 , 10 3 , and 10 7 , from left top to right bottom. The temperature T /A = 0.35. The lattice sites which have different particles in the configuration at finite time t > 0 from the original configuration at t = 0 are filled with cubes. Here, at this temperature, the autocorrelation function has a plateau when 10 1 t 10 3 . For the configuration at t = 10 7 , the lattice sites belonging to the largest connected cluster are filled with blue cubes, and those for the other clusters are filled with red cubes.

1 α

 1 FIG. S3. The inverse relaxation time τ -1 α as a function of (T -T d )/T d and a power-law curve ∼ (T -T d ) γ . The dynamical transition temperature Tc and the exponent γ are estimated by fitting using the data at 0.4 ≤ T /A ≤ 0.6, as T d = 0.348(4) and γ = 4.18(12), respectively.

FIG

  FIG. S4. The relaxation time as a function of inverse temperature for L = 5, 10, and 20 for the local and non-local dynamics. The inverse dynamical transition temperature A/T d is estimated in Supplemental Item 2.

  35, the autocorrelation function has a plateau from time t 10 to t 10 3 , see Fig 2. In this time regime, the size of connected clusters n c

The authors thank A. Ikeda for many useful discussions. Y.N. is grateful to J. Takahashi, T. Takahashi, S. Takabe, H. Yoshino, M. Ozawa, L. Berthier and D. Coslovich for critical comments and useful discussions. This work was supported by JSPS KAKENHI Grant Numbers 17J10496, 17H02923 and 19H04125. Numerical simulation in this work has mainly been performed by using the facility of the Supercomputer Center, Institute for Solid State Physics, the University of Tokyo.


Mirko Fiacchini

Yet another computation-oriented necessary and sufficient condition for stabilizability of switched linear systems

Keywords: Switched linear systems, stabilizability, convex analysis

This paper presents a computational method to test the stabilizability of discrete-time switched linear systems. The existence of a conic cover of the space on whose elements a convex condition holds is proved to be necessary and sufficient for stabilizability. An algorithm for computing a conic partition that satisfies the new necessary and sufficient condition is given. The algorithm, that allows also to determine bounds on the exponential convergence rate, is proved to overcome the conservatism of conditions equivalent to periodic stabilizability and is applied to a four dimensional system.

I. INTRODUCTION

Switched systems are characterized by a finite set of possible dynamics, among which the system evolves [START_REF] Liberzon | Switching in Systems and Control[END_REF], [START_REF] Sun | Stability Theory of Switched Dynamical Systems[END_REF]. In the last decades, this class of systems attracted an increasing attention since it allows to model complex dynamics such as networked and interconnected systems and the interaction between physical systems and digital devices. Classical issues in automatic control, such as stability analysis and control design, result to be rather involved even for switched linear systems, though, yielding to the necessity of a tailored theory [START_REF] Lin | Stability and stabilizability of switched linear systems: a survey of recent results[END_REF], [START_REF] Sun | Stability Theory of Switched Dynamical Systems[END_REF]. Concerning the problem of stability analysis, several works appeared presenting sufficient and necessary and sufficient conditions. For instance, [START_REF] Margaliot | Stability analysis of switched systems using variational principles: An introduction[END_REF] provides stability analysis based on a variational approach, whereas the work [START_REF] Daafouz | Stability analysis and control synthesis for switched systems : A switched Lyapunov function approach[END_REF] considers mode-dependent Lyapunov functions to provide a sufficient condition for stability.

Concerning the problem of stabilizability of switched linear systems, the necessity of using nonconvex or time-varying functions have been recognized, see [START_REF] Blanchini | Stabilizability of switched linear systems does not imply the existence of convex Lyapunov functions[END_REF], and often employed to reduce the conservatism in constructive methods, for instance those based on Lyapunov-Metzler conditions [START_REF] Geromel | Stability and stabilization of continuoustime switched linear systems[END_REF], [START_REF] Geromel | Stability and stabilization of discretetime switched systems[END_REF], on convex conditions [START_REF] Fiacchini | On the stabilizability of discrete-time switched linear systems: novel conditions and comparisons[END_REF], [START_REF] Sun | Stability Theory of Switched Dynamical Systems[END_REF], and on quadratic time-varying Lyapunov functions [START_REF] Deaecto | Stability and performance of discretetime switched linear systems[END_REF]. Also necessary and sufficient conditions for stabilizability appeared, [START_REF] Fiacchini | Necessary and sufficient condition for stabilizability of discrete-time linear switched systems: A set-theory approach[END_REF], [START_REF] Sun | Stability Theory of Switched Dynamical Systems[END_REF], highlighting the inherent complexity of the problem of determining whether a switched linear systems is stabilizable. As a matter of fact, testing whether a necessary and sufficient condition for stabilizability holds is equivalent to check if the unit ball is contained in the interior of the union of a potentially arbitrarily big number of ellipsoids, given by the ball preimages. The stabilizability problem has been proved, indeed, to be undecidable in [START_REF] Jungers | On feedback stabilization of linear switched systems via switching signal control[END_REF], that addresses the problem in terms of joint spectral subradius, see also [START_REF] Jungers | The Joint Spectral Radius: Theory and Applications[END_REF]. M. Fiacchini is with Univ. Grenoble Alpes, CNRS, Grenoble INP, GIPSAlab, 38000 Grenoble, France. mirko.fiacchini@gipsa-lab.fr.

Research funded in part by ANR via project HANDY, number ANR-18-CE40-0010.

A strictly weaker condition, referred to as uniform convergence [START_REF] Stanford | Some convergence properties of matrix sets[END_REF], consistent [START_REF] Sun | Stabilizability and insensitivity of switched linear systems[END_REF], [START_REF] Sun | Stability Theory of Switched Dynamical Systems[END_REF] or periodic [START_REF] Fiacchini | On the stabilizability of discrete-time switched linear systems: novel conditions and comparisons[END_REF] stabilizability, has been also characterized and its relation with general stabilizability analysed. While [START_REF] Stanford | Some convergence properties of matrix sets[END_REF]- [START_REF] Sun | Stability Theory of Switched Dynamical Systems[END_REF] already proved that periodic stabilizability is only sufficient for stabilizability, [START_REF] Fiacchini | On the stabilizability of discrete-time switched linear systems: novel conditions and comparisons[END_REF] provides a necessary and sufficient convex condition for periodic stabilizability and its relation with Lyapunov-Metzler conditions. A hint on the conservatism of conditions for periodic stabilizability has been given in [START_REF] Heemels | On Lyapunov-Metzler inequalities and S-procedure characterizations for the stabilization of switched linear systems[END_REF], expressing it in terms of S-procedure application. Thus, although necessary and sufficient conditions are available, their complexity led to inspect alternative conditions, just sufficient but more computationally affordable. The notable exception [START_REF] Jungers | On feedback stabilization of linear switched systems via switching signal control[END_REF] proposes algorithms to test necessary and sufficient conditions for stabilizability and to compute tight bounds on the convergence rate.

This paper aims at providing an algorithmic method to test whether a necessary and sufficient condition for stabilizability is satisfied for a switched linear system. The method is based on a novel necessary and sufficient stabilizability condition, consisting in the existence of a conic cover of the space on whose cones a convex condition holds. For switched systems, in fact, conditions based on homogeneous functions and defined on conic partions permit to characterize global properties, see for instance [START_REF] Iervolino | Cone-copositive piecewise quadratic Lyapunov functions for conewise linear systems[END_REF], [START_REF] Lin | Synthesis of uniformly ultimate boundedness switching laws for discrete-time uncertain switched linear systems[END_REF], [START_REF] Lin | Switching stabilizability for continuoustime uncertain switched linear systems[END_REF]. A major contribution of this paper is the formulation and application of a computational method for generating a sequence of conic covers of the state space for which the new necessary and sufficient condition for stabilizability is satisfied. The method allows also to determine bounds on the exponential convergence rate for a given switching horizon. Although no guarantee for the termination of the procedure can be given, being the problem itself undecidable [START_REF] Jungers | On feedback stabilization of linear switched systems via switching signal control[END_REF], rather tight numerical results have been obtained. In particular, the method is applied to certify the stabilizability of a system for which the sufficient conditions have been proved not to hold [START_REF] Fiacchini | On the stabilizability of discrete-time switched linear systems: novel conditions and comparisons[END_REF]. Moreover, it has been applied to a two dimensional classical example for which tight bounds on the convergence rate are available in the literature [START_REF] Jungers | On feedback stabilization of linear switched systems via switching signal control[END_REF], showing that analogous estimations are obtained. Finally, a specific four dimensional switched system has been built for which the convergence rate bound can be inferred but hardly computed with alternative necessary and sufficient conditions. The numerical method proved to be able to certify the bound on the convergence rate even for the four dimensional system.

Notation: Given n ∈ N, define N n = {x ∈ N : 1 ≤ x ≤ n}. The subindex denotes, with slight abuse of notation, both the i-th row of a matrix (or a vector) and the i-th element of a countable set of elements. Given Ω ⊆ R n , int(Ω) is the interior of Ω. Given P ∈ R n×n with P > 0, define E (P) =

x ∈ R n : x T Px ≤ 1 . The q modes of the switched systems are denoted by I , that is I = N q , then I k = ∏ k j=1 I are all the possible sequences of modes of length k and I [N] = N k=1 I k , the sequences of modes of length up to N. Given N ∈ N, the number of elements i ∈ I

[N] is N = ∑ N k=1 q k . Given i ∈ I [N] define A i = A i N . . . A i 1 and B i = {x ∈ R n : x T A T i A i x ≤ 1}, with B unit ball: B = {x ∈ R n : x T x ≤ 1}. Define J = {-1, 1} and given β ∈ J n , define the matrix T (β ) ∈ R n×n as T (β ) = diag(β ). The simplex in R n is denoted as ∆ n , that is ∆ n = {λ ∈ R n : λ ≥ 0, ∑ n i=1 λ i = 1}. The sets {Γ i } i∈N p are a cover of R n if i∈N p Γ i = R n .

II. PROBLEM STATEMENT

Consider the switched linear system

x k+1 = A σ (k) x k , (1)
where x k ∈ R n is the state at time k ∈ N, σ : N → I is the switching law and {A i } i∈I , with A i ∈ R n×n for all i ∈ I , the modes transition matrices. Denote with x σ k (x) the state of (1) at time k for a switching law σ (•) and an initial condition x. The objective is to provide a constructive condition for stabilizability, defined below for the case under analysis.

Definition 1: The system (1) is globally exponentially stabilizable if there are c ≥ 0 and λ ∈ [0, 1) and, for all x ∈ R n , there exists a switching law σ : N → I , such that

x σ k (x) ≤ cλ k x , ∀k ∈ N. (2)
Since asymptotic and exponential stabilizability are equivalent for this class of systems, there is no loss of generality in considering exponential stabilizability. Necessary and sufficient conditions exist for system (1) to be stabilizable, given first in [START_REF] Sun | Stability Theory of Switched Dynamical Systems[END_REF] and in [START_REF] Fiacchini | Necessary and sufficient condition for stabilizability of discrete-time linear switched systems: A set-theory approach[END_REF] in a set-theoretic version, recalled here.

Theorem 1 ( [START_REF] Fiacchini | Necessary and sufficient condition for stabilizability of discrete-time linear switched systems: A set-theory approach[END_REF]): The switched system (1) is stabilizable if and only if there exists N ∈ N such that

B ⊆ int i∈I [N] B i . (3)
These results claim substantially that system (1) is stabilizable if and only if a state-dependent switching law exists such that the Euclidean norm of the state can be reduced within the horizon of length N. The main problem, that is the objective of the present study, is that checking whether (3) holds might be hardly manageable, for nontrivial cases. Indeed, the problem has been proved to be inherently undecidable, [START_REF] Jungers | On feedback stabilization of linear switched systems via switching signal control[END_REF].

Other conditions, that are more computationally affordable but just sufficient, have been given in the literature. A first condition is the one based on Lyapunov-Metzler inequalities, [START_REF] Geromel | Stability and stabilization of continuoustime switched linear systems[END_REF], [START_REF] Geromel | Stability and stabilization of discretetime switched systems[END_REF], posed in terms of bilinear matrix inequalities. Another sufficient condition for stabilizability, in form of convex optimization problem, is given in [START_REF] Fiacchini | On the stabilizability of discrete-time switched linear systems: novel conditions and comparisons[END_REF] and recalled hereafter.

Theorem 2 ( [START_REF] Fiacchini | On the stabilizability of discrete-time switched linear systems: novel conditions and comparisons[END_REF]):

The switched system (1) is stabilizable if there exist N ∈ N and η ∈ R N such that η ∈ ∆ N and ∑ i∈I [N] η i A T i A i < I. (4)
Condition ( 4) has been proved in [START_REF] Fiacchini | On the stabilizability of discrete-time switched linear systems: novel conditions and comparisons[END_REF] to be equivalent to periodic stabilizability, also referred to as uniform convergence [START_REF] Stanford | Some convergence properties of matrix sets[END_REF] or consistent stabilizability [START_REF] Sun | Stabilizability and insensitivity of switched linear systems[END_REF], [START_REF] Sun | Stability Theory of Switched Dynamical Systems[END_REF], and less conservative than the Lyapunov-Metzler condition. Moreover, it has been proved in [START_REF] Stanford | Some convergence properties of matrix sets[END_REF], [START_REF] Sun | Stability Theory of Switched Dynamical Systems[END_REF] that periodic stabilizability is only sufficient, not necessary for stabilizability. This means that switched linear systems exist such that condition (3) is satisfied with an appropriate N ∈ N whereas (4) cannot hold, for every N ∈ N, nor Lyapunov-Metzler conditions can be satisfied, see the Example 17 in [START_REF] Fiacchini | On the stabilizability of discrete-time switched linear systems: novel conditions and comparisons[END_REF] that will be considered in this paper. A hint on the source of conservatism of the condition for periodic stabilizability is given in [START_REF] Heemels | On Lyapunov-Metzler inequalities and S-procedure characterizations for the stabilization of switched linear systems[END_REF], where condition ( 4) is proved to be obtainable by resorting to the S-procedure.

The objective of this work is to provide a computational method for obtaining a test for stabilizability which is less conservative than the sufficient ones, in the line of [START_REF] Jungers | On feedback stabilization of linear switched systems via switching signal control[END_REF], and to show its efficiency even in four dimensions.

III. NECESSARY AND SUFFICIENT CONDITIONS FOR STABILIZABILITY

A useful result, whose proof is in [START_REF] Fiacchini | On the stabilizability of discrete-time switched linear systems: novel conditions and comparisons[END_REF], is recalled. Lemma 1 ([6]): Given P i ≥ 0, i ∈ N m , the following equality holds

π∈∆ m E ∑ i∈N m π i P i = i∈N m E (P i ).
A first novel necessary and sufficient condition for stabilizability is presented below. Based on this condition, another one, more computation oriented, is given afterward.

Theorem 3: The switched system (1) is stabilizable if and only if there exists N ∈ N and, for every

x ∈ R n , there is η(x) ∈ R N such that η(x) ∈ ∆ N and x T ∑ i∈I [N] η i (x)A T i A i x < x T x. (5)
Proof: To prove sufficiency, suppose η(x) ∈ ∆ N exists for every x ∈ R n such that (5) holds. Condition ( 5) is equivalent to the fact that if x is in B then it also belongs to the interior of the set

E ∑ i∈I [N] η i (x)A T i A i = {z ∈ R n : z T ∑ i∈I [N] η i (x)A T i A i z ≤ 1}. (6) Denoting Λ = {η(x) : x ∈ R n }, that is a subset of ∆ N , and from Lemma 1, it follows that B ⊆ x∈R n int E ∑ i∈I [N] η i (x)A T i A i = η∈Λ int E ∑ i∈I [N] η i A T i A i ⊆ η∈∆ N int E ∑ i∈I [N] η i A T i A i ⊆ int   η∈∆ N E ∑ i∈I [N] η i A T i A i   = int i∈I [N] B i ,
and then (3) holds. For necessity, suppose that (3) is satisfied. Then every x ∈ B, i.e. such that x T x ≤ 1, belongs also to the interior of i∈I [N] B i which implies, from Lemma 1 and ( 6), that η(x) ∈ ∆ N exists such that

x T ∑ i∈I [N] η i (x)A T i A i x < 1.
Hence, (5) holds. Theorem 3 substantially claims that the existence of η(x) ∈ ∆ N such that (5) is satisfied, for every x ∈ R n , is necessary and sufficient for the stabilizability of the switched linear system. On the other hand, the determination of such η(x) defined on the whole R n might be computationally intractable. Note that Theorem 3 could have been alternatively proved by noticing that ( 5) is equivalent to the existence, for every x ∈ R n different from the origin, of i(x) ∈ I [N] , such that A i(x) x 2 < x 2 , for appropriate N ∈ N, and then also equivalent to [START_REF] Boyd | Linear Matrix Inequalities in system and control theory[END_REF].

Remark 1: The sufficient condition (4) can be seen as the particular case of ( 5) with η = η(x) for all x ∈ R n .

The alternative computation-oriented necessary and sufficient condition can now be given.

Theorem 4: The switched system (1) is stabilizable if and only if there exist: N ∈ N; a conic cover {C ( j) } j∈N p of R n ; and {µ ( j) } j∈N p with µ ( j) ∈ ∆ N such that

x T ∑ i∈I [N] µ ( j) i A T i A i x < x T x, ∀x ∈ C ( j) . (7)
Proof: Sufficiency follows directly from Theorem 3 with

η i (x) = µ ( j) i
for all i ∈ I [N] , j ∈ N p , and x ∈ C ( j) . To prove necessity, suppose the system is stabilizable, that is equivalent to satisfaction of (3) for

N ∈ N. Then ρ > 0 exists such that 1 √ 1 -ρ B ⊆ int i∈I [N] B i . (8)
Condition [START_REF] Geromel | Stability and stabilization of continuoustime switched linear systems[END_REF] implies that for every x ∈ R n there exists η(x) ∈ R N such that x T P(x)x < -ρx T x,

with

P(x) = ∑ i∈I [N] η i (x)A T i A i -I, (10)
from Theorem 3. The aim is to prove that (9) implies the existence of the conic covering of R n and µ ( j) for every cone such that (7) holds. From homogeneity, only x ∈ ∂ B, i.e. with x 2 = 1, are considered in the following. Given x ∈ ∂ B, the function defined for y ∈ R n as f x (y) = (x + y) T P(x)(x + y)x T P(x)x is continuous in y and such that f x (0) = 0, and thus is smaller than ρ for y in a neighborhood of the origin. Then for every x ∈ ∂ B, there exists ε(x) positive such that f x (y) < ρ for all y ∈ ε(x)B. This implies that, taking ε * = min x∈∂ B ε(x), which exists and is positive, then [START_REF] Geromel | Stability and stabilization of discretetime switched systems[END_REF]. Hence, for all x on the boundary of the unit ball, the condition

(x + y) T P(x)(x + y) = f x (y) + x T P(x)x < 0, ∀y ∈ ε * B for all x ∈ ∂ B, from
z T ∑ i∈I [N] η i (x)A T i A i z < z T z,
holds for all z such that xz 2 ≤ ε * , and also for z in the cone generated by x + εB, i.e.

C x = {y ∈ R n : ∃z ∈ R n , ∃λ ∈R s.t. x -z 2 ≤ ε * , y = λ z, λ ≥ 0} ( 11
)
from homogeneity. Finally, defining: a grid of x ( j) ∈ ∂ B, with j ∈ N p , such that for all x ∈ ∂ B there is x ( j) such that x ( j)x 2 ≤ ε * ; the conic covering {C x ( j) } j∈N p from [START_REF] Iervolino | Cone-copositive piecewise quadratic Lyapunov functions for conewise linear systems[END_REF]; and µ ( j) = η(x (i) ), condition (7) holds. Theorem 4 means that the existence of a fine enough conic cover of the space such that (7) holds is a necessary and sufficient condition for stabilizability. Moreover, condition [START_REF] Fiacchini | Necessary and sufficient condition for stabilizability of discrete-time linear switched systems: A set-theory approach[END_REF] boils down to a linear condition if the conic sets C ( j) are defined as the negative sublevel sets of quadratic forms, as will be illustrated below. More details on the computational implication will be given in the following section.

Note that, being the condition given in Theorem 4 necessary and sufficient for stabilizability, it is equivalent to the other ones, [START_REF] Fiacchini | Necessary and sufficient condition for stabilizability of discrete-time linear switched systems: A set-theory approach[END_REF], [START_REF] Sun | Stability Theory of Switched Dynamical Systems[END_REF], and strictly less conservative than sufficient conditions, for instance those based on convex conditions [START_REF] Fiacchini | On the stabilizability of discrete-time switched linear systems: novel conditions and comparisons[END_REF], [START_REF] Sun | Stability Theory of Switched Dynamical Systems[END_REF], on Lyapunov-Mezler ones [START_REF] Geromel | Stability and stabilization of continuoustime switched linear systems[END_REF], [START_REF] Geromel | Stability and stabilization of discretetime switched systems[END_REF] and on time-varying Lyapunov functions [START_REF] Deaecto | Stability and performance of discretetime switched linear systems[END_REF], see [START_REF] Fiacchini | On the stabilizability of discrete-time switched linear systems: novel conditions and comparisons[END_REF] for a relation between some of them. Moreover, the satisfaction of the condition in Theorem 4, as for the equivalent ones, implies that the classical min-switch control strategies stabilize the switched system.

IV. CONIC PARTITION AND NONNEGATIVE QUADRATIC FORMS

The first issue is how to compute a set of quadratic forms determining a cover of the whole space, on whose cones condition [START_REF] Fiacchini | Necessary and sufficient condition for stabilizability of discrete-time linear switched systems: A set-theory approach[END_REF] has to be tested. That is, the design of a set of matrices {M j } j∈N p such that

j∈N p {x ∈ R n : x T M j x ≥ 0} = R n . (12)
Moreover, it might be necessary to iteratively generate finer covers of the state space. The presented method is based on generating covers of the state space composed by the closed convex cones determined by a basis matrix. Hence, for every cone, a quadratic form that is nonnegative on the cone is defined. Finally a method to generate finer covers of the state space, by replacing a cone with a conic cover of it, is given.

A. Initial conic partition generation

Any basis of the state space is the linear transformation of the standard basis of R n , denoted {e i } i∈N n . Indeed, clearly, the elements of any given basis

{b i } i∈N n with b i ∈ R n are Be i = [b 1 b 2 . . . b n ]e i = b i , for all i ∈ N n . Moreover, given a nonsingular matrix B ∈ R n×n the set C(B) = {x ∈ R n | ∃λ ∈ R n : Bλ = x, λ ≥ 0} ( 13
)
is a closed convex cone. Then given the basis {b i } i∈R n of R n , the closed convex cones C(BT (β )) defined for all β ∈ J n , are such that β ∈I n C(BT (β )) = R n . Thus, every nonsingular matrix B defines a conic partition of R n , that will be used in the algorithm as initial partition.

B. Nonnegative quadratic form

Given β ∈ J n and the cone C(BT (β )), a quadratic form that is nonnegative on C(BT (β )) must be defined. Given β ∈ J n , denote V = BT (β ) in this section, to ease to notation.

Lemma 2: Given {v i } i∈N n basis of R n and V = [v 1 v 2 . . . v n ], if M(V ) ∈ R n×n symmetric satisfies v T i M(V )v j ≥ 0, ∀i, j ∈ N n s.t. i ≤ j (14)
then x T M(V )x ≥ 0 for all x ∈ C(V ). Proof: By definition of the convex cone C(V ), for every x ∈ C(V ) there exists λ ∈ R n with λ ≥ 0 such that x = ∑ n i=1 λ i v i . Then x T M(V )x is given by a finite sum of terms λ i λ j v i M(V )v j , that are nonnegative from [START_REF] Liberzon | Switching in Systems and Control[END_REF]. Hence x T M(V )x ≥ 0 for all x ∈ C(V ).

Lemma 2 can be used to determine a quadratic form that is nonnegative on the cone generated by every {v i } i∈N n basis of R n . Recall that M(V ) might not be positive definite to determine nontrivial cones.

Proposition 1:

Given {v i } i∈N n basis of R n and V = [v 1 v 2 . . . v n ] the following inclusion holds C(V ) ∪C(-V ) ⊆ {x ∈ R n : x T M(V )x ≥ 0} ( 15
)
for all M(V ) satisfying ( 14).

Proof: The fact that C(V ) ⊆ {x ∈ R n : x T M(V )x ≥ 0} follows directly from Lemma 2. Moreover, if M(V ) satisfies the conditions in [START_REF] Liberzon | Switching in Systems and Control[END_REF], than it also satisfies them for {-v i } i∈N n , which implies that also

C(-V ) ⊆ {x ∈ R n : x T M(V )x ≥ 0}.
An implication of Proposition 1 is that some of the sets in the cover β ∈I n C(BT (β )) could be disregarded from the analysis, since redundant.

C. Finer partition generation

As it will be clearer from the following section, a method is necessary to split a closed convex cone C(V ) in smaller cones whose union contains C(V ). In particular, given the nonsingular matrix V , a first method is provided for determining n nonsingular matrices U j , with j ∈ N n such that

C(V ) = j∈N n C(U j ). ( 16
)
The simpler approach proposed consists basically in using the barycenter vector v(V ) = 1/n ∑ j∈N n v j to split the cone:

U j = [v 1 . . . v j-1 v(V ) v j+1 . . . v n ] ∀ j ∈ N n . (17)
As v(V ) belongs to the convex hull of {v i } i∈N n and matrices U j are nonsingular, for every j ∈ N n , then their columns are basis of R n and satisfy [START_REF] Lin | Switching stabilizability for continuoustime uncertain switched linear systems[END_REF]. This means that, given a conic partition of R n , another partition of R n is obtained by replacing the cone C(V ) with the cones C(U j ) as in [START_REF] Lin | Stability and stabilizability of switched linear systems: a survey of recent results[END_REF]. Alternatively, a cone can be splitted along the longer edge of the set co({v j } j∈N n ). This consists in defining

(i * , j * )(V ) ∈ arg max (i, j)∈N 2 n v i -v j , v * (V ) = 0.5(v i * + v j * ) (18)
to split the cone in two

U i * = [v 1 . . . v i * -1 v * (V ) v i * +1 . . . v n ], U j * = [v 1 . . . v j * -1 v * (V ) v j * +1 . . . v n ]. (19)
Also in this case matrices U i * and U j * are nonsingular and satisfy C(V ) = C(U i * ) ∪ C(U j * ). Moreover, the maximal distances between the vectors v j determining a partition would converge to zero by iterating the splitting procedure.

V. LOCAL CONTRACTION CONSTRAINT AND ALGORITHM

Now that a conic partition is given and a quadratic form that is nonnegative on every cone can be determined, the contraction conditions can be posed over every cone. Suppose that a set of basis V = {V ( j) } j∈N p determining the cover of R n as j∈N p C(V ( j) ) has been obtained as illustrated above, i.e. starting with an initial cover and then refining it if necessary. Also the following set of cones

j∈N p {x ∈ R n : x T M(V ( j) )x ≥ 0} is a cover of R n , from Proposition 1. Proposition 2: Given {v i } i∈N n basis of R n and V = [v 1 v 2 . . . v n ], the condition (7) is satisfied with C ( j) = C(V ) ∪C(-V ) and µ ( j) = µ(V ) if ∑ i∈I [N] µ i (V )A T i A i -I + M(V ) < 0, (20)
holds with M(V ) satisfying ( 14).

Proof: From Proposition 1, condition [START_REF] Fiacchini | Necessary and sufficient condition for stabilizability of discrete-time linear switched systems: A set-theory approach[END_REF] with

C ( j) = C(V ) ∪C(-V ) and µ ( j) = µ(V ) holds if x T ∑ i∈I [N] µ i (V )A T i A i x < x T x, ∀x s.t. x T M(V )x ≥ 0 (21)
is satisfied, and then also if [START_REF] Sun | Stabilizability and insensitivity of switched linear systems[END_REF] holds, from S-procedure [START_REF] Boyd | Linear Matrix Inequalities in system and control theory[END_REF]. Thus, condition [START_REF] Sun | Stabilizability and insensitivity of switched linear systems[END_REF] implies satisfaction of [START_REF] Fiacchini | Necessary and sufficient condition for stabilizability of discrete-time linear switched systems: A set-theory approach[END_REF] with

µ ( j) = µ(V ) and C ( j) = C(V ) ∪C(-V ).
The interest of this result lies in the fact that ( 20) is more suitable than [START_REF] Fiacchini | Necessary and sufficient condition for stabilizability of discrete-time linear switched systems: A set-theory approach[END_REF] from the computational point of view.

Finally, a sketch of algorithm is given for testing whether the condition in Theorem 4 is satisfied and hence if the system (1) is stabilizable, see Algorithm 1. If Algorithm 1 terminates, then the system is stabilizable.

Algorithm 1 Testing stabilizability condition in Theorem 4. Input:

Horizon N ∈ N, matrices {A i } i∈I [N] , basis matrix B. 1: V = {BT (β ) : β ∈ J n }; Initial conic partition 2: while V = / 0 do 3: extract V from V ; 4: compute α(V ) such that α(V ) = min α,µ(V ), M(V ) α s.t. ∑ i∈I [N] µ i (V )A T i A i + M(V ) ≤ αI µ(V ) ≥ 0, ∑ i∈I [N] µ i (V ) = 1, M(V ) = M(V ) T v T i M(V )v j ≥ 0, ∀i ≤ j, 5: if α(V ) ≥ 1 then 6:
add U i * and U j * to V ; Partition refinement 7:

end if 8: end while Remark 2: The conditions presented in this paper as well as the Algorithm 1 can be adapted to determine, besides the property of stabilizability, bounds on the exponential convergence rate λ , sometimes referred to as stabilization radius [START_REF] Jungers | On feedback stabilization of linear switched systems via switching signal control[END_REF]. It is sufficient, indeed, to multiply all the matrices A convex problem, in LMI form, has to be solved at every iteration. On the other hand, there is no guarantee, in general, of finite termination of Algorithm 1, as the existence of a conic partition cannot be ensured unless the system is stabilizable. This is reasonable since the problem is undecidable by its nature. Moreover, the number of elements in an eventual conic partition satisfying the necessary and sufficient condition could be arbitrarily big. Nonetheless, as illustrated in the next section, the application of Algorithm 1 to examples from the literature and to a four dimensional system shows the efficiency of the algorithm in testing stabilizability and computing tight bounds on the exponential convergence rate.

VI. NUMERICAL EXAMPLES

Three examples will be considered to illustrate the results and the application of the algorithm. The set of matrices M(V ), determining the conic cover, and the parameters α(V ) and µ(V ) resulting from the application of Algorithm 1 for each example are available at [5], in form of matlab data files, together with simple scripts to verify the conic covering of the space and the stabilizability condition (20) satisfaction.

1: Consider Example 17 in [START_REF] Fiacchini | On the stabilizability of discrete-time switched linear systems: novel conditions and comparisons[END_REF], employed to prove that condition (4) is just sufficient for stabilizability. The system is given by three modes with matrices

A 1 = AR(0), A 2 = AR 2π 3 , A 3 = AR -2π 3
,
where

A = 0.6 0 0 0.6 -1 , R(θ ) = cos(θ ) -sin(θ ) sin(θ ) cos(θ ) .
The system is proved in [START_REF] Fiacchini | On the stabilizability of discrete-time switched linear systems: novel conditions and comparisons[END_REF] to be stabilizable since the necessary and sufficient condition (3) is satisfied with N = 1.

Fig. 2. Set B, black, ellipsoidal preimages B i for i ∈ I [9] , in red, and basis vectors of {V ( j) } j∈N 23 , in black lines.

Then, the Euclidean norm decreases at every step, but (4) does not hold, for any N ∈ N. This implies also that Lyapunov-Metzler conditions cannot validate stabilizability either. Algorithm 1 with N = 1 and the identity as initial basis matrix, i.e. B = I, has been applied to this example, validating the stabilizability with a final cover of 8 cones {V ( j) } j∈N 8 , depicted in Figure 1. Figure 1 represents also the sets B and its preimages B 1 , B 2 and B 3 that graphically validate condition (3) with N = 1. Note that from Proposition 1, it is sufficient to generate a cover of half of the space R 2 .

Example 2: The second switched system, with modes

A 1 = R π 4 , A 2 = 2 0 0 1/2 (22)
is based on an example from [START_REF] Stanford | Some convergence properties of matrix sets[END_REF] and used in [START_REF] Jungers | On feedback stabilization of linear switched systems via switching signal control[END_REF] to illustrate methods for computing bounds of the stabilizability radius. These methods, based on conditions evaluated on a grid of the unit circle in R 2 , hence on a one dimensional space, provides two bound estimations, of 0.886 and 0.88.

To infer the value of the convergence rate λ , consider the system whose matrices are given by Ā1 = λ -1 A 1 and Ā2 = λ -1 A 2 . For N = 9 and λ = 0.8855 the set B and the preimages B i for all I [9] , whose number is N = 1022, are drawn, see Figure 2. From geometric inspection it can be proved that B is contained in the union of ellipsoids and λ = 0.8855 seems a rather tight bound for N = 9, see the Figure 3 Algorithm 1 is then applied with B = I, validating the condition of stabilizability with convergence rate of λ = 0.8855 by testing the convex condition over 23 cones covering the half space, also depicted in Figure 2.

Comparing this result with those obtained in [START_REF] Jungers | On feedback stabilization of linear switched systems via switching signal control[END_REF] it is worth noting that, besides providing a tighter bound with respect to one proposed, the method presented here provides guarantee of stabilizability, that cannot be ensured in general when evaluating conditions on a grid of points. Another benefit is the fact that, as the proposed method does not resort to a grid Fig. 3. Zoom of set B, black, ellipsoidal preimages B i for i ∈ I [9] , in red, and basis vectors of {V ( j) } j∈N 23 .

of the space, it can be applied to higher dimensional systems, as illustrated in the next example. For this, the algorithm has also been applied to this example with N = 5, leading to the validation of a bound of λ = 0.96 with 10 cones. Example 3: This example is specifically built such that a tight estimate of the convergence rate is known but other methods for its sharp estimation, based on the space gridding, might be not applicable due to the space dimension. To illustrate the possibility of applying the proposed method in higher dimension, in fact, a system is considered that is composed by two subsystems equal to the one of Example 2. Any method based on gridding the space is hardly applicable in R 4 . Moreover, an accurate estimation of the convergence rate is known by construction, and then the precision of the estimation obtained with Algorithm 1 can be objectively evaluated. Consider the four dimensional system with 4 modes

Â1 = λ -1 A 1 0 0 δ I , Â2 = λ -1 A 2 0 0 δ I , Â3 = δ I 0 0 λ -1 A 1 , Â4 = δ I 0 0 λ -1 A 2 ,
with λ the known estimation of the convergence rate, see below, and A 1 and A 2 as in (22). Parameter δ is posed equal to 0.999, to guarantee a decreasing of the function x T x within the horizon N of the four dimensional system. The objective is to check whether Algorithm 1 is able to validate the convergence rate for the four dimensional system. Since N = 9, with q = 4, leads to more than a half million of sequences, the smaller horizon of N = 5 has been employed. Note that N = 5 gives N = 1364 switching elements in I [5] . Since the bound of λ = 0.96 has been validated in the previous example for the single subsystem, it must hold also for the four dimensional system. The objective is to test whether the tight bound λ = 0.96 for N = 5 can be validated also in this case.

Indeed, the application of Algorithm 1, with a randomly generated non-singular matrix B, provides a validation of the bound λ = 0.96 for the four dimensional system, the semidefinite conditions [START_REF] Sun | Stabilizability and insensitivity of switched linear systems[END_REF] holding in each of the 1827 cones composing the generated cover of R 4 .

VII. CONCLUSIONS

In this paper a novel necessary and sufficient condition for stabilizability has been proposed. The condition consists in the existence of an appropriate conic cover of the space and a convex condition on its elements. This condition led to a procedure for iteratively generating a sequence of conic covers and linear problems allowing to test the stabilizability of systems for which no periodic stabilizability can hold. The method has been also applied to estimate the convergence rate for switched linear systems, even in dimension four.

Fig. 1 .

 1 Fig. 1. Set B, black, ellipsoidal preimages B 1 , B 2 and B 3 , in red, and basis vectors of {V ( j) } j∈N 8 , in thin lines.


Wei-Ting Yeh

Misaki Ozawa

Kunimasa Miyazaki

Takeshi Kawasaki

Ludovic Berthier

Glass stability changes the nature of yielding under oscillatory shear

We perform molecular dynamics simulations to investigate the effect of a glass preparation on its yielding transition under oscillatory shear. We use swap Monte Carlo to investigate a broad range of glass stabilities from poorly annealed to highly stable systems. We observe a qualitative change in the nature of yielding, which evolves from ductile to brittle as glass stability increases. Our results disentangle the relative role of mechanical and thermal annealing on the mechanical properties of amorphous solids, which is relevant for various experimental situations from the rheology of soft materials to fatigue failure in metallic glasses.

Amorphous solids are generated by decreasing the temperature or by increasing the pressure of supercooled liquids, colloidal suspensions, and granular materials [1][2][3]. When deformed, their mechanical response is initially elastic but plastic deformations appear at larger strains, which eventually trigger material flow. This elastic-toplastic transformation is a yielding transition, which is currently an active research topic [4,5]. Understanding the microscopic nature of yielding not only helps in understanding the physics of glasses [6,7], but also guides the design of amorphous materials for industrial applications [8,9].

Yielding in amorphous solids is broadly classified into two classes: ductile and brittle yielding. Ductile materials accumulate plastic activity before yielding and show a smooth crossover from elastic-to-plastic regimes. Many soft glassy materials, such as foams [10], emulsions [11], and colloidal glasses [12] belong to this class. In brittle materials, yielding is characterized by a sharp stress overshoot accompanied by strain localization, which triggers material failure at large deformation. Metallic glasses [9,13], window glasses, and concretes [14] are typical brittle materials.

Brittleness is not necessarily an intrinsic material property. For example, brittleness can gradually increase by changing the cooling rate to prepare the system [15][16][17][18][19][20] or by physical aging [21][22][23][24][25][26]. In athermal conditions, this evolution was recently described as a phase transition and confirmed by atomic glass simulations under uniform shear [27]. Evidence was provided that a critical point separates ductile and brittle behaviors, demonstrating the importance of glass stability to understand yielding.

In many studies of the yielding transition, a large amplitude oscillatory shear (LAOS) protocol is used instead of a uniform shear flow. This is useful to relate macroscopic properties with microscopic trajectories of the constituting atoms [28]. LAOS is also the first step toward understanding the response of amorphous solids under dynamical loading conditions [29], relevant to the rheology of soft materials like emulsions and colloids [30,31], and the mechanism of fatigue failure in metallic glasses [32,33]. Previous numerical work under oscillatory shear found that yielding in athermal conditions is associated with a microscopic irreversible transition at the particle level [34][35][36][37]. It was also shown that small strain amplitude oscillatory shear can increase glass stability [37][38][39][40][41]. Such "mechanical annealing" [42] was proposed as a route distinct from thermal annealing to prepare stable glasses [43]. When the strain amplitude is increased above yielding, the system cannot find a stable energy state, and the plastic activity associated with energy dissipation is observed at each cycle [34,37,38]. Large-scale simulations [37,39,44] and a mesoscopic model [45,46] in this regime showed that macroscopic shear bands can form at long times, via a mechanism proposed to be similar to shear band formation under uniform shear at large deformation [47][48][49]. However, the role of glass preparation on the nature of yielding under oscillatory shear has not been tested over a broad range of glass stabilities. Therefore, the relative role of thermal and mechanical annealings on mechanical properties remains unclear [34,[37][38][39][40][42][43][44].

Our goal is to unify the yielding transition of thermally and mechanically annealed systems and to study how yielding changes when glass preparation is varied over a broad enough range to mimic the physics of both colloidal and metallic glasses under oscillatory shear. We use atomistic computer simulations of a simple glass former over a wide range of preparation temperatures [50]. We carefully disentangle the relative roles of mechanical and thermal annealings, and the various physical regimes accessible in the (strain amplitude, preparation temperature) phase diagram, regarding both the nature of the yielding transition and the distinct mechanisms of shear band formation in different glassy materials.

We simulate a three-dimensional size polydisperse system [50] with a pairwise soft potential given by

U (r ij ) = ǫ 0 d ij r ij 12 +c 0 +c 1 r ij d ij 2 +c 2 r ij d ij 4 , (1)
where r ij = |r ir j | is the distance between particles i and j, and

d ij ≡ (d i + d j )(1 -0.2|d i -d j |
)/2 controls the nonadditive interaction between particles of diameters d i and d j . The cutoff distance of the potential is set to r (ij) cut = 1.25d ij , and the constants c 0 , c 1 , and c 2 are chosen so that the first and second derivatives of the potential vanish at the cutoff. We, respectively, use the average diameter d and energy scale ǫ 0 as our length and energy units. To generate amorphous solids with a broad range of stabilities, we first produce equilibrium configurations of the supercooled liquid over a wide range of temperatures T init , using the swap Monte Carlo method [50]. We study the range T init ∈ [0.062, 5] and work at fixed number density ρ = 1.02. These equilibrium configurations are then instantly quenched to zero temperature (using the FIRE algorithm [51]) to form the initial glass configurations to be sheared. The model and physical properties of these glasses were documented in Refs. [27,50,52,53].

In particular, the energy of the inherent structures, E IS , is almost constant at high T init , and starts decreasing as T init drops below the onset temperature T onset ≈ 0.18. Several physical properties change qualitatively when T init decreases below the mode-coupling crossover temperature, T MCT ≈ 0.108 [53]. The estimated experimental glass temperature is T g ≈ 0.072. We deform these zero-temperature amorphous solids using a simple oscillatory shear with a strain amplitude γ 0 at constant volume. We employ both the athermal quasistatic (AQS) shear protocol [54] and the finite (but small) strain rate protocol. A homogeneous shear strain is realized using Lees-Edwards boundary conditions [55]. In the AQS protocol, the system is strained at each step by a small incremental strain |∆γ| = 2×10 -4 , after which the configuration is again relaxed into the nearest energy minimum. In the finite strain rate protocol, the system is driven by an overdamped dynamical equation [34]

ξ dr i dt -γ(t)r i (t) • ŷx + j =i ∂U (r ij ) ∂r i = 0, (2)
where x and ŷ are the unit vectors in the x and y directions, r i is the position of particle i, and ξ is the friction coefficient that sets the time unit dξ/ǫ 0 . The applied strain is γ(t) = γ 0 sin(2πt/T cyc ), with an oscillation period T cyc . Most reported results are for the finite strain rate protocol with a small strain rate γ0 ≡ 2πγ 0 /T cyc = 6.2832 × 10 -4 , and total number of particles N = 12000. In the Supplemental Material, we show that this finite strain rate is small enough to reproduce results obtained from the AQS limit [56]. Where needed, we also discuss the results of AQS simulations and larger systems N = 48000. A complete system size dependence is presented in the Supplemental Material [56]. For each data point, we average over at least 50 cycles for at least three different initial conditions after the system reaches a steady state. Unless specified otherwise, error bars in each figure represent the standard deviation of the measurements.

We first present the behavior of macroscopic rheological quantities in steady state after many cycles. We calculate the linear description of the nonlinear stressstrain relation [28] by fitting the time-dependent shear stress σ(t) (see the Supplemental Material for the definition [56]) to σ 0 sin(2πt/T cyc + δ 0 ), where σ 0 is the stress amplitude and δ 0 is the phase lag. We then obtain the storage modulus as G ′ = σ 0 cos(δ 0 )/γ 0 . The γ 0 dependence of these quantities is shown in Figs. 1(a-c) for several preparation temperatures T init . Data for the less annealed systems, 0.12 ≤ T init ≤ 5, whose stabilities would correspond to typical soft materials [27], essentially collapse on curves that are continuous with a cusp at the same yield strain amplitude γ Y ≈ 0.08 irrespective of T init ; In Fig. 1(a) σ 0 displays a smooth overshoot at γ Y , consistent with Refs. [34,37,38]; δ 0 starts to increase sharply in Fig. 1(b), but almost continuously at γ Y . In Fig. 1(c) G ′ decreases with a kink at γ Y , as also shown in Ref. [34]. A qualitatively distinct behavior is observed when T init is below a critical value T c ≈ 0.1. The data at small γ 0 no longer collapse and are distinct for each T init . One observes that more stable glasses have larger shear and storage moduli, and are less dissipative than poorly annealed ones. At the yielding transition point, all data display a sharp discontinuity, which depends on the glass stability. To our knowledge, such large discontinuities were not observed before in computer simulations under oscillatory shear, because it was impossible to produce stable enough glasses [37][38][39][40][41]. Beyond yielding, all curves merge again, as the plastic activity drives the glass away from its stable initial conditions, and memory is eventually lost [37,38]. Interestingly, the number of cycles required to reach the steady state remains finite below T c in the vicinity of γ Y , whereas it diverges for poorly annealed systems [34,35,37] (see the Supplemental Material [56]). In Fig. 1(d), we show the magnitude of the jumps of the various rheological quantities at γ Y as a function of T init (See the Supplemental Material for the precise definitions [56]). This representation eloquently demonstrates that the nature of yielding changes qualitatively at the critical value of T c ≈ 0.1. This value is close to the critical temperature discussed in the context of uniform shear [27], but more precise measurements and a finite size scaling analysis would be needed to establish a more precise connection.

To analyze the relative role of thermal and mechanical annealings, we study the dependence of the steadystate inherent structure energy, E IS , on T init and γ 0 , see Fig. 2(a). We again find two qualitatively distinct behaviors, depending on the value of T init before yielding. For γ 0 = 0, E IS decreases with decreasing T init , reflecting the increasing stability of the initial conditions by thermal annealing [53]. For T init > T c , E IS is a decreasing function of γ 0 , which confirms that poorly annealed glasses can access deeper energy states when submitted to oscillatory shear of modest amplitude. This is mechanical annealing (MA), as reported in previous works [42,43]. By contrast, for T init < T c [START_REF]We have not tried to mix thermal and mechanical annealings to access glassy states with lower energies[END_REF], mechanical annealing is not observed. This implies that thermal annealing (TA) is so efficient for these glasses that mechanical annealing is no longer able to drive them toward lower energy states. The recent results of Ref. [41] can be interpreted as representative of glasses aged near or slightly below T c .

Driving poorly annealed glasses above yielding produces instead higher energy configurations, and the yielding transition for those materials appears therefore as a cusp in E IS , whereas it again appears as a sharp discontinuity when T init < T c , which emerges at T c . We conclude that mechanical annealing operates for small γ 0 and high enough T init , but becomes inefficient when the effect of thermal annealing increases at low T init . The emergence of the critical temperature T c can be physically understood as a direct consequence of this competition. It corresponds to the temperature below which mechanical annealing is no longer useful. Since T c corresponds to a sharp change in the nature of yielding, we see no conceptual reason to relate it to the mode-coupling temperature, which describes a smooth physical crossover in finite dimensions and equilibrium conditions [START_REF] Bhaumik | [END_REF]. Physically, lowering T c decreases the quenched disorder and presumably the density of weak regions in the initial glass configurations [27].

These findings can be summarized in the nonequilibrium phase diagram in the (γ 0 , T -1 init ) plane shown in Fig. 2(b). The yield strain amplitude γ Y separates the "elastic" region at low strain amplitude where the plastic activity is suppressed in steady state, from the "yielded" region where plastic events occur during each cycle at large strain amplitude. γ Y is constant in the regime dominated by mechanical annealing (MA regime) for T init > T c [38], but it increases with decreasing T init in the regime dominated by thermal annealing (TA regime). The value of T c can be determined by a power-law fit of γ Y [56], which gives T c ≈ 0.101. Although the nature of yielding changes dramatically with T init , it is accompanied for all stabilities by a similar discontinuous irre- versible transition at the particle scale (see, for example, Fig. S9 in the Supplemental Material [56]).

The existence of two distinct regimes for yielding under oscillatory shear also influences the nature of shear banding above yielding. Recent works on the poorly annealed glasses in the AQS limit and small strain rate suggested that persistent shear bands may emerge under slow oscillatory shear for large systems at long times [37,39,44]. We confirm that a persistent shear band forms in our simulations at long times in similar conditions, but it does not appear when the strain rate is large or γ 0 is away from γ Y [56]. In Fig. 3(a), we illustrate the slow emergence of a shear band after many cycles using AQS simulations and larger system size, N = 48000, starting from a poorly annealed glass with T init = 0.5. In steady shear, the stress-strain rate relation is monotonic and no shear band is formed [START_REF] Singh | [END_REF] and thus the observed shear band purely originates from the oscillatory nature of the shearing. A color map of the one-cycle nonaffine deformation d na,i for each particle i [60] is used to visualize strain localization after many cycles. To quantify the gradual formation of the shear band, we introduce an order parameter |H|, which quantifies the extent of spatial inhomogeneities of d na,i . See the Supplemental Material for the precise definition of d na,i and |H| [56]. In Fig. 3(c), we report the corresponding slow increase of |H|. The finite value of |H| at steady state implies the presence of a persistent shear band.

Shear band formation is very different for T init < T c , as demonstrated in Fig. 3(b). In the first few cycles, a very small amount of localized plastic activity is observed. This is followed by a sudden, macroscopic failure of the material accompanied by the instantaneous formation of a shear band at n = 29 cycles. The shear band is formed within a single cycle, which is very different from the gradual emergence reported for T init > T c at long times [37,39,44]. The observed behavior is reminiscent of fatigue failure in metallic glasses under cyclic loading in the sense that, after some loading cycles, the system fails abruptly [61]. The time dependence of |H| in Fig. 3(c) confirms the sudden formation of the band, which gradually broadens as the number of cycles increases. At very long times, the memory of the initial condition is lost and both types of systems display similar dynamics.

In this work, we systematically analyzed the mechanical properties under oscillatory shear of amorphous solids prepared at different depth in the energy landscape and the interplay between thermal and mechanical annealings. We have identified a critical temperature T c separating two distinct regimes for yielding. For systems prepared at T init > T c in the MA regime (corresponding to typical soft materials), mechanical annealing can drive the system toward lower energy states at small strain amplitude, but memory is gradually lost at large strain amplitude. In the TA regime, T init < T c (corresponding to typical metallic glasses), no mechanical annealing occurs at small strain amplitudes, and yielding corresponds to a macroscopic failure of the material accompanied by a discontinuous jump of macroscopic rheological quantities and energy, associated with the sudden appearance of a shear band. Therefore, the MA and TA regimes are separated by a ductile-to-brittle transition. The present work dealing with oscillatory shear extends qualitatively similar findings obtained for uniform shear in AQS conditions [27] to time-dependent flows.

We observed distinct mechanisms of shear band formation in the MA and TA regimes. A shear band forms gradually at long times in the MA regime, becoming persistent in the steady state after yielding. In the TA regime, however, an abrupt emergence of a shear band in a single cycle is found. In particular, before shear band formation, a small number of localized precursor events are observed (see Fig. 3(b)), which would initiate a macroscopic shear band [62]. Future work should address in more detail the distinct mechanisms for shear banding in these two regimes. Some preliminary results on this aspect are shown in the Supplemental Material [56].

We investigated the important role of glass stability on the yielding under oscillatory shear relevant for various experimental situations from the rheology of soft materials to fatigue failure in metallic glasses. Our work elucidates the interplay between mechanical and thermal annealing processes, providing useful information for material design in industrial applications [63], as well as theoretical descriptions of amorphous materials under cyclic deformation [45,46,64].

  . (a)-(c) Various rheological quantities as a function of the strain amplitude γ0 for different Tinit's [see the legend shown in (a)]. (a) Stress amplitude σ0. (b) Phase lag δ0. (c) Storage modulus G ′ . Dashed vertical lines in these figures represent the discontinuous jumps of these quantities. (d) The magnitude of the discontinuities of these quantities at yielding as a function of 1/Tinit vanishes at Tc ≈ 0.1.

FIG. 2 .

 2 FIG. 2. (a) Steady-state inherent structure energy EIS as a function of the strain amplitude γ0 for several Tinit. The dashed vertical lines represent the discontinuous jumps. (b) Nonequilibrium phase diagram in the (γ0, T -1 init ) plane. The critical temperature Tc (red long dashed line) divides the phase space into the mechanical annealing (MA) and thermal annealing (TA) regimes. The black solid line separates the elastic and yielded regions (see the text).

FIG. 3 .

 3 FIG. 3. The snapshots show a color map of the one-cycle nonaffine deformation dna,i for (a) Tinit = 0.5 and (b) 0.075 after n = 10, 30, and 200 cycles (from top to bottom). (c) Evolution of the shear band order parameter |H|. The arrows indicate the times shown in (a, b). Yielding in (a) is accompanied by the slow emergence of a diffusive shear band, whereas in (b) a sharp shear band is suddenly formed when the material fails macroscopically. Data obtained using the AQS protocol and N = 48000 with γ0 = 0.11.

We thank Srikanth Sastry for discussions. We also thank the anonymous referee for bringing our attention to Ref. [41]. This research is supported by Japan Society for the Promotion of Science (JSPS) KAKENHI (Nos. 16H06018, 16H04034, 18H01188, 19K03767, and 19H01812) and a grant from the Simons Foundation (No. 454933, L. B.).