Tom Vincent-Dospital
email: vincentdospitalt@unistra.fr

Renaud Toussaint
email: renaud.toussaint@unistra.fr

Alain Cochard

Eirik G Flekkøy

Knut Jørgen Måløy

Thermal dissipation as both the strength and weakness of matter. A material failure prediction by monitoring creep

Keywords:

In any domain involving some stressed solids, that is, from seismology to general engineering, the strength of matter is a paramount feature to understand. We here discuss the ability of a simple thermally activated sub-critical model, which includes the auto-induced thermal evolution of cracks tips, to predict the catastrophic failure of a vast range of materials. It is in particular shown that the intrinsic surface energy barrier, for breaking the atomic bonds of many solids, can be easily deduced from the slow creeping dynamics of a crack. This intrinsic barrier is however higher than the macroscopic load threshold at which brittle matter brutally fails, possibly as a result of thermal activation and of a thermal weakening mechanism. We propose a novel method to compute this macroscopic energy release rate of rupture, Ga, solely from monitoring slow creep, and we show that this reproduces the experimental values within 50% accuracy over twenty different materials, and over more than four decades of fracture energy.

I.

Although seminal, the early theoretical descriptions of crack dynamics, such as Griffith's [1] or Slepyan's [2,3], was somewhat binary: beyond a critical mechanical load, matter suddenly breaks. It is however acknowledged that, at load levels below the critical one, a far slower crack propagation already occurs, that will here be referred to as 'creep'. This phenomenon was successfully modelled with Arrhenius-like sub-critical growth laws [4,5], and is hence sometimes called 'stress corrosion'. With the increasing number of experimental work, the description of such a slow dynamics was quickly refined, and five propagation stages were notably distinguished [5]. Let us start this manuscript by summarising them. Figure 1 illustrates these stages in a V -G plot, where V is the crack velocity for a given load G, which is the 'energy release rate', that is, the energy that the fracture consumes to advance by unit surface [1]. At stage 0, while under only a mild mechanical input, cracks do not actually propagate forward. This was notably explained by the existence of some healing processes, that there efficiently compete with the failure ones [4]. From this state, when the load is increased above a given threshold, some slow fracture growth starts to be observed (stage I). The propagation velocity V increases exponentially with the crack's energy release rate G. In a sub-critical (i.e., Arrhenius-like) description, it implies that V is to first order explained by an activation mechanism dependent on G, in a chemical-like rupture reaction [6]. Logically, this regime was observed to also depend on the surrounding temperature and on the fluid which is present in the fracture [7], which affects the chemical reaction involved in molecular bond break-ing. When reaching a faster propagation, some velocity plateau might then hold (stage II), possibly as the transport of fluid corrosive elements toward the tip cannot efficiently cope with the crack advance. Such plateau is, in this case, only a transition to a sub-critical growth 'invacuum-condition', where the dynamics becomes notably insensitive to the fracture fluid (stage III). Finally, when a particular threshold is reached for the energy release rate, the velocity jumps to a far quicker regime: the material fails (stage IV). We will denote 1 this threshold G a in J m -2 , with 'a' standing for 'avalanche'. In this work, we will show how studying the slow creep regime allows to predict this particular failure load. This can lead to methods to characterise natural or lowly controlled materials, where the critical energy release rate G a is not well known a priori, but where the monitoring of creep allows to infer it. In a previous study [8], we indeed proposed a unifying model of the slow creep and the fast regime, holding a precise quantification of the energy budget and the heating of the crack tip, which is coupled with an Arrhenius-type activation law. We have shown how it accounts, in some polymers [9], for seven decades of propagation velocities and for the transition, at the avalanche load, from creep to sudden failure. Here, we present how well this thermodynamics based model can predict the threshold G a for a broad range of materials, by comparing its forecasts to actual experimental failure thresholds from twenty data sets from the literature. By FIG. 1. Summary of the different forward crack velocity regions observed in experimental velocity curves. After Fracture of Brittle Solids, Lawn [5] .

doing so, one can actually identify the microscopic rupture energy of the breaking bonds, G c , and show how this quantity is related to, yet different from, the macroscopic G a . The agreement between the predictions and the realisation is obtained for materials spanning more than 4 orders of magnitude in energy release rate, indicating the robustness of this description among different types of materials and the versatility of the theoretical framework.

II. THE THERMAL WEAKENING MODEL

We consider that the propagation of a crack follows an Arrhenius sub-critical growth law, in which the temperature term accounts for the induced heat generated at the plastic crack tip [10,11]. Such a model, introduced in Refs. [8] and [9], writes as

V = V 0 min exp - d 3 0 (G c -G) 2lk B (T 0 + ∆T ) , 1 (1)
∂(∆T ) ∂t = λ C ∇ 2 (∆T ) + φGV Cπl 2 f, (2)
where the first equation describes the Arrhenius growth (i.e., the term in brackets is a probability for the thermal bath to overcome an energy barrier), and the second one in the diffusion equation governing the thermal evolution around the crack front. The heat conductivity and volumetric heat capacity of the solid matrix are respectively denoted λ and C. V 0 is a nominal atomic speed related to the collision frequency in the thermal bath, and should typically be comparable to the mechanical wave velocity of the studied media [2,12]. The activation energy is modelled proportional to (G c -G), where G c is the surface energy barrier to overcome in order to break atomic bonds. d 3 0 is the characteristic volume for the bonds (d 0 ∼ 1 Å), k B is the Boltzmann constant, T 0 the ambient temperature and ∆T any variation away from it at the crack tip. A percentage φ of the power consumed per unit of crack length GV is uniformly dissipated as heat over a zone of support function f and of radius l. This heating zone is a subset of the process zone (that is, the full extent of plasticity around the tip), and we assume that it also constrains the stress level σ at the tip, as verified in Vincent-Dospital et al. [9]: σ ∼ GE/l (e.g., see Ref. [5]), where E is the materials Young's modulus. This assumption is the reason why l also intervenes in Eq. ( 1), where the elastic energy stored in the rupturing link typically writes as d 3 0 σ 2 /(2E), and is thus similar to d 3 0 G/(2l). In this framework, an acceleration of a fracture is thus both related to an increase of stress at the tip, which reduces the rupture energy barrier d 3 0 (G c -G)/(2l), and to a related increase in volumetric internal energy (C∆T ) at the fracture's head. Note that it was shown [9,13] that at low velocities (i.e., the creep velocities we are interested in), ∆T computed from Eq. ( 2) can, more simply, be approximated to

∆T ∼ (φGV h)(l/V ) C(πδ 2 h) = φGV λ , (3)
where φGV h is the thermal power deposited along a portion of length h of the crack front, l/V is the characteristic time for the front heating (i.e., the time a given location of the material stays in the advancing heating zone) and δ ∼ λl/(πCV ) is the skin depth of heat diffusion upon the same time. According to Eq. ( 3), ∆T does not depend on C or l, notably because if the crack advances slowly enough, the temperature elevation is constrained by the heat diffusion skin depth δ rather than by the size of the heat production zone, as the former is, in this case, big compared to the latter. Approximating Eqs. ( 1) and ( 2) by their steady state solutions, two stable propagation branches are derived from this model [8], as shown in Fig. 2: a fast phase, which is obtained for a hot crack tip and corresponds to the catastrophic failure of matter, and a slow one corresponding to the creep regime, when ∆T T 0 . In between these two branches, a hysteresis situation holds with a third unstable phase. In this study, we are mainly interested in the slow to fast regime transition (i.e., that leads to quick material failure). When approaching this transition, the velocity deviates from its negligible heating asymptotic expression, which is a simple exponential increase with the load G:

ln V V 0 ∼ (G -G c ) d 3 0 2lk B T 0 , (4)
as the rise in temperature ∆T in Eq. ( 1) becomes comparable to the room temperature T 0 . The particular energy release rate G a is then reached, at which ∂V /∂G → +∞, and beyond which the crack can only avalanche to a velocity which is orders of magnitude higher (see Fig. 2). Matter suddenly breaks. As a result of thermal activation, G a is actually less than the actual surface energy barrier for breaking bonds G c . Although rarely regarded today, such an importance of the auto-induced heat to explain brittleness was early developed [14][15][16]. These studies reckon that the dissipated energy favours failure by locally softening the material at the tip (i.e., a decrease of its elastic moduli with temperature), in particular in the case of soft rubber-like materials. In the formalism of Eq. ( 1), this view would correspond to a brutal reduction of the energy barrier d 3 0 (G c -G)/(2l) at the onset of the instability, rather to (or maybe in addition to) a dramatic increase in temperature, due to chemicophysical phase changes of the matrix. Our model neglects such a softening effect and instead considers that the reaction rate for rupture is increased from the elevated temperature, only as understood by statistical physics. Of course, both views are not mutually exclusive, and one of the two mechanisms may prevail depending on the considered material, or depending on the crack velocity. Indeed, for a softening effect to be at stake, the typical time for a material temperature-related change in phase must be less than the typical warming time l/V of a given location of the crack trajectory. By contrast, a change in velocity, as understood by statistical physics, should be as quick as a few atomic vibrations [6] (∼ 10 13 Hz) and, hence, should be more likely to explain fast cracks in hard solids. Such a qualitative discussion underlines the importance of knowing the actual rheology close to crack tips in given materials, and its sensibility to eventual temperature bursts. Note that in both approaches of thermal weakening here discussed, the G value of interest (i.e., G a ) remains similar in its concept: the threshold for which ∆T is high enough so that a quick avalanche can be generated by a given mechanism. In the next sections and the rest of this manuscript, we will focus on the thermally activated, statistical, model that we have developed above, where the accurate rheology at the tip is actually considered to be of second order on the crack dynamics, and only the increase of ∆T is considered.

III. MODEL PREDICTIONS VERSUS REPORTED FAILURES

Extensive fracturing experiments on numerous materials can be found in the literature. Hence, we can compare the model predictions of G a to some experimentally reported avalanche thresholds, that are often referred to as 'critical energy release rate' or 'material toughness', although it does not correspond to what is here denoted G c , which is an intrinsic (microscopic) medium property not directly measurable at lab scale. In our framework, G c is called the critical energy release rate because it differentiate between an actual subcritical propagation (at G < G c ), where thermal activation is required for the crack to advance, and an overcrital propagation (at G > G c ), where enough mechanical energy is given to the crack so that thermal activation is not strictly needed. In this latter case, Eq. ( 1) simplifies to V = V 0 . A fast 'dynamical' crack can be subcritical (see Fig. 2), in particular when G a < G < G c , and we thus refrained to call G a a critical energy release rate. Note that Eq. ( 1) does not account for all of the creep regimes summarized in Fig. 1, that one can meet with an experimental test, but displays a unique low velocity slope (i.e., from Eq. ( 4)). We have indeed discarded any healing processes, needed to explain stage 0, as they are beyond the topic of the current study. Such processes can however be included in the model [9]. We have also assumed no rate-limiting environmental factor, that is, no significant chemical interaction of the matrix with the fracture fluid (i.e., no stage I or II). We have hence restricted our comparison to experimental data to such a non-environmental creep case (stage III), although distinguishing it with certitude is not always straightforward. When available, we have notably preferred data sets of in-vacuum or dry (in air) experiments or with lowly corrosive (e.g., neutral pH) fracture fluids (see the creep plots in the ESI † for an exhaustive list of the considered experiments, and their environmental conditions). Indeed, high fluid pH or moisture tend to increase the creep dynamics of given materials (e.g., [17]), presumably by altering the chemistry of the rupture process at the tip (i.e., stages I and II), and such effects are not included in our model. Note however that, when some fluid-matrix interaction does take place, the model could still be somewhat applied, if failure is preceded by a unique slope (i.e., if it occurs before the slope break between stages I and II), or after it, once clearly having entered in regime III. In the former case, the definition of the surface energy barrier G c may slightly change: from an intrinsic strength of the solid to an equivalent (lower) strength under a given chemical environment. To predict G a , it is of course needed to know, for each material, the values of the model constitutive parameters, that is, the parameters of Eqs. ( 1) and ( 2), which describe the evolution of the crack velocity as a function of the applied energy release rate. Although they are not many, most of these parameters are not usually considered, and are hence unknown. It is however possible to estimate their order of magnitude from known material properties, or to assess them from the slow (creep) part of the loading curve. We have first considered that V 0 is of the order the mechanical wave velocity. It could ideally be that of the Rayleigh waves [12], but it is often simpler to instead estimate the -similar-shear wave velocity of solids, V S ∼ µ/ρ, as the shear modulus µ and the density ρ of most materials are easily available. The heat conductivity λ is also known in most cases, and T 0 is nothing but the room temperature at which a given reported experiment took place. We assume the interatomic space d 0 to be 1 Å. While it could be two or three times bigger depending on the materials, which would have an order of magnitude effect on the term d 3 0 , this uncertainty would only impact the estimation of l, as the ratio d 3 0 /l is here of importance. We indeed have to deduce l and G c from the slope and intercept of the slow sub-critical growth, that is, from the two terms of Eq. ( 4) fitted to the experimental curves with the fit parameters V null and b: ln(V /V null ) = bG, V null being the velocity a crack would have at a null load (G = 0), in the absence of fluid corrosion or healing effects and b the slope of the creep curve of dimension m 3 s -1 J -1 . This gives

l = d 3 0 2bk B T 0 (5)
for the size of the thermal zone, and, for the intrinsic strength of the solid:

G c = 2lk B T 0 d 3 0 ln V 0 V null . ( 6
)
This implies that we can predict G a if relying on some creep observations, that can yet be at loads far below the failure threshold. The only remaining model parameter, the percentage φ of energy converted into heat is mostly unknown. While qualitative statements, such as larger φ in metals rather than, say, polymers, are tempting, we have here arbitrarily fixed this percentage to 50 % in all materials, except for a couple of instances where we could estimate it [9,13]. For different materials, we however show how φ affects the prediction of our model in the ESI † . Note that, while the velocity is often reported in relation to the stress intensity factor K rather than the energy release rate G, we have here converted from one to the other with the following relation [5]: G ∼ K 2 /E to derive a and b, and then l and G c . Backwardly, with the here proposed method, we will thus predict the toughness, K avalanche ∼ √ EG a , based on the creep measurement. All the introduced parameters can now be estimated, and we did so for twenty materials for which the creeping behaviour was studied in the literature [9,[17][18][19][20][START_REF] Wiederhorn | Influence of water vapor on crack propagation in soda-lime glass[END_REF][START_REF] Wiederhorn | Crack growth in sapphire[END_REF][START_REF] Meredith | Fracture toughness and subcritical crack growth during high-temperature tensile deformation of westerly granite and black gab-bro[END_REF][START_REF] Nara | Study of subcritical crack growth in andesite using the double torsion test[END_REF][START_REF] Holder | Experimental determination of subcritical crack growth parameters in sedimentary rock[END_REF][START_REF] Wang | Subcritical crack growth induced by stress corrosion in hardened cement paste[END_REF][START_REF] Nara | Subcritical crack growth and long-term strength in rock and cementitious material[END_REF][START_REF] Oates | Subcritical crack growth in lead zirconate titanate[END_REF][START_REF] Chevalier | Subcritical crack propagation in 3Y-TZP ceramics: Static and cyclic fatigue[END_REF][START_REF] Evans | Crack propagation and failure prediction in silicon nitride at elevated temperatures[END_REF][START_REF] Hénaff | Creep and creepfatigue crack growth in aluminium alloys[END_REF][START_REF] Huang | Internal hydrogeninduced subcritical crack growth in austenitic stainless steels[END_REF][START_REF] Sastry | Subcritical crack-growth under sustained load in Ti-6AI-6V-2Sn[END_REF][START_REF] Nadeau | Subcritical crack growth in vitreous carbon at room temperature[END_REF]. The corresponding G to V curves and a table of all the inferred parameters values are shown in the ESI † . We can then solve numerically the two non linear equations ( 1) and ( 2), now taking into account the temperature rise ∆T . In other words, for any value of the energy release rate G, we can compute the possible crack propagation velocities according to the model, as for instance shown in Fig. 2, where, depending on G, one to three velocities are possible. The inflection of the obtained curve, where ∂V /∂G → +∞ (see Fig. 2), can be identified as G a and compared to the reported experimental thresholds. The detailed procedure is discussed and applied to two materials (PMMA and Pressure Sensitive Adhesive) in Ref. [9]. This comparison is summarised for all the media in Fig. 3, and our model displays there a good general description of catastrophic failure. In the same figure, the surface energy barrier G c is also displayed for comparison, as well as the relative error made in the estimation of G a .

IV. ANALYTICAL APPROXIMATION

While, to derive the modelled G a , one should compute the full crack dynamics (i.e., as displayed in Fig. 2), and search for the points where ∂G/∂V = 0, we explain in the ESI † how Eqs. ( 1) and ( 2) also approximately lead to

G a ∼ λT 0 φV 0 exp(R a ) R a , (7)
where R a is the activation energy at the avalanche threshold counted in thermal energy units: R a = d 3 0 (G c -G a )/(2lk B T 0 ). As this ratio notably depends on G a , Eq. ( 7) only implicitly defines the threshold. Although a numerical solver is there still required to compute the threshold, it is simpler and far quicker than finding the accurate solution, and potentially easy to use in potential engineering applications. We show, in the ESI † , how the approximation of Eq. ( 7) is a slight overestimation of the real solution, by about 0 to 10%. It however gives further insight on the influence of each parameter. The remaining dependence on (G c -G a ) shows that, to obtain thermal weakening, one must already be close to the actual microscopic energy barrier of rupture G c . The stress at the tip thus remains an important driving mechanism of the crack propagation.

The brutal failure threshold G a increases with the ambient temperature T 0 , as a higher temperature elevation ∆T is then required at the tip to significantly overcome a higher thermal bath. Similarly, G a increases with the λ/φ ratio, as a high λ/φ indicates a rather cool crack tip, either because little energy is converted into heat (i.e., a small φ) or because this energy is efficiently evacuated away from the crack (i.e., a high λ). The failure threshold also decreases with V 0 , as, at a given load, a higher V 0 means a faster (and then hotter) creeping crack from a higher atomic collision frequency V 0 /d 0 .

V. MICROSCOPIC VS MACROSCOPIC FRACTURE ENERGY, AND LOCAL VS BULK ENERGY DISSIPATION

In Fig. 3, one can notice that the surface energy barrier G c is always similar in order of magnitude to the rupture threshold G a . Yet, the rupture always occurs at a load less than G c , with G a being about twice lower in average for all the displayed solids. We have here explained how a weakening mechanism, as the thermal view that we have here developed, allows to account for this discrepancy. Having gathered various exponential creep data, and derived l and G c from their slope and intercept in their ln V -G representations, we can notably infer the intrinsic crack energy barrier in each material:

U c = d 3 0 G c /(2l) = k B T 0 ln(V 0 /V null ).
As shown in Fig. 4, this quantity is always in the order of 10 -19 J ∼ 1 eV, logically comparable to the energy level necessary to unbind single atomic covalent bounds [START_REF] Kaufman | Inorganic chemistry: principles of structure and reactivity, 4th ed[END_REF], which confirms the relevance of a simple thermally activated model for the description of stage III (in-vacuum-like) creep. The actual values of U c are yet often slightly smaller than the typical covalent strength. This could derive from an averaging effect. Indeed, the Arrhenius law of Eq. ( 1) is a mesoscopic statistical law, which we have fitted to some macroscopic measurements of crack propagation (i.e., the creep experiments). With cracks that are prone to follow the weakest paths, that possibly includes weak inter-molecular bonds (such as Van der Waals and hydrogen links) and dislocations or atomic voids (i.e., when the distance between two consecutive breaking bonds is more than a few ångströms), U c is likely representative of the average rupture energy in a disordered landscape. For instance, in polymers, part of the rupture shall be inter-molecular, and, in rock-type materials, the crack dynamics might benefit from the intrinsic porosity. However, due to the simplicity of the model, care should be taken when interpreting U c beyond its order of magnitude. It is clear however that the value of G c varies by a factor 10 4 for different materials, while its counterpart U c does not. As most materials have the same U c and d 0 , in order of magnitude, the large variability in G c (and hence in G a ) which is observed is, in this description, attributable entirely to the variability in the scale for the release of heat. We indeed infer that l varies from the radius of a single atom, for the weakest materials, up to 1 µm, for the ones with the highest G c (see Fig. 5). The wider the plastic area that shields crack tips, the stronger matter is. But backwardly, we have discussed how the heat dissipation might be the root cause for dramatic ruptures in brittle solids, if the heat is not efficiently evacuated away from the rupture front. Overall, Eq. ( 1) should be understood as:

V = V 0 exp - U c -U (Load [+]
, Thermal radius [-] k B T Dissipation [+] , Diffusion [-] ,

(8) where U is the mechanical energy corresponding to the stress actually transmitted to the crack tip covalent bond of average strength U c , and where [+] and [-] indicate whether the T and U functions are increasing or decreasing with the specified concepts (i.e., with the mechanical load, the core (thermal) dissipation radius l, the dissipation of heat or its diffusion). This equation emphasises that, in our formalism, thermal dissipation is both the strength of matter (from the shielding of the mechanical energy U actually transmitted to the crack tip) and the weakness of matter (from the amplification in internal energy around the tip). We can compare the values of l with the more typical plastic radius predicted by a Dugdale view [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] of the process zone, l macro ∼ G c E/σ y 2 , where σ y is the tensile yield stress, beyond which macroscopic samples lose their elasticity. As shown in Fig. 5, the latter is consistently five to seven orders of magnitude higher than what we predict for l. This likely translates to the fact that plasticity (here understood as the dissipation of mechanical energy in any form) ought to be a rather heterogeneous phenomenon, with a greater density of energy dissipation close to the front than away from it. A higher density of energy is indeed likely to be dissipated close to the tip singularity than at the edge of the l macro plastic radius. Thus, the scale of a process zone can be characterised either by its core radius l, where most of the heating due to the dissipation takes place, or by its full extent l macro , where the rheology becomes non elastic. While the former is to include significant thermal losses (quantified by φG), the latter can also encompass various other dissipation mechanisms for another portion of G, namely, the nucleation of dislocations, the release of extra heat over a greater volume around the tip, the emission of phonons and photons, or even some material change in phase (e.g., softening). In our formalism, l is also the main parameter that describes the shielding of the stress near the crack tip, and this stress is indeed likely to be higher at the centre of the process zone than at its periphery, calling for l < l macro . Large scale bulk energy losses (e.g., [START_REF] Slootman | Quantifying rate-and temperature-dependent molecular damage in elastomer fracture[END_REF][START_REF] Samaca Martinez | Heat and strain measurements at the crack tip of filled rubber under cyclic loadings using full-field techniques[END_REF]) are hence permitted by our model, and we only state that it is of second order effect on the temperature elevation and the stress level at the crack front and, hence, on the crack dynamics.

Besides its straight derivation from a canonical kinetics law (i.e., Eq. ( 1)), there are various indicators sug-FIG. 4. Microscopic fracture energy Uc = d 3 0 Gc/(2l) as a function of the macroscopic energy barrier Gc, for the same materials as in Fig. 3. See the ESI † for the exhaustive list. Note that the accuracy of Uc is not better than an order of magnitude. The horizontal lines show some typical covalent cohesion energies for comparison [START_REF] Kaufman | Inorganic chemistry: principles of structure and reactivity, 4th ed[END_REF]. gesting that, although smaller than the more common l macro , the sizes l that we have inverted bear a strong physical significance. For materials where l is here computed to be relatively big (i.e., in the micrometer range) infrared emissions can easily be recorded around (hot) crack tips by standard resolution infrared cameras, such as in paper [13] or steel [START_REF] Palumbo | Damage monitoring in fracture mechanics by evaluation of the heat dissipated in the cyclic plastic zone ahead of the crack tip with thermal measurements[END_REF]. In the case of paper, we inferred l close to the micrometric domain, which is about the dimension of the fibril forming this media. For running cracks in materials with a smaller l, such as PMMA or even glass, a hot temperature could also be derived by the characterisation of light emission [START_REF] Fuller | The temperature rise at the tip of fast-moving cracks in glassy polymers[END_REF][START_REF] Pallares | Fractoluminescence characterization of the energy dissipated during fast fracture of glass[END_REF], although its observation is not as straightforward. In glass [START_REF] Pallares | Fractoluminescence characterization of the energy dissipated during fast fracture of glass[END_REF], in particular, the light emission lies in the visible domain, calling for blackbody [START_REF] Jain | The physics of blackbody radiation: A review[END_REF] temperatures elevations of thousands of kelvins, and the light emitting zone was shown to have a nanometric radius. In PMMA, secondary microcracks tend to nucleate ahead of fast fronts from cavities/bubbles of about 100 nm in radius [9], possibly initiated by sublimation at hot temperature of the solid matrix. In rocky materials, nanometric damages can be observed [START_REF] Siman-Tov | Nanograins form carbonate fault mirrors[END_REF] on fault planes. These planes are sometimes referred to as fault mirrors, due to their glossy appearance arising from their very low roughness at the visible light wavelengths. Intense thermal effects, often referred to as 'flash heating', are notably suspected for the instability of some seismic faults [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF][START_REF] án | Seismic evidence for thermal runaway during intermediate-depth earthquake rupture[END_REF]. Overall, some intense, extremely localised (i.e., l), dissipation processes close to rupture fronts are thus likely to explain these various observations, and, in our framework, could also explain the brittleness of matter.

VI. DISCUSSION AND CONCLUSION

We have thus presented a model that gives reasonable predictions of the rupture load, over a broad range of materials. We did this with a full expression (Eqs. ( 1) and ( 2)), or in simplified forms (Eqs. ( 1) and (3) or Eq. ( 7)). This predicted load is still, however, overestimated by about 25% (in average for all media, see the errors in Fig. 3). This could derive from numerous causes. First, most of our parameters were only broadly estimated, when not arbitrarily fixed. We have in particular assumed that G c is a homogeneously distributed constant, whereas it is likely to hold some level of quenched disorder [START_REF] Tallakstad | Local dynamics of a randomly pinned crack front during creep and forced propagation: An experimental study[END_REF][START_REF] Cochard | Thermally activated crack fronts propagating in pinning disorder: simultaneous brittle/creep behavior depending on scale[END_REF]. In this case, the overall creep dynamics (i.e. the slow branch of Fig. 2, described by Eq. ( 4)) would not be strongly affected, as it shall mainly depend on an average value of G c . The failure, however, would be prone to occur on weaker locations [8], that are controlled by a lower G c , which would explain our overestimation on G a . It corresponds to the common idea that the overall strength of a material is highly dependent on its heterogeneities. As discussed in section V, our computation of the microscopic energy barrier U c likely represents the average of a disordered landscape in dissociation energy, as its value lies in between that of an actual covalence energy and those of weaker interactions. Additionally, we derived G c from the mean fit to experimental creep curves that, in some materials, hold some significant data dispersion (see the ESI † ). Some of this dispersion could arise from mesoscopic fluctuations of macroscopic fracture energy G c in these materials, and the avalanche threshold G a could be highly dependent on the standard deviation of these fluctuations. Another explanation for the overestimation of G a could be that the experimental error on the measurement of this parameter may in practice be important, as the avalanches occur in a regime where the crack velocity diverges with G, just before test samples snap at a velocity comparable to that of the mechanical waves. Hence, the last mechanical load accurately measured before rupture is, by essence, to be slightly below the actual physical threshold. Note also that, sometimes, the actual creep stage (i.e., 0 to III in Fig. 1) that we fit to derive our parameters is not unambiguously identifiable on the experimental curves, while our theory does not encompass fluid-to-matrix interactions. Besides these considerations, the model is extremely simple, applying mesoscopic laws (i.e., Fourier conductivity and Arrhenius growth) at atomic scales. For instance, a propagative description [START_REF] Tang | On non-fourier temperature wave and thermal relaxation time[END_REF] of the heat transport (i.e., not assuming, such as Fourier diffusion, an infinite transport velocity) could be needed, due to the small time and space scales that are here considered. Overall, a transposition of the model into a, more complicated, atomistic solver [START_REF] Rountree | Atomistic aspects of crack propagation in brittle materials: Multimillion atom molecular dynamics simulations[END_REF] would be beneficial. Of course, the fact that the model reproduces an instability does not necessarily mean that it is the only explanation for britleness. Other models (e.g., the thermal softening of the matrix around the tip [16], or the perturbation of the stress field ahead of the crack by the rupture-induced emission of high frequency phonons [2]) have been proposed for this instability, and, in practice, depending on the materials, more than one physical process could here be at stake. Still, the model we propose gave, in some instances [9], a comprehensive explanation of the full dynamics of failure. Additionally, we have here showed how G c , the intrinsic surface energy barrier of materials, shall only depend on a heat dissipation scale around the crack tip, and that the accumulation of this induced heat is effectively reducing the mechanical resistance of matter (G a < G c ). Countering this latter effect could be a key to design advanced strong materials, in particular as some intriguingly tough solids such as graphene [START_REF] Hwangbo | Fracture characteristics of monolayer cvdgraphene[END_REF][START_REF] Balandin | Superior thermal conductivity of single-layer graphene[END_REF] or arachnid silks [START_REF] Huang | New secrets of spider silk: Exceptionally high thermal conductivity and its abnormal change under stretching[END_REF], are indeed very conductive. Interestingly, the conductivity of spider threads even increases with deformation [START_REF] Huang | New secrets of spider silk: Exceptionally high thermal conductivity and its abnormal change under stretching[END_REF], which could be a nature made adaptive defence mechanism for the stability of nets, whenever they are pressurized. Replicating such a behaviour with a man-made material would then be an important achievement that could lead to high performance cables or bulk materials. For instance, a first step could be the engineering of highly conductive atomic networks, integrated into strong solid matrices, thus limiting any local rise in temperature that could weaken matter. A more down-to-earth application of the model could be the monitoring of structures and infrastructures, as we have shown how their creep rate can be used to predict their failure. This would be of particular interest for bodies that have aged in uncontrolled conditions, in which the change in mechanical properties becomes uncertain with time, but could be inverted from their creep. Finally, and although we have only treated about fracture in mode I, we suggest that most of the effects that we have discussed shall be valid for mixed-mode fracturing and solid friction. The latter is also suspected to hold some non negligible, thermal related, weakening mechanisms [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF][START_REF] án | Seismic evidence for thermal runaway during intermediate-depth earthquake rupture[END_REF], which could notably be a key in geophysics and in understanding the stability of seismic faults. In particular, when increasing the background temperature T 0 , it was shown that the model holds a critical point, beyond which not enough heat can be generated to trigger instabilities in the dynamics of cracks [8], which may physically explain the brittleductile transition in the Earth's crust [START_REF] Scholz | The brittle-plastic transition and the depth of seismic faulting[END_REF][START_REF] Aharonov | The brittle-ductile transition predicted by a physics-based friction law[END_REF], below which rocks tend to flow rather than break. Now that some straightforward expressions for ∆T are known, we can move on to infer G a . Our model, the Arrhenius law as considered in the main manuscript, defines a function S(V, G) such that S(V, G) = V :

S(V, G) = V 0 min exp - α 2 [G c -G] k B [T 0 + ∆T (V, G)] , 1 . (4)
To lighten the equations that will follow, we have here denoted α 2 the ratio d 3 0 /(2l). We have discussed, in the main manuscript, how this relation might have one to three solutions depending on G (see Fig. 1). Two particular energy release rates mark the passages from a singular to multiple solutions: the avalanche threshold G a , of interest in this study, and another threshold, G s , which is the load at which an avalanche has to stop. All functions being continuously smooth, the switch from one solution to three solutions implies that S(V, G) is tangent to the identity function for these two particular G, as illustrated in Fig 1 . G a and the corresponding velocity V a must therefore verify the following system of equations:

S(V, G) = V (5) ∂S ∂V (V, G) = 1. (6)
To solve this system, we assume that the transition towards the fast phase happens in a regime where the temperature elevation still increases linearly with the crack velocity (i.e., ∆T = ∆T slow (V, G) (1)). Equation ( 6) then becomes:

φGλα 2 (G c -G) k b (λT 0 + φGV ) 2 S(V, G) = 1. (7)
Inserting Eq. ( 5) back into (7) leads to the following quadratic equation in V :

φGV λT 0 2 + 2 + α 2 (G -G c ) k b T 0 φGV λT 0 + 1 = 0. ( 8
)
While it might of course hold two solutions, only the lower one is of interest to derive the avalanche threshold G a . The upper solution would indeed correspond to the 'arrest' of the crack avalanche, but the initial hypothesis of ∆T = ∆T slow would there be wrong anyway, as this 'arrest' occurs while on the quick (hot) propagation branch. Focusing therefore on the lower solution of (8), we have:

V a = T 0 λ 2φG a (R a -2 -R a 1 -4/R a ), (9)
with R a = α 2 (G c -G a )/(k B T 0 ). This equation indicates at which slow velocity a crack avalanches, given the corresponding G a threshold. Substituting ( 9) in ( 5), one finally derives the equality that defines the avalanche threshold:

G a ∼ λT 0 2φV 0 R a -2 -R a 1 -4/R a exp -2 1 -1 -4/R a . ( 10
)
Such an expression gives a fairly good approximation of G a as predicted by the model. The only hypothesis was indeed the validity of Eq. ( 1), that is λl/(V a πC)/l 1 and, for the materials that we have studied in our manuscript, this ratio ranges from 300 to 1500. While Eq. ( 10) is easy to solve for G a with any numerical method, it can however be further simplified by grossly assuming that R a 4 and by developing the term 1 -4/R a . We thus obtain the equation presented in the manuscript:

G a ∼ λT 0 φV 0 exp(R a ) R a . ( 11
)
Figure 2 shows the quality of the approximation for G a , off by a few percents as, as shown in Tab. I, the R a 4 hypothesis is not strictly valid.

II. THE ARREST THRESHOLD (FOR COMPLETENESS)

Similarly, one can solve ( 5) and ( 6) at the 'arrest' point: the transition from a quick regime back to the low velocity phase, occurring at the particular load G s . While G a is vastly reported for a lot of materials, making it the topic of this manuscript, G s is more rarely reported, so that the following computation is given for completeness. We here assume that the transition arises when the crack cools down from the plateau temperature ∆T = ∆T fast (G) (2), along the intermediate slope defined by ∂∆T /∂V = ∂∆T mid (V, G))/∂V (3). We thus turn the system into a quadratic equation of √

V : φG √ V 4πλClT 0 2 + 2 + α 2 (G -G c ) 2k b T 0 φG √ V 4πλClT 0 +1 = 0, (12)
the upper solution of which, together with Eq. ( 5), leads to:

V s = πλClT 0 2 4(φG s ) 2 R s -4 + R s 1 -8/R s 2 , (13)
where R s = α 2 (G c -G s )/(k B T 0 ). When inserting (13) back into (5), one gets:

4(φG s ) 2 V 0 πλClT 0 2 = R s -4 + R s 1 -8/R s 2 exp α 2 (G s -G c ) k b [T 0 + φG s /(πCl)] . (14)
Assuming that R s 8 and ∆T fast T 0 , Eq. ( 14) further simplifies to:

G s ∼ T 0 φ πλCl V 0 exp πClT 0 2φG s R s , (15)
which gives a relatively simple expression to invert for G s .

III. SENSITIVITY OF THE φ PARAMETER

The percentage φ of energy that is converted into heat in the heat zone of radius l is a parameter that we have broadly assumed to be 0.5. In practice, this parameter is unknown, and shall likely be material dependent. In Fig. 3, we show the variation in our model prediction of the failure threshold G a when varying φ, for three materials spanning the whole range of considered toughnesses. The higher φ, the less is G a (e.g., see Eq.( 11)), as it makes the crack tip hotter at a lower load. As shown in the figure, when considering small values of this heat efficiency, one can also obtain an all ductile behavior 2 , as the crack becomes too cold to hold an instability in its dynamics. In this case, the modelled crack velocity only increases exponentially to V 0 , as a function of the energy release rate G, and no velocity jump (defining G a ) is obtained.

FIG. 3. Variation of the model predicted failure threshold Ga as a function of the heat efficiency φ, relative to the experimentally reported Ga, for a Titanium alloy (Gc ∼ 190 kJ m -2 ), TZP ceramic (Gc ∼ 1900 J m -2 ), and soda-lime glass (Gc ∼ 12 J m -2 ). The considered model parameters for these materials are shown in table I.

IV. MATERIALS CREEP CROSSPLOTS AND PARAMETERS TABLE

A summary of the model parameters considered for each media is also provided in Tab. I. These parameters are deduced, as explained in the main manuscript, from the V to G creep data of these materials, shown in Fig 4 to Fig 21. One can notably notice the variability in fit quality for these datasets, that of course impacts our inversion work, but also how it is not always straightforward to know to which subcritical phase the data correspond (i.e., phase I to III, from environmental induced corrosion to void-like conditions).

FIG. 4. Creep data of dry soda-lime glass, from Wiederhorn [3] , figure 3. A rather complex creep law holds there so that we only roughly fitted the last part (i.e., stage III). FIG. 5. Creep data of dry sapphire (r-plane), from Wiederhorn and Krause [4] . λ (SI) φ (-) V0 (m/s) l ( Å) T0 (K) Gc (J/m 2 ) Ga real (J/m 2 ) Ga model (J/m FIG. 6. Creep data of quartz in vacuum, from Dove [5] , figure 4.

FIG. 7. Creep data of Scioto sandstone, from Holder et al. [6] , figure 3.

FIG. 8. Creep data of Kumamoto andesite in moist air at 67 °C, from Nara and Kaneko [7] , figure 9.

FIG. 9. Creep data of Westerly granite in moist air at 20 °C, from Meredith and Atkinson [8] , figure 7.

FIG. 10. Creep data of vitreous carbon, from Nadeau [9] , figure 4.

FIG. 11. Creep data of high strength ultra low porosity concrete in moist air, from Nara et al. [10] , figure 9.

FIG. 12. Creep data cement in water, from Wang et al. [11] , figure 4a.

FIG. 13. Creep data of paper in air, from Santucci [12] , figure 3.32.

FIG. 14. Creep data of hot silicon nitride at 1200 °C, from Evans and Wiederhorn [13] , figure 5.

FIG. 15. Creep data of Lead Zirconate Titanate at ambient conditions, from Oates et al. [14] , figure 2 (open circuit).

FIG. 16. Creep data tetragonal zirconia (TZP) in vacuum, from Chevalier et al. [15] , figure 5.

FIG. 17. Interfacial creep data in sintered PMMA plates in air, from Lengliné et al. [16] , figure 5.

FIG. 18. Creep data of high density polyethylene, from Yoda et al. [17] , figure 4.

FIG. 19. Creep data of aluminium 2650 T6 alloy in vacuum at 175 °C, from Hénaff et al. [18] , figure 6.

FIG. 20. Creep data in Ti-6A1-6V-2Sn titanium alloy in moist air, from Sastry et al. [19] , figure 6a (beta annealed).

FIG. 21. Creep data of AISI 310S austenitic stainless steel in air, from Huang and Altstetter [20] , figure 1 (uncharged plot).

FIG. 2 .

 2 FIG. 2. Modelled crack velocity as a function of energy release rate, as per Eqs. (1) and (2). Stages III and IV correspond to the one labelled in Fig.1. As explained in the text, stages 0, I and II are here not covered. In our model, the failure occurs when the cracks becomes hot enough, that is, when ∆T ∼ T0. The dashed line corresponds to a cold case ∆T T0 in Eq. (1).

FIG. 3 .

 3 FIG. 3. (Bottom): Modelled Ga thresholds (squares) and modelled surface energy barrier Gc (triangles) compared to the experimental thresholds from the literature. The black line is the identity. The labels locate different materials. The unlabelled rock materials are quartz, sapphire, granite and andesite. See the ESI † for an exhaustive list. (Top): Relative error on the avalanche threshold.

FIG. 5 .

 5 FIG. 5. (Left):Core size of the process zone l, as understood by our model, versus macroscopic plasticity scale lmacro, as derived from reported tensile yield stresses. The straight lines mark a factor 10 5 , 10 6 and 10 7 between both views. The two unlabelled points in the vicinity of Glass represent Quartz and Concrete. (Right): Simplified spacial distribution of the intensity (arbitrary unit) at which energy is dissipated (i.e., plasticity) around the crack tip, as a possible explanation for the difference in scale. The graph is not to scale, as lmacro l.

FIG. 1 .

 1 FIG.1. Representation of V = S(V, G) for three values of G: Gs, Ga and a mid-value between Gs and Ga (plain plot). The intersections of SG with the identity plot (straight line) give the possible crack velocities for a given energy release rate, as per Eq. (4). The axes are not annotated for the sake of generality. See Ref.2 for further information.

FIG. 2 .

 2 FIG. 2. (Bottom): Ga threshold, as approximated by Eq. (11) versus the accurate numerical solution of the model. The black line is the identity. (Top): Relative error from the approximation

TABLE I .

 I Model parameters for various materials of the literature. The real and modelled Ga thresholds are compared in the two former last columns. The cells colour help to highlight standing out values for λ and T0.

										2 ) Ra (-)
	Acrylic adhesive	0.4	1	30	10	296	150	90	97	5.7
	Paper	0.035 0.12	1300	1000 296	25000	14000	9500	15.6
	Bulk PMMA	0.18 0.2	880	80	296	1300	700	580	10.9
	Interfacial PMMA	0.18 0.2	880	8	298	275	140	190	13.5
	HD Polyethylene	0.4	0.5	900	8500 293	200000	70000	87000	16.6
	Soda lime glass	1	0.5	3400	0.3	296	12	8	10	8.3
	Sapphire	24	0.5	6000	0.8	296	36	20	32	6
	Quartz	8	0.5	3400	0.6	293	21	13	18	5.7
	Westerly Granite (ambient)	2	0.5	3000	4	293	120	68	92	8.5
	Westerly Granite (hot)	2	0.5	3000	0.7	573	43	24	35	6.8
	Kumamoto Andesite	1	0.5	2200	3	330	120	80	97	8.8
	Scioto Sandstone	2	0.5	2000	2	296	55	37	44	7.3
	Cement paste	1	0.5	2200	3	298	310	250	280	10.7
	HSULP Concrete	0.8	0.5	3000	1	293	44	38	40	9.9
	Vitreous carbon	5	0.5	2600	0.2	296	15	13	14	7.2
	Lead Zirconate Titanate (PZT)	1	0.5	2000	1	296	40	24	33	11.3
	Tetragonal zirconia (TZP)	2	0.5	1600	40	298	1900	1500	1530	10.9
	Silicon nitride	30	0.5	5500	45	1573	510	260	400	8.9
	2650 T6 Aluminium alloy	150 0.5	3100	1000 448	54500	27000	39000	10.1
	AISI 310S Stainless Steel alloy	14	0.5	3000	9000 298	265000	102000	158000	13.4
	Ti-6A1-6V-2Sn Titanium	7	0.5	3100	8000 298	190000	72000	93000	14.9

ACKNOWLEDGEMENTS

The authors declare no competing financial interests in the publishing of this work and acknowledge the support of the IRP France-Norway D-FFRACT, of the Universities of Strasbourg and Oslo and of the CNRS INSU ALEAS program. We thank the Research Council of Norway through its Centres of Excellence funding scheme, project number 262644.

Thermal dissipation as both the strength and weakness of matter.

A material failure prediction by monitoring creep ------------------------Electronic Supplementary Information (ESI)

Tom Vincent-Dospital, Let us start this additional material with the analytical approximations of the temperature at a running crack tip. Assuming a quasi-constant velocity and energy release rate, some simplified expressions can indeed be derived 1 for ∆T . At low velocity, the typical diffusion skin depth is large compared to the radius of the heat production zone ( λl/(V πC)/l 1) and the heat diffusion is hence the ruling process:

At high velocity, however, the rise in temperature is limited by the scale over which heat is produced and:

Between these two cases, and typically for V ∼ λ/(πCl), an intermediate regime holds:

We invite the reader to a more in-depth derivation of these equations in Toussaint et al. [1] or Vincent-Dospital et al. [2] .


Giovanni Azzolina

Roman Bertoni

Eric Collet
email: eric.collet@univ-rennes1.fr

General Landau Theory of Non-Symmetry-Breaking and Symmetry-Breaking Spin Transition Materials

Spin-transition materials, including the families of spin-crossover and charge-transfer systems, and more generally molecular-based materials exhibiting electronic and/or structural bistability, may undergo various types of phase transitions. The change of electronic state is stabilized by molecular reorganizations and both phenomena, which are usually non-symmetry breaking, can be described through the evolution of an order parameter q. Due to symmetry, q linearly couples to volume change. It is known that such elastic interactions are responsible for cooperative phenomena in non-symmetry-breaking spin-transitions. However, spin-transition materials may also exhibit symmetry-breaking phenomena related to various types of orders such as structural order as well as spin-state concentration waves. The universal framework of the Landau theory of phase transition is relevant for describing such ordering processes through the evolution of a symmetry-breaking order parameter . The simultaneous or sequential occurrence of spin transition and symmetry breaking phenomena are reported for numerous spin-transition materials and the coupling between these two types of instabilities is responsible for the emergence of various types of functions. In this work, we use the Landau approach to describe both symmetry-breaking phenomena and non-symmetry-breaking spin transition. We discuss how their coupling can generate sequences of phase transitions, from simple spin crossover to spin transition, continuous or discontinuous symmetry breaking, including ferroelasticity or stepwise spin transitions.

I. Introduction

Spin-transition materials, including the families of spin-crossover, [START_REF] Halcrow | Spin-crossover materials : properties and applications[END_REF][START_REF] Gutlich | [END_REF][3] charge-transfer (CT) systems [4][5][6][7] , may undergo various types of phase transitions. The concept of molecular bistability expands to various types of systems exhibiting electronic and/or structural bistability between states of different entropy and volume. [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24] Phase transitions may be driven by external stimuli including temperature, pressure, magnetic field, electric field or light, allowing for the emergence of functions. In many cases, both the change of electronic state and the coupled structural reorganizations are non-symmetry breaking, which can be described through the evolution of an order parameter q. This order parameter q transforms as the identity representation of the symmetry group and therefore linearly couples to the volume strain. That is the reason why cooperative phenomena in non-symmetry-breaking spin-transitions can occur, due to the relative change of the bonding or antibonding nature of the electronic distribution that couples to lattice expansion.

Materials may also exhibit symmetry-breaking phenomena, related to diverse types of electronic and/or structural orders, including ferromagnetism, ferroelectricity or ferroelasticity. The universal Landau theory of phase transitions allows describing these symmetry-breaking phenomena through the evolution of one (or several) symmetry-breaking order parameter(s) , the amplitude of which measures the deviation from the high symmetry phase.

The non-symmetry-breaking change of electronic state is an instability, which may couple to a symmetry-breaking phenomenon as recently highlighted. 25 In this paper, we extend the Landau approach [START_REF] Halcrow | Spin-crossover materials : properties and applications[END_REF][START_REF] Gutlich | [END_REF][25][26][27][28][29] to study the key role of the coupling between the non-symmetry-breaking order parameter q, describing the electronic instability and a symmetry-breaking instability , describing the structural order. The richness of behaviors extracted from the models explains the diversity of spin transition curves reported so far, with non-symmetry-breaking and symmetry-breaking instabilities that may occur simultaneously or sequentially. We discuss how the coupling between parameters may generate sequences of phase transitions, from simple spin crossover to spin transition, with continuous or discontinuous symmetry breaking. The coupling terms enlarge the hysteresis domain of bistability or generate stepwise spin transitions with or without spin-state ordering.

II. Spin transition or crossover without symmetry breaking

The bistable and non-symmetry-breaking nature of molecular spin states can be represented by an Ising variable (qi=±1) corresponding to one of the two states. This model has been intensively used in the case of spin-crossover molecules: qi= +1corresponds to the high spin (HS) state and qi=-1 to the low spin (LS) state. 16,[30][31][32][33][34] What we discuss hereafter is more general than the description of spin transition materials only as the states +1 and -1 may refer to spin state for spin-crossover materials and more generally to one of the bistable molecular states (CT, isomerization…) for other systems.

Different types of thermal spin state conversions occur, depending on the strength of the interactions, of elastic nature, among molecules, from continuous spin-crossover for weak interactions to discontinuous spin transition with thermal hysteresis for strong interactions. The spin-state transition is described through the fraction  of molecules in the high spin state or through the order parameter q:  = 𝑁 𝐻𝑆 𝑁 𝐻𝑆 +𝑁 𝐿𝑆 and 𝑞 = 𝑁 𝐻𝑆 -𝑁 𝐿𝑆 𝑁 𝐻𝑆 +𝑁 𝐿𝑆 (1) with  = 𝑞+1 2

. 𝑁 𝐻𝑆 and 𝑁 𝐿𝑆 denote the number of sites in HS or LS states and 𝑞 = +1 in the full HS state while in the fully LS state 𝑞 = -1. The order parameter q does not break symmetry and the spin transition (ST) is described by expanding the thermodynamic potential of the system in powers of q. Since q transforms as the identity representation of the symmetry group, the thermodynamic potential includes all powers of q. As discussed in previous works, 25,35 truncating the potential at fourth order allows for renormalizing q to eliminate the third-order term for simplifying the potential to:

𝐹 = 𝐴′(𝑇)𝑞 + 1 2 𝐵′𝑞 2 + 1 4 𝐶𝑞 4 + 𝜆 𝑞 𝑣 𝑠 ( 1-𝑞 2 ) + 1 2 𝐶 𝑠 0 𝑣 𝑠 2 (2)
𝐴′(𝑇) changes sign with temperature. C>0 is required for stability. ) is the elastic coupling term of q to the volume strain 𝑣 𝑠 that scales as (

Therefore, the elastic coupling shifts the equilibrium line between q<0 and q>0 and decreases the q 2 coefficient. We then rewrite the potential as:

𝐹(𝑞, 𝑇) = 𝐴(𝑇)𝑞 + 1 2 𝐵𝑞 2 + 1 4 𝐶𝑞 4 (4)
Since the goal of the present paper is to discuss the different types of thermal behaviors, we will use potential (4) and consider for pedagogical purpose that B is constant and that 𝐴 changes linearly with temperature: A=-a0(T-TST). We will discuss the temperature dependence of q for different values of B. As A and B also depend on pressure, our approach is analogous to crossing the phase transition line of the (𝑃, 𝑇) phase diagram not at constant pressure but in an oblique way. The stability condition of q, 𝑑𝐹 𝑑𝑞 = 0 and

𝑑𝐹 2
𝑑𝑞 2 > 0, provides the spin transition curve q(T), also given as (T), from predominantly LS (q<0) below 𝑇 𝑆𝑇 to predominantly HS (q>0) above. In this potential, B>0 corresponds to a gradual spin state conversion, the so-called spin crossover behavior, with a continuous evolution from q=1 at high temperature, q=0 at T=TST and q= -1 at low temperature (Figure 1a and 2a). B<0 promotes cooperativity (Figure 1b) and the resulting spin transition curve (T), or q(T), has a the characteristic "S shape" (Figure 2b). Indeed, F(q) exhibits two symmetric This model does not provide stable phases at half conversion (q=0) over a broad temperature range, as it is the case for stepwise spin transition discussed in sections IV and V. This phenomenological model is relevant for describing spin conversion. It is clear that the degree of cooperativity, described through the parameter B, is the key for explaining spin conversion of spin-crossover type (B>0) or spin-transition type (B<0). Equations ( 2) and (3) highlight that it is the elastic coupling of q to the volume strain 𝑣 𝑠 that makes the q 2 coefficient more negative, which is responsible for the appearance of spin transition or spin-crossover.

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. The appearance of spin-crossover or spin transition was discussed for the molecular system Fe(PM-BIA)2(NCS)2, which may crystallize in two polymorphs that show drastically different spin conversions. [36][37][38][39] The orthorhombic polymorph (Pccn) exhibits an abrupt spin transition with a thermal hysteresis (width ≈5 K), which corresponds to the cooperative case in Figure 2b (B<0), while the monoclinic polymorph (P21/c) shows only a gradual spin-crossover type conversion, which corresponds to the crossover case in Figure 2a (B>0). It was shown that the different intermolecular couplings of both polymorphs play an important role for the degree of cooperativity, 40 which translates in our model in different coupling strength 𝜆 𝑞 to the volume strain 𝑣 𝑠 in Eq (2). This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0041453

PLEASE CITE THIS ARTICLE AS

Various techniques allow monitoring the spin state switching, including Mossbauer spectroscopy, magnetic susceptibility, optical, Raman or X-ray spectroscopies to name a few. 41 X-ray diffraction can monitor the structural changes at the molecular scale. The average metal-ligand bond length < 𝑀 -𝐿 > obtained from structural refinement is the weighted contributions of HS (< 𝑀 -𝐿 > 𝐻𝑆 ) and LS (< 𝑀 -𝐿 > 𝐿𝑆 ) states: < 𝑀 -𝐿 > (𝑇) = 𝛾(𝑇) < 𝑀 -𝐿 > 𝐻𝑆 + (1 -𝛾(𝑇)) < 𝑀 -𝐿 > 𝐻𝑆 . 28 < 𝑀 -𝐿 > 𝐻𝑆 and < 𝑀 -𝐿 > 𝐿𝑆 are measured in the full HS and LS states and the thermal dependence of the HS fraction, 𝛾(𝑇), is experimentally extracted from:

𝛾(𝑇) = | <𝑀-𝐿>(𝑇)-<𝑀-𝐿> 𝐿𝑆 <𝑀-𝐿> 𝐻𝑆 -<𝑀-𝐿> 𝐿𝑆 | = 𝑞(𝑇)+1 2 (5)
This simplest model for describing non-symmetry-breaking spin conversion does not allow describing multi-step spin conversion and/or symmetry-breaking phenomena that may couple to the change of spin state. We discuss hereafter how to take into account symmetry breaking and how the coupling of symmetry breaking and spin transition can lead to various types of phase transitions and phase diagrams.

III. Symmetry breaking phase transition

Symmetry breaking (SB) may be associated with structural ordering phenomena, including ferroelectric phase transition, ligand, anion/cations motions, or ferroelastic phase transitions 8 for which the crystal system changes (from cubic to tetragonal, from orthorhombic to monoclinic…). [START_REF] Salje | Phase Transitions in Ferroelastic and Co-elastic Crystals[END_REF] In many cases, when a symmetry breaking occurs, the low symmetry phase is a subgroup of the high symmetry phase, i.e. some symmetry operators of the high symmetry phase are lost in the low symmetry phase, while the others are preserved. Then, the Landau theory of phase transition can be applied. [START_REF] Landau | [END_REF] For measuring the deviation from the high symmetry phase, the theory uses the amplitude of , the symmetry-breaking order parameter (OP), which obeys the transformation properties of the irreducible representation of the symmetry group that is responsible for the symmetry change. The symmetry-breaking OP may have one, two or three dimensions depending of the representation. In the case of ferroelastic phase transitions, the symmetry breaking corresponds to a change from one crystalline system to a lower symmetry one. For example, the cubic-tetragonal ferroelastic distortion is characterized by the difference between the a and c crystalline axis, which are equal in the high symmetry cubic phase:  ∝ (𝑐-𝑎)

𝑎

. 25 The monoclinictriclinic distortion is characterized by the deviation of lattice angles  and  from 90°... Symmetry breaking may also occur without change of lattice, due to the relative motions ∆𝑟 ⃗ of atoms from the high symmetry position and  ∝ |∆𝑟|.

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0041453

The Landau theory describes the evolution of the thermodynamical potential F of the system around the phase transition temperature, through the expansion of F in power series of . The orders allowed depend on the representation of the order parameter. Hereafter we discuss the simplest case where  is one-dimensional and scalar, which allows two symmetry-equivalent domains (±) in the low symmetry phase Since both domains are equivalent by symmetry, 𝐹(𝜂) = 𝐹(-𝜂), only even orders in  are allowed. The simplest temperature dependence of F with 𝜂 near the phase transition is given by:

𝐹(𝜂, 𝑇) = 1 2 𝑎(𝑇)𝜂 2 + 1 4 𝑏𝜂 4 (6)
b>0 is required for stability. a(T)=a0(T-TSB) changes sign at TSB (a0>0), where the symmetrybreaking phase transition occurs. This potential describes a second-order phase transition between the high and low symmetry phases (Figure 3a), as  deviates continuously from 0 below TSB (Figure 4a). We use the same a0 in Eq ( 4) and Eq ( 6) in order to reduce the number of parameters in the model.

We may also consider a sixth order term:

𝐹(𝑇) = 1 2 𝑎(𝑇)𝜂 2 + 1 4 𝑏𝜂 4 + 1 6 𝑐𝜂 6 (7)
with c>0 for stability. b>0 gives the same scenario of continuous phase transition described above.

On the contrary, b<0 gives rise to a first-order phase transition, for which the high symmetry phase (=0) is stable down to TSB and the low symmetry phase (≠0) is stable up to 𝑇 1 = 𝑏 2

8𝑎𝑐

+ 𝑇 𝑆𝐵 (Figure 3b). The symmetry-breaking curve   (T) exhibits then a thermal hysteresis where  changes discontinuously at 𝑇 𝑆𝐵 and 𝑇 1 (Figure 4b).

A similar first-order symmetry-breaking phase transition was found in the case of the cubictetragonal ferroelastic transition occurring in RbMnFe charge-transfer systems. In this case, the ferroelastic symmetry-breaking order parameter  belongs to a bidimensional representation, 25 and F includes an  3 term, limiting symmetry-breaking phase transition to first-order types. Therefore  changes discontinuously from the high symmetry phase (=0) to the low symmetry phase (≠0)

and the thermal evolution of the order parameter  is then similar to what is shown in Figure 4a.

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. The volume change at the molecular scale, due to a redistribution of electrons between more or less bonding states, is an important feature common to many spin-transition materials. In crystals, the average change of electronic state q induces an important volume strain that may reach 10%, as observed in spin transition materials as well as CoFe or MnFe Prussian blue analogues. 17,25,44 Both the non-symmetry-breaking order parameter q and the symmetry-breaking order parameter  can couple to the volume strain, 25 which results in an effective coupling between q and  2 .

PLEASE CITE THIS ARTICLE AS

Hereafter we consider the total free energy including the contributions of the non-symmetrybreaking order parameter q and the symmetry-breaking order parameter 𝜂, as well as various symmetry-allowed coupling terms.

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0041453

IV. Linear-quadratic coupling of spin transition and symmetry-breaking order parameters

Since q transforms as the identity representation and 𝜂 as a symmetry-breaking representation, the symmetry-allowed coupling term of lowest order is 𝐷𝑞𝜂 2 . The potential reaches:

𝐹(𝜂, 𝑞, 𝑇) = 1 2 𝑎𝜂 2 + 1 4 𝑏𝜂 4 + 1 6 𝑐𝜂 6 + 𝐴𝑞 + 1 2 𝐵𝑞 2 + 1 4 𝐶𝑞 4 + 𝐷𝑞𝜂 2 (8)
We studied previously a similar potential, and highlighted the elastic nature of the 𝐷𝑞𝜂 2 coupling term. 25 We use 𝑎 0 (𝑇 𝑆𝐵 -𝑇 𝑆𝑇 ) = -𝐴 -𝑎, to monitor the difference of temperature instability between the symmetry-breaking (𝑇 𝑆𝐵 ) and the spin (𝑇 𝑆𝑇 ) transitions. The different phases that may appear in the (𝑇, 𝑇 𝑆𝐵 -𝑇 𝑆𝑇 ) space, are characterized by the equilibrium values of the OP corresponding to a minimum of the potential in the (𝑞, 𝜂) space with:

𝑑𝐹 𝑑𝜂 = 0, 𝑑𝐹 𝑑𝑞 = 0, 𝑑𝐹 2 𝑑𝜂 2 > 0, 𝑑𝐹 2
𝑑𝑞 2 > 0 and

𝑑𝐹 2 𝑑𝑞𝑑𝜂 > 0 (9)
Here we do not explore the regions of the phase diagram were only the spin transition or only the symmetry breaking transition occurs, as this is already discussed in previous sections. Instead, we focus our attention on the spin transition curve q(T), also scaled in (T), and symmetry breaking curve (T). The two transitions may occur sequentially or simultaneously, depending on 𝑇 𝑆𝐵 -𝑇 𝑆𝑇 and the coupling strength D. The different phases are characterized by their spin state ( > Therefore, we use D>0 for stabilizing the low temperature LSls phase (𝑞 < 0, 𝜂 ≠ 0). The relative stability of the different phases depends on 𝑇 𝑆𝐵 -𝑇 𝑆𝑇 , the temperature T and the coupling strength D. The order of symmetry-breaking and spin transitions depends on 𝑇 𝑆𝐵 -𝑇 𝑆𝑇 . For the calculations we used b=±2, B=±2, c=6, C=4, D=2-4.

Figure 5 shows the phase diagram in the (𝑇 𝑆𝐵 -𝑇 𝑆𝑇 , T) space, when the spin transition (B<0) and the discontinuous symmetry-breaking (b<0), are both associated with thermal hysteresis. Figure 6 shows different sequences of symmetry-breaking and/or spin transitions of continuous or discontinuous nature, calculated for the different combinations (b<0, B<0), (b>0, B<0), (b<0, B<0)

or (b<0, B>0). Each panel shows the curves for 𝑇 𝑆𝐵 < 𝑇 𝑆𝑇 , 𝑇 𝑆𝐵 = 𝑇 𝑆𝑇 and 𝑇 𝑆𝐵 > 𝑇 𝑆𝑇 .

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. (b<0, B<0). The results are similar to the ones of the potential including a symmetry-allowed  3 , restricting symmetry-breaking phase transition to first-order types only. 25 When TST and TSB are different, the symmetry-breaking and spin transitions occur sequentially (top and bottom panels), with intrinsic hysteresis widths. Due to the coupling, the symmetry-breaking phase transition is accompanied by a weak discontinuous change of q, and when TST and TSB are close (middle panel), the symmetry-breaking and spin transition occur simultaneously with a broader hysteresis loop.

PLEASE CITE THIS ARTICLE AS

This case maps several scenarios, including non-symmetry breaking CT-based phase transition in the Rb0.73MnFe compound, 45 or purely symmetry-breaking phase transitions in Prussian blue analogues 17 . The broadening of the hysteresis width with the coupling strength, is similar to the broadening observed under chemical pressure, 46 or the thermal shift of the hysteresis under hydrostatic pressure. 47 temperature. This last sequence of phase transitions corresponds to the ones observed for the [FeH2L 2-Me ](ClO4)2 complex, 49 exhibiting a cooperative and non-symmetry-breaking spin transition at high temperature followed by a continuous symmetry-breaking at lower temperature.

For this system, the symmetry breaking is associated with a cell doubling, resulting from a longrange ordering of bonds on the ligand and motions of the anions. Figure 7 is a schematic representation of the sequence of phases with a spin transition at high temperature and a symmetry breaking at low temperature that may involve cell doubling or ferroelastic distortion. ST coupled to a tetragonal to monoclinic ferroelastic transition(I41/a to P21/c) 51,52 and the potential in equation ( 8) including only even orders in  must be used. A discontinuous symmetry breaking phase transition coupled to a spin-crossover corresponds also to the case of the [Fe(ptz)6](BF4)2 system. [53][54][55][56][57][58] At atmospheric pressure this system exhibits a spin conversion coupled to a first-order ferroelastic phase transition and both phenomena occur simultaneously, 54 while both phenomena occurs separately under pressure. 59,60 V Non-symmetry-breaking stepwise phase transitions Non-symmetry breaking stepwise phase transitions occur in systems for which the crystalline lattice is made of sub-lattices with active molecular sites independent by symmetry. In the case of two sub-lattices, the spin transition of each sub-lattice is monitored by 𝑞 1 and 𝑞 2 . We may approximate the total thermodynamical potential F as a sum of the thermodynamical potentials 𝐹 1

and 𝐹 2 of the sub-lattices and the elastic coupling term 𝐹 𝑐 between 𝑞 1 and 𝑞 2 :

𝐹(𝑞 1 , 𝑞 1 ) = 𝐹 1 + 𝐹 2 + 𝐹 𝑐 = 𝐴 1 𝑞 1 + 1 2 𝐵 1 𝑞 1 2 + 1 4 𝐶 1 𝑞 1 4 + 𝐴 2 𝑞 2 + 1 2 𝐵 2 𝑞 2 2 + 1 4
𝐶 2 𝑞 2 4 + 𝐷 12 𝑞 1 𝑞 2 (10) When the elastic coupling between the sub lattices is negligible (𝐷 12 = 0), the global spin conversion is then described through the weighted contribution from each sub-lattice 𝑞 = . In this case, the intermediate state at half conversion corresponds to a crystalline structure with molecules mainly in HS state on sub-lattice 1 and LS state on sub-lattice 2.

We discuss now the effect of the elastic coupling 𝐷 12 between two sub-lattices undergoing cooperative spin transitions (𝐵 1 < 0 and 𝐵 2 < 0). The coupling may be of anti-ferroelastic nature (𝐷 12 > 0), which stabilizes sub-lattices in different spin states. This anti-ferroelastic coupling favors therefore the half conversion state at 𝛾 = 1 2

, for which a single sub-lattice undergoes the spin transition. This state with one sub-lattice mainly HS and one mainly LS is stabilized over a broader temperature range (Figure 8d) compared to the case without coupling (Figure 8c). On the contrary, a ferroelastic coupling 𝐷 12 < 0 favors the same spin state for both lattices. This ferroelastic coupling destabilizes the half-conversion step and when the coupling is strong enough there is a single step spin transition (Figure 8e), where both sub-lattices switch from LS to HS states simultaneously.

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0041453

Figure 8 Non-symmetry-breaking stepwise spin transition curves (T) for different potentials (10) including spin crossover and/or spin transition for the two sub-lattices at TST1 and TST2 (a-c). The anti-ferroelastic coupling (D12>0) stabilizes the half conversion step (d), while the ferroelastic coupling (D12>0) favors a single step spin transition (e).

The compound [Fe(btr)3](ClO4)2, 61 exhibits such a LS ↔ HS spin conversion occurring in two steps: the low-temperature step is very abrupt and occurs with a thermal hysteresis (3 K width), while the high-temperature step is rather gradual. This situation is described in our model with 𝐵 1 > 0 and 𝐵 2 < 0 (Figure 8a). In this system, there are two different Fe 2+ sub-lattices, which interact differently with the anions, and in the intermediate phase the Fe1 sites remain HS while the Fe2 sites are LS (Figure 9a). Other systems exhibit similar behaviours. 62,63 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. A similar half spin transition curve was also reported in [Fe(3PhOH-trz)2Pt(CN)4]•2H2O with an hidden hysteresis for the LS-LS phase. 65 However, the process is different since it is associated with a symmetry change, forming two independent Fe sites in the intermediate phase from a single site in the fully HS and high symmetry phase. We discuss this relation between stepwise transition and symmetry breaking hereafter.

PLEASE CITE THIS ARTICLE AS

VI Stepwise and symmetry breaking spin transition:

The competition between short-range elastic interaction, favoring alternation of molecules in different spin states with different volume, and long-range elastic interactions due to the global crystalline volume change, favoring molecules in the same spin state, can generate stepwise spin transition curves, 30,31,66 as reported in many Fe II , Fe III , Co II , and Mn III complexes. 28,30,31,48,61,63,[67][68][69][70][71][72][73][74][75][76][77][78] Variety of stepwise of processes were studied theoretically. 79,80 More and more systems are reported to undergo sequences of phases with a fully HS state at high temperature, an intermediate (𝑞 1 -𝑞 2 ). The symmetry change can also involve atomic motions, ligand or lattice torsion and rotation... As it is the case for non-symmetry-breaking spin transition, the HShs phase is characterized by (0 < 𝑞 ≤ 1,  = 0) and the LShs phase is characterized by (-1 ≤ 𝑞 < 0,  = 0). In the INT phase, if 𝑞 = 0 and the sub-lattice 1 is fully HS and 2 is fully LS ( = 1) the order is complete. The other domain, or the antiphase depending on the nature of the symmetry-breaking, corresponds to sub-lattice 1 fully LS and 2 fully HS ( = -1). These two domains can equally form, and in most cases the order is partial in the symmetry-broken INT phase (0 < || < 1). Here again, a phase transition where spin transition and symmetry breaking occur simultaneously requires coupling between q and . In addition to the symmetry-allowed coupling term 𝐷𝑞𝜂 2 discussed above, it is the bi-quadratic coupling term 𝐸𝑞 2 𝜂 2 , with 𝐸>0, that stabilizes the half conversion step, by minimizing the thermodynamical potential around q=0 for  ≠ 0. For limiting the number of terms in the potential, we neglect the 𝐷𝑞𝜂 2 term, which shifts the stability region of the INT phase as explained by Chernyshov, 35 and we use the most simple potential:

𝐹(𝜂, 𝑞) = 1 2 𝑎𝜂 2 + 1 4 𝑏𝜂 4 + 𝐴𝑞 + 1 2 𝐵𝑞 2 + 1 4 𝐶𝑞 4 + 𝐸𝑞 2 𝜂 2 (11)
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0041453

Here again the stability condition is found with:

𝑑𝐹 𝑑𝜂 = 0 = (𝑎 + 2𝐸𝑞 2 + 𝑏𝜂 2 )𝜂 = 0.
𝜂 = 0 is a trivial solution and corresponds to the non-symmetry-breaking spin transition case.

𝜂 ≠ 0 is found for 𝜂 2 = - 𝑎+2𝐸𝑞 2 𝑏
, with a stability condition

𝑑𝐹 2 𝑑𝜂 2 = 𝑎 + 3𝑏𝜂 2 + 2𝐸𝑞 2 > 0.
Substituting 𝜂 2 leads to 𝑎 + 2𝐸𝑞 2 < 0, which can also be written 𝑇 -𝑇 𝑆𝐵 < 2 𝐸 𝑎 0 𝑞 2 . The low symmetry INT phase with HS-LS order and 𝜂 ≠ 0 is then stabilized in a temperature range

[𝑇 𝑆𝐵2 -𝑇 𝑆𝐵1 ] around q=0. Figure 10 shows the spin transition curve q(T) and the symmetry breaking curve  2 (T) corresponding to cooperative symmetry-breaking stepwise transitions, with thermal hysteresis.

In many cases reported so far, the symmetry breaking is associated with a change of translation symmetry (cell doubling, change of lattice type or crystalline system) or the loss of nonsymmorphic symmetry operators (screw axis or glide planes), making molecular sites inequivalent in the crystal. The appearance of new Bragg peaks related to such symmetry-breaking is then often a direct characterization of the symmetry change, and the intensity of the Bragg peaks is proportional to  2 . 30 A more direct way to probe the amplitude of the symmetry breaking due to the ordering of HS and LS molecules on different sub-lattices, is to measure by diffraction techniques the difference of the average metal-ligand bond lengths <M-L> between sub-lattices 1 and 2. 28  scales with the difference between the sites normalized to the difference between fully HS and fully LS states:

 = | <𝑀-𝐿> 1 -<𝑀-𝐿> 2 <𝑀-𝐿> 𝐻𝑆 -<𝑀-𝐿> 𝐿𝑆 | = |(𝑞 1 -𝑞 2 )| (12)
Figure 10 Symmetry-breaking stepwise spin transition curves (T) and  2 (T) for potential (11).

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0041453

Here again, the order is complete when one sub-lattice is fully HS and the other one fully LS, i.e.

when the bond lengths <M-L> of one sub-lattice site corresponds to the LS state and the one of the other sub-lattice to the HS state. For partial order the symmetry is broken (0 < || < 1) and the bond lengths are intermediate between the fully HS and LS states. When the phase transition is continuous, the difference between the sites is weak and reaches zero at the phase transition, and the signal related to symmetry breaking is therefore more difficult to observe. 30 As shown in Figs. 10 and11, 𝛾 may evolve with temperature on approaching the INT phase, above 𝑇 𝑆𝐵1 or below 𝑇 𝑆𝐵2 , as well as in the INT phase. The discontinuous change of  at 𝑇 𝑆𝐵1 and 𝑇 𝑆𝐵2 is accompanied by a discontinuous change of 𝛾 due to the coupling. The high temperature step is less discontinuous than the low temperature step due to entropy. Various scenarios may appear with i) two strongly discontinuous phase transitions, ii) a half spin crossover at high temperature with a continuous symmetry change and a discontinuous phase transition at low temperature iii) two continuous symmetry-breaking phase transitions with a continuous change of q. Chernyshov described a similar case, with two discontinuous symmetry-breaking phase transitions, without thermal hysteresis. 35 Since the HS fraction 𝛾 differs on the sub-lattices, the long-range HS-LS ordered structure forms the so-called spin-state concentration waves (SSCW), 27,28,82 which is a meaningful way to describe the spin state order. As temperature decreases, the average HS fraction 𝛾 decreases from 1 to 0 (q evolves from 1 to -1). In the symmetry-broken INT phase there is a spatial modulation of the HS fraction over the different sites around the average value 𝛾 (related to an average q) and the amplitude of the SSCW is related to the symmetry-breaking order parameter . The symmetry breaking due to HS-LS ordering relaxes structural constrains and may therefore be accompanied by rotation, torsion or motion of the molecules, counter-ions or solvent.

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS

VII Stepwise and reconstructive spin transition:

The [Fe(3PhOH-trz)2Pt(CN)4]•2H2O system exhibits an unusual behavior associated with a half spin conversion from the fully HS phase (Pn) to a long-range ordered HS-LS phase (Pc). 65 However, there is no group-subgroup relationship between the two space groups, as the HS and HS-LS phases share the same a, b, c lattice. The symmetry change comes from the fact that the four sites in the unit cell split in two sub-lattices (Figure 12). In the Pn HS phase, for each sub- identify thirty six different spin-state ordered phases. 87 Boukheddaden used an elastic model, taking into account both volume change upon spin transition and elastic frustration to explain the existence of organized spatial modulations of HS state inside the plateau of two-step spin transitions. 31 The stabilization of the steps is balanced by the competition between short-range and long-range interactions. 30 The [Fe II (bipydz)(Au I (CN)2)2]⋅4(EtOH) systems illustrates perfectly this concept of multi-step SSCW. This 3D Hofmann-like material exhibits an hysteretic four-step spin conversion, between the high symmetry HS phase ( = 1) and LS phase ( = 0) with three intermediate phases corresponding to the steps  = , with a long-range periodic ordering of spin states. 86 This situation corresponds to the steps shown in red in Figure 13. However, each step has its own lattice periodicity R. Compared to the HS phase with periodicity r, the periodicity is 3r for the steps at  = . Therefore there is no group-subgroup relationship between the different phases on the steps and the phase transition between each step are reconstructive. Consequently, in this case the phase transition between the steps are limited to first-order, except for the first and last steps with the high symmetry phases. However, the concept of order parameter also applies for measuring the deviation from the HShs phase. 27,30 The description of aperiodic SSCW, with a spatial modulation of the HS fraction incommensurate with the lattice periodicity, remains more challenging, even though the concept of order parameter also applies to the super space approach. 27,88 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS

Conclusion

The Landau theory, with proper identification and use of symmetry-breaking and non-symmetrybreaking order parameters is very useful to describe phase transitions in multi-stable materials.

Indeed, non-symmetry-breaking instabilities related to the OP q, monitoring change of spin state or CT, and a symmetry-breaking instability related to an OP , monitoring ferroelastic, ferromagnetic, ferroelastic or structural orders, may occur simultaneously, or sequentially. The phenomenological approach presented here, discussing the coupling between the order parameters, provides the necessary formalism for describing and disentangling spin transition and symmetry breaking. The diverse spin transition curves (T), obtained for the different types of potentials, map most of the scenarios found in spin-transition materials. These include spin-crossover and spin transition, coupled or not to continuous or discontinuous symmetry-breaking, with or without intermediate step.

We can also found in the literature situations that are more complex than the ones listed in the different sections above. Some of them correspond to a mixing of the scenarios discussed in this symmetry one (LSls) with a symmetry breaking different from the one occurring on the step. The HShs and LSls phases correspond then to the linear-quadratic coupling scenario between spin transition and symmetry-breaking order parameters. The potential to be used to map these mixed scenarios is more complex, since different symmetry-breaking order parameters have to be used: one for describing the HS-LS order and one for the structural distortion in the LSls phase. The [Fe(salpm)2]ClO4•0.5EtOH complex corresponds to another mixed case. It is made of two sublattices in the HShs phase. On cooling, one sub-lattice remains HS and the other one undergoes a half spin conversion with HS-LS symmetry-breaking order. At lower temperature, both sublattices switch to the LS state. This scenario is a mixture of a stepwise symmetry breaking spin transition for one sub-lattice and a spin-transition on the other. Here again, the potential required to map this mixed scenario requires including a coupling between the sub-lattices to explain their simultaneous spin transition at lower temperature.

The phenomenological approach discussed in this paper allows for a rigorous description and understanding of spin transition materials, by explaining how to disentangle the non-symmetrybreaking (spin state, charge-transfer…) and the symmetry-breaking (ferromagnetic, ferroelectric, ferroelastic, structural order…) instabilities. The elastic coupling between both types of instabilities can enhance the hysteresis width, drive simultaneous spin transition and symmetrybreaking and generate stepwise spin transitions. The coupling between symmetry-breaking and non-symmetry-breaking instabilities is also an important route for enhancing the (multi-) functionalities. For example, charge transfer can enhance the electric dipole, which, combined to ferroelectric order, enhances the electric polarization. Spin transition and charge transfer can enhance the magnetic dipole, which combined to ferromagnetic order enhances the magnetization.

Mixing symmetry-breaking and non-symmetry-breaking instabilities is therefore an important route for understanding, controlling and developing multi-functional materials.

1 -𝑞 2 )

 12 to be zero in the high temperature HS phase.The contribution of the volume strain 𝑣 𝑠 , or volume contraction, to the potential is included by finding its equilibrium value that minimizes (2):𝑣 𝑠 = -

𝑇 2 .

 2 stable solutions at T=TST ; 𝑞 = ± 𝐵 𝐶 , while q=0 is unstable. The limits of stability of q>0 and q<0 correspond to temperatures were the first and the second derivatives of the potential are 0: 𝑇 𝑢 = The thermal hysteresis, inherent to first order spin transitions, is characterized by its width ∆𝑇 𝐶𝑇 = 𝑇 𝑢 -𝑇 𝑑  4𝐶(

Figure 1

 1 Figure 1 Evolution with temperature of the thermodynamical potential F(q,T) (4) and the equilibrium value of q for B>0 (a) and B<0 (b).

Figure 2

 2 Figure 2 Spin conversion curves q(T) describing a first-order spin transition for B>0 (a) and a spin crossover for B<0 (b).

Figure 3

 3 Figure 3 Evolution with temperature of the thermodynamical potentials (7) F and the equilibrium value of  for b>0 (a) and b<0 (b).

Figure 4

 4 Figure 4 Symmetry breaking curve  2 (T) extracted from potential (7) describing a second order phase transition for b>0 (a) and a first-order phase transition for b<0 (b).

  , with high symmetry (hs, 𝜂 = 0) or low symmetry (ls, 𝜂 ≠ 0): HShs, HSls, LShs or LSls. Usually the high temperature phase is High Spin high symmetry (HShs), because high temperature favors both the high entropy HS state and the disordered state (high symmetry).

Figure 5

 5 Figure 5 Phase diagrams in the in (𝑇 𝑆𝐵 -𝑇 𝑆𝑇 , T) space exhibiting a spin transition (B<0) and the discontinuous symmetry-breaking (b<0) for D=2. The spin transition occurs around T=TST between low and high spin (LS and HS) states, while the symmetry-breaking temperature occurs at TSB (diagonal) between low and high symmetry (ls and hs) states. The dark thick lines represent the thermal hysteresis and the doted lines are the temperature path shown in Figure 6.

Figure

  Figure 6a corresponds to discontinuous spin transition (ST) and symmetry breaking (SB) with

Figure 6

 6 Figure 6 Spin transition curves (T), also scaled in q(T), and symmetry breaking curves  2 (T). (a), (b), (c) and (d) correspond to the coupling of discontinuous (b<0, B<0) or continuous (b>0, B>0) spin conversion and symmetry breaking. Each panel shows the curves for 𝑇 𝑆𝐵 < 𝑇 𝑆𝑇 , 𝑇 𝑆𝐵 = 𝑇 𝑆𝑇 and 𝑇 𝑆𝐵 > 𝑇 𝑆𝑇 .

Figure 7

 7 Figure 7 Representation of the sequence of phases with a non-symmetry-breaking spin transition at TST from HShs to LShs and a symmetry-breaking phase transition at TSB from LShs to LSls. A spin transition occurs with mainly HS molecules above and LS below. The symmetry breaking may have various nature, like cell doubling or ferroelastic distortion, and results in atomic motions ∆𝑟 ⃗, and motions or torsions of ligand, solvent or counter ions.

Figure 6c corresponds to

  Figure 6c corresponds to a continuous symmetry breaking and a continuous spin crossover (b>0, B>0) and in this case the effect of the coupling is not pronounced.

Figure

  Figure6dcorresponds to a case, where a discontinuous symmetry breaking phase transition couples to a spin-crossover (b<0, B>0). In the middle panel both phenomena occur simultaneously due to the coupling, as observed for [Co(dpzca)2].50 This system exhibits an abrupt and complete

  Since the sub-lattices may be associated with different thermodynamical parameters (ligand fields, entropy …), 𝑞 1 may undergo a spin crossover (𝐵 1 > 0) or a spin transition (𝐵 1 < 0) around 𝑇 𝑆𝑇1 and 𝑞 2 may undergo a spin crossover (𝐵 2 > 0) or a spin transition (𝐵 2 < 0) around 𝑇 𝑆𝑇2 . The different scenarios shown in Figures 8a-c provide more or less pronounced plateaus at half spinconversion 𝛾 = 1 2

Figure 9

 9 Figure 9 Sequence of phases for stepwise spin transitions. a) Non-symmetry-breaking stepwise spin transition, found in [Fe(btr)3](ClO4)2 for example, with two sub-lattices (1 square -2 circle) in the unit cell (rectangle) coupled elastically (D12). b) Symmetry-breaking stepwise phase transitions, as observed for [Fe(2-pic)3]Cl2•EtOH or [FeH2L 2-Me ][PF6]2 for example, where the symmetry-equivalent sites in the HShs phase are no more equivalent in the INT ls phase, with a doubling of the unit cell for example.

( 2 ,

 2 INT) phase corresponding to a step at half conversion (q=0) and a fully LS phase at low temperature. The INT phase is associated with a long-range and periodic order of molecules in LS and HS states as observed in [Fe(2-pic)3]Cl2•EtOH 71 or [FeH2L 2-Me ][PF6]2 for example.30 However, on the contrary to the case of non-symmetry-breaking stepwise spin transition discussed in section IV with different sub-lattices, in the case of symmetry-breaking stepwise spin transition the mainly HS and LS sites appearing in the INT phase are initially equivalent in the fully HShs phase, and also sometimes in the LS phase as shown in Figure9b.71,81 Therefore, in the INT phase, where the HS fraction is 𝛾 ≈ 1 a symmetry breaking occurs, forming two sub-lattices, one being mainly HS the other mainly LS. The symmetry-breaking order parameter 𝜂 belongs to the active representation of the symmetry group driving the symmetry change. Its amplitude describes the amplitude of the difference between the molecular sites forming the sub-lattices, including the different HS fraction:  = 1 2

Figure 11

 11 Figure 11 Schematic representation of the average HS fraction 𝛾 and symmetry-breaking order  at different temperatures (a-g) during the symmetry-breaking stepwise phase transition.

  lattice one site is equivalent to another one by the n glide plane. Both sub-lattices are HS. In the Pc HS-LS phase, the sites are no more equivalent within the same sub-lattice of the HS phase, one site being mainly HS and the other one mainly LS, as the n plane is lost. However, the two HS sites from the different sub-lattices of the Pn phase, initially inequivalent in the HS phase, become equivalent by the c glide plane in the Pc phase. Therefore, the two LS sites initially inequivalent in the HS phase also become equivalent by the c glide plane. Such a structural reorganization is reconstructive and the Landau theory does not apply. Indeed, in this case, there is no symmetry This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0041453 deviation from a high symmetry phase, there are just two different phases with different symmetries. This HS-LS phase is stable down to 10 K or less. Reverse-LIESST allows reaching a fully LS state, which as the same Pc space group of the [HS-LS] phase, but both Fe sites are now LS. This reverse-LIESST, without symmetry change, from a phase with HS and LS sub-lattices, is then analogous to the reverse-LIESST from the HS-LS phase found in [Fe(trzpy)2Pt(CN)4]•3H2O.64

Figure 12

 12 Figure 12 The reconstructive phase transition in [Fe(3PhOH-trz)2Pt(CN)4]•2H2O. In the HS phase (150 K) there are two Fe1 and Fe2 sub-lattices and the Fe centers of each sub lattice are equivalent by the n glide plane. Different sub lattices form in the HS-LS phase (100 K) one LS and HS and the Fe centers of each sub lattice are equivalent by the c glide plane. Reproduced with permission from Angew. Chem. Int. Ed. 58, 11811 (2019). Copyright 2019 John Wiley and Sons.

2 3( 3 ( 2 the 2 3

 2322 + +sequence) and  = 1 + -sequence), while on the step at  = 1 periodicity is 2r (+ -+sequence). From the point of view of translation symmetry, the phase formed by the SSCW at  = groups of the HS phase. However the translation symmetry 3r allowed on the steps =

Figure 13

 13 Figure 13 Devil's staircase for SCO materials, the steps with T correspond to fractions 𝛾 = 𝑛 𝑚 with long-range ordering of LS (-)and HS (+) sites. The modulation vector R=nr and there are n HS sites among m sites. Structure of the simplest (thick segments) steps are represented. Reproduced with permission from Angew. Chem. Int. Ed. 55, 8675 (2016). Copyright 2016 John Wiley and Sons.

  paper. For example, in the case of the [FeH2L 2-Me ][PF6]2 and [FeH2L2-Me ][AsF6]2 systems,67,68 an INT phase associated with HS-LS order forms on cooling from the HShs phase, which correspond to the stepwise and symmetry breaking spin transition scenario. However, the LS phase is a low This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0041453

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0041453

Acknowledgement

The authors gratefully acknowledge Agence Nationale de la Recherche for financial support undergrant ANR-16-CE30-0018, ANR-19-CE30-0004 and University Rennes 1 for funding. We thank C. Ecolivet, H. Cailleau and S. Miyashita for scientific discussions.

Availability of data

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0041453