Elastic Origin of the Unsymmetrical Thermal Hysteresis in Spin Crossover Materials: Evidence of Symmetry Breaking

The jungle of experimental behaviors of spin-crossover materials contains a tremendous number of unexpected behaviors, among which, the unsymmetrical hysteresis loops having different shapes on heating and cooling, that we often encounter in literature. Excluding an extra effect of crystallographic phase transitions, we study here these phenomena from the point of view of elastic modeling and we demonstrate that a simple model accounting for the bond lengths misfits between the high-spin and low-spin states is sufficient to describe the situation of unsymmetrical hysteresis showing plateaus at the transition only on cooling or on heating branches. The idea behind this effect relates to the existence of a discriminant elastic frustration in the lattice, which expresses only along the high-spin to low-spin transition or in the opposite side. The obtained two-step transitions showed characteristics of self-organization of the spin states under the form of stripes, which we explain as an emergence process of antagonist directional elastic interactions inside the lattice. The analysis of the spin state transformation inside the plateau on cooling in terms of two sublattices demonstrated that the elastic-driven self-organization of the spin states is accompanied with a symmetry breaking.


Introduction
The spin-crossover transition (SCO) is a phenomenon that occurs in transition metal coordination compounds based on Fe (II), Fe (III), Mn (II), Mn (III), Co (III), Cr (II) [1][2][3][4][5] which change their spin state due to external stimuli such as a variation of temperature [6][7][8], pressure [9][10][11], light irradiation [12][13][14][15][16], electrical or magnetic fields [17,18] etc. In the literature, this phenomenon is most commonly studied on iron (II), with 3d 6 electronic configuration. From the point view of the ligand field theory [19], the central Fe(II) ion surrounded by nitrogen atoms in an octahedral environment may have two spin configurations, namely low-spin (LS) or high-spin (HS). According to the strength, strong or weak, of the ligand field, compared to the electronic pairing energy of the electrons in the metal, one can stabilize respectively, the LS or the HS electronic configuration. From this theory, the change in the internal distribution of electrons inside the atomic orbitals of the central metal causes the conversion from a paramagnetic high-spin ( , , = 2) stable at high temperature to a diamagnetic low-spin ( , , = 0) which becomes the ground state at low-temperature. In general, the transition from LS to HS is accompanied by an increase of the Fe-N bond lengths between ~2.00 Å in the LS state and ~2.2 Å in the HS state [20][21][22], leading to sudden expansion of the molecular volume, that is the volume of the coordination sphere, by around 30% [22].
Historically, the thermally induced change of spin state was first observed by Cambi and co-workers [23] in a spin transition compound based on iron (III) and has since been extensively studied in iron (II) complexes. It was found that the change of spin state induced in the spin-crossover compounds by the temperature occurs when the difference in the Gibb's free energies between the two spin states involved is on the order of thermal energy, [24]. Indeed, when the temperature increases, the entropy (magnetic and vibrational) of the HS state increases as a result of its nonzero spin value and its soft character, which leads to its thermodynamic stabilization, while on lowering the temperature, the entropy decreases and the diamagnetic LS state which has a rigid character is stabilized by the ligand field which becomes stronger. Overall, the thermodynamic stabilization of these two states is also due to the lattice distortion which accompanies the spin transition. On the other hand, spin-crossover materials are known for their specific property of showing long-lived photo-excited states at low-temperature, known as LIESST (Light-Induced Spin State Trapping) effect, which is obtained after photoexcitation of the LS state with a particular wavelength (usually green light) [25]. When trapped in this metastable state, the SCO materials may return to their stable LS through a thermally-induced relaxation mechanism subsequent to an increase of temperature.
In cooperative spin-crossover systems, these local volume expansions take place at several points in the lattice and their interference cause a global volume expansion of the whole lattice by ~3-5% [26] which also corresponds to the volume expansion of the unit cell. This value is very small compared to the 30% of the volume expansion of the molecular coordination sphere, which means that a large part of the molecular volume expansion is absorbed inside the lattice in the form of reorientations of the ligands and local degrees of freedom which do not affect the unit cell volume. Many experimental studies [27][28][29][30][31][32] were devoted to this phenomenon, which appears in organometallic complexes, due to the rich variety of collective thermodynamic behaviors such as: (i) gradual spintransitions [33] of independent SCO molecules, described by a simple Boltzmann distribution over all energy levels involved in weak cooperative systems, (ii) first-order transitions [21] with large thermal hysteresis, (iii) incomplete spin transitions with significant fractions of molecules blocked in the HS state at low-temperature [34] (iv) as well as twostep or multi-step transitions [35,36] which often express the existence of antagonist interactions inside the lattice. Furthermore, several experimental investigations conducted in some spin-crossover complexes have shown the presence of two-or multi-steps transition states which appeared only on cooling or on heating leading to an unsymmetrical transition [37][38][39][40]. From the chemistry point of view, the richness of these behaviors is obtained by several chemical changes, like monitoring the nature, the size or the coordination mode of the ligand or the anions, which affect the crystal packing as well as the mechanical strains, resulting in huge changes of the macroscopic magnetic and mechanical properties of the bistable molecular solids. Many models have been proposed to describe the thermal properties of SCO materials. Some of them were based on classical thermodynamic approaches [41,42], or continuum mechanics modeling where the interactions have considered in mean-field approximation [43]. More recent models used first two-state Ising model to mimic the spin transition, in which the interaction although elastic in nature was written under a form of an exchange term [44,45]. The effect of the interaction between the spin state of the molecule and the lattice was introduced latter by considering spin phonons interactions [46,47] or vibronic approaches where the spin and local molecular distortion are coupled in a quantum way [48]. Later, new type of models taking into account for the lattice volume change at the transition appeared. They are called mechanoelastic [49][50][51][52][53], anharmonic [54,55] or electro-elastic [56] and have been proposed, for some of them, to explain the experimental observations of optical microscopy investigations [57] which showed that on single crystals, the spin transition proceeds via a single domain nucleation which appears at the crystal edge/corner and then propagates over the whole lattice, preventing any other nucleation for opposite sides. In most of these elastic models, the spin-crossover molecules are modeled as atoms with a stable LS state at low-temperature and a stable HS state at high-temperature, coupled via harmonic or anharmonic springs whose elastic constant may depend on the spin states of the connected spins [56]. In particular, the case of multi-step conversion has attracted a lot of the attention due to the possible application of the materials showing this property for multi-byte information storage. From the theoretical point of view, they are described using Ising-like models [44,45,56,[58][59][60][61][62] combing short-range antiferromagnetic-like interactions and long-range ferromagnetic ones; both of these interactions are introduced in the phenomenological way.
The present theoretical work, based on the extension of our electro-elastic model to investigate the case of unsymmetrical thermal hysteresis, focusses on the effect of the lattice parameter misfit on the thermal and bistability properties of SCO solids. Thus, we here consider the general case of the genuine electro-elastic model [58], and investigate the effects of the lattice misfit parameter between LS and HS states on the thermal behavior of the HS fraction and that of the average lattice bond lengths. The model is solved in its general form, with elastic constants depending on the spin states of the linked SCO sites. We observed that the tuning of the lattice parameters and elastic constants between HS − HS, HS − LS and LS − LS configurations leads to a variety of thermal behaviors of the HS fraction and average lattice parameters. Among them, hysteretic first-order transitions are recovered [58], as well as the appearance of an intermediate plateau on cooling or on heating with the presence of structural and spin state self-organizations in the plateau regions. Moreover, unsymmetrical thermal hysteresis loops and incomplete transitions are also obtained, which allows to propose a physical mechanism for the appearance of these behaviors in the experiment.
The manuscript is organized as follows: in Section 2, introduces the presently used electro-elastic model to describe spin-crossover materials and demonstrates its equivalence with an Ising-like Hamiltonian including an infinite long-range interactions term. Section 3 is dedicated to the Monte Carlo resolution of the Hamiltonian, the presentation and the discussion of the obtained results. In Section 4, we conclude and explore some extensions of the current investigations.

The Model
The electro-elastic model describing the spin transition phenomenon is considered here for a 2D square lattice of size, × , schematized in Figure 1, in which each site, representing a SCO molecule, is associated with a fictitious spin state, = ±1, where = +1 refers to the high-spin (HS) state and = −1 to the low-spin (LS) state.
To take into account for the local variation of the volume accompanying the spin change, the sites are linked by harmonic springs [56], whose elastic constants and equilibrium distances between the nearest neighbors (nn) and next-nearest neighbors (nnn) sites depend on their connected spin states. The electro-elastic Hamiltonian of such system is given as The first term of Equation (1) represents the effective temperature-dependent energy gap between the HS and LS states. It contains the ligand field energy contribution, Δ, arising from the nitrogen atoms surrounding the SCO metal center, and the entropic contribution, ln , originating from the electronic and vibrational degeneracy ratio, = , between the HS and LS states. The second and the third terms of (1) are the elastic interaction energies between the nn and nnn of spin-crossover units, respectively. Here, and are the corresponding elastic constants and, and are the instantaneous distances, while, and are their associated equilibrium distances.
In the similar way, we consider different rigidities for the LS and HS states, which leads to elastic constants between nn (resp., nnn) atoms depending on their spin states. They will be written, , and (resp., , and ) for the HS − HS, HS − LS and LS − LS configurations, respectively. For simplicity reasons, and according to the 2D square topology of the lattice, the equilibrium nnn distances in the three configurations are taken as = √2, = √2 and = √2. It is straightforward to prove that the general expressions of the nn and nnn equilibrium bond lengths, ( , ) and (S , S )), and those of the elastic constants, ( , ) and ( , ) satisfying the previous constraints, can be written as function of the spin states, and under the following forms: where: where, = − , is the lattice misfit between the HS and LS phases. Similarly, the nn ( ) and nnn ( ) elastic constants are written as follow: where, the parameters ( = 0, 1, 2) and ( = 0, 1, 2) write under the general form: It is interesting to notice that in the general expressions (2a, 2b) of the nn equilibrium bond length, , , represents the average lattice parameter, which appears as a center of mass between the average lattice spacing of the four spin configurations (HH, HL, LH, LL) of the connected nn sites. On the other hand, , relates to the bond length misfit between HS and LS states, and stands for the gap between the lattice parameter of HL configuration and that of the average lattice parameter, , between HS and LS phases.
Inserting the expressions of the lattice bond lengths (2a) and the elastic constants (3a) in Hamiltonian (1), and after some simple algebra, it is easy to demonstrate that (1) can be mapped in the form of the following Ising-like model, where, denotes the cohesion elastic energy of the lattice, whose expression is: where, = ( + 2 ) + ( + 2 ) + ( + 2 ) + 2( + 2 ) is a constant contribution. In Equation (4a), and are the local nn and nnn exchange-like interactions, ℎ is a local effective field. All these quantities include the contributions of nn and nnn elastic interactions, and so they couple the spin states to the local distortions of the lattice. In contrast, expresses the pure elastic energy of the system which does not depend on the spin states. This contribution is identified here as the cohesion energy.
The expressions of the effective nn and nnn local exchange-like parameters, and , are given in Equations (5a) and (5b). They show linear and quadratic contributions of the lattice elastic field, and, where, Overall, one can conclude that the local exchange-like interactions, and , contain constant short-range terms, and , and long-range elastic contributions, , and , , recognized through their dependences on the lattice positions, whose sign depends on several parameters, namely , , , and . The local field-like term, ℎ , writes as follows, where, Here, (resp., ) runs over the nn (resp., nnn) of a given site and = 4 is the lattice coordination number of the 2D square model considered in this study. Equations (7a)-(7d) also show that the initial electronic ligand field, (∆ − ln ), acting on each site is renormalized by the elastic contributions of the neighbors in a complex way whose signs depend on the model parameters.
It is interesting to remark in the mapped Ising-like Hamiltonian (4a), the presence of intrinsic frustration resulting from the exchange-like contributions, and , along the nn and nnn lattice parameters. This intrinsic elastic frustration is easily evidenced in the special case, where the elastic constants are independent on the spin states, i.e., = = 0 and = = 0. Here, the nn and nnn exchange interactions become = and = 2 and both of them are positive. In this case, the short-range couplings and , stabilize an antiferroelastic order in the lattice, while the local effective field, ℎ , which contains a strong long-range ferro-type interaction (in particular for = 0 or = ), wants to impose a ferroelastic order. It results a situation of elastic frustration in the lattice which is not produced by geometrical effects but rather due to competing elastic interactions that may lead to special self-organization of the spin states, driven by the lattice effects. On the other hand, it is quite well known that an Ising model with antiferro-like nn and nnn interactions already leads to frustration due to the antagonism of these two types of interactions [63][64][65][66]. For a better understanding of the role of the local interaction parameters on the possible emergence of self-organization of spin states, so as to predict in which conditions they stabilize or destabilize the LS and HS phases, for the benefit of other spatially modulated states, we perform in Section 3 MC simulations for several values. The latter are varied between 1.0 and 1.2 nm, which are the values of the lattice spacing in LS and HS states, respectively.

Model Parameters and Monte Carlo Procedure
In this study of the thermal properties of the spin crossover materials, we consider a 2D square lattice of size × = 50 × 50, with free boundary conditions in order to account for the global lattice deformation arising from local volume expansions/contractions consequent to spin flips. We assume for simplicity that the lattice deformations remain inside the plane, and the simulations do not explicitly involve external pressure effects. Table 1 summarizes the used values of the equilibrium lattice parameters and elastic constants in the numerical simulations, except that of which is monitored along the numerical investigations. Considering the following expressions: R = for the average equilibrium distance between HS and LS phases, δ = R − R , for the lattice misfit, the previous ρ , ρ and ρ parameters given in Equation (2b) can be written as, A schematic configuration of these different bond lengths is given in Figure 2. The present work, is dedicated to the investigations of the thermal properties of the electro-elastic 2D lattice by tuning the values of lattice spacing, , between and . Three situations of values, will be studied. They correspond to (i) = , which means that the HS-LS configuration plays a neutral role with respect to the HS and LS ones; (ii) located in the interval , and (iii) located in the interval [ , [.
For that, Monte Carlo technique running on spin and lattice position variables are implemented in a sequential way. For both spin and lattice position degrees of freedom, a Metropolis algorithm is implemented to drive the spin flips from toand the lattice displacements from the initial position ⃗ to ⃗ + ⃗ , where ⃗ is the elementary displacement chosen in a random direction with a random amplitude ‖ ⃗ ‖ chosen so that ‖ ⃗ ‖ ≪ , which is the smallest lattice parameter value used in the simulations. Here, we used the value, ‖ ⃗ ‖ = 0.05 nm. In this problem, where each site, , has two degrees of freedom, and ⃗ , the following double MC procedure is performed. First, we call randomly a site and update its spin state following the Metropolis algorithm. Whatever the result, accepted or rejected spin flip, we perform another MC process by displacing each lattice position randomly in order to minimize the system's elastic energy. We update all lattice positions by MC several times (here 10 times) in order to reach the mechanical equilibrium. Then, we call randomly another spin site and we repeat the procedure until visiting all spin sites. When all spin sites have been visited once, we call this 1 MC step (MCS). Thus, for a lattice of × (here 50 × 50) sites, in 1 MCS, each lattice position is updated 10 × times, which makes the procedure highly time consuming. Along the simulation, and depending of the studied situation, the system is cooled down from 200 to 1 K and then warmed up until 200 K, or inversely. In this range of temperatures, the HS fraction, , the average bond length distance, < >, as well as the spatial spin and position configurations are recorded every 1 K after a waiting time corresponding to 5000 Monte Carlo steps in order to reach the thermodynamic equilibrium. Among these 5000 MC steps, 2000 MCS are used to evaluate the average quantities of interest. So at the end of the process, and when the total elastic energy reaches its stationary state, we perform all averages of physical quantities of interest. In all simulations, the ligand field energy is maintained constant, = 450 K and the value of degeneracy ratio between the HS and the LS states is taken as, = 150 (ln = 5). This effective degeneracy value contains the electronic degeneracy 15 of the Fe 2+ SCO complexes in the HS state to which is included the intra-molecular vibrational entropy change at the transition as well as that of the lattice. Taking = 150 leads to an entropy change Δ ~ 40 J//K/Mol which is in fair agreement with available experimental data [63,66]. As a consequence, the equilibrium transition temperature, which corresponds to the temperature which zeros the effective ligand field is expected to be, = = 90 K, when the elastic constant are spinstates independent [58]. In the present case, since all constants (nn and nnn) depend on the spin state, the equilibrium transition temperature is obtained by putting ℎ = 0 in Equation (7b), which leads to , only in the case where = = = = 0.

Results and Discussions
This section is devoted to the presentation of the thermodynamic properties of the studied electro-elastic lattice. We will focus on the shape of the thermal hysteresis and its dependence on the lattice parameter which is considered here as a variable. The nucleation and growth mode of the spin states as well as spatial distribution or self-organization along the thermal hysteresis will be analyzed. Two cases, corresponding to two situations of elastic frustrations, will be considered according to values: (i) ≥ and (ii) ≤ .   Figure 3a show, for = 1.1 and 1.11 nm, the presence of a first-order transition between the HS ( = 1) and the LS ( = 0) states accompanied with a thermal hysteresis. The origin of this thermal hysteresis is clearly due to the elastic interactions acting between the spin states through local volume changes accompanying the spin state switching. These changes of the local volumes deploy over the whole lattice causing longrange effects.

Results of
When increasing the value of beyond, = 1.11 nm, an intermediate plateau appears on cooling from HS to LS while it is absent on heating from LS to HS. For larger values, 1.14 ≤ ≤ 1.2, an incomplete transition takes place on cooling due to the stabilization of intermediate states, and the obtained plateau persists until 1 K, due to the vanishing of thermal fluctuations, leading to a residual HS fraction whose amplitude increases with value. In our simulations, where the MC waiting time is fixed over the whole study, this behavior appears beyond the threshold value, = 1.14 nm.
When inspecting Figure 3a, we remark that the cooling branches of the hysteresis are shifted upward downward when we increase from 1.1 to 1.1336 and then split into two branches with the appearance of a plateau. The upper one continue to shift upward while the lower one shifts downward and disappears for situated in the interval 1.1336 and 1.14. Beyond this value, incomplete spin transitions are obtained because the temperature at which the plateau region starts, reaches zero. Indeed, for = 1.14 nm, the HS fraction do not totally convert at very low temperature and remain blocked at = 0.55 leading to a mixed phase where HS and LS states coexist.
As a result, on heating with the used MC kinetics, the system stays in the plateau and follows a succession of metastable states before to reach the HS state at temperatures (~70 K) well below those of the heating branches of the thermal hysteresis (~150 K). This discontinuity in the change of the value of the heating branch of the thermal hysteresis can be left by starting the "measurements" from the LS state. For that, one can prepare the LS state from the intermediate plateau state at low-temperature by applying pressure [67] which stabilizes the LS state or by light-irradiation through reverse-LIESST (Light-Induced-Spin-State Trapping) effect [68].
In modeling, if one starts the simulation by warming and then cooling, we expect the same results as those of Figure 3a,b for values ranging in the interval 1.0-1.1336, while different trends are expected for ≥ 1.14 nm. Indeed, in this case, the heating process from the LS state leads to an increasing upper transition temperature of the thermal hysteresis (following the same behaviour as those of < 1.14), and then we recover the "right" heating branch, , which is in the continuity of the other heating branches obtained for < 1.14, as well indicated in the phase diagram vs. of Figure 4. In contrast, when the system reaches the HS state, the cooling branch will be the same as that already represented in Figure 3a,b for the HS fraction and the average nn bond length, respectively. Now we turn to the explanation of the global shift of the thermal hysteresis for increasing values. To understand quantitatively this point, one has to consider the energy barrier that faces a spin state during the spin flip and to express it as a function of . However, a qualitative approach provides a first explanation. Starting from the HS state, whose lattice parameter is and cooling down, the appearance of a LS state changes the equilibrium lattice parameter to . If increases towards , then the transition from HS to LS will take place at lower temperature, because the "distance" -decreases. On the other hand, starting from LS with lattice spacing equal , the increase of will increase the "distance" -which will result in an increase of the LS to HS transition temperature. Overall, the total thermal hysteresis shifts to higher temperature regions, as found in Figure 3a,b. In the following section, we will discuss these behaviors observed in relation with elastic properties. , the lower branch of the thermal hysteresis, splits, for ≃ 1.12, into two branches with inflexion points and . While, depicts an increasing behavior with , falls down around, ≃ 1.14, beyond which the system stays quenched in a metastable state, leading to an incomplete spin transition. For, ≥ 1.12, the system depicts self-organized spin states.

Long Range Effect in the Competition between Ferro-and Antiferro-Elastic Interactions on the Thermal Dependence of the HS Transition and Generation of Elastic Frustration
To discuss quantitatively the behavior of the thermal dependence of the HS fraction and the average distance parameter observed in the previous section, we focus on the signs of the various interaction parameters involved in the effective Ising-like Hamiltonian (4a). According to values, the exchange-like interaction terms and , as well as the local effective field contribution, ℎ , combine short-and long-range interactions which can be ferroelastic or antiferroelastic. So the, competition between local ferroelastic and antiferroelastic interactions can be generated by monitoring variable. It is therefore expected that the thermal behavior of the HS fraction will be affected by this parameter, by stabilizing elastically new phases that would never appear for = = .
Let us start with the case, = = 1.1 nm, for which the parameter, = , is systematically equal to zero which means that the local exchange-like interactions, and , do not depend on the atomic positions. Consequently, these two parameters will do not include any long-range effect in this case. In contrast, the ligand field term, ℎ , does contain long-range interaction terms. Thus, in this case, nn and nnn short-range antiferroelastic interactions are identified in local exchange-like interaction terms, which induces an intrinsic frustration in the system. During the thermal transition, we notice a competition between short-range antiferroelastic interactions arising from the exchangelike interaction and long-range ferroelastic interactions resulting from the effective fieldlike contribution, ℎ . However, the strength of these short-range antiferroelastic interactions in this case is very weak, and consequently the thermal dependence of the HS fraction is mainly governed by the effective field-like interaction, ℎ , which combines an effective ligand field energy, ℎ = (∆ − ln ) + z( + 2 ) competing with the additional long-ranged elastic field energy [see Equation (7a)]. The obtained thermal-dependence of the HS fraction and lattice parameter are given in Figure 5a1 When we inspect the spatial mode of the nucleation of the spin states (Figure 5c1) we remark that the LS (resp., HS) states on cooling (resp., heating) grow from the corners in a single domain fashion and then coalesce at the center of the lattice. This behavior is very different from that of the usual Ising model which shows multi-droplet nucleation with ramified spin structures. The present behavior is attributed to the existence of long-range elastic interactions between the spin states, deployed as through the strain field which prevents the nucleation from other sites thanks to the internal pressure produced by the converted phase.
For > , becomes negative and the antiferroelastic long-range interactions in the effective exchange-like interactions, and , and their values will increase with and will participate in the thermal behavior of the HS fraction. As we have already seen, for = 1.1 nm, the strength of these antiferroelastic long-range interactions is null, while for = 1.11 nm, it becomes different from zero but its contribution in the thermal behavior of the HS fraction and average lattice parameter is hardly seen (compare Figure 5a1,b1 and Figure 5a2,b2) since the shape of the thermal hysteresis remains almost unchanged. Despite of the absence of any indication on ( ) and < > ( ), the analysis of the spatiotemporal behavior of the HS fraction along the spin transitions shows important differences compared to the case = 1.1. Indeed, the comparison between Figure 5c1 and Figure 5c2 shows clear differences in the spatial organization of the spin states at the transition. Although, HS stripes appear on cooling in both cases, the latter are more dense in the case = 1.11, which stabilizes this configuration. It is then expected to observe an enhancement of this type of self-organization by increasing, . This point will be discussed in the next section. For the values, = 1.12, 1.123, (Figure 3a,b) the thermal evolution of the HS fraction does show a shift of the thermal hysteresis to higher temperature regions, as already mentioned, but without any signature of plateau. The latter appears beyond the threshold value ∼ 1.126 nm where the lower branch of the thermal hysteresis exhibits two gradual transitions between the HS state and the intermediate state at = 54 K and then to the LS state at = 37 K. When we increase the plateau width is enhanced which pushes the temperature to higher temperature regions and to lower values. At = 1.1336 nm, reaches 0 and one gets an incomplete spin transition. In this region well self-organized structures are expected as we will see in the next section.

Spatial Organization of the Spin States in the Course of the Spin Transition
The thermal evolution of the HS fraction and average lattice parameter are discussed here for the case ≤ ≤ including the spatial organization of the spin states. Figure 5a1 depicts the case = , where we do not expect any self-organization since the average lattice parameter of the HS-LS configuration does not favor neither the HS-HS nor the LS-LS configuration. In this case, usually the expected transition temperature is = = 90 K, and it is situated in the middle of the thermal hysteresis [58].
However, in the present study the elastic constants ( and ) depend on the spin states, which results in a different expression of which depends on and . As a consequence, the transition temperature is shifted with respect to , which is no more in the middle of the thermal hysteresis. However, if one takes = = = = 0 [See Equation (3a)], we recover the behaviors of Ref. [58]. As it clearly appears in the snapshots of Figure 5c1,c2, for = 1.1 and 1.11 nm, the spin transition along the cooling (resp., heating) branch starts from the four corners following the usual scenario of nucleation and propagation of single domains of the stable LS (resp., HS) phase inside the metastable HS (resp., LS) state. After the formation of homogeneous LS domains at the corners, inhomogeneous "multidroplets" nucleation grows around the interfaces. This phenomenon is ascribed to a macroscopic nucleation [69] and is attributed to the propagation of the ferroelastic stresses that prevent nucleation from the center due the high elastic energy cost. Although snapshots Figure 5c1,c2 are very similar, one can remark that there are some significant differences in the organization of the spin states on cooling where the HS strings appear more clearly in Figure 5c2 for which = 1.11. For all parameter values, = 1.126, 1.13, and 1.1336 nm, the nucleation on heating occurs as macroscopic single HS domains starting from the corners and propagating towards the center of the lattice. This transformation is accompanied with a significant lattice deformation. This is the usual domain propagation already reported in previous elastic models devoted to SCO phenomenon. In contrast, on cooling, an emergence of a plateau region, with self-organized spin state structure (labyrinths), takes place. These structures, observed for all above mentioned values are formed through long entangled and ramified HS and LS strings. These self-organizations of the spin states are induced by the elastic-field and are clearly attributed to the competition between long-range ferroelastic interactions [field-like term of Equation (7a)] and antiferro-elastic frustrated short-range interactions emerging from the contribution of the parameter, (< 0) which appears in the local effective exchange-like interactions, and , given in Equations (5a) and (6a). The occurrence of the plateau at the transition clearly results from the presence of antagonist ferro-and antiferro-elastic interactions which produce an inhomogeneous internal pressure between the SCO sites. This effect combines with the lattice shape [70,71] to produce consecutively 1D HS patterns and 1D LS strings.
From the value = 1.14 ( Figure 3), the LS state is no more "stable" at 0 K, within the MC kinetic used in the simulations, and a residual HS fraction appears at low temperature, as illustrated in Figure 3a. This behavior may have a kinetic origin, since all simulations have been launched from the HS state where all distance are set equal to . However, this tendency clearly confirms that increasing stabilizes the HS-LS configurations leading to large plateaus on cooling. Interestingly, we do not obtain any plateau on heating for ≤ ≤ , which then leads to non-symmetric hysteresis (plateau on cooling and sharp transition on heating). This type of behavior was frequently observed in experiment [37][38][39][40], but was never reproduced or modeled. So, here we bring a possible explanation of this unusual behavior, based on the existence of an elastic frustration due to lattice parameter misfit, which manifests itself only during the transformation from HS to LS, while it is silent in the opposite process.
Let us now discuss the self-organization of the spin states in the plateau regions of Figure 6c1-c3. In all these figures, we first remark in panels B of all snapshots that, on cooling, the nucleation of the LS state starts from everywhere in the center of the lattice for weak LS concentrations. A brief inspection of the snapshots B indicates that for weak LS concentration, the LS sites prefer to have HS sites as nn, which means that the system maximizes the numbers of HS-LS bonds. This behavior can be understood from the point of view of elastic energy. Here, we simply consider a system with uniform elastic constants for simplicity.
When starting from the HS state, the spin flip of a HS site to LS costs the elastic energy 4( + 2 )( − ) which decreases as gest closer to .
Let us now increase the number of LS sites, and take for example 2 LS sites in the lattice. We have three possible configurations: (i) they can be far from each other, (ii) also they can be nn or (iii) nnn sites. If the two LS sites are isolated then the total energy is just the double of that of one isolated LS site: 8( + 2 )( − ) . On the other hand, if the two sites are nn, the total elastic energy is 8 × 2 ( − ) + 6 ( − ) + ( − ) . And, finally if the two LS sites are nnn, the elastic energy cost before relaxation is 8 ( − ) + 2 × 6( − ) + 2 ( − ) . Since the quantity − is always bigger than − , it is clear that for weak LS concentrations, the system will prefer isolated LS species. Thus a simple numerical evaluation of these energies within the model parameter values used in the simulations gives, 1290 K for the isolated LS site and 1962 K (resp., 1977 K) for the configuration of two nn (resp., nnn) LS sites.
However, as soon as the concentration of LS species increases (panels C) small LS clusters will appear. According to the previous calculations, the energies ( and ) of LS clusters made of two nn and nnn LS sites are such as, − = (4B − 2A) ( − ) − < 0. Consequently, the system prefers to grow LS nano-strings surrounded by 2 HS nn. This, forces the HS phase to self-organize under the form of ramified stripes.

Evidence of Symmetry Breaking
We could calculate the number of nn HS sites belonging to stripes by scanning the lattice along − (resp., −) direction and counting all HS sites having two HS nn along the horizontal (resp., vertical) direction and two nn LS along the vertical (resp., horizontal) one. So doing, we determine the size distribution of the chains as well as their density. We see in all cooling snapshot panels of Figure 6 that due to the system compactness, the density of HS chains among the total HS fraction, is maximum when the HS phase becomes minority, as in all panels E where also the average chain length becomes maximum. To confirm this behavior we plot in Figure 7 the temperature-dependence of the density of HS chains (whatever their length) along with the HS fraction, . To characterize this specific behavior, where minority species HS (resp., LS) tend to self-organize into chains. We first remark that every HS (resp., LS) site belonging to a HS (resp., LS) chain (except those at the extremities) has two HS and two LS nearest neighbors and two HS neighbors. If we denote by the sum of neighboring spins of the HS (resp., LS) central site in question, then every HS (resp., LS) site belonging to HS (resp., LS) chain fulfills the condition, = 0 , for its surrounding spins, except as said, for sites located at the extremities of the chain. Then, we evaluate numerically the quantities, which represents the respective numbers of HS and LS sites that fulfill the condition = 0. Here, , is the Kronecher delta function, while is the sum of nearest neighbors spins ( = ∑ ) of site .
Moreover, we determine, at each temperature, the densities of HS (resp., LS) sites belonging to HS (resp., LS) chains, defined as, In Equation (10), (resp., ) is the number of bulky SCO molecules or sites in the HS (resp., LS) state. It is worth noticing that when all HS (resp., LS) atoms belong to HS (resp., LS) chains, the quantity, (resp. ), becomes equal to 1. To furthermore enhance the relevance of these results, we introduced an additional constraint on the occupation of HS-HS and LS-LS pairs in the lattice. To do so, we calculate the respective quantities, Here again, when all HS (resp., LS) sites belong to HS (resp., LS) chains, we expect the following relation, = − (resp., = − ) and therefore, for large chains one may consider, ≃ (resp., ≃ ). Thus, the quantities, − (resp., − ) allow to evaluate the probability for the HS (resp., LS) sites to not belong to HS (resp., LS) chains. Figure 7 summarizes the density, = × (resp., = × ) of HS and LS sites fulfilling both conditions (9) and (11). We see in Figure 7, that on heating, and in agreement with the snapshots of Figure 6c3, the lattice does not show any significant amount of HS chains around the transition temperature, = 125 K. In contrast, in the cooling process, the density (red curve) shows a very characteristic behavior, since it becomes maximum around the entire plateau region and falls down in other regions. Indeed, it suddenly increases below 75 K until 50 K, due to the increase of the LS fraction, which induces the appearance of HS-LS species. As a result, the quantity decreases (that is the number of HS molecules) while raises up due to the transformation of HS sites to LS sites. Consequently, = is a strongly increasing function in this region. On the other hand, decreases below 75 K, following the trend of the HS fraction (black curve). Overall, the density, , is mainly driven by the change of , to which it is slaved, and therefore increases in the temperature range 75 − 50 K and reaches a maximum value, equal to , around 25 K. Below 15 K, becomes equal to 1, and -which drastically decreases-controls the behavior of the density . The density, , of LS sites belonging to LS chains, represented by the blue curve in Figure 7, shows very similar behavior as that of . However, we remark two different aspects: (i) starts to increase at a lower temperature (~ 60 K), because in the region 75 K-60 K, the system grows an antiferro-like phase (see panel B of Figure 6c3). Below 60 K, the HS state is already self-organized in the form of stripes, which forces the emergence of inclusions of LS stripes. Here, in all this region is driven by the change of , which reaches its maximum value around = 25 K. Then, below 15 K, the stripped LS inclusions grow up and form LS compact domains, which causes the increase of (whose maximum value is equal to 1) and vanishing of . Overall, also decreases in this second part of the plateau.
On the other hand, to demonstrate the occurrence of a symmetry breaking-like behavior in the plateau region, we divided the lattice into alternate horizontal and vertical A and B sublattices to which we associate the respective couples of HS fractions order parameters ( , ) and ( , ). Thus, for a total HS fraction, = 0.5, and perfectly ordered horizontal (resp., vertical) alternate HS, LS chains, one expects the values = 1 and = 0 (resp. = 1, = 0) and = = 0.5 (resp. = = 0.5). Figure 8, displays the thermal dependence of the order parameters associated with sublattices A, B, together with the total HS fraction. While in the HS and LS phases the sublattices A and B are equivalent, we clearly see the occurrence of a bifurcation around the plateau region, with the presence of partially-ordered phase, where and departure from the curve of the HS fraction. This behavior is a signature of the existence of symmetry breaking. To well emphasize this character, we plotted in the inset of Figure 8, the thermal behavior of ( − + 0.5) and ( − + 0.5) around the plateau region, which shows clearly features of symmetry breaking, through the equivalence of the two sublattices on both sides of the plateau region. The main reasons for obtaining a partial order is here attributed to "kinetic character" of the plateau region, in which the derivative is different from zero, which means that the obtained structures do not have a long lifetime. Next, we investigate the kinetic nature of the obtained plateaus in the simulations, by performing isothermal MC simulations at = 27 K for = 1.1336, corresponding to point D of Figure 6c3. For that, we performed several runs of MC simulations, relaxing the spin states and the lattice positions by starting from the HS state (all spins equal to +1 and all bond lengths equal to ). The obtained results at = 27 K are first summarized in Figure 9a which displays the temporal evolution of the HS fraction, = (blue curve) and the two points nn correlation = + (red curve) along the relaxation process. First of all, we see that the system switches from HS to a full LS state, which means that at 27 K, the LS state is the stable phase. In addition, the shape of both relaxation curves depicts the existence of three regimes along the transition process, with the presence of a clear plateau (between points B and E) and finally a sigmoidal relaxation to the LS state. The first rapid regime, between points A and B on the HS and HS-LS fraction curves, corresponds to the relaxation from the unstable HS state to the "metastable" intermediate state (the plateau). During this regime the HS fraction switches from 1 to ~0.6, while the parameter behaves in an explosive way and almost reaches its maximum value in a very short time, indicating that this first regime is mainly controlled by the short-range correlations. Then follows a slow relaxation regime from this states towards a plateau regime for both order parameters. This intermediate regime is found for other relaxations performed at different temperatures, as displayed in Figure 9b. This transient phase shows a very characteristic self-organization of the spin states under the form of interpenetrated stripes. This fact is supported by the corresponding snapshots of Figure 9c depicting the spins organization along the relaxation process. The self-organization concerns first the LS phase and starts around ~0.6 leading to highly ordered small LS chains completely surrounded by HS sites. The fact that these configurations are less disordered than those of Figure 6c3 is due to low value of temperature (27 K) at which they are formed, and for which less fluctuations perturb the system, in comparison with those of Figure 6c3. are considered here as the relevant system's order parameters. Then we plotted in Figure 10 the system's behavior in the phase space ( ( ), ( )). It is remarkable that all these flow diagrams lead to the same trajectory of the system in its phase space, leading to a single universal behavior, indicating that all these relaxations go through similar self-organized structures.
Moreover, we reported in Figure 10 by the red curve the equilibrium path ( ( ), ( )), obtained from the thermal-dependence of the cooling branch of Figure  6a3. We see that except the case = 10 K, which departures in the first regime, all other curves of the flow diagram ( ( ), ( )) started form the HS state at different temperatures clearly overlap and match with that of the thermal transition. This means that during the relaxation, the systems follows adiabatically the local equilibrium behavior, thus confirming again the strong attractor of the striped configuration.
On the other hand, if one considers a random distribution (or mean-field) of the spin states, then the simple relation connecting ( ) and ( ) is = (1 − ). In the phase space, ( , ) this leads to the parabola given by the black dashed line in Figure   10 which undoubtedly results in a different trajectory. A meticulous inspection of the mean-field curve shows that the latter matches quite well with that of the various system's trajectories obtained with MC simulations. The agreement is quite good in the region ∈ [1, 0.85], which corresponds to low values of LS fractions, which means that in this region the LS species appear randomly in the lattice. However, beyond this point, one sees clearly the presence of a departure between the two sets of curves due to the existence of short-range correlations between the LS species which start to self-organize in the form of strings. This disagreement persists and becomes more and more evident as the relaxation process advances. It reaches a maximum in the plateau region, corresponding to interval values ∈ [0.2, 0.6 , where the system self-organizes. In the second part of these investigations, we check the robustness of the attractor found in Figure 10. For that, we prepare the system at = 27 K in different initial spin configurations, corresponding to different values of ( = 0) and ( = 0) , with bond length distances equal to those of HS and we relax it from these states to the LS state. Several situations were considered, like mixed phases of HS and LS domains with for each situation, Monte Carlo simulations are performed on both spin states and positions until reaching the unique stable LS state ( = 0 , = 0). Again, as summarized in Figure   11, we found that whatever the chosen initial state in the phase space, all trajectories converge to the system's trajectory with a full HS initial state (black curve) which was already found (in Figure 10) that it matches that of the local equilibrium path. It is worth emphasizing the remarkable feature of trajectories of Figure 11, which occur in two steps: (i) first, the system reaches the "parabola", with a short range parameter ( ) which sizably varies compared to which adapts to the constraints, (ii) and then it reaches the nearest attractor by joining the parabolic path. At this stage, the correlation parameter ( ) and the HS faction, , simultaneously vary. In addition, below = 0.6, a remarkable linear relation is established between and , which means that the evolution of the local probability density follows adiabatically the long-range parameter .

Incomplete Spin Transition
In the extreme case where = = 1.2 ( Figure 12) a very high plateau of the HS fraction, is obtained on cooling, from the HS state. Figure 12a shows that in the plateau region, the HS fraction reaches the value, = 0.7, while the corresponding average bond length (Figure 12b) remains unchanged, < > ≃ 1.2 nm, overall the temperature's interval of the study. It is interesting to notice that the analysis of the distribution of unit cell configurations (H , H L , H L , H L , L ), whose histogram is presented in the inset of Figure 12b, indicates that at low temperature, 80% of the unit-cells contain 3HS lattices and 1 LS lattice, and almost 20% of the remaining unit cells contain 2HS and 2LS atoms. According to these macroscopic information, the HS fraction can be evaluated to 0.8 × ¾ + 0.2 × which gives exactly, = 0.7, in excellent agreement with data of Figure 12a. It is worth mentioning that starting the simulation from the LS state leads (see inset of Figure 12a) to observe a "stable" LS state over a wide temperature range and subsequent sharp LS to HS transition on heating at ~ 180 K, accompanied with a formation of single domains (not shown here), propagating from the corner towards the center, without any self-organization of the spin states.
The analysis of the associated snapshots (self-organization) has shown here a different type of self-organization, made of alternate HS strings and antiferro-like strings in both directions. In a perfectly ordered case, this structure leads to = 3/4, which is quite close to the value = 0.7, found in the simulations, indicating the existence of 5% of additional LS species which also aggregates in the form of 1D clusters as, we can see in panel E of snapshots of Figure 12c. Analysis of the Snapshots: Evidence of Self-Organization and Symmetry Breaking The corresponding snapshots giving the spatial distribution of the spin states along this process are summarized in Figure 12c. They clearly evidence the existence of a selforganization of the spin states on cooling, when approaching the incomplete state at lowtemperature of HS fraction, = 3/4. In the simulations, we found that for ≥ 1.14, the snapshots show the presence of strings in the and directions containing alternate ferro and anti-ferro phases mixed with small LS domains, forming a macroscopic labyrinth structure. This spatial organization of the spin states is accompanied by slight deformations of the lattice which vanishes when the equilibrium distance parameter, , becomes close to . When we look into the details of the macroscopic distribution, we remark that these alternate F and AF strings tend to form special configurations of the unit cells.

Residual HS Fraction
To understand the origin of the residual HS fraction appearing at low-temperature, we perform MC simulations at 0 K, using the equilibrium HS − LS distance parameter, , as a control parameter. We also perform analytical investigations on the energy barriers between the HS and LS in order to explain the behavior of the obtained MC data. The simulations as realized as follows: we initially prepare the lattice in the HS state from both spin and bond lengths distances ( = 1 ∀ = 1, and = ). Then MC simulations on the spin and lattice degrees of freedom are realized by varying values, in the range 1 to 1.2 nm, with steps Δ ≈ 0.01. For each value, we realize 3000 Monte Carlo steps to reach the equilibrium followed with 2000 MCS for the statistics, from which we calculate the average HS fraction and nn lattice bond length, minimizing the total energy of the lattice. The results are summarized in Figure 13a,b.  decreases following a similar trend as that of and reaches the value, < > = 1, confirming that system switches from HS to LS as a function of . The HS fraction (resp., average nn bond lengths) stays equal to 0 (resp., 1.0) in the interval 1.095 < < 1.134 and starts to increase again in a singular way to reach the value, = 0.7 and < > ≃ 1.2 for = 1.2 as we found in Figure 12a,b. The insets of Figure 13a,b show that in the region 1 < < 1.2, both HS fraction and lattice parameter follow a power-law function, ( − * ) as a function of , as clearly confirmed by the linear behavior of ln( ) and ln(< > ).

Case ≤ : Appearance of an Intermediate Plateau on Heating Only
Now, we investigate the situation where the equilibrium distance decreases towards .  Let us first consider the case where = 1.0, for which Figure 14a1 shows that on heating from LS state, the HS fraction stays equal to zero (LS state) until, ~60 . Beyond this value, the HS fraction increases slowly indicating the occurrence of a plateau region. This plateau persists until ∼ 180 K. It is worth mentioning that in all this region, the average bond length, < > (see Figure 14b1) stays almost constant, due to the choice, = . Indeed, in all this thermal interval, the switching of LS sites mainly produces HS-LS pairs and only few HS-HS pairs, as well indicated in the inset of Figure 14a1 showing the thermal dependence of the occupancy of pair probabilities. Here, the HS-HS pair probability remains negligible until 100 K, while that of HS-LS pairs start to increase beyond 60 K, in phase with the HS fraction. This explains the delay in the thermal behavior of < > compared to that of , which is also displayed in the inset of Figure 14b1, reporting the system "trajectory" in its phase space ( , < >), where we can identify the existence of two regimes in the LS to HS transition. Figure 14b1 indicates that beyond = 180 K, the average bond length suddenly increases to reach 1.2 nm and the HS fraction similarly shows a sharp increase. This behavior is fully consistent with the thermal-dependence of HS-HS pairs. On cooling from the HS state, the system stays trapped in the HS state due to the large elastic energy barrier resulting from the creation of a LS site inside the HS phase. Indeed, the total energy barrier is given by Δ = 4( + 2 )( − ) = 4( + 2 ) which is maximum. As clearly evidenced in the snapshots, in all plateau region, the spin state transformation is accompanied with a negligible lattice distortion.
In the corresponding snapshots, we see clearly that, in the plateau region, the HS transformation starts from the center under the form of multi-droplets or small HS clusters, which rapidly transform to labyrinth networks (snapshots C and D) which subsist over all the plateau region, before to disappear when the system reaches the HS state. Incidentally, it is interesting to notice that the transition from the plateau to the final HS state starts from the corner as evidenced in snapshot D.
When we increase the value of , the plateau disappears progressively due to the long-range nature of the local field interactions. For the case, = 1.09, the nucleation starts from the corners (snapshot B) and the transformation is accompanied with a strong lattice deformation (snapshot C).
On the other hand, on cooling from HS state, the lattice remains blocked in the HS state ( (0 )~1) for = 1.0, except for = 1.09 where an emerging HS to LS transformation starts from the corners (see snapshots H and I) leading to ( (0 )~0.9) but the temperature is already too low, and the transformation cannot continue due to the MC kinetics. The analytical calculations of the energy barrier when the system goes from HS to LS at = 50 K for = 1.0 gives Δ → ( = 50 )~5.5 10 which is higher than the electronic energy gain −2Δ = −900 K corresponding to spin flip from HS to LS at low-temperature.

Conclusions
In conclusion, we have investigated the physical conditions allowing the emergence of unsymmetrical thermal hysteresis in spin-crossover materials, using a simple elastic model including the interactions between the molecules' spin states and the lattice. The model developed here is based on the general genuine electro-elastic model designed for 2D square lattices for which nn and nnn and equilibrium distances and elastic constants parameters depend on the spin states of the connected sites. The interaction between the sites goes through the local volume changes accompanying the spin state changes between LS and HS. This volume change naturally produces competing short-(antiferro) and long-range ferroelastic interactions, which may lead in particular situations to the emergence of self-organization of the spin state, driven by the structure. In most of the previous studies, the equilibrium distance (denoted, ) between a HS and a LS nearestneighbors is simply taken as equal to the average distance ( = ) between the equilibrium distances of HS-HS ( ) and LS-LS ( ) configurations. Here, we demonstrate that by tuning the value of between these two quantities leads to unexpected results.
Indeed, while for, = = , usual symmetric thermal hysteresis with widths depending on the misfit Δ = − are recovered. In contrast, when becomes smaller or bigger than, , asymmetrical thermal hysteresis, with a plateau only in one branch, emerge from the simulations. It is interesting to notice that in most of the models developed for SCO materials, when a plateau appears on heating, another plateau also exists on cooling, while here we could discriminate between both processes thanks to an adapted elastic frustration that plays only on cooling or on heating. Indeed, tuning the lattice bond length for the HS-LS configuration, generates a selective antiferro-elastic interaction which acts on the LS to HS transition when < or on the HS to LS transition when, > . This effect was possible, because this model considers that the instantaneous distance between two SCO sites ( , ) does not depend only of the sum ( + ) of the connected fictitious spins but also on their product through the relation, , = + + + , where the factor, = , changes its sign from both sides of .
The analysis of the spin states configurations of the lattice along the cooling and heating of the same thermal hysteresis branches showed different type of nucleation and growth modes related to the different of shapes of the two branches. Thus in the case of < , the thermal hysteresis branch on cooling, which does not display any plateau, shows a single domain nucleation process starting from the corners. In contrast, the nucleation of the HS phase on heating shows a very different behavior: the HS domains appears everywhere in the lattice as large domains surrounded by ramified and interpenetrated LS strings, forming closed loops. A similar behavior is also observed in the cooling branch for, > . The asymmetrical character of the thermal hysteresis is then much deeper than the simple fact of different shapes of the heating and cooling branches of the curves of the HS fraction but also concerns the difference of microscopic organization of the spin states and their interactions along these two branches belonging to the same thermal hysteresis. Indeed, this organization in the form of ramified stripes has been largely discussed in the case of = 1.1336, which showed a wide plateau on cooling. The analysis of the lattice into alternate 1D horizontal and vertical A and B sublattices (stripes) with their own order parameters demonstrated the occurrence of a symmetry breaking in this region. Indeed, while sublattices A and B are equivalent in HS and LS phases, their order parameters clearly show a bifurcation in the plateau region, with the presence of partially-ordered phase. In addition, the calculations of the probability density for a site to belong to a chain has shown that in the plateau region, minority species tend to be in a chain. Finally, the study of the time and temperature-dependences of the relevant parameter = + and its behavior in the flow diagram highlights the strong attractor nature of the striped configuration as well as the "local equilibrium" character of the lowtemperature HS to LS relaxation processes. Overall, these facts are interesting from the fundamental point of view since the energy landscape of the system depends on the transformation pathway, which leads to bistable materials having different physical properties according to the branch of the hysteresis, which is visited.