# Ising-Like Model Hasn't Yet Said Its Last Word: Exact and Mean-Field Investigations of 2, 3 and 4-Body Interactions in 1D Ising Chain of Binuclear Spin-Crossover Solids

Kamel Boukheddaden, Nour El Islam Belmouri, Nicolas di Scala

► To cite this version:

Kamel Boukheddaden, Nour El Islam Belmouri, Nicolas di Scala. Ising-Like Model Hasn't Yet Said Its Last Word: Exact and Mean-Field Investigations of 2, 3 and 4-Body Interactions in 1D Ising Chain of Binuclear Spin-Crossover Solids. ChemPhysChem, 2024, 25 (20), 10.1002/cphc.202400238. [https://arxiv.org/abs/2002.06317](https://arxiv.org/abs/2002.06317)

HAL Id: hal-04850631

https://hal.science/hal-04850631v1

Submitted on 19 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Ising-Like Model Hasn't Yet Said Its Last Word: Exact and Mean-Field Investigations of 2, 3 and 4-Body Interactions in 1D Ising Chain of Binuclear Spin-Crossover Solids

Kamel Boukheddaden,*[a] Nour El Islam Belmouri,[a] and Nicolas di Scala[a]

We investigate the static properties of a new class of 1D Ising-like Hamiltonian for binuclear spin-crossover materials accounting for two-, three-, and four-body short-range interactions between binuclear units of spins ($s_1^A, s_1^B$) and ($s_2^A, s_2^B$). The following 2-, 3-, and 4-body $J_1(s_1^A + s_1^B)(s_2^A + s_2^B)$, $K_1 s_1^A s_1^B (s_2^A + s_2^B)$, and $K_2(s_1^A s_1^B)(s_2^A s_2^B)$ terms are considered, in addition to intra-binuclear interactions, such as effective ligand-field energy and exchange-like coupling. An exact treatment is carried out within the frame of the transfer matrix method, leading to a 4x4 matrix from which, we obtained the thermal evolution of the thermodynamic quantities. Several situations of model parameter values were tested, among which that of competing intra-

and inter-molecular interactions, leading to the occurrence of (i) one-step spin transition, (ii) two-, three-, and four-step transitions, obtained with a reasonable number of parameters. To reproduce first-order phase transitions, we accounted for inter-chains interactions, treated in the mean-field approach. Hysteretic multi-step transitions, recalling experimental observations, are then achieved. Overall, the present model not only suggests new landscapes of interaction configurations between SCO molecules but also opens new avenues to tackle the complex behaviors often observed in the properties of SCO materials.

## 1. Introduction

Spin-crossover (SCO) materials are extensively studied due to their potential applications as sensor devices, molecular memories, displays etc.[1-6] as well as for their fascinating fundamental aspects related to the coupling between their electronic and elastic structures. Such spin transition materials have been commonly investigated in iron(II) complexes with $3d^6$ configuration in octahedral symmetry. From the magnetic properties point view, they exhibit two kinds of electronic configuration, high-spin state (HS) at high-temperature and low-spin state (LS) at low-temperature with their total spin, $S=2$ (paramagnetic state) and $S=0$ (diamagnetic state), respectively. From an optical and elastic point of view, the switching between HS and LS states, which is achieved by light,[7,11] temperature,[7,8,10,12,13] magnetic[14] and electrical field,[15] pressure[16-18] etc., is accompanied by a change of color and a significant variation of molecular volume related to the variations of the bond lengths between Fe and N atoms. So, the Fe-N bond lengths increase between ~2.00 Å in the LS state and ~2.2 Å in the HS state,[7,8,10] leading to an abrupt expansion by around 30% of the molecular volume of the coordination sphere.[19] At the same moment, the expansion of the unit cell is ~3-5%,[20] very small compared to that of the coordination sphere. Indeed, this important difference of volume result from the crystal packing whose a large part of the molecular volume expansion is

absorbed inside the lattice in the form of reorientations of the ligands and others local degrees of freedom which do not affect the unit cell volume, but induce huge distortions of the crystal. Indeed, in cooperative spin-crossover systems, the nucleation and the growth of the spin states during the spin transition arise with local volume expansions which take place at several regions in the lattice (short-range nature of the interactions) and their propagation (coalescence) and interference in the whole lattice (due to the long-range character of the interactions) cause a global volume expansion accompanied by a large deformation of the crystal lattice, which produces inhomogeneous mechanical stresses inside the system. The interplay between the spin-crossover and the lattice properties of the spin transition compounds has been studied experimentally as well as theoretically in the way to understand the main mechanisms that lead to the spin-crossover phenomenon. From the experimental point of view, the investigations on the molecular complexes relate to the experiments of x-ray diffraction, optical microscopy, calorimetry, Mössbauer and uv-Vis spectroscopies, diffuse reflectivity, photoluminescence[19,21-31] etc. and show many features of these materials such as gradual transition which follows the Boltzmann statistics, first-order transition with hysteresis, incomplete spin transition with residual HS fraction at low-temperature, two- and multi-step transition characterized by intermediate plateaus in which the spin states are self-organized.[12,32-39] Theoretically, the mechanism behind all of these behaviors are investigated using macroscopic or microscopic descriptions such as the regular solutions model[40] based on a thermodynamical approach where the interaction parameters related to the weak intermolecular interactions in iron(III) compounds are introduced in a phenomenological way, continuous medium model[41] which

<footnote>[a] K. Boukheddaden, N. E. I. Belmouri, N. di Scala
Université Paris-Saclay, CNRS-Université de Versailles Saint-Quentin-en-Yvelines Groupe d'Études de la Matière Condensée, UMR 8635,
45 Avenue des Etats Unis, 78035 Versailles, France
E-mail: kamel.boukheddaden@uvsq.fr</footnote>
<figure><img src="image_5.png" /></figure>

does not explicitly describe the particular micro-organization of
the HS and the LS domains, and also Ising-like models[42,44-48]
which has a microscopic origin. Among the Ising-like models,
we distinguish the vibrational models such as the spin-phonon
model[49] accounting for the coupling between spin and phonon
in the lattice, which discard the deformation of the lattice, and
the models with anharmonic potentials such as Lennards-Jones
pair potential[50] which contains a repulsive short-range and
attractive long-range contributions or Morse potential,[51] me-
chanoelastic model,[52] harmonic electroelastic model[53-55] ac-
counting for the lattice deformation. All these models are based
on the interactions (short- and long-range) between the SCO
neighboring atoms whose amplitude of forces (weak or strong)
provides information about the degree of cooperativity of these
systems. Indeed, in the electroelastic model[53,56] the spins of the
SCO sites as well as their positions are combined to study the
spatiotemporal and thermodynamic features of spin-crossover
solids. This model reproduces several SCO behaviors and offers
a lot of possible extensions. For example, in recent works,[57-59]
we used it to describe the photoexcitation of the long-lived
metastable HS state and also the light-induced thermal
transition effects as well as multi-step transitions with symmetry
breaking in the plateau region for different modeled
situations.[60-63]

Binuclear spin-crossover complexes represent an important
step in the study of polynuclear spin transition systems. They
constitute good examples to investigate the intimate relation
between the spin transition and the intra- and inter-molecular
interactions. Since 1992, where the first deep studies emerged
in the literature[64-76] a large number of binuclear SCO systems
have been synthesized. Most of them present a double-step
spin transition with temperature, with the presence of a plateau
during the transition. The nature of this plateau has been
deeply studied and in most of the cases[77] it was found that the
two atoms within the molecule occupy the configuration [HS-
LS]. It was concluded, that steric anti-ferro elastic interactions
are responsible for this behavior. However, since the last ten
years, several new binuclear complexes have shown different
types of thermal-dependence of the HS fraction, like a one-step
sharp first-order transition where the binuclear unit switches
from the LS-LS to the HS-HS configuration without any
intermediate state and others where one or sometimes several
intermediate states are got through the transition from the LS
to the HS state.

From the theoretical point of view, Ising-like models,[43,78-86]
due to their high adaptability and despite their phenomeno-
logical character, remain a source of an important activity in the
description of SCO materials, helping in the understanding of
several of their physical properties. In particular, 1D Ising-like
models benefited in the past in several investigations concern-
ing the static properties of spin-crossover chains,[87-89] by
combining nearest-neighbor short-range interactions with infin-
ite long-range interactions.[90,91] One of the advantages of these
models is the possibility of accessing their exact thermodynam-
ic properties using the analytical transfer matrix method.
Extensions including next-nearest neighbor interactions have
also been investigated using the same method and the main

results of these approaches[87,89,92] can be summarized as follows:
(i) thermal hysteresis may occur due to the competition of
short- and long-range interactions; (ii) the shape of the
hysteresis loop changes from “S” to “Z” on increasing the ratio
of the short-range interaction; two kind of microscopic organ-
izations have been found for the plateau region, HS-LS-HS-LS-
HS-LS or HS-HS-LS-LS-HS-HS-LS-LS depending on the sign of
the interaction parameters involved in the model. These studies
concerned only mononuclear SCO materials and the extension
to binuclear SCO chains is still an open problem.

In the present study, we deal with the case of 1D SCO
chains made of binuclear SCO units. Then, each molecule
contains two atoms, denoted A and B, which may be high-spin
(HS) or low-spin (LS). This work is inspired by experimental
results which showed that in some two-step transitions, in the
intermediate plateau, the spin states configuration is of the
form (HS-HS)-(LS-LS)-(HS-HS)-(LS-LS), while in most of the cases,
the antiferro configuration (HS-LS) which can be ordered like
(HS-LS)-(HS-LS)-(HS-LS)[93] or disordered such as (HS-LS)-(LS-HS)-
(LS-HS)-(HS-LS) is found. In the case of the presence of an
ordered intermediate state, one speaks about the presence of
symmetry breaking, a concept that will be explicitly discussed
in this work. On the other hand, the stabilization of config-
urations like (HS-HS)-(LS-LS)-(HS-HS)-(LS-LS) indicates that the
binuclear units react as a single object due to the presence of a
strong intra-molecular correlation. To reproduce such behavior,
it is then necessary to go beyond the classical Ising-like model
and to include interactions between strongly correlated
binuclear units, that is the origin of the presence of four- and
three- body interactions in the present model. As we will see,
according to the sign of these interactions, several interesting
behaviors can be obtained, allowing to reproduce a large panel
of experimental spin transitions, going from incomplete, to first-
order and multi-step transitions.

Historically, the two-level microscopic models developed for
cooperative spin-crossover solids are based on the Ising-like
Hamiltonian, following the pioneering work of Wajnsflasz and
Pick.[42,43] Such two-state model can be viewed as a simple Ising
model under an applied field (a temperature-dependant field to
account for the different degeneracies of the levels).[94] Its
dynamical extension[79] led exactly to the phenomenological
formula given by Hauser, Spiering and Gütlich[95] for the
sigmoidal relaxation. Here we extend this model to study
exactly 1D binuclear SCO systems.

The manuscript is organized as follows: Section 2 introduces
the model's Hamiltonian and the method of resolution;
Section 3 summarizes the obtained results and the discussion of
the thermal behavior of the system's order parameters accord-
ing to the sign of the interactions (ferro-, antiferro-type).
Section 4 is devoted to summary.
## 2. 1D Hamiltonian of Binuclear SCO System and Transfer Matrix Method

We present an exact examination of thermodynamic properties of one-dimensional binuclear Ising-like model,[42,43,47] in which nearest and next-nearest neighbor interactions are included. Let's denote by A and B the two atoms of a dimer, by $J_{AB}$ the intra-dimer interaction, and by $K_{AA}$, $K_{BB}$, $K_{AB}$ the inter-dimer interaction between the respective AA, BB, and AB atoms, located on different molecules. The ligand field energy is denoted by $\Delta$ and the degeneracy of the high spin (HS), g, plays a role of an additional entropic term $-k_B T \ln g$, equivalent to a temperature-dependent field term, $h = \Delta - k_B T \ln g$. It is necessary to recall that all the interaction parameters of the model are elastic in nature.[96] A schematic view of the chain is given in Figure 1. The Hamiltonian of a binuclear chain of SCO system writes

$$H = \sum_{i}^{N} H_{i}^{\text{intra}} + \sum_{i,j}^{N} H_{ij}^{\text{inter}}. \qquad (1)$$

Here $H^{\text{intra}}$ and $H^{\text{inter}}$, respectively express the intra-binuclear and the inter-binuclear interactions.

The intramolecular interactions contain the ligand-field energies of both sites belonging to the binuclear i and the intramolecular interaction between the two sites A and B inside the binuclear. At site i, $H_{i}^{\text{intra}}$ writes:

$$H_{i}^{\text{intra}} = -J_{AB} s_{i}^{A} s_{i}^{B} + h(s_{i}^{A} + s_{i}^{B}), \qquad (2)$$

with $h = \Delta - k_B T \ln g$.

The interactions between the binuclear units involve several types of contributions, between A–A, B–B, and A–B sites belonging to nearest-neighbors (nn) binuclear sites. In addition, we introduce in this work a new type of interaction which considers that the binuclear unit may also respond as a single and compact object, which then leads to three body interactions of type $(s_i^A + s_i^B)s_j^A s_j^B$ as well as to four body interactions of type, $s_i^A s_j^B \cdot s_j^A s_i^B$.

A schematic view of the various types of couplings between the binuclear molecules are sketched in Figure 1 which summarizes the intra- and inter-binuclear interactions considered in the present study.

In the general case, the Hamiltonian including the interactions between the binuclear species can be written as follows

$$H_{ij}^{\text{inter}} = -(K_{AA}s_{i}^{A}s_{j}^{A} + K_{BB}s_{i}^{B}s_{j}^{B}) - K_{AB}(s_{i}^{A}s_{j}^{B} + s_{i}^{B}s_{j}^{A}) \\ -K_{1}((s_{i}^{A} + s_{i}^{B})s_{j}^{A}s_{j}^{B} + s_{i}^{A}s_{j}^{B}(s_{j}^{A} + s_{j}^{B})) - K_{2}s_{i}^{A}s_{j}^{B}s_{j}^{A}s_{j}^{B}. \qquad (3)$$

In addition, we assume periodic boundary conditions, $s_{N+1}^A = s_1^A$ and $s_{N+1}^B = s_1^B$, where N is the number of A (B) spins in the chain and so the chain contains $2N$ atoms.

With the notation $\beta = 1/k_B T$, where $k_B$ is the Boltzmann constant, the partition function is given by

$$Z = \text{Tr}[e^{-\beta H}] = \sum_{s_1^A, s_1^B} \sum_{s_2^A, s_2^B} \dots \sum_{s_N^A, s_N^B} T(s_1^A, s_1^B; s_2^A, s_2^B) \dots T(s_N^A, s_N^B; s_1^A, s_1^B). \qquad (4)$$

The kernel $T$ is defined by

$$\begin{aligned} T(s_1^A, s_1^B; s_2^A, s_2^B) \\ = \exp \beta \left[ \frac{J_{AB}}{2} (s_1^A s_1^B + s_2^A s_2^B) \right] \exp \beta \left[ K_{AA} s_1^A s_2^A + K_{BB} s_1^B s_2^B \right] \\ \times \exp \left[ \beta K_{AB} (s_1^A s_2^B + s_1^B s_2^A) \right] \exp \left[ -\frac{\beta h}{2} (s_1^A + s_1^B + s_2^A + s_2^B) \right] \\ \times \exp \left[ \beta K_1 (s_1^A s_1^B (s_2^A + s_2^B) + s_2^A s_2^B (s_1^A + s_1^B)) \right] \\ \times \exp \left[ \beta K_2 s_1^A s_1^B s_2^A s_2^B \right]. \end{aligned} \qquad (5)$$

The kernel $T$ is symmetric by construction, i.e. $T(s_1^A, s_1^B; s_2^A, s_2^B) = T(s_2^A, s_2^B; s_1^A, s_1^B)$ and has a dimension $4 \times 4$. Following standard studies, like those of ANNNI (axial next-nearest neighbour Ising) model,[97,98] where similar four-point kernel $T(s_1, s_2, s_3, s_4)$ are obtained, the partition function can be written as $Z = \text{Tr}[T]^N$, where the matrix $T$ contains 16 matrix elements whose expressions are given in the Appendix.

To find the eigenfunctions $\psi_i$ of the transfer matrix $T$ we only solve one integral equation. The matrix $T$ writes

$$\sum_{s_1^A} \sum_{s_1^B} T(s_1^A, s_1^B; s_2^A, s_2^B) \psi_i(s_2^A, s_2^B) = \lambda_i \psi_i(s_2^A, s_2^B), \qquad (6)$$

Here, $\lambda_i$ corresponds to the eigenvalue associated with the eigenfunction $\psi_i(s_2^A, s_2^B)$. It is known from Perron-Frobenius

<figure><img src="image_6.png" /><figcaption>Figure 1. Schematic views of the chain and interaction configuration between the binuclear SCO molecules whose thermal properties are exactly studied using the transfer matrix method. Each metallic center A (resp. B) of the dimer molecule is associated with a fictitious spin $s^A$ (resp. $s^B$). Intra-binuclear two-bodies interactions coexist with inter-molecular two-, three- and four-bodies couplings.</figcaption></figure>
theorem[99] that the biggest eigenvalue is real and non-degenerate. The transfer matrix of Eq. (6) can be expressed as an expansion in terms of eigenvalues and eigenfunctions by

$$T(s_1^A, s_1^B; s_2^A, s_2^B) = \sum_{i=1}^{4} \lambda_i \psi_i(s_1^A, s_1^B) \psi_i(s_2^A, s_2^B). \quad (7)$$

Let us denote by $\lambda_0$ the largest eigenvalue. Using the Frobenius theorem, we know that $\lambda_0$ is real and unique. We have also verified numerically that all the other eigenvalues are also non-degenerate. The partition function $Z_N$ is the trace of $T^N$, so that in the thermodynamic limit it is simply given by

$$Z_N = \lambda_0^N. \quad (8)$$

The free energy is

$$F = -k_B T \ln Z_N = -N k_B T \ln \lambda_0. \quad (9)$$

Hence, by simple numerical derivation, we can obtain the heat-capacity $C_v$ through the relation, $C_v = -T \frac{\partial^2 F}{\partial T^2}$.

The fraction of molecules in the HS state, $n_{HS}$, is related to the average magnetization $\langle s \rangle$ per site through $n_{HS} = (1 + \langle s \rangle)/2$. In the calculations, we will mainly plot $\langle s \rangle$ which is equal to -1 in the LS state and +1 in the HS state. The average fictitious magnetization $\langle s \rangle$ is obtained from the knowledge of the eigenfunction of the fundamental state as:

$$\langle s \rangle = \sum_{s_1^A} \sum_{s_1^B} \psi_0(s_1^A, s_1^B) \left( \frac{s_1^A + s_1^B}{2} \right) \psi_0(s_1^A, s_1^B). \quad (10)$$

Despite the rather simple procedure for obtaining the partition function, the average magnetization in the previous calculations, the two-spin correlation functions require the knowledge of both eigenvalues and eigenfunctions of the whole energy spectrum. For the present model, we are interested in the pair correlation function, $\langle s_1^A s_1^B \rangle$ which gives information about the intra-binuclear organization of the spin states. Following the general method introduced in Ref. [99], we obtain in the thermodynamical limit, i.e. $N \to \infty$

$$\lim_{N \to \infty} \langle s_1^A s_1^B \rangle_N = \sum_{i=1}^{4} \left(\frac{\lambda_i}{\lambda_0}\right) \sum_{s_1} \sum_{s_2} \left( \begin{aligned} & \psi_0(s_1^A, s_1^B) s_1^A \psi_1(s_1^A, s_1^B) \\ & (\psi_0(s_1^A, s_1^B) s_1^B \psi_1(s_1^A, s_1^B)) \end{aligned} \right) \quad (11)$$

It is worth mentioning that the higher order spin-spin correlations, like $\langle s_1^A s_2^B s_1^A \rangle$ or even $\langle s_1^A s_2^B s_2^A s_2^B \rangle$ can be calculated using the same procedure described above or by differentiating the total free energy with respect to the interaction term associated with the quantity of interest.

## 3. Transfer Matrix Results

We have performed the numerical calculations of the eigenvalues and the eigenfunctions at each temperature using the routine F02EBF of the NAG library,[100] where all the eigenvalues and eigenvectors have been obtained with a precision of less than $10^{-6}$. We calculate the thermal behavior of the HS fraction and its corresponding heat capacity for different situations using the exact formulation of the transfer matrix method. In the present contribution, we aim to stress on several aspects related to the effect of the nearest (nn) and next-nearest (nnn) neighbors' interactions on the emerging thermal behavior of the HS fraction and that of the average state of the binuclear units $\langle s^A s^B \rangle$, which will be developed below.

### 3.1. On the Role of Each Model Parameter

The total Hamiltonian $H_{intra} + H_{inter}$ given in Eqs. (2) and (3) contains several contributions whose energetic effects merit to be dissected in a detailed way. Obviously, in the intra-molecular part (Eq. (2)), the ligand-field energy, $\Delta$ stabilizes the LS state while the intra-binuclear interaction $J_{AB} < 0$ favors the configurations $++$ and $--$ inside the binuclear, where $+$ and $-$ stand for the HS and LS configurations, respectively. Already at this stage, we can see the existence of a competition between $\Delta$ and $J_{AB}$. To consider the general case including the interbinuclear interactions, $K_{AAr}$, $K_{BBr}$, $K_{ABr}$, $K_1$ and $K_2$, one may examine the energy spectrum of two interacting binuclear units at 0 Kelvin for several situations of interaction parameters. However, since there are 16 energy levels, it is a heavy task to handle all situations by this way. In contrast, it is possible to discuss qualitatively the effect of all these parameters on the electronic configurations. Thus, for example, positive $K_{AAr}$, $K_{BBr}$ parameters favors the appearance of ferro-type interactions between A-A and B-B sites of neighboring binuclear units. Within these signs, the following configurations of neighboring binuclear units, $(++)-(+)$, $--)(-)$, $(+-)(-+)$, $--)$ have the same energy. However, with $K_{AB} > 0$, one stabilizes the configurations $(++)-(+)$, $--)(-)$, $(+-)(-+)$, $--)$, $(-+)(-+)$ which have the same energy and destabilizes the configuration $(+-)(-+)$ $(-+)(-+)$ as well as $(\pm \pm)(-\mp)$ which are favoured when $K_{AB} < 0$. The coupling $K_1 > 0$ (resp. $K_1 < 0$) favors the configuration $(++)-(+)$ (resp. $--)(-+)$) even at low-temperature which then clearly indicates that strong positive $K_1$ values may lead to recovering the HS state at low-temperature for $K_1 > 0$, thus producing re-entrant phases. The final coupling parameter, $K_2 > 0$ (resp. $K_2 < 0$) stabilizes all configurations leading to a positive (resp. negative) product $s_1^A s_1^B s_2^A s_2^B$, like $(++)-(+)$, $--)(-)$, $(+-)(-+)$ etc. So, without even reviewing all possible configurations, we see clearly that the present model opens a large possibility of thermal behaviors for the HS fraction and the spin-spin correlations. In the following, we examine some of these thermodynamic behaviors. We note in passing that in calculations that follow, we will consider, as is customary in statistical physics, the
<figure><img src="image_7.png" /></figure>

Boltzmann's constant as equal to 1, so all energies will be given in Kelvin.

## 3.2. Case of Short-Range Interactions

The first part of the study is devoted to the case of short-range interactions only, which is exactly the case of the total Hamiltonian (Eq. (1)), combining the intramolecular contribution (Eq. (2)) and the intermolecular one (Eq. (3)). Several configurations of interactions, including ferro- and antiferro types in the intra- and inter-binuclear units will be considered in this section. Although the absence of long-range interactions prevents the emergence of hysteretic transitions in this case, the competition between ferro- and antiferro short-range interactions is at the origin of a rich variety of thermal behaviors and spin organization, as will be shown in this section.

### 3.2.1. Case of Non-Interacting Binuclear Species

The case of non-interacting binuclear species is quite obvious to investigate, however, it has the merit to help to understand the system behavior when all intra- and inter-binuclear interactions are suppressed. The results are summarized in Figure 2 which shows a simple gradual and continuous transition of the average fictitious spin state, $\langle s \rangle$ and the free energy $F$, which appears as a decreasing function of temperature, as expected by basic thermodynamics. On one hand the transition temperature, obtained for $\langle s \rangle = 0$, has the expression: $T_{eq}^0 = \frac{\Delta}{k_B \ln g} = 120$ K. On the other hand, the average correlation function inside the dimers, $\langle s^A s^B \rangle$, shows a peak around the transition temperature, with a value $\langle s^A s^B \rangle = 0$ indicating a disordered state in the configurations of the dimers at this transition.

<figure><img src="image_8.png" /><figcaption>**Figure 2.** Thermal evolution of the average fictitious magnetization, $\langle s \rangle$ (black curve) and intramolecular spin-spin correlation pair $\langle s^A s^B \rangle$ (red curve) and corresponding free energy per site (blue curve) for the case of non-interacting binuclear units. The parameter values are:</figcaption></figure>

$J_{AB} = K_{AA} = K_{BB} = K_{AB} = K_1 = K_2 = 0$ K, $\Delta = 600$ K and $\ln g = 5$.

### 3.2.2. Sharp One Step Transition: Case $K_{AB} > 0$, $K_{AA} > 0$ and $K_{BB} > 0$

In the restricted case of ferro intra-binuclear interactions $J_{AB} > 0$, and when $K_{AA} = K_{BB} = K_{AB} > 0$ and $K_1 = K_2 = 0$, our problem can be solved analytically.[99] The expected effect of these ferro-like short-range interactions is to steep the slope of the HS fraction $n_{HS}$ (related to the average fictitious magnetization through $n_{HS} = \frac{1+\langle s \rangle}{2}$) curve at the transition temperature, i.e. for $n_{HS} = \frac{1}{2}$ ($\langle s \rangle = 0$).

This indicates that the chain approaches a "first-order" phase transition-like behavior, even if we know that there is no true first-order transition with (nn) and (nnn) interactions in 1-D systems. However, including later a weak long-range ferro- or antiferro-like interaction will lead to producing the genuine first-order transitions. In the present case, the ferro-type interaction, considered inside the dimer means that the molecular "cage" is not rigid enough to prevent the breathing of each dimer and to cause steric effects between them. As a result, the sharp transition, observed in Figure 3 indicates a direct transition from LS-LS ($S^A = -1$, $S^B = -1$) dimer state to a HS-HS ($S^A = +1$, $S^B = +1$) state, as proved by the thermal behavior of the average $\langle s^A s^B \rangle$ order parameter which stays almost equal to +1 for all temperatures, except around the transition temperature where a sharp peak (with a low amplitude) appears.

### 3.2.3. Case of Antiferro-Like Intra-Binuclear Interaction: $J_{AB} < 0$

Here, we investigate the case of non-interacting binuclear entities with intra-molecular coupling, $J_{AB} < 0$. Due to this antiferro-like interaction, which favors the appearance of HS-LS species, along the thermal transition, a two-step transition is obtained, as shown in Figure 4a. Although not shown here the width of the plateau depends on the strength of the intra-antiferro interaction, which originates from steric effects between the two constituents of the dimer. To confirm the self-organization of the chain in the plateau region, we calculate the thermal dependence of the average two-spins correlations function, $\langle s^A s^B \rangle$ (red curve in Figure 4a) which shows peaks in this region around the value $\langle s^A s^B \rangle \approx -0.5$, which means that 75% of the dimers occupy the state HS-LS. The corresponding heat capacity (blue curve in Figure 4a) displays a double Schottky anomaly whose peaks correspond to the splitting of the previous transition temperature of Figure 2 to two transition temperatures, estimated as $T_{eq}^{\pm} \approx \frac{\Delta}{k_B \ln g} \pm J/2$ ($\sim 70$ and $\sim 170$ K). To predict this behavior, we depict in Figure 4b, the field $h$ ($h = \Delta - k_B T \ln g$) dependence of the energy of the sixteen configurations of two neighboring SCO binuclear units. Here, strongly negative (resp. positive) $h$-values correspond to high- (resp. low) temperatures. According to this energy diagram, we see that starting from highly positive values of $h$ and following the minimum energy curves leads to ($--$; $--$) as a fundamental state. When decreasing $h$, a crossover occurs for $h^\pm = J_{AB} = -100$ K with the state ($-+$; $-+$) which is four times degenerate $[(\pm \mp; -+)$ and $(\pm \mp; + -)]$. Later, this state ceases being the
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 3.** Thermal evolution of the average fictitious magnetization, 〈 *s* 〉 showing a sharp spin transition between the LS and HS states in the case of ferro-type interactions inside and between the binuclear units. Left inset is the intramolecular spin-spin correlation pair 〈 *s*<sup>A</sup>*s*<sup>B</sup> 〉 which remains almost equal to +1 and the right inset is the corresponding heat capacity, given in *k*<sub>B</sub> units. The parameter values are: *J*<sub>AB</sub> = *K*<sub>AA</sub> = *K*<sub>BB</sub> = *K*<sub>AB</sub> = 100 K and *K*<sub>1</sub> = *K*<sub>2</sub> = 0 K, Δ = 600 K and ln *g* = 5.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 4.** a) Thermal evolution of the average fictitious magnetization, 〈 *s* 〉 (black curve) and intramolecular spin-spin correlation pair 〈 *s*<sup>A</sup>*s*<sup>B</sup> 〉 (red curve) and corresponding heat capacity per site (blue curve), given in *k*<sub>B</sub> units, showing two peaks for the case of non-interacting binuclear units. b) Evolution of the energy spectrum as a function of the field *h* = Δ−*k*<sub>B</sub>T ln *g* indicating the evolution of the configuration energies of two (nn) binuclear units. The parameter values are: *K*<sub>AA</sub> = *K*<sub>BB</sub> = *K*<sub>AB</sub> = *K*<sub>1</sub> = *K*<sub>2</sub> = 0 K − *J*<sub>AB</sub> = −100 K.</figcaption></figure>

fundamental one below $h^+ = -J_{AB} = 100$ K for the benefit of the non-degenerate state (++; + +). From the thermodynamic point of view, these two values, $h^+$, of the effective field correspond the temperatures $T^\pm = T_{eq}^0 \pm \frac{2J_{AB}}{k_B \ln g'}$ which are valued to 80 and 160 K, in excellent agreement with the maximums of the heat capacity of Figure 4, obtained at $T^\pm = 75$ and $T^\pm = 160$ K.

It is interesting to mention that the plateau observed in Figure 4 does not consist in an ordered state of the HS and LS spins species, and so no symmetry breaking occurs during this double transition. Indeed, keeping nonzero only the intra-molecular coupling parameters, $J_{AB}$ and $\Delta$, leads to Hamiltonian (Eq. (2)), which is invariant by the permutation of the spins $s^A$ and $s^B$, which makes them equivalent. On the other hand, one always has in this case, $\langle s^A \rangle = \langle s^B \rangle$, where $\langle s^A \rangle$ and $\langle s^B \rangle$ are the

average thermodynamic of $s^A$ and $s^B$. Consequently, the plateau region in Figure 4 is constituted of a disordered state mixing, (+,−) and (−,+), binuclear states. Later on, the case of symmetry breaking in the plateau will be discussed when we introduce the long-range interactions between the binuclear units.

### 3.2.4. Case of Antiferro Inter-Binuclear Short-Range Interaction: $K_2 < 0$

The special case of antiferro interactions ($K_2 < 0$) between binuclear entities is considered here, where only the interaction parameter $K_2$ is nonzero ($K_{AA} = K_{BB} = K_{AB} = K_1 = 0$). According to the previous discussions, and the Hamiltonian Eq. (1), this
<figure><img src="image_12.png" /></figure>

situation favors the emergence of (+;±), (−;±) configurations and their symmetric forms (±;+), (±;−) between (nn) binuclear molecules. At the same time, the effective ligand, h = Δ−k$_{B}$TIng stabilizes the (−;−) and (+, +) configurations in the low temperature (LT) and high temperature (HT) regions, respectively. First of all, it is important to notice that in the present situation where only h and K$_{2}$ are non-zero, the total Hamiltonian becomes invariant by exchanging s$_{i}^{A}$ and s$_{i}^{B}$. As a result, we have $\langle s^{A} \rangle = \langle s^{B} \rangle = \langle s \rangle$.

### 3.2.4.1. Thermal Dependence of the Configuration Probabilities

To characterize the self-organization of the spin states along the spin transition, let's derive the probability of each of the 16 configurations of two neighboring binuclear units. Let's denote by $P(s_1, s_2, ..., s_N)$, the global probability of finding the chain at temperature, T, with the spin values, s$_{1}$ for the 1st site, s$_{2}$ for the 2nd site etc. Remarking that the quantity $\frac{1}{2}(1 + s_{i}^{\prime})$ is equal to +1 for $s'_{i} = s_{i}$ and to 0 for $s'_{i} = -s_{i}$ (doing this we run over all possible values of $s'_{i} = \pm 1$), one can write the following identity for the general case of N spins:

$$P(s_1, s_2, ..., s_N) = \frac{1}{2^N} \sum_{\{s_i\} = (\pm s_i)} (1 + s_1 s'_1)(1 + s_1 s'_1) \dots (1 + s_N s'_N) P(s'_1, s'_2, ..., s'_N) \quad (12)$$

where $\{s_j\} = \{s_1, s_2, ..., s_N\}$ and $\{s'_j\} = \{s'_1, s'_2, ..., s'_N\}$. The previous expression can be written as follows:

$$P(s_1, s_2, ..., s_N) = \frac{1}{2^N} \prod_{i=1}^{N} \sum_{s_i = \pm s_i} (1 + s_1 s'_1)(1 + s_2 s'_2) \dots (1 + s_N s'_N) P(s'_1, s'_2, ..., s'_N) \quad (13)$$

Although these developments can be continued in the general version, we limit ourselves from now to the case of two binuclear units ($s'_1, s'_2, s'_3, s'_4$) for the sake of simplicity. In the present case, the probability of each configuration for four neighboring sites (N = 4) leads to the following equation:

$$P(s_1, s_2, s_3, s_4) = \frac{1}{2^4} \prod_{i=1}^{4} \sum_{s'_i = \pm s_i} (1 + s_1 s'_1)(1 + s_2 s'_2) (1 + s_3 s'_3)(1 + s_4 s'_4) P(s'_1, s'_2, s'_3, s'_4) \quad (14)$$

After simple developments, Eq. (14) can be re-written as a function of the macro-variables of the system:

$$P(s_1, s_2, s_3, s_4) = \frac{1}{2^4} \left[ 1 + \sum_{i=1}^{4} s_i \cdot m_i + \sum_{i \neq j} s_i s_j \cdot r_{ij} + \sum_{i \neq j \neq k} s_i s_j s_k \cdot R_{ijk} + \sum_{i \neq j \neq k \neq \ell} s_i s_j s_k s_\ell \cdot Q_{ijkl} \right] \quad (15)$$

Where $m_i = \langle s_i \rangle = \sum_{\{s'\}} s'_i P(s'_i, s'_j, s'_k, s'_\ell)$, $r_{ij} = \langle s_i s_j \rangle = \sum_{\{s'\}} s'_i s'_j P(s'_i, s'_j, s'_k, s'_\ell)$, $R_{ijk} = \langle s_i s_j s_k \rangle = \sum_{\{s'\}} s'_i s'_j s'_k P(s'_i, s'_j, s'_k, s'_\ell)$ and $Q_{ijkl} = \langle s_i s_j s_k s_\ell \rangle = \sum_{\{s'\}} s'_i s'_j s'_k s'_\ell P(s'_i, s'_j, s'_k, s'_\ell)$.

In the case where the total Hamiltonian is invariant through permutations of spin s$_{i}$ and s$_{j}$ ($\forall i \neq j$), one can write $\langle s_i \rangle = m$, $\langle s_i s_j \rangle = R$, $\langle s_i s_j s_k \rangle = Q$ leading to the following expression of the configuration probabilities:

$$P(s_1, s_2, s_3, s_4) = \frac{1}{2^4} [1 + (s_1 + s_2 + s_3 + s_4)m + (s_1 s_2 + s_1 s_3 + s_1 s_4 + s_2 s_3 + s_2 s_4 + s_3 s_4)r + (s_1 s_2 s_3 + s_1 s_2 s_4 + s_1 s_3 s_4 + s_2 s_3 s_4)R + s_1 s_2 s_3 s_4 Q] \quad (16)$$

The derivation of the previous probabilities (Eq. (12)) is a strong tool allowing the determination of the dominant contributions in the system while accounting for all spin-spin correlations. In the following, we will use (Eq. (16)) to derive the thermal dependence of the spin configuration probabilities of two neighboring dimers to understand the microscopic mechanism of the evolution of microscopic populations of states and their competition as a function of temperature. In the particular case of antiferro interactions between the binuclear units, the access to the occupation probabilities is of great help in identifying the microscopic self-organization of the system, which can not be easily disentangled from the macroscopic variables, due to the highly degenerate nature of the energy levels. In the end, one can easily write down the possible states with degeneracies $\Omega$ as follows:

$$\begin{align}
P(+, +, +, +) &= \frac{1}{16} (1 + 4m + 6r + 4R + Q); \quad \Omega = 1 \nonumber \\
P(+, -, +, +) &= \frac{1}{16} (1 + 2m - 2R - Q); \quad \Omega = 4 \nonumber \\
P(+, -, -, -) &= \frac{1}{16} (1 - 2m + 2R - Q); \quad \Omega = 4 \nonumber \\
P(+, -, +, -) &= \frac{1}{16} (1 - 2r + Q); \quad \Omega = 4 \nonumber \\
P(+, +, -, -) &= \frac{1}{16} (1 - 2r + Q); \quad \Omega = 2 \nonumber \\
P(-, -, -, -) &= \frac{1}{16} (1 - 4m + 6r - 4R + Q); \quad \Omega = 1 \nonumber
\end{align} \quad (17)$$

### 3.2.4.2. Results of the Numerical Calculations

Let's now go back to discussing the results of the simulation presented in Figure 5. First of all, from Figure 5a depicting the thermal dependence of the total magnetization, we identify a
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 5.** a) Thermal evolution of the average fictitious magnetization, $\langle s \rangle$ showing multi-step transitions caused by the antiferro short-range intermolecular coupling, $K_2$. b) Evolution of the energy spectrum of two coupled binuclear units as a function of the field $h = \Delta - k_B T \ln g$, indicating several crossovers corresponding to the multistep transitions. Large negative [resp. positive] $h$ values correspond to the high temperature (HT) [resp. low temperature (LT)] region. c) Corresponding thermal evolution of the intramolecular spin-spin correlation pair $\langle s_A^A s_B^B \rangle$ and corresponding d) heat capacity per site (in $k_B$ units) showing several peaks located at the transition temperatures. e) Thermal evolution of the four spin correlation $\langle s_A^A s_B^B s_A^B s_B^B \rangle$ in blue showing a minimum at 108 K with a value $\langle s_A^A s_B^B s_A^B s_B^B \rangle = -0.72$. f) Thermal dependence of four spin configuration probabilities where the black, magenta, cyan, red, blue, green curves correspond to the configurations $(+; +; +)$, $(+; -; -)$, $(+; -; +)$, $(+; +; -)$, $(-; -; +)$, $(-; -; -)$ respectively. The parameter values used for the simulation are: $\Delta = 600$ K, $\ln g = 5$, $J_{AB} = 0$ K, $K_{AA} = K_{AB} = K_{BB} = 0$ K, $K_1 = 0$ K and $K_2 = -100$ K.</figcaption></figure>

three-step spin transition which is obtained by adjusting only the parameter $K_2$. This multi-step transition reflects in the energetic behavior of the phase diagram given in Figure 5b, as well as in the intramolecular nn correlation $\langle s_A s_B \rangle$ (Figure 5c), in the heat capacity (Figure 5d) and in the other order parameters (Figure 5e) and probability occupations of the energy levels (Figure 5f). Thus, starting from the HT phase (~200 K), where the HS state is dominant, a first partial and very smooth transition occurs around 170 K corresponding to the first “peak” in Figure 5d representing the thermal dependence of the heat capacity. In this region, the order parameter, $\langle s \rangle$, it changes from 1 to 1/2 ($n_{HS} = 3/4$) as shown in Figure 5a. This behavior announces the existence of spin configurations of type $(\pm\mp, + +)$ as well as their symmetric forms $(+ +, \pm\mp)$ for consecutive binuclear units. To predict the consecutive transitions observed in Figure 5a, we calculate the energetic phase diagram of a couple of two interacting binuclear units related to the present case. The following energy levels are obtained

$$E(s_1, s_2, s_3, s_4) = h(s_1 + s_2 + s_3 + s_4) - K_2 \cdot s_1 \cdot s_2 \cdot s_3 \cdot s_4, \quad (18)$$

where $h = \Delta - k_B T \ln g$. Thus, the negative (resp. positive) value of $h$ corresponds to the high- (resp. low-) temperature phase. Due to the highly symmetric form of the expression of the energy in

Eq. (18), the system presents only five levels whose three are degenerate, as depicted in Figure 5b, which depicts the effective field-dependence of the energies of the accessible spin state configurations. Based on the crossing lines of the energy diagram of Figure 5b, the configurations $(\pm\mp, + +)$ are the most favorable to appear first on cooling as a crossing between the energy of the fourth degenerate $(\pm\mp, + +)$ configuration and the $(+ +; + +)$ configuration. The intersection occurs for $h \approx -200$ K, which corresponds to $T = 180$ K). This is further supported by the increase in the probability of the $(+ +; - +)$ configuration represented by the red curve in Figure 5f, obtained from Eq. (16). Indeed, the value of $P(+, +, +, -)$ in the curve representing the total probability of the four equivalent forms $(+ +; - +)$, $(+ +; - +)$, $(+ -; + +)$, and $(- +; + +)$ increases as the temperature decreases and the species transform from HS to LS state, while the value $P(+, +, +, +)$ decreases and the probability of all the other configurations stay almost equal to zero indicating that at the microscopic level, the appearance of the $(+ +; - +)$ self-organization and its symmetric forms are favored upon the first gradual conversion. In addition, to ensure the predicted dominant configuration is obtained in the simulation, we plot the thermal dependence of the mean value of the correlation $\langle s_i^A s_{i+1}^B s_{i+1}^A s_{i+1}^B \rangle$ in Figure 5e, where we also see a decrease of the
<figure><img src="image_15.png" /></figure>

four body correlation from $\langle s_i^A s_j^B s_{i+1}^A s_{j+1}^B \rangle \sim 1$ at 250 K, to $\langle s_i^A s_j^B s_{i+1}^A s_{j+1}^B \rangle = -0.13 = -0.13$ at 170 K, thus indicating a change in the order of the spin states towards the predicted configuration (+, +, +).

By lowering the temperature, a second transition occurs at 118 K (Figure 5a) which corresponds to the second sharp peak (maximum) in the heat capacity curve, where this time the order parameter in Figure 5a goes to 0, indicating a new phase transition with the appearance of a new ordering of the binuclear units in the form of (±, +, −) and its symmetric form (−, +, ±). Interestingly, at T=118 K, the correlation function $\langle s_i^A s_j^B \rangle$ drops to zero abruptly, as a result of the coexistence of the above-mentioned configurations and the (+, +, ±) and (±, +, +) ones. This behavior is anticipated by the phase diagram plot in Figure 5b, which shows that at h=0 corresponding to T=T<sub>eq</sub>, a crossing between the (+, +, +) and the (−, −, +) configurations takes place, indicating that at T=T<sub>eq</sub> = 117 K, these configurations would be the dominant ones.

This is further established by the thermal dependence of the configuration probabilities curves of Figure 5f, as the red curve representing the probability of the four times degenerate state (+, +, +, −) which was at a maximum value after the first transition decreases, while the blue curve representing P(−, −, −, +), the probability of the four equivalent forms (−, −, +, +), (−, −, +, +), (+, −, −, +), and (+, −, −, +) increases from zero to its maximum value at the end of the second transition, with a crossing between the two curves at T=T<sub>eq</sub>, temperature at which although being dominant, these eight configurations coexist with the two times degenerated state (+, +, −, +) and the four times degenerated state (+, −, +, +), whose contributions appear as marginal. The results yielded by the simulation are in agreement with these predictions as the value of the four body correlation stays at its minimum as shown in Figure 5e with $\langle s_i^A s_j^B s_{i+1}^A s_{j+1}^B \rangle \sim -0.5$, which can be observed only if the dominant configurations are the (+, +, +, −) and (−, −, +, +) ones.

The third and last transition occurs at 67 K as seen in the heat capacity plot of Figure 5d which peaks at low temperature, as well as from the mean value of the magnetization curve Figure 5a where the order parameter changes from −1/2 to −1, indicating the appearance of (−, −, −) short-range binuclear units. This transition is also predicted by the phase diagram of Figure 5b which depicts a crossing point between the (−, −, +) and the (−, −, −) energetic lines for h ≃ 200 K, corresponding to T ≃ 80 K, a value which is in fair agreement with simulations findings. After this transition, the system goes to the LS state with only the (−, −, −) configuration existing as confirmed by the thermal dependence of the configuration probabilities in Figure 5f.

### 3.3. Case of Interacting Chains and Thermal Hysteresis

Let's consider now the case where the SCO chains interact with each other. The purpose of this section is to get closer to experience by considering inter-chain interactions in order to produce first-order transitions. Indeed, in real 1D SCO solids, the chains are not isolated and steric or/and electrostatic inter-

actions play an important role in driving abrupt and hysteretic transitions.[101] This is the case in SCO chain compounds[102,103] in which non-coordinated water molecules and counter-anions connect the polymeric chains, leading to the existence of true first-order phase transitions, accompanied with wide thermal hysteresis.

To extend the present model, in order to account for inter-chain interactions,[55,104,105] we use the previous Hamiltonian of the isolated chain as a starting point. First of all, to simplify the presentation of the equations, we adopt from now the following new variables to describe the spin states or spin configurations of the binuclear unit

$$t_i = s_i^A s_i^B \quad \text{and} \quad S_i = s_i^A + s_i^B. \qquad (19)$$

In Eq. (19), the knowledge of $t_i$ and $S_i$ for the ith-binuclear site, uniquely defines its state. Within these notations, the previous total Hamiltonian (Eq. (1)) of an isolated chain of interacting dinuclear SCO units, can be simply written in a compact form as,

$$H^c = \sum_i (hS_i - J_{AB}t_i) - \sum_j (KS_i S_{i+1} - K_1(S_i t_{i+1} + t_i S_{i+1}) - K_2 t_i t_{i+1}). \qquad (20)$$

Within this new formulation, the interaction constant K does not depend on sites A and B, we implicitly assume that $K_{AA} = K_{AB} = K_{BB} = K$. However, we'll come back to the general case later when we give the expression for the transfer matrix.

The total Hamiltonian including inter-chain effects writes

$$H_{tot} = H^c + \sum_{i \neq j} V_{ij}^{c-c}, \qquad (21)$$

where $H^c$ is the Hamiltonian of the chain, given by Eq. (20), and $V_{ij}^{c-c}$ is the interchain interaction potential energy. Let us denote by i and j two neighbouring chains as depicted in Figure 6, and by (i, k) and (j, l) the coordinates of two SCO dinuclear units belonging respectively to chains i and j and labelled by the indices k and l, respectively.

Let us denote by $S_{i,k}$, $S_{j,l}$, $t_{i,k}$ and $t_{j,l}$, their associated fictitious "spins". The inter-chain interaction potential between two-chains i and j writes,

$$V_{ij}^{c-c} = -J' \sum_{k,l} S_{i,k} S_{j,l} - K_1 \sum_{k,l} (S_{i,k} t_{j,l} + t_{i,k} S_{j,l}) - K_2 \sum_{k,l} t_{i,k} t_{j,l}, \qquad (22)$$

where $J'$, $K_1'$ and $K_2'$ are the inter-chain interactions, which are schematized in Figure 6. Obviously, one may solve this new Hamiltonian (21) using Monte Carlo simulations (exact calculations). Here, we choose to keep the exact treatment for the intra-chain Hamiltonian and treat the inter-chain contribution in the mean-field approach, which is known to be equivalent to considering infinite homogeneous long-range interactions between the chains, where the coupling has the same strength, whatever the distance between the sites. Within this approach, the inter-chain potential (22) becomes,
<figure><img src="image_16.png" /><figcaption>Figure 6. Schematic view of two SCO chains, some interchain interaction configurations between the binuclear units to enlighten the connections between the chains. See text for more explanations; in particular the structure of the $V_{ij}^{C-C}$ interchain potential (22) helps to identify the set of considered inter-chain interactions.</figcaption></figure>

$$V_{LR}^{C-C} = -J_{LR} \sum_{i} S_{i} \langle S \rangle - K_{1LR} \sum_{i} (S_{i} \langle t \rangle + t_{i} \langle S \rangle) - K_{2LR} \sum_{i} t_{i} \langle t \rangle, \quad (23)$$

where, the previous interchain interaction constants are now noted: $J_{LR}$, $K_{1LR}$, and $K_{2LR}$. It is worth noticing that the present infinite long-range interactions are well-adapted with their elastic nature in SCO materials. These elastic interactions between mononuclear units was found to belong to the mean-field universality class, as proven in Ref. [106].

Finally, the resulting total mean-field Hamiltonian, accounting for intra- and inter-chain interactions can be written under the following form

$$H_{tot}^{MF} = (h - J_{LR}\langle S \rangle - K_{1LR}\langle t \rangle) \sum_{i} S_{i} - (J_{AB} + K_{1LR}\langle S \rangle + K_{2LR}\langle t \rangle) \sum_{i} t_{i} - K \sum_{i} S_{i}S_{i+1} - K_{1} \sum_{i} (S_{i}t_{i+1} + t_{i}S_{i+1}) - K_{2} \sum_{i} t_{i}t_{i+1}, \quad (24)$$

where the index i runs along the chain.

It is useful at this point to recall the following notations:

$$S_i = (s_i^A + s_i^B), t_i = s_i^A s_i^B, \langle S \rangle = \frac{\langle s_i^A \rangle + \langle s_i^B \rangle}{2} \text{ and } \langle t \rangle = \langle s_i^A s_i^B \rangle.$$

Thus, the total mean-field Hamiltonian (24) expresses a self-consistent problem in the basis of the 1D chain Hamiltonian (1). The thermodynamic properties of the system are then obtained self-consistently using the Transfer Matrix method. The new kernel T, given here in the general case, using the intra-chain interaction constants $K_{AA}$, $K_{AB}$ and $K_{BB}=K$, is defined by:

$$\begin{aligned} T(s_1^A, s_1^B; s_2^A, s_2^B) &= \exp \beta \left[ \frac{J_{AB} + K_{1LR} \langle S \rangle + K_{2LR} \langle t \rangle}{2} (s_1^A s_1^B + s_2^A s_2^B) \right] \\ &\quad \times \exp \beta [K_{AA} s_1^A s_2^A + K_{BB} s_1^B s_2^B + K_{AB} (s_1^A s_2^B + s_1^B s_2^A)] \\ &\quad \times \exp \beta \left[ \frac{J_{LR} \langle S \rangle + K_{1LR} \langle t \rangle - h}{2} (s_1^A + s_1^B + s_2^A + s_2^B) \right] \\ &\quad \times \exp \beta [K_1 (s_1^A s_1^B (s_2^A + s_2^B) + s_2^A s_2^B (s_1^A + s_1^B))] \\ &\quad \times \exp \beta [K_2 s_1^A s_1^B s_2^A s_2^B] \end{aligned} \quad (25)$$

where $J_{LR}$ ($J_{AB}$), $K_{1LR}$ ($K_1$), and $K_{2LR}$ ($K_2$) are respectively the long-range (resp. short-range) 2-, 3-, and 4-body interaction parameters.

The calculations of the order parameters $\langle s^A \rangle$, $\langle s^B \rangle$ and $\langle s^A s^B \rangle$ from the transfer matrix lead to three self-consistent equations which are solved numerically using Newton-Raphson method. The procedure is very simple: we first start the calculations from the HS state, by fixing all average “order parameters” to their expected values ($\langle S \rangle = 1$, $\langle t \rangle = 1$) and then we diagonalize numerically the transfer matrix, from which we derive the greatest eigenvalue, noted $\lambda_0$ and therefore the free energy per site $f = -k_B T \ln \lambda_0$.

There are two ways of calculating $\langle S \rangle$ and $\langle t \rangle$ at each iteration. One may simply calculate numerically the following derivatives $\langle S \rangle = -\frac{\partial f}{\partial \langle S \rangle}$ and $\langle t \rangle = \frac{\partial f}{\partial \langle t \rangle}$ or determine these average values by using the eigenvectors of the transfer matrix through Eqs. (10) and (11). We have checked that both procedures lead to the same results. In addition, high-order moments, like $\langle s_1^A s_1^B (s_2^A + s_2^B) \rangle = 1/2 \frac{\partial f}{\partial \langle s_1^A \rangle}$ and $\langle s_1^A s_1^B s_2^A s_2^B \rangle = \langle t_1 t_2 \rangle = \frac{\partial f}{\partial \langle s_1^A s_1^B \rangle}$, can be achieved using the first method, without having to go through the eigenvectors of the transfer matrix.

Finally, the variational free energy per site of this self-consistent system writes
<figure><img src="image_17.png" /></figure>

$$f = -k_B T \ln \lambda_0 - k_B T \ln g + \frac{1}{2} J_{LR} \langle S \rangle^2 + \frac{1}{2} K_{2LR} \langle t \rangle^2 + K_{1LR} \langle S \rangle \langle t \rangle. \quad (26)$$

According to the values and the sign of the interaction parameters $J_{LR}$, $K_{1LR}$ and $K_{2LR}$, one may obtain simple first-order transitions with hysteresis or hysteretic multi-step transitions. Several cases can be distinguished in the long-range interactions' effects, which we report in the following sections.

### 3.3.1. Effect of Two Bodies $J_{LR}$ Long-Range Interaction and Symmetry Breaking

Here we analyze the simple case where only two body long-range interactions are acting between the chains, considering a general interaction scheme, where A–A, A–B and B–B sites belonging to different chains may interact differently. The interaction potential writes in this specific case,

$$V^{C-C} = - \sum_{i \neq j} \sum_{k} \sum_{\ell} \left[ J_{LR}^{AA} s_{i,k}^{A} s_{j,\ell}^{A} + J_{LR}^{BB} s_{i,k}^{B} s_{j,\ell}^{B} + J_{LR}^{AB} (s_{i,k}^{A} s_{j,\ell}^{B} + s_{i,k}^{B} s_{j,\ell}^{A}) \right]. \quad (27)$$

In the mean-field approximation, and after absorbing the coordination number in the interaction parameters, for which we keep, for convenience, the same notation, the previous contribution becomes

$$V_{LR}^{C-C} = - \sum_{i} \left[ J_{LR}^{AA} s_{i}^{A} \langle s^{A} \rangle + J_{LR}^{BB} s_{i}^{B} \langle s^{B} \rangle + J_{LR}^{AB} (s_{i}^{A} \langle s^{B} \rangle + s_{i}^{B} \langle s^{A} \rangle) \right], \quad (28)$$

where i runs now along the chain sites. It is quite easy to demonstrate that the interaction potential (28) renormalizes the previous effective ligand field energies ($h = \Delta - k_B T \ln g$) of sites A and B which now become,

$$\left\{ \begin{aligned} h_{\text{eff}}^{A} &= \Delta - k_{B} T \ln g - J_{LR}^{AA} \langle s^{A} \rangle - J_{LR}^{AB} \langle s^{B} \rangle \\ h_{\text{eff}}^{B} &= \Delta - k_{B} T \ln g - J_{LR}^{BB} \langle s^{B} \rangle - J_{LR}^{AB} \langle s^{A} \rangle. \end{aligned} \right. \quad (29)$$

In the general case, the coupling parameters, $J_{LR}^{AA}$, $J_{LR}^{BB}$ and $J_{LR}^{AB}$ can be different, which automatically breaks the symmetry between the sublattices A and B. A similar study was done in detail a long-time ago using a two-bodies interaction within the Ising-like model by Bousseksou et al.[43,78] and demonstrated a rich variety of the phase diagram of this model with the occurrence of a symmetry breaking between the sub-lattices A and B in the plateau region.

For simplicity and to show the case of symmetry breaking in the plateau region, we considered the situation where the long-range coupling $J_{LR}^{AA} \neq J_{LR}^{BB}$ and $J_{LR}^{AB} = 0$. In this situation, exchanging $s_i^A$ by $s_i^B$ in the Hamiltonian does not keep it invariant. The obtained results are summarized in Figure 7a where we display the thermal-dependence of the two-order parameters $\frac{\langle s^A + s^B \rangle}{2}$ (black curve) and $\frac{\langle s^A - s^B \rangle}{2}$ (blue dotted curve). Contrasting with the upcoming case of Figure 8, the present two-step transition shows a symmetry breaking in the plateau region where the parameter $\frac{\langle s^A - s^B \rangle}{2} \neq 0$, which means that the A and B sublattices become inequivalent in this region. This result is confirmed by the thermal behavior of $\langle s^A \rangle$ and $\langle s^B \rangle$, shown in Figure 7b, depicting the existence of two bifurcation points on heating and cooling between the two order parameters.

a)

<figure><img src="image_18.png" /><figcaption>Figure 7. a) Thermal evolution of the average fictitious magnetization, $\langle s \rangle = \frac{\langle s^A + s^B \rangle}{2}$ (bold black curve) showing a double step transition for the long-range interaction parameter values $J_{LR}^{AA} = 300$ K and $J_{LR}^{BB} = 100$ K. A symmetry breaking behavior is revealed through the thermal-dependence of $\frac{\langle s^A + s^B \rangle}{2}$ (blue dotted curve) which shows a characteristic behavior indicating the existence of a partial order HS=LS-HS=LS-HS=LS in cooling and LS=HS-LS=HS-LS=HS in heating between the binuclear units. b) Thermal-dependence of $\langle s^A \rangle$ and $\langle s^B \rangle$ showing the different behavior of the two sites ($\langle s^A \rangle$ in red and $\langle s^B \rangle$ in blue) around the plateau region. On both heating and cooling processes, B sites switch first then followed by A sites. The other parameter values are: $J_{AB} = -100$ K, $K_{AA} = K_{BB} = K_{AB} = K_1 = K_2 = 0$ K, $J_{LR}^{AB} = 0$ K, $\Delta = 600$ K and $\ln g = 5$.</figcaption></figure>

<figure><img src="image_19.png" /></figure>
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 8.** a) Thermal evolution of the average fictitious magnetization, $\langle s \rangle$ in black showing a double step transition for $J_{LR}^{AA} = J_{LR}^{BB} = J_{LR}^{AB} = 250$ K. Both order parameter $\langle s^A \rangle$ and $\langle s^B \rangle$ remain equal at all temperatures. The intermediate state is disordered as confirmed by the behavior of $\langle s^A s^B \rangle$ in blue and $\frac{\langle s^A - s^B \rangle}{2}$ in the inset and so no symmetry breaking occurs in the plateau region. b) Temperature dependence of the free energy showing a decreasing function with temperature allowing to see the true equilibrium transition temperatures, $T_{eq1}$ and $T_{eq2}$. Inset: enlargement of the transition region showing the location of the metastable states on heating and cooling. The other parameter values are: $J_{AB} = -100$ K, $K_{AA} = K_{BB} = K_{AB} = K_1 = K_2 = 0$ K, $\Delta = 600$ K and $\ln g = 5$.</figcaption></figure>

To convince the reader about the existence of two-step transitions without symmetry breaking, we performed similar calculations by taking $J_{LR}^{AA} = J_{LR}^{BB} = J_{LR}^{AB} = 250$ K for the same order parameter values as those of Figure 7. In the present case, the Hamiltonian becomes invariant with respect to the exchange of $s_i^A$ and $s_i^B$ sites which guaranties the relation $\langle s^A \rangle = \langle s^B \rangle$, thus excluding the occurrence of a symmetry breaking, despite the existence of a plateau at the transition, caused by the antiferrotype intra-molecular interaction $J_{AB}$.

The obtained results, summarized in Figure 8a, show that the double step transition takes place on both sublattices concomitantly, with order parameters $\langle s^A \rangle = \langle s^B \rangle$ which proves that the symmetry breaking is absent in this situation. Thus, the plateau region, results in the existence of a disordered state in both A and B sublattices, leading to $\langle s^A \rangle = \langle s^B \rangle = 0$ to which corresponds a value $\langle s^A s^B \rangle = -0.5$. This is in very good agreement with the behavior of the order parameter $\frac{\langle s^A - s^B \rangle}{2}$ which is equal to zero at all temperatures, proving definitely that there is no self-organization in the plateau of Figure 8a. In addition, to check the stability of the solutions in Figure 8a, we plot in Figure 8b the corresponding temperature-dependence of the variational free energy, given by Eq. (26), where we have identified the equilibrium transition temperatures, which are the temperatures at which the free energies of the ordered (LS or HS) and disordered phases are equal. According to basic thermodynamics, the free energy should be a decreasing function with temperature, which is the case here, except for metastable and unstable states, well characterized by the presence of the cusps in the inset of Figure 8b. In addition, the identified values, $T_{eq1} \approx 60$ K and $T_{eq2} \approx 187.3$ K, of the transition temperatures, from the free energy panel of Figure 8b are in excellent agreement with those deduced from the thermal evolution of the fictitious magnetization, plotted in Figure 8a.

### 3.3.2. Three and Four Steps Transition Caused by Negative Four Bodies Inter-Dimers Interactions

When negative short-range intra- ($J_{AB} = -200$ K) and inter-dimers interactions, $K_2 < 0$ K, coexist with ferro long-range couplings $J_{LR}^{AA} = J_{LR}^{BB} = J_{LR}^{AB} = 250$ K and ferro-type long-range of four bodies inter-dimer interactions $K_{2LR} = 300$ K, three and even more steps transitions are obtained and accompanied with several hysteresis'. Figures 9a and 9b report the thermal-dependence of the average fictitious magnetization in two cases corresponding to weak and strong anti-ferro short-range 4-bodies interactions (parameter $K_2 = -100, -200$ K) in the presence of ferro long-range 4-bodies coupling ($K_{2LR} = 300$ K). The thermal-dependencies of their associated variational free energies are given respectively in Figures 9c and 9d.

In the case of strong antiferro ($K_2 < 0$) interactions, 3-step transitions (Figure 9a) are obtained with intermediate plateaus at $\langle s \rangle = 0.5$ and $-0.5$, corresponding to HS fractions values equal to 3/4 and 1/4 which has been several times obtained in Ref. [107]. The case of weaker $|K_2|$ values presented in Figure 9b is even more spectacular since it leads to several steps combining hysteretic and gradual transitions, denoting the prominent role of competing interactions. Here also, similar curves are often observed in the experiments[36] and sometimes without thermal hysteresis as in the case of the so-called "devil staircase" reported in Ref. [108]. In most cases, these behaviors relate to the presence of structural phase transitions and the occurrence of symmetry breaking in the plateau regions, which accompanies the spin transition process.

For both previously discussed situations, we have calculated the thermal dependence of the free energy function, given by Eq. (26). In the case of 3-step transitions of Figure 9a, the corresponding graph of the free energy (Figure 9c) shows three accidents at temperatures $T_{eq1} \approx 62.3$ K, $T_{eq2} \approx 119.2$ K and $T_{eq3} = 194.2$ K which fall exactly in the regions of thermal transitions
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>b)</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>c)</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>d)</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure 9.** Thermal evolution of the average fictitious magnetization, 〈 *s* 〉 showing in a) a 3-step hysteretic transitions separated between them by large plateaus and in b) a 4-step transition, for K₂ = −200 K and K₂ = −100 K respectively. c, d) plots are the respective thermal dependencies of the free energies showing the three and four transition temperatures. The inset in panel c) is the enlargement of the transition regions, where the metastable phases, presenting cusps in free energy behavior, are highlighted by the dashed ellipses. Similar metastable states are also identified in panel d). The other parameter values are:</figcaption></figure>

$$
J_{LR}^{AA} = J_{LR}^{BB} = J_{LR}^{AB} = 250 \text{ K}, K_{2LR} = 300 \text{ K}, J_{AB} = -200 \text{ K}, K_{AA} = K_{BB} = K_{AB} = K_1 = 0 \text{ K}, J_{LR}^{AB} = 0 \text{ K}, \Delta = 600 \text{ K and } \ln g = 5.
$$

of the corresponding effective magnetization. Similarly, Fig-
ure 9d which presents the free energy of the 4-step transitions
case (corresponding to Figure 9b) leads to transition temper-
ature, T<sub>eq1</sub> &asymp; 97.7 K, T<sub>eq2</sub> &asymp; 112.0 K, T<sub>eq3</sub> = 126.7 K and T<sub>eq4</sub> =
147.8 K, in excellent agreement with those of panel b).

**4. Conclusions**

In conclusion, we have examined the case of binuclear SCO
materials in the frame of an extended version of the Ising-like
Hamiltonian that includes, in addition to the usual two bodies
interactions, new binuclear-binuclear couplings, where the
binuclear units which are considered as a single object (the
molecule). Indeed, when the two iron centers belonging to the
same dimer are strongly correlated, interactions between the
binuclear units, reacting as a single object, cannot be neglected
anymore. Previous Ising-like models, designed to describe the
thermodynamic properties of coupled binuclear systems,

accounted only for two-bodies intra- and inter-binuclear
interactions, where the intra-molecular one is considered as
antiferro-like and the inter-molecular interaction as ferro. These
competing interactions, lead naturally to two-step thermal
transitions of the HS fraction, which can be both of first-order
or gradual or one gradual and the other of first-order depend-
ing on the model parameters. Thus, in the former models, the
interaction terms, between two binuclear units labelled *i* and *j*
and containing spins (*s*<sub>*i*</sub><sup>*A*</sup>, *s*<sub>*i*</sub><sup>*B*</sup>) (binuclear *i*) and (*s*<sub>*j*</sub><sup>*A*</sup>, *s*<sub>*j*</sub><sup>*B*</sup>) (binuclear *j*),
can only contain the mixing product (*s*<sub>*i*</sub><sup>*A*</sup> + *s*<sub>*i*</sub><sup>*B*</sup>) · (*s*<sub>*j*</sub><sup>*A*</sup> + *s*<sub>*j*</sub><sup>*B*</sup>), which
does not take into account for the total electronic information
contained in each binuclear. Indeed, each binuclear unit *i* is
totally defined by the knowledge of the sum (*s*<sub>*i*</sub><sup>*A*</sup> + *s*<sub>*i*</sub><sup>*B*</sup>) and the
correlation *s*<sub>*i*</sub><sup>*A*</sup> · *s*<sub>*i*</sub><sup>*B*</sup>. Thus, in the general case, interactions of type
(*s*<sub>*i*</sub><sup>*A*</sup> + *s*<sub>*i*</sub><sup>*B*</sup>) · (*s*<sub>*j*</sub><sup>*A*</sup> + *s*<sub>*j*</sub><sup>*B*</sup>), also (*s*<sub>*i*</sub><sup>*A*</sup> + *s*<sub>*i*</sub><sup>*B*</sup>)(*s*<sub>*j*</sub><sup>*A*</sup>*s*<sub>*j*</sub><sup>*B*</sup>) (and its symmetric form)
for the three-body as well as (*s*<sub>*i*</sub><sup>*A*</sup>*s*<sub>*i*</sub><sup>*B*</sup>)(*s*<sub>*j*</sub><sup>*A*</sup>*s*<sub>*j*</sub><sup>*B*</sup>) for the four-body,
should also be taken into account when treating the interaction
scheme between two binuclear units *i* and *j*. Thus the present
model can be viewed as a general extension of the Ising-like
model for poly-nuclear SCO materials since similar issues can be considered in strongly-coupled trinuclear and tetra-nuclear systems. According to this new extended version in which, in addition to the short-range interactions, we included long-range interactions representing the stress-driven elastic nature of the SCO phenomenon, new order parameters appear in the model describing the additional degrees of freedom conferred to the system. Here, the model is exactly solved in its 1D version using the well-known transfer matrix method which leads to diagonalizing an 8×8 matrix, similarly as for the 1D Ising model with nearest and next-nearest neighbors interactions. Among the various new results obtained with this extended model, one can quote (i) the three steps gradual transitions, (ii) two steps transitions with and without symmetry breaking, (iii) three steps transitions with three hysteresis' and (iv) finally four steps transitions with open hysteresis'. All these results have been obtained by analysing the phase diagrams of the evolution of the energy spectrum of two coupled binuclear species as a function of the ligand-field Δ. Moreover, one should notice that by considering two (nn) binuclear units, one has 16 energetic configuration levels. Thus by considering different interaction parameters between A-A sites and B-B, for example, one may lift additional degeneracies and access to more steps along the thermal transition between the LS and the HS states. It is also important to notice that the extension of the model to 2D including 3- and 4-bodies interactions is quite obvious and is planned as a next target. We expect in this case to realize unprecedented self-organization of the spin states in the plateau regions. Moreover, the extension of this study to include polynuclear SCO materials, among which the case of 1D trinuclear SCO chains, as the one reported by Pittala et al.,[109,110] is in progress.

APPENDIX

The matrix elements of the transfer matrix $T(s_1, s_2, s_3, s_4)$ are given in the basis $(s_A, s_B)$ connecting two spins from the first dimer $(s_1, s_2)$ and two spins of the neighbouring dimer $(s_3, s_4)$. Due to the symmetric nature of the kernel, we only give the diagonal and upper diagonal elements of the matrix; while the other elements are deduced by symmetry.

$$
\begin{align}
T(++;++) &= T(+1, +1; +1, +1) \nonumber \\
&= e^{\beta(J_{AB}+K_{AA}+K_{BB}+2K_{AB}+4K_1+K_2-2h)} \nonumber \\
T(++;-) &= T(+1, +1; +1, -1) \nonumber \\
&= e^{\beta(K_{AA}-K_{BB}-2K_1-K_2-h)} \nonumber \\
T(++;-) &= T(+1, +1; -1, +1) \nonumber \\
&= e^{\beta(-K_{AA}+K_{BB}-2K_1-K_2-h)} \nonumber \\
T(++;-) &= T(+1, +1; -1, -1) \nonumber \\
&= e^{\beta(J_{AB}-K_{AA}-K_{BB}-2K_{AB}+K_2)} \nonumber \\
T(--;--) &= T(+1, -1; +1, -1) \nonumber \\
&= e^{\beta(-J_{AB}+K_{AA}+K_{BB}-2K_{AB}+K_2)} \nonumber \\
T(--;--) &= T(+1, -1; -1, -1) \nonumber \\
&= e^{\beta(K_{AA}-K_{BB}+2K_1-K_2+h)} \nonumber \\
T(--;--) &= T(-1, -1; -1, -1) \nonumber \\
&= e^{\beta(J_{AB}+K_{AA}+K_{BB}+2K_{AB}-4K_1+K_2+2h)} \tag{30}
\end{align}
$$

The expression of the transfer matrix is then given by

$$
\begin{equation}
\begin{split}
& +1, +1 \quad +1, -1 \quad -1, +1 \quad -1, -1 \\
& +1, +1 \quad \left\{
\begin{array}{@{}l@{}}
  T(++;++) \quad T(++;+) \quad T(++;-) \quad T(++;-) \quad T(++;-) \\
  T(+;++) \quad T(+;+) \quad T(+;-) \quad T(+;-) \quad T(+;-)
\end{array}
\right. \\
& -1, +1 \quad \left\{
\begin{array}{@{}l@{}}
  T(+-;++) \quad T(+-;+) \quad T(+-;-) \quad T(+-;-) \quad T(+-;-) \\
  T(---;++) \quad T(---;+) \quad T(---;-) \quad T(---;-)
\end{array}
\right. \\
& -1, +1 \quad \left\{
\begin{array}{@{}l@{}}
  T(---;++) \quad T(---;+) \quad T(---;-) \quad T(---;-) \quad T(---;-)
\end{array}
\right.
\end{split}
\tag{31}
\end{equation}
$$

The transfer matrix $T$ has eigenvectors $|\psi_i\rangle$, corresponding to the eigenvalues $\lambda_i, i=1,2,3,4$. Each eigenvector, writes in the basis $|s_A, s_B\rangle$ as follows,

$$ |\psi_i\rangle = \alpha_{1,i}|+, +\rangle + \alpha_{2,i}|+, -\rangle + \alpha_{1,j}|-+, +\rangle + \alpha_{1,j}|-,-\rangle, \quad (32) $$

where the coefficients $\alpha_i$ are obtained by the numerical diagonalization of $T$, which then writes as

$$ T = \sum_i |\psi_i\rangle \lambda_i \langle \psi_i| \quad (33) $$

Acknowledgements

This research was funded by ANR (Agence Nationale de la Recherche Scientifique), grant number Mol-CoSM ANR-20-CE07-0028-02. This project has received financial support from the CNRS through the MITI interdisciplinary programs through its exploratory research program. The Universities of Versailles and Paris-Saclay are also acknowledged for their support.
<figure><img src="image_27.png" /></figure>

Conflict of Interests

The authors declare no conflict of interest.

Data Availability Statement

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

Keywords: Ising-like model · Spin crossover · Multi-spin interactions · Transfer matrix method · Multistep transitions · First-order phase transitions · Mean-field theory

[1] K. Boukheddaden, M. H. Ritti, G. Bouchez, M. Sy, M. M. Dirtu, M. Parlier, J. Linares, Y. Garcia, *J. Phys. Chem. C* **2018**, *122* (14), 7597–7604.

[2] J. Linares, E. Codjovi, Y. Garcia, *Sensors* **2012**, *12* (4), 4479–4492.

[3] O. Kahn, C. J. Martinez, *Science* **1998**, *279* (5347), 44–48.

[4] S. P. Vallone, A. N. Tantillo, A. M. dos Santos, J. J. Molaison, R. Kulmaczewski, A. Chapoy, P. Ahmadi, M. A. Halcrow, K. G. Sandeman, *Adv. Mater.* **2019**, *31*, 1807334.

[5] P. J. von Ranke, *Appl. Phys. Lett.* **2017**, *110*, 181909.

[6] B. Benaicha, K. Van Do, A. Yangui, N. Pittala, A. Lusson, M. Sy, G. Bouchez, H. Fourati, C. J. Gómez-García, S. Triki, K. Boukheddaden, *Chem. Sci.* **2019**, *10*, 6791–6798.

[7] P. Gütlich, A. Hauser, H. Spiering, *Angew. Chem. Int. Ed.* **1994**, *33*, 2024.

[8] K. Boukheddaden, H. Fourati, Y. Singh, G. Chastanet, *Magnetochemistry* **2019**, *5* (2), 21.

[9] G. Chastanet, M. Lorenc, R. Bertoni, C. Desplanches, C. R. Chim. **2018**, *21* (12), 1075–1094.

[10] P. Gütlich, A. Hauser, *Coord. Chem. Rev.* **1990**, *97*, 1–22.

[11] S. Decurtins, P. Gütlich, K. M. Hasselbach, A. Hauser, H. Spiering, *Inorg. Chem.* **1985**, *24*, 2174–2178.

[12] J. Krober, E. Codjovi, O. Kahn, F. Grolifcre, C. Jay, *J. Am. Chem. Soc.* **1993**, *115*, 9810–9811.

[13] S. W. Biernacki, B. Clerjaud, *Phys. Rev. B* **2005**, *72*, 024406.

[14] A. Bousseksou, N. Negre, M. Goiran, L. Salmon, J.-P. Tuchagues, M.-L. Boillot, K. Boukheddaden, F. Varret, *Eur. Phys. J. B* **2000**, *13*, 451–456.

[15] T. Mahfoud, G. Molnar, S. Bonhommeau, S. Cobo, L. Salmon, P. Demont, H. Tokoro, S.-I. Ohkoshi, K. Boukheddaden, A. Bousseksou, *J. Am. Chem. Soc.* **2009**, *131* (41), 15049–15054.

[16] D. Pinkowicz, M. Rams, M. Mišek, K. V. Kamenev, H. Tomkowiak, A. Katrusiak, B. Sieklucka, *J. Am. Chem. Soc.* **2015**, *137* (27), 8795–8802.

[17] V. Ksenofontov, G. Levchenko, H. Spiering, P. Gütlich, J. F. Létard, Y. Bouhedja, O. Kahn, *Chem. Phys. Lett.* **1998**, *294* (6), 545–553.

[18] I. Maurin, M. Itoi, J. M. Cain, D. R. Talham, T. Gacoin, K. Boukheddaden, J. P. Itié, *J. Appl. Phys.* **2021**, *129*, 23.

[19] H. Oyanagi, T. Tayagaki, K. Tanaka, *J. Lumin.* **2006**, *119*, 361–369.

[20] Y. Garcia, V. Niel, M. C. Munoz, J. A. Real, *Spin Crossover Transition Met. Comp. J* **2004**, *233*, 229–257.

[21] J. A. Wolny, V. Schünemann, Z. Németh, G. Vankó, C. R. Chim. **2018**, *21* (12), 1152–1169.

[22] K. Boukheddaden, M. Sy, M. Paez-Espejo, A. Slimani, F. Varret, *Physica B: Condensed Matter* **2016**, *486*, 187–191.

[23] M. Castro, O. Roubeau, L. Piñeiro-López, J. A. Real, J. A. Rodríguez-Velamazán, *J. Phys. Chem. C* **2015**, *119* (30), 17334–17343.

[24] M. Sorai, J. Ensling, K. M. Hasselbach, P. Gütlich, *Chem. Phys.* **1977**, *20* (2), 197–208.

[25] N. N. Adarsh, M. M. Dirtu, A. Rotaru, Y. Garcia, *Hyperfine Interact.* **2017**, *238*, 1–6.

[26] M. Sy, F. Varret, K. Boukheddaden, G. Bouchez, J. Marrot, S. Kawata, S. Kaizaki, *Angew. Chem. Int. Ed.* **2014**, *53* (29), 7539–7542.

[27] S. Pillet, J. Hubsch, C. Lecomte, *Eur. Phys. J. B* **2004**, *38* (4), 541–552.

[28] Y. Garcia, F. Robert, A. D. Naik, G. Zhou, B. Tinant, K. Robeyns, L. Piraux, *J. Am. Chem. Soc.* **2011**, *133* (40), 15850–15853.

[29] C. Lochenie, K. Schotz, F. Panzer, H. Kurz, B. Maier, F. Puchtler, B. Weber, *J. Am. Chem. Soc.* **2018**, *140* (2), 700–709.

[30] M. Estrader, J. Salinas Uber, L. A. Barrios, J. Garcia, P. Lloyd-Williams, O. Roubeau, G. Aromí, *Angew. Chem. Int. Ed.* **2017**, *56* (49), 15622–15627.

[31] F. Varret, C. Chong, A. Slimani, D. Garrot, Y. Garcia, A. D. Naik, *Spin Crossover Mater.* **2013**, *42*, 452–441 (John Wiley & Sons Ltd).

[32] T. Togo, S. A. Amolegbe, R. Yamaguchi, T. Kuroda-Sowa, M. Nakaya, K. Shimayama, S. Hayami, *Chem. Lett.* **2013**, *42* (12), 1542–1544.

[33] O. Roubeau, *Chem. Eur. J.* **2012**, *18* (48), 15230–15244.

[34] N. Moliner, C. Muñoz, S. Létard, X. Solans, N. Menéndez, A. Goujon, J. A. Real, *Inorg. Chem.* **2000**, *39* (23), 5390–5393.

[35] P. J. Van Koningsbruggen, Y. Garcia, O. Kahn, L. Fournès, H. Kooijman, A. L. Spek, P. Gütlich, *Inorg. Chem.* **2000**, *39* (9), 1891–1900.

[36] Z. Y. Li, J. W. Dai, Y. Shiota, K. Yoshizawa, S. Kanegawa, O. Sato, *Chem. Eur. J.* **2013**, *19* (39), 12948–12952.

[37] C. J. Adams, M. C. Muñoz, R. E. Waddington, J. A. Real, *Inorg. Chem.* **2011**, *50* (21), 10633–10642.

[38] C. Zheng, S. Jia, Y. Dong, J. Xu, H. Sui, F. Wang, D. Li, *Inorg. Chem.* **2019**, *58* (21), 14316–14324.

[39] Z. Y. Li, H. Ohtsu, T. Kojima, J. W. Dai, T. Yoshida, B. K. Breedlove, M. Yamashita, *Angew. Chem. Int. Ed.* **2016**, *128* (17), 5270–5275.

[40] C. P. Slichter, H. G. Drickamer, *J. Chem. Phys.* **1972**, *56* (5), 2142–2160.

[41] H. T. Spiering, N. Willenbacher, *J. Phys. Condens. Matter* **1989**, *1* (50), 10089.

[42] J. Wajnflasz, R. Pick, *J. Phys. Colloq.* **1971**, *32*, 1–91.

[43] A. Bousseksou, J. Nasser, J. Linares, K. Boukheddaden, F. Varret, *J. Phys. I* **1992**, *2*, 1381–1403.

[44] K. Boukheddaden, J. Linares, E. Codjovi, F. Varret, V. Niel, J. A. Real, *J. Appl. Phys.* **2003**, *93*, 7103–7105.

[45] S. Miyashita, H. Takano, *Prog. Theor. Phys.* **1985**, *73* (5), 1122–1140.

[46] A. Bousseksou, H. Constant-Machado, F. Varret, *J. Phys. I* **1995**, *5* (6), 747–760.

[47] J. Linares, J. Nasser, A. Bousseksou, K. Boukheddaden, F. Varret, *J. Magn. Magn. Mater.* **1995**, *1503*, 140–144.

[48] J. Pavlik, W. Nicolazzi, G. Molnár, R. Boca, *Eur. Phys. J. B* **2013**, *86*, 292.

[49] J. A. Nasser, K. Boukheddaden, J. Linares, *Eur. Phys. J. B* **2004**, *39*, 219–227.

[50] W. Nicolazzi, S. Pillet, C. Lecomte, *Phys. Rev. B* **2008**, *78*, 174401.

[51] C. Fourmental, S. Mondal, R. Banerjee, A. Bellec, Y. Garreau, A. Coati, V. Repain, *J. Phys. Chem. Lett.* **2019**, *10* (14), 4103–4109.

[52] C. Enachescu, M. Nishino, S. Miyashita, L. Stoleriu, A. Stancu, *Phys. Rev. B* **2012**, *86*, 054114.

[53] M. Ndiaye, K. Boukheddaden, *J. Phys. Soc. Jpn.* **2020**, *89* (1), 014004.

[54] R. Traiche, H. Oubouchou, K. Boukheddaden, *Symmetry* **2021**, *13*(10), 1836.

[55] R. Traiche, H. Oubouchou, K. Boukheddaden, *Crystals* **2023**, *13* (6), 937.

[56] A. Slimani, K. Boukheddaden, K. Yamashita, *Phys. Rev. B* **2015**, *92*, 014111.

[57] H. Fourati, M. Ndiaye, M. Sy, S. Triki, G. Chastanet, S. Pillet, K. Boukheddaden, *Phys. Rev. B* **2022**, *105*, 174436.

[58] A. C. Felts, J. M. Cain, M. J. Andrus, A. R. Ahir, K. A. Abboud, A. Slimani, M. W. Meisel, K. Boukheddaden, D. R. Talham, *J. Am. Chem. Soc.* **2018**, *140*, 5814–5824.

[59] A. Slimani, K. Boukheddaden, *Phys. Chem. Chem. Phys.* **2018**, *20*, 28583–2859.

[60] M. Paez-Espejo, M. Sy, K. Boukheddaden, *J. Am. Chem. Soc.* **2016**, *138* (9), 3202–3210.

[61] Y. Singh, H. Oubouchou, M. Nishino, S. Miyashita, K. Boukheddaden, *Phys. Rev. B* **2020**, *101*, 054105.

[62] M. Ndiaye, N. E. I. Belmouri, J. Linares, K. Boukheddaden, *Symmetry* **2021**, *13* (5), 828.

[63] N. E. I. Belmouri, N. di Scala, K. Boukheddaden, *J. Phys. Chem. Solids* **2024**, *190*, 111985.

[64] J. A. Real, H. Bolvin, A. Bousseksou, A. Dworkin, O. Kahn, F. Varret, J. Zarembowitch, *J. Am. Chem. Soc.* **1992**, *114*, 4650–4658.

[65] J. A. Real, A. B. Gaspar, M. C. Muñoz, P. Gütlich, V. Ksenofontov, H. Spiering, *Top. Curr. Chem.* **2004**, *233*, 167.

[66] D. Fedaoui, Y. Bouhadja, K. Kaiba, P. Guionneau, J. Létard, *Eur. J. Inorg. Chem.* **2008**, *7*, 1022–1023.

[67] J. J. M. Amoore, C. J. Kepert, J. D. Cashion, B. Moubaraki, S. M. Neville, K. S. Murray, *Chem. Eur. J.* **2006**, *12*, 8220–8227.

[68] D. Pelleteter, R. Clérac, C. Mathonière, E. Harté, W. Schmitt, P. E. Kruger, *Chem. Commun.* **2009**, *2*, 221–223.

[69] R. J. Archer, H. S. Scott, M. I. J. Polson, B. E. Williamson, C. Mathonière, M. Rouzières, R. Clérac, P. E. Kruger, *Dalton Trans.* **2018**, *47*, 7965–7974.

[70] L. Li, A. R. Craze, R. Akiyoshi, A. Tsukiashi, S. Hayami, O. Mustonen, M. Bhadhade, S. Bhattacharyya, C. E. Marjo, Y. Wang, et al., *Dalton Trans.* **2018**, *47*, 2543–2548.
<figure><img src="image_28.png" /></figure>

[71] M. Koike, K. Murakami, T. Fujinami, K. Nishi, N. Matsumoto, Y. Sunatsuki, *Inorg. Chim. Acta* **2013**, *399*, 185-192.

[72] H. J. Shepherd, P. Rosa, L. Vendier, N. Casati, J. F. Létard, A. Bousseksou, P. Guionneau, G. Molnár, *Phys. Chem.* **2012**, *14*, 5265-5271.

[73] O. Roubeau, P. Gamez, S. J. Teat, *Eur. J. Inorg. Chem.* **2013**, *2013*, 934-942.

[74] R. W. Hogue, S. Singh, S. Brooker, *Chem. Soc. Rev.* **2018**, *47*, 7303.

[75] K. S. Murray, *Eur. J. Inorg. Chem.* **2008**, *2008*, 3101-3121.

[76] E. Milin, S. Belaid, V. Patinec, S. Triki, G. Chastanet, M. Marchivie, *Inorg. Chem.* **2016**, *55*, 9038.

[77] V. Ksenofontov, A. B. Gaspar, V. Niel, S. Reiman, J. A. Real, P. Gütlich, *Chem. Eur. J.* **2004**, *10* (5), 1291-1298.

[78] A. Bousseksou, F. Varret, J. Nasser, *J. Phys. I* **1993**, *3*, 1463-1473.

[79] K. Boukheddaden, I. Shteto, B. Hôo, F. Varret, *Phys. Rev. B* **2000**, *62*, 14796.

[80] K. Boukheddaden, I. Shteto, B. Hôo, F. Varret, *Phys. Rev. B* **2000**, *62*, 14806.

[81] C. Parreira, C. Enachescu, J. Linares, K. Boukheddaden, F. Varret, *J. Phys. Condens. Matter* **2000**, *12*, 9395.

[82] B. Hôo, K. Boukheddaden, F. Varret, *Eur. Phys. J. B* **2000**, *17*, 449.

[83] M. Nishino, K. Boukheddaden, S. Miyashita, F. Varret, *Phys. Rev. B* **2003**, *68*, 224402-224412.

[84] L. Turban, *J. Phys. A* **2016**, *49*, 355002.

[85] L. Turban, *J. Phys. C* **1982**, *15*, L65.

[86] J. M. Debierre, L. Turban, *J. Phys. A* **1983**, *16*, 3571.

[87] J. Linares, H. Spiering, F. Varret, *Eur. J. Phys. B* **1999**, *10*, 271.

[88] K. Boukheddaden, J. Linares, H. Spiering, F. Varret, *Eur. Phys. J. B* **2000**, *15*, 317.

[89] K. Boukheddaden, J. Linares, R. Tanasa, C. Chong, *J. Phys. Condens. Matter* **2007**, *19*, 106201.

[90] A. De Masi, *Spin Systems with Long Range Interactions*. In: P. Picco, J. San Martin, (eds) From Classical to Modern Probability. Progress in Probability, vol. 54. Birkhäuser, Basel **2003**.

[91] K. Fukui, S. Todo, *J. Comput. Phys.* **2009**, *228* (7), 2629-2642.

[92] A. Gudyma, I. Gudyma, *J. Appl. Phys.* **2021**, *129*, 12.

[93] Y. Garcia, C. M. Grunert, S. Reiman, O. van Campenhoudt, P. Gütlich, *Eur. J. Inorg. Chem.* **2006**, *2006*(17), 3333-3339.

[94] S. Doniach, *J. Chem. Phys.* **1978**, *68*, 11.

[95] A. Hauser, P. Gütlich, H. Spiering, *Inorg. Chem.* **1986**, *25* (23), 4245-4248.

[96] H. Spiering, N. Willembacher, *J. Phys. Condens. Matter* **1989**, *1* (50), 10089.

[97] B. Carazza, E. Rastelli, A. Tassi, *Z. Phys. B* **1991**, *84* (2), 301-308.

[98] M. G. Pini, A. Rettori, *Phys. Rev. B* **1993**, *48*, 3240.

[99] J. M. Yeomans, *Statistical mechanics of phase transitions*, Clarendon Press, **1992**.

[100] G. H. Golub, C. F. van Loan, *Matrix Computations*, Johns Hopkins University Press (3rd Edition), Baltimore, **1996**.

[101] F. Thétiot, S. Triki, K. Boukheddaden, G. Chastanet M Marchivie, *Chem. Mater.* **2017**, *29* (2), 490-494.

[102] Y. Garcia, V. Niel, M. C. Munoz, J. A. Real, *Top. Curr. Chem.* **2004**, *233*, 229-257.

[103] O. Roubeau, M. Castro, R. Burriel, J. G. Haasnoot, J. Reedijk, *Phys. Chem. B.* **2011**, *115*, 3003-3012.

[104] K. Boukheddaden, S. Miyashita, M. Nishino, *Phys. Rev. B* **2008**, *75*, 094112.

[105] R. Traiche, M. Sy, K. Boukheddaden, *J. Phys. Chem. C* **2018**, *122*, 4083-4096.

[106] S. Miyashita, Y. Konishi, M. Nishino, H. Tokoro, P. A. Rikvold, *Phys. Rev. B* **2008**, *77*, 014105.

[107] C. Zhang, L. Lian, H. Huang, S. Bala, Z. Ni, M. Tong, *Chem. Commun.* **2019**, *55*, 11033-11036.

[108] E. Trzop, D. Zhang, L. Piñeiro-Lopez, F. J. Valverde-Muñoz, M. Carmen Muñoz, L. Palatinus, L. Guerin, H. Cailleau, J. A. Real, E. Collet, *Angew. Chem. Int. Ed.* **2016**, *55*, 8675-8679.

[109] N. Pittala, F. Thétiot, C. Charles, S. Triki, K. Boukheddaden, G. Chastanet, M. Marchivie, *Chem. Commun.* **2017**, *53*, 8356-8359.

[110] N. Pittala, E. Cuza, D. Pinkowicz, M. Magott, M. Marchivie, K. Boukheddaden, S. Triki, *Inorg. Chem. Front.* **2022**, *9*, 6468.

Manuscript received: March 4, 2024

Revised manuscript received: May 16, 2024

Accepted manuscript online: June 5, 2024

Version of record online: October 7, 2024