Elastic Modeling of Two-Step Transitions in Sterically Frustrated 1D Binuclear Spin-Crossover Chains

Among the large family of spin-crossover materials, binuclear systems play an important role due to their specific molecular configurations, allowing the presence of multi-step transitions and elastic frustration. Although this issue benefited from a significant number of spin-based theories, there is almost no elastic description of the spin transition phenomenon in binuclear systems. To overcome this deficiency, in this work we develop the first elastic modeling of thermal properties of binuclear spin-crossover solids. At this end, we investigated a finite spin-crossover open chain constituted of elastically coupled binuclear (A = B) blocks, ⋯A=B−A=B−A=B⋯, in which the considered equivalent A and B sites may occupy two configurations, namely low-spin (LS) and high-spin (HS) states. The sites of the binuclear unit interact via an intramolecular spring and couple to the neighboring binuclear units via other springs. The model also includes the change of length inside and between the binuclear units subsequent to the spin state changes. When injecting an elastic frustration inside the binuclear unit in the LS state, competing interactions between the intra- and the inter-binuclear couplings emerge. The latter shows that according to the intra- and inter-binuclear elastic constants and the strength of the frustration, multi-step transitions are derived, for which a specific self-organization of type (HS = HS)-(LS-LS)-(HS = HS)⋯ is revealed and discussed. Finally, we have also studied the relaxation of the metastable photoinduced HS states at low temperature, in which two relaxation regimes with transient self-organized states were identified when monitoring the elastic frustration rate or the ratio of intra- and intermolecular elastic interactions. These behaviors are reminiscent of the thermal dependence of the order parameters of the system. The present model opens several possibilities of extensions of elastic frustrations acting in polynuclear spin-crossover systems, which may lead to other types of spin-state self-organizations and relaxation dynamics.


Introduction
Spin-crossover materials (SCO) emerged recently as one of the new paradigms of switchable molecular solids with serious potential applications in reversible high density memories fast switching molecular electronics as pressure and thermal sensors, in addition to their possible applications in spintronic [1][2][3][4][5][6][7][8][9][10][11]. Among the large families of spincrossover materials, iron (II)-based SCO solids [1,3,[12][13][14][15][16][17] constitute very interesting objects, due to their bistable vibronic nature which allows them to switch between a nonmagnetic low-spin (LS) state (LS, 0 23], electric current, static electric field [24][25][26], high magnetic or pulsed field, etc. [27]. From the microscopic point of view, and at the molecular scale, the SCO phenomenon can be considered as originating from a vibronic interaction which couples the electronic and the vibrational structures of the molecule. This results in the electronic switching of the iron center between the LS and the HS states, and is accompanied with a change of the iron-ligands distances. Earlier crystallographic studies [5,28] demonstrated that the volume of the coordination sphere of the molecule expands by almost 30% during the switching from the HS to LS, while the unit cell volume increases by only 3-5%. Thus, locally, each conversion of a SCO molecule causes an avalanche of acoustic waves that travel in the crystal and reflect at the surface, which in turn are felt by the other molecules of the lattice. For the thermally-induced SCO transitions, several molecules may switch at the same time in the lattice, each one producing acoustic waves reflecting at the surface and interfering, leading to the generation of complex compressive and tensile stresses, depending on the direction, inside the lattice. This is the physical elastic origin of the longrange interactions between the SCO units, which is clearly directional and may lead to competing interactions. From the experimental point of view, a large variety of thermallyinduced spin transition behaviors is observed, among which are the case of (i) continuous spin transitions [29], (ii) discontinuous first-order transitions [30], (iii) incomplete transitions [31] with residual HS fractions and (iv) multi-step spin transitions [32][33][34][35]. It is also well known that the spin transition is sensitive to the nature of the ligands [36] as well as to the solvent molecules [37].
Despite these different experimental thermal behaviors, these materials share the same physics. Indeed, it was demonstrated experimentally that changing the ligand or the anion in the material transforms some transitions from first-order to double or three-step transitions [38]. In addition, under applied external pressure, some SCO transitions change from first-order transitions to continuous or incomplete transitions [39]. In particular, the case of two-step transitions attracted a lot of attention, which is clearly the result of the microscopic organization of the HS and LS molecules in the plateau region, and is not due to the existence of some intermediate spin state at the molecular level between the HS and the LS state in Fe(II)-based SCO solids. Interestingly, spin transition materials made of interacting polynuclear SCO, among which the most famous are the dinuclear SCO systems [40,41] usually lead to two-step SCO transitions where the transitions can be continuous, such as for [Fe(2-pic)3]-Cl2.EtOH [42,43], or both are of first-order [44] or even one is gradual and the other is discontinuous [45]. In all these cases, the two transitions are separated by a plateau region, in which the thermal transition strongly "slows down". The thermal width of this plateau usually depends on the elastic rigidity of the chemical bond linking the two iron centers. The presence of the plateau at the transition in the case of dinuclear systems can be explained by steric effects between the two metal centers that hinder the concomitant transition from LS to HS, which involves the increase in the bondlengths of both transition metals. As a result, the transition of one metal center strongly affects the ligand field of the other, which then delays its transition to a higher temperature. It is also worth noticing that several experimental examples found that the two-step transition can be accompanied by a symmetry breaking [46][47][48][49][50][51][52], which results in the emergence of self-organized spin states in the plateau region. In addition, recent reports have shown the existence of multi-steps (three, four and sometimes even more) along the spin transition [49,50,53]. Among these was the famous "Devil's Staircase" [50], in which the number of steps was also enhanced by the presence of inequivalent iron centers. Interestingly, similar behaviors have been also observed in one-dimensional SCO solids [54][55][56], which means that the elastic frustration from which these phenomena originate presents a common trend.
From the theoretical point of view, spin-crossover materials have been modeled by various type of models. Some of them were first based on Bragg-Williams descriptions of regular solutions leading to the so-called domain model of Sorai and Seki [57], while others were based on the Ising-type description initiated by Wajnflasz and Pick [58,59]. Later, it was demonstrated that these two models are equivalent [60].
On the other hand, the modeling of multi-step spin transitions started first with Isinglike models with competing ferro-and antiferro-like interactions [59,61] solved in meanfield approximation or using Monte Carlo method. Very recently, considerable progress has been made towards an elastic description of multi-step transitions phenomena in SCO systems based on the concept of elastic frustration resulting from the presence of antagonist steric effects, which was introduced in ref. [47]. In this paper, the authors have studied the case of a 2D lattice with elastic frustration along the diagonal direction, which shrinks along the LS to HS transition while the unit cell concomitantly expands along the lateral direction.
Later on, in ref. [62], we investigated the case of 1D frustrated elastic SCO mononuclear chains, in which the frustration was introduced through conflicting interactions between the nn and nn equilibrium distances. We have shown in this study that according to the strength of the frustration, several thermodynamic behaviors can be reproduced, among which multi-step [30,[49][50][51][52][63][64][65][66][67][68][69] and incomplete spin transitions are retrieved.
In the present work, we extend the elastic modeling of SCO systems to that of binuclear SCO solids. For this, we first investigate the thermodynamic properties of an isolated open 1D chain made of elastically coupled SCO binuclear units (A = B), where both A and B can switch from the LS to HS with temperature. The sites A and B are coupled via an intramolecular spring, and coupled to the other sites belonging to neighboring binuclear units through other springs. To broaden the horizons of the model, an elastic frustration is injected in the intramolecular part, which results in competing interactions between the intra-and the inter-binuclear couplings.
This electro-elastic Hamiltonian, which includes both spins and positions variables, is solved using Monte Carlo (MC) simulations preformed on both degrees of freedom. Analytical calculations based on a homogeneous elastic description of the chain are also developed to validate the MC results. The results show that the misfit of intra-and interbinuclear elastic interactions is already sufficient to create an elastic frustration between the binuclear units, leading to stabilization of intermediate states when moving from the LS to the HS state. These intermediate states display the emergence of spatial self-organization of the spin states, in which the binuclear units are coupled antiferro while the organization inside the binuclear is of ferro type. Similar results are obtained by monitoring the strength of the elastic frustration.
The manuscript is organized as follows. Section 2 introduces the model Hamiltonian and the details of the simulations; in Section 3, we present the results related to the equilibrium properties and derive the phase diagrams; in Section 4, we present the behavior of the isothermal relaxation of the metastable HS for various intra and intermolecular elastic constants and frustration rates; and in Section 5 we conclude and outline some possible extensions of the model.

The Model Hamiltonian
The present SCO model for binuclear lattices is inspired from our previous elastic modeling [70][71][72] of the spin-crossover material. Let us first recall the basic concepts of the elastic modeling of SCO materials, accounting for their spin state and volume degrees of freedom. In the case of 1D SCO systems, the Hamiltonian contains electronic (spin state) and lattice (elastic) contributions.
In this work, we consider a spin-crossover chain ( Figure 1) constituted of elastically coupled binuclear units, = , in which A and B sites may occupy two configurations, namely, low-spin (LS) and high-spin (HS). Using the same ideas of Ising model, each spin state of the molecule is associated with a fictitious two-state spin operator, . Thus, = +1 will correspond to the HS state, while the value = −1 will be associated with the LS state. The sites are connected by springs and are supposed to move only along the chain, i.e., in the direction. Thus, the binuclear at site contains two atoms A and B with positions and , and their associated spin states are denoted with and . Each site is connected elastically to its nearest neighbors (nn) by a spring of elastic constant 1 inside the binuclear unit, by 1 between the binuclear units and to its next nearest neighbors by another spring of constant elastic 2 , as schematized in Fig  is the intramolecular elastic constant, while 1 and 2 are the nn intermolecular and nnn elastic constants. The bond lengths distances between the sites depend on the spin states of the connected sites.
The equilibrium bond length connecting two nn sites and is noted as � , �. This quantity can be written as a function of the spin states, so as to have a bigger length of the chain in the HS state in comparison with that of the LS state. In the same way, the equilibrium distance between two nnn sites and , denoted ( , ), is also a function of the spin states, and . The Hamiltonian describing this 1D systems, taking into account the electronic and elastic degrees of freedom of the lattice, is derived from our previous 2D electro-elastic model [71], which has been investigated for several situations of crystal shapes and sizes [72][73][74][75][76][77][78][79], including core-shell [80] and elastically-frustrated lattices [47,62]. The Hamiltonian of this 1D version designed for 1D binuclear chain of Figure 1 writes, The first term of Hamiltonian (1) represents the effective ligand-field energy, which here is temperature-dependent due to the presence of the degeneracy ratio, g, between the HS and the LS states. Thus, the total contribution, writes, ∆ 0 − ln( ), where 0 , is 0 K ligand field energy and ln( ), is the entropy of the HS states, which originates from electronic and vibrational degrees of freedom. The molar entropy change at the transition, ∆ = ln , easily derived from calorimetric data [34], allows the determination of values, which are found as ranging from 100 to 2000. The second term of the Hamiltonian refers to the intramolecular elastic interactions inside the binuclear units, while the third term accounts for the elastic inter-binuclear A-B closest atoms (see Figure 1), which can be considered here as nearest neighbors (nn). The fourth and the fifth terms of Equation (1) describe the elastic interactions between next nearest neighbors (nnn) atoms of type (A-A and B-B) belonging to different binuclear sites. Here, all elastic constants, 1 , 1 and 2 are assumed, for simplicity, as independent on spins and lattice positions. Furthermore, for simplicity, the ligand fields, Δ 0 , and degeneracies, , of A and B sites are assumed to be equal.

The Equilibrium Distances
The molecules (here, the sites) interact via elastic springs, as depicted in Figure 1, and the equilibrium bond lengths between two nn and ± (resp. nnn. and ± ) sites with spins and (resp. and , and ) are noted 0 ( , ) (resp. 0 � , � and 0 � , �). We denote by 0 , 0 and 0 , the equilibrium distances between nn HS-HS, LS-LS and HS-LS sites. Then, we have 0 = 0 (+1, +1), 0 = 0 (−1, −1) and 0 = 0 (+1, −1) . Moreover, we assume the following relation, 0 = 1 2 ( 0 + 0 ), which states that the equilibrium bond length between HS-LS configuration is equal to the average of those of LS-LS and HS-HS configurations. Using the previous parameters, 0 , 0 and 0 , it is straightforward to demonstrate that general relations connecting the equilibrium bond lengths and the spin states are written as, where = ( 0 − 0 ) represents the lattice parameter misfit between the LS and HS states.

Elastic Frustration inside the Binuclear Units
The concept of elastic frustration supposes the existence of a lattice mismatch between the nn and nnn bond lengths, leading to antagonist equilibrium distances. From the experimental point of view, these behaviors may result from the existence of steric effects imposed by the interactions between ligands of neighboring molecules, due to their extended spatial occupations [34,60,81]. For example, a strong − stacking interaction may contribute to cooperative ferro-elastic interactions, while other electrostatic interactions such as "weak" hydrogen bonding may prevent the molecular volume change along the transition between LS and HS, leading to antiferro-elastic interactions.
For this purpose, we chose the HS state free from any elastic frustration, while we inject an elastic frustration between the nn and nnn for the other HS = LS and LS = LS configurations. In particular, in the present study, the elastic frustration is considered only inside the binuclear units as follows, where > 0 is the frustration strength.
Using expression (5), the HS state ( = +1, = +1 ∀ ) is indeed non-frustrated, since (+1, +1) = 0 ∀ . When > 0, the equilibrium distances between nn pairs for L=L and H=L configurations are respectively, while those of intermolecular nnn sites are maintained to 2 0 and 2 0 , respectively. Table 1 summarizes the sixth spin configurations of a pair of neighboring binuclear molecules one and two. It is clear from Table 1 that when the system is in the LS phase, the elastic frustration strength leads to antagonist distances between nn and nnn distances. Indeed, since the sites are restricted to move only along the direction of the chain, the nnn distance for example between A1 and A2 should be equal to the sum of the distances between A1-B1 and B1-A2. This is not the case when the frustration is injected as we can see for the LS case, where the total nnn equilibrium distance (LS1-LS2, see Table 1) is equal to 2 0 (no frustration), while the nn distances impose a total value equal to 2 0 + . As a result, the system will try to share the elastic energy excess with the other bonds in order to optimize the spatial distribution of elastic energy, avoiding large compressive and tensile stresses. Obviously, the resulting distribution of bond length distance will depend on the strengths of the nn and nnn elastic constants. Table 1. Nearest neighbors and next nearest neighbors spin configurations of a pair of binuclear units (1 and 2) and their corresponding equilibrium distances showing their antagonist character, except for = 0, where nn and nnn distances fit each other.

Spin Configurations nn Equilibrium Distances nnn Equilibrium Distances
We start by studying the frustration for a "homogeneous" binuclear chain, i.e., 1 = 1 . In the simulations, we considered the following values for the non-frustrated equilibrium distances, 0 = 1 nm, 0 = 1.2 nm, 0 = 1.1 nm, leading to the lattice misfit value, = 0.2 nm . The nn elastic constants have been fixed to 1 = 1 = 34,000 K · nm −2 = 34 meV · A −2 , and the nnn elastic constant is taken as 2 = 0.1 1 while Δ 0 = 450 K and ln = 5. This value leads to the entropy change, ∆ = 41 . . −1 , a value which is in fair agreement with experimental finding of calorimetry [32].

Results and Discussion
In the electro-elastic model for SCO materials, the spin transition is described using two macroscopic order parameters, which are the HS fraction, , connecting to the average magnetization < >, and the average lattice bond length, < >. Another quantity , which measures the degree of order between A (B) sites belonging to consecutive binuclear units, is also used. The intramolecular information about the binuclear units can be also studied by monitoring the temperature dependence of the quantities < > and < >, the average of intra and inter binuclear distances. These quantities can be expressed simply as, , and < > = where < > is the average spin value and � − � is the distance between neighboring sites A and B inside the same binuclear unit at position . This is while � +1 − � is the distance between neighboring B (in binuclear ) and A (inside binuclear + 1) sites belonging to consecutive binuclear units. The Hamiltonian (1), which includes spin and distortion variables, cannot be solved analytically in view of its complex structure. Therefore, we use the Monte Carlo (MC) technique to study its thermal properties with frustration, the strength of which is controlled by the parameter, .
The investigations of the physical properties of the adapted electro-elastic 1D model for binuclear solids is implemented on a chain with = 60 atoms, consisting in 30 binuclear units. Here, each SCO atom interacts with its surrounding nearest and next nearest neighbors. When starting the simulations, we first prepare the system in the HS (non-frustrated by construction) phase by fixing all the spins to = +1 and all nn lattice distances equal to 0 . The Monte Carlo simulations, based on the Metropolis algorithm, concerned both spin and lattice position variables. The stochastic algorithm is performed in the following way: first we select randomly a site, , with spin and position , and we change its spin state to ′ = − without position change. This new state is accepted or rejected by the usual Metropolis criterion. In the second step, whatever the result of the first MC process, we mechanically relax the lattice by a slight motion of all the nodes randomly selected, with a quantity = 0.001 nm, which is rather small compared to the distance between two consecutive sites (average distance ~1.1 nm). This procedure is repeated 10 times for each spin flip. When all nodes of the lattice are visited for the spin change, we define the unit Monte Carlo step, denoted "MCS". The same procedure is performed at each temperature. To determine the thermal properties of the system, we first cool down from 140 to 1 K and then warm up to 140 K, with a 1 K increment. At each temperature, we perform 10 5 MCS to reach the equilibrium state and we use 10 5 other MCS for the statistics. Using the parameter values of the ligand field energy, Δ 0 , and the degeneracy ratio, , and according to the general law that at the transition the total effective ligand field becomes equal to zero, we arrive immediately to the expression of the transition temperature (0) , through the relations: < > = 0 ( = 0.5) and < > = 0 , Figure 2 summarizes the thermal behavior of the HS fraction values for the elastic constant's values, 1 = 1 = 34,000 K · nm −2 = 34 meV · A −2 , for the frustration rate, = 0. We see clearly that the chain presents a first-order transition between the LS and HS states with a large thermal hysteresis, which is characteristic of a strong cooperative SCO system. The MC equilibrium temperature is evaluated by taking = + + − 2 ≈ 90 , where + and − are the upper and lower transition temperatures at = 0.5. The results of Figure 2 indicate the presence of a thermal hysteresis curve while we deal with a 1D system. This effect can be attributed to the long-range character of the elastic interactions, which belong to mean field universality class, as demonstrated by Miyashita et al. [81] and as confirmed in many recent works [82] related to the study of the dynamic phase transition in 1D. These investigations demonstrated that the hysteresis widths increase with the length of the chain.
The strength of the cooperative interactions is proportional to the misfit of elastic energy per site, � 1 + 2 2 � 2 . It is worth noticing here that this behavior contrasts with that of the usual Ising 1D chains with short-range interactions, the kinetic hysteresis of which are significantly less wide. One may think that the presence of elastic interactions, which are known to lead to long-range interactions, clearly lead to larger energy barriers of metastable phases, which may then enhance the existence of these larger kinetic thermal hysteresis. In experimental systems of 1D spin-crossover systems [66], the SCO chains are rarely isolated and several types of interactions, such as weak hydrogen bonding or strong stacking − , may connect them, which then lead to 3D topology of interactions with a strong 1D character, due to this structural anisotropy. In such situations, true first-order transitions or very sharp transitions are allowed.
The spatial distribution of the HS and LS species along the cooling and heating processes of the phase transition are summarized in Figure 2b,c. Here, we see that the transformation starts from both extremities of the chain without showing the presence of any special self-organization.  Figure 3a presents the thermal evolution of the HS fraction for different values of the frustration rate . We noticed that in the absence of frustration, the SCO transition is of first-order with hysteresis (black curves of Figure 3a), with a width of ∆ = 21 . In contrast, when the elastic frustration between nn and nnn sites comes into play, the HS fraction curve shifts to a lower temperature region, and the previous first-order transition transforms to a non-hysteretic two-step transition with a large plateau around = 0.5 (pink curves of Figure 3a). A further increase of the frustration strength enlarges the plateau and stabilizes the intermediate state. Similar information is also obtained from the thermal dependence of the average lattice parameter, < >, given in Figure 3b, which exhibits a sizable and continuous increase in < > in the LS state. In particular, it is remarked that the plateau does emerge in the simulations from the threshold value = 0.6, for which an average distance corresponds, < > ≃ 1.09 nm (~ 0 ). The analysis of the thermal behavior of the average intra-and inter-dinuclear distances, < > and < >, given in Figure 3c,d, leads to the same behavior as that of Figure 3b, due to the homogenous character of all equilibrium HS-HS, HS-LS and LS-LS distances considered inside the binuclear and between the binuclear species.
In Figure 3e, we plotted the corresponding thermal dependence of the order parameter , given in Equation (7), which quantifies the degree of ordering between consecutive binuclear units. In the HS and LS states, one has = 0, which means that sublattices A and B are equivalent. In contrast, in the plateau region where < + > = 0 (i.e., = 0), the curves ( ) reach a maximum indicating the existence of LS = LS-HS = HS ordering. This order is partial for = 0.2 and becomes almost complete in the whole plateau region for = 1, indicating the existence of a symmetry breaking between the sublattices which couple in an antiferro-elastic fashion. Figure 3f summarizes the spatial configurations of the spin states (red dots = HS and blue dots = LS) inside the lattice along the thermal transition of the two-step behavior obtained for = 1. Here, in the temperature interval 34 ≤ ≤ 45 , we see the emergence of self-organization of spin states inside the plateau region under the form of spin configurations LS = LS-HS = HS-LS = LS-HS = HS, noted for simplicity L2H2…, with the presence of antiphase boundaries. The latter take place because the symmetric configuration H2L2 has the same energy. It is interesting to remark that while these configurations occur in spin crossover molecules coupled ferro elastically inside the binuclear units, the binuclear units themselves are coupled in an antiferro fashion.

Analytical Expressions of the Relaxed Lattice Bond Lengths
This section is devoted to the analytical investigations of the dependence of the relaxed nn distance, , for the three ordered electronic configurations of Figure 3, namely H2H2 (HS = HS-HS = HS…), L2L2 (LS = LS-LS = LS…) and L2H2 (LS = LS-HS = HS…), corresponding respectively to a system fully HS, fully LS and antiferro-magnetical between binuclear units. When the system is fully HS, its elastic energy is always zero because the nn and nnn distances are not frustrated. In contrast, when the system adopts the LS and the ordered antiferro-like structures, the frustration appears, and their total relaxed elastic energies remain non-zero, due to the intrinsic incompatibility between their nn and nnn distances.
For this, we consider a homogenous elastic system, in which we replace the local nn (resp. nnn) distances , +1 (resp. , +2 ) by an invariant distance (resp. ) in the lattice [83]. We denote by 2 2 and 2 2 (resp. 2 2 and 2 2 ) the relaxed nn (resp. nnn) distances corresponding to L2L2 and L2H2 macroscopic states. The total non-relaxed energies of the system (including electronic and elastic contributions) for the configurations H2H2, L2L2 and L2H2 become, For a chain with homogenous lattice interspacing ( = ) and homogenous nn elastic constant 1 = 1 = 1 , the total non-relaxed energy density after substituting � , � and ( , ) by their expressions given in Equations (3)- (5), is given in the case of a long chain ( ≫ 1) by The relaxed distances are obtained by minimizing these expressions with respect to , which gives, is the ratio between the nnn and nn elastic constants.
In Equations (11b) and (11c), we see that the relaxed nn distances ( 2 2 ) in the LS and in the ordered antiferro-like L2H2 state ( 2 2 ), are bigger than their respective nonfrustrated values, 0 and 0 , and are simply linear increasing functions of the frustration parameter > 0 and the lattice parameter misfit , while they are decreasing functions of the elastic constant ratio . Figure 4 displays the -dependence of the analytical values of ( ), given by Equation (11b) together with the results of MC simulations of Figure 3, at 0 K. An excellent agreement is obtained between these two sets of data, confirming the relevance of the present analytical approach.

Effect of the Elastic Frustration on the Transition Temperatures and Phase Diagram
In this section, we aim to find the analytical dependence of the transition temperature on the frustration parameter . Two situations corresponding to the case of first-order transition and double step transitions, reported in Figure 3, are analyzed. In the latter case, we will express the dependence of the transition temperatures between the LS and the intermediate H2L2 state on the one hand, as well as that of the second transition temperature between the H2L2 state and the full HS state on the other hand.
We notice that in the case of the first-order transition with a thermal hysteresis, the analytical method will not give the limiting temperatures ( − , + ) of the hysteresis, which depend on the kinetic aspects of the process (lifetime of metastable states), but will rather lead to find the equilibrium temperature , while the MC simulations do exactly the contrary. In this last case, the equilibrium temperature is derived from MC data as ~( − + + )/2.
To determine analytically the transition temperatures, , as a function of the strength of frustration , we study the case of the hysteretic transition of Figure 3a. This cooperative transformation between full LS and HS states takes place between two ordered states, namely LS and HS, in which their corresponding configurational entropies are null. Therefore, the condition of equal free energies at the transition temperature transforms to that of equal energies, and we then have � = , < >= +1, = � = � = , < >= −1, = �. Using the total energies calculted in Equation (7), we arrive after some simplifications to the following expression of the transition temperature, In Equation (12), the first term corresponds to the transition temperature of the usual Ising-like model and the second corresponds to the contribution of the elastic frustration which appears as an additional effective ligand field energy, − 2 ln 2 (1+4 ) 2 2 , to ∆ 0 . It is also observed that the global transition temperature of the system shows a quadratic decrease with the frustration parameter.
The results are shown in the region (i) of the phase diagram of Figure 5, where we see that the analytical curve ( ) is in good agreement with the data of MC simulations.

Figure 5.
Phase diagram ( ) summarizing the evolution of the equilibrium transition temperatures with the frustration parameter ξ. Squares represent the MC data and the continuous lines correspond to the analytical calculations. Region (i) with 0.2 < < 0.5 relates to the single-step firstorder spin transition between HS and LS states, while region (ii) with 0.5 < < 1 delimits the complete two-step transitions behavior made of two consecutive gradual transitions. Black (resp. blue) and blue (resp. violet) squares (resp. lines) are MC (resp. analytical) data of − ( ) and + ( ).

Condition of Appearance of the Two-Step Spin Transitions
According to the results of Figure 3a, for large values of , the thermal hysteresis transforms from first-order to a two-step transition with the emergence of a plateau around ~0.5. In this region, the spin states show an antiferromagnetic-like self-organization (L2H2L2…) between the binuclear units.
As before, we calculate the analytical expressions of the transition temperatures, − ( ) and + ( ), corresponding to the transition between the full LS state and the intermediate L2H2 state, and between the L2H2 state and the full HS state. At these temperatures, the total energies of the elastically relaxed LL, L2H2 and HH states must obey the following conditions, Using the expressions of energies found in Equations (10a)-(10c), we arrive immediately after some simplifications to, The two-step transitions exist when the two following conditions are fulfilled: + ( ) > − ( ) and − ( ) > 0. In Figure 5, which represents the system's phase diagram of the system, we report both analytical ( ) curves and MC results. It is found that there is very good agreement between the analytical predictions and the MC data. In particular, two important regions of values ( < 0.5 and > 0.5) corresponding to firststep and two-step transitions are well reproduced as well as the "critical" value, ≅ 0.5, beyond which a "bifurcation-like" behavior of the transition temperature is evidenced.

Effect of the Intramolecular Elastic Constant
We now turn to the analysis of the effect of intramolecular elastic interactions inside the binuclear units on the global behavior of the chain, and for that we consider the last frustrated SCO chain with the parameter of frustration = 1. We fix the values of the nn and nnn intermolecular elastic constant to 1 = 34,000 K · nm −2 and 2 = 0.1 1 and we vary the intramolecular elastic constant, 1 . We denote by = 1 1 the ratio between the intra and inter nn elastic constants. The values of the parameters taken in the simulations are the same as the previous ones.

a/Thermal dependence of the HS fraction and lattice bond lengths:
In Figure 6a-e we show the behavior of ( ), < > ( ), < > ( ), < > ( ) and curves for different values of the ratio between intra and intermolecular con- In the first figure, we see that increasing (ie. 1 ) leads to an increase in the width of the plateau for the two-step transition, which reaches a saturation for the high values of . It is also observed (see Figure 6) that at a low temperature, the nn intramolecular distances increase and saturate rapidly to the value of the HS equilibrium distance while the intermolecular nn distances decrease when increases. This is a result of the counter-reaction of the increase in the intramolecular distance by maintaining the nnn distance free from frustration. Figure 6f presents the spatial configurations of the spin states (red dots = HS and blue dots = LS) inside the lattice of the chain along the thermal transition for = 1 and = 40, giving rise to the two-step behavior. We see again the presence of the self-organized spin state structure L2H2 in the plateau region, i.e., for 26 ≤ ≤ 45 . In addition, antiphase boundaries are also evidenced because the symmetric configuration (H2L2) has the same energy. We denote by and 2 2 the relaxed nn distances corresponding to L2L2 and L2H2 macroscopic states. We denote by the nn intramolecular distances and by the nn intermolecular distances. The total non-relaxed energy of the system is still given by Equations (9a)-(9c) for all these configurations.
The analytical dependence of 〈 〉, 〈 〉 and 〈 〉 = 〈 〉+〈 〉 2 for the three studied configurations are given in Figure 7 together with the MC results, showing a very good agreement between these two sets of data.

b/Transition temperatures and phase diagram:
The thermal dependence of the HS fraction of Figure 6a shows that for large values of , the width of the plateau increases. In the following, we determine analytically the expressions of the transition temperatures, − ( ) and + ( ) , associated with the L2L2↔L2H2 and the H2H2↔L2H2 transitions, respectively. At these temperatures, the total energies of the elastically relaxed L2L2, L2H2 and H2H2 states must respect the following conditions, Using the expressions of energies found in Equations (15a)-(15c), we arrive immediately after some simplifications to, Here again, the comparison between the analytical results and MC data, summarized in Figure 8, shows a very good agreement, confirming the relevance of the homogeneous elastic approach to reproduce the MC results.

Relaxation of the HS Fraction at Low Temperature
This last section is devoted to the analysis of the dynamical properties of the present binuclear chain through the analysis of the time evolution of the HS fraction and the average nn distances along the relaxation at low temperature (10 K) of the metastable HS state. For this, the chain is initially prepared in the HS state by setting all spin state values equal to +1 and all nn distances equal to 0 . In this state, by construction the chain does not experience any elastic frustration, which means that its elastic energy is equal to zero. Since at 10 K, the LS state is the stable one (see Figure 3a), the HS chain will relax towards the LS state. Here, we are interested in the effect of the frustration parameter, , and elastic interactions ratio, , on the shape of the relaxation curves.

Elastic Frustration Effects on the Low Temperature Relaxation of the Metastable HS State
We start by analyzing the effect of the elastic frustration on a homogeneous binuclear chain in which all involved (intra-and intermolecular) elastic constants are equal. Here, the values of the nn elastic constants have been fixed to 1 = 1 = 51,000 K · nm −2 , and those of nnn are taken as 2 = 0.1 1 . The obtained results are summarized in Figure 9a,b, where we observe that for the values of frustration parameter, 0 < < 0.6, the lifetime of the HS metastable state decreases and the relaxation process accelerates with increasing values. This behavior is ascribed to the decrease in the effective elastic energy barrier facing the flip of an HS site to LS. Indeed, for = 0, the elastic energy cost for the HS to LS spin flip is easily found as Δ 0 = � This energy barrier has a non-monotonous be-havior with , since it goes through a minimum for = 1 + 4 . It has a decreasing function in the interval 0< < 1 + 4 and its derivative becomes positive for > 1 + 4 . In the situation in Figure 9, the ratio, = 2 1 = 1 10 , and the barrier, Δ , is then expected to completely vanish for = 0.7. Δ ( ), is then clearly lower than Δ 0 in the region < 0.7 and has a decreasing behavior as increases. As a result, the relaxation of the metastable HS state is expected to be faster for higher values in this region. This analytical prediction is in excellent agreement with the MC data in Figure 9a In this second slow regime, taking place after 1000 MCS, the increase in results in an increase in the lifetime of the metastable state. This behavior is in very good agreement with the previous analysis based on the non-monotonous behavior of the elastic energy barrier with the frustration parameter, .

Effect of Intramolecular Elastic Constant on the Relaxation Curves
Here we study the effect of intramolecular elastic interactions inside the binuclear units on the relaxation process of the previous chain. For this, we monitor the ratio, = 1 1 , between the intra and intermolecular elastic constants in the range 1-1.7. The values of the other parameters are the same as those in Figure 9. Figure 10a,b displays the time dependence of the HS fraction, , as well as that of the average nn bond length distance, < >, along the relaxation process from HS to LS at T = 10K for different values of . Both figures show that an increasing leads to an increase in the lifetime of the HS metastable, which extend the relaxation process with an f) emergence of a slight plateau for = 1.7. The corresponding snapshots for the case = 1.7 are given in Figure 11.

Conclusions
1D binuclear SCO chains are widely explored in the literature due to a rich variety of behaviors, among them being first-order, gradual, two-step as well as incomplete or partial spin transitions. However, the elastic description of 1D binuclear SCO chains has never been investigated in the theoretical literature of SCO materials. To address this short-coming, we have extended the previous electroelastic model and adapted it to mimic the thermal behavior of a 1D elastic chain made of elastically-coupled binuclear spin-crossover blocks, ••• = − = − = •••. To account for the steric effects inside and between the binuclear units, we included an elastic frustration leading to antagonist equilibrium distances between nearest and next nearest neighbors. The sites of the binuclear unit interact via a harmonic intramolecular spring and the binuclear unit itself couples to the neighboring binuclear units via other springs. We have solved this model by combining MC simulations performed on the spin states and another MC procedure for the lattice positions. To understand the relevant parameters acting on the outputs of the model, we have also adapted an analytical method, developed in the context of recent studies of frustrated 1D chains [62]. By doing so, we could predict the dependence of the HS fraction on several model parameters, such as the elastic frustration strength and the intramolecular elastic constant.
In the present work, we have first demonstrated that the injection of an elastic frustration between nn and nnn equilibrium distances, i.e., inside the binuclear units, as a parameter leading to expand the intramolecular distance in the LS state, leads to stabilize intermediate spin states along the thermal transition between the LS and the HS phases. We found that according to the strength of this elastic frustration, well organized antiferro-elastic structures (LS = LS-HS = HS-LS = LS…) are stabilized. This result is contrasting with the usual antiferro HS-LS-HS-LS… configurations which appeared in the previous frustrated 1D mononuclear model [62]. However, if one considers the binuclear unit as a single component, then it becomes clear that the binuclear unit organizes itself in the plateau under the usual HS-LS-HS-LS antiferro configurations. We find that the spin state organization in the plateau is made of self-organized structures. In the second step of this work, we identified that the intramolecular elastic interaction linking the atoms inside the binuclear units plays a leading role in the thermal properties of the studied spin transition between the LS and the HS phases. Thus, an increase in this elastic constant with respect to the intermolecular elastic constant results in two-step transition behaviors with an increase in the width of the plateau which reaches saturation. On the other hand, the stabilization of the previous LS = LS-HS = HS-LS = LS… configuration expresses that in this situation the two sites inside the binuclear are so strongly correlated that each binuclear unit reacts as a single molecular object.
We have also studied the time dependence of the HS fraction and lattice parameter of the present frustrated chain, thus simulating the relaxation of a metastable photoinduced HS state. The obtained results clearly demonstrate that the shape of the relaxation curves transform from continuous to two-step behaviors, which is reminiscent of the thermal dependence of the order parameters of the system at equilibrium thermodynamics.
Lastly, the present work opens several opportunities of extensions related first with the examination of other situations of frustrations. Indeed, here we considered the HS state free from frustration, but one may consider the LS state free from frustration and the HS and HS-LS states as frustrated. This will surely lead to other types of thermal behavior of the binuclear chain. Moreover, this model deserves to be extended to two-dimensions to well analyze the dimensionality effects on the modes of self-organization of the spin states.