# A comparative study between Monte Carlo entropic sampling method and local mean field investigations of thermal properties of spin-crossover nanoparticles based on Ising-like model

Jorge Linares, Catherine Cazelles, Pierre-Richard Dahoo, Nour-El-Islam Belmouri, Kamel Boukheddaden

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

## ► To cite this version:

Jorge Linares, Catherine Cazelles, Pierre-Richard Dahoo, Nour-El-Islam Belmouri, Kamel Boukheddaden. A comparative study between Monte Carlo entropic sampling method and local mean field investigations of thermal properties of spin-crossover nanoparticles based on Ising-like model. Chemical Synthesis, 2023, 3 (4), 10.20517/cs.2023.23. hal-04447558

HAL Id: hal-04447558

https://hal.uvsq.fr/hal-04447558

Submitted on 8 Feb 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>
# A comparative study between Monte Carlo entropic sampling method and local mean field investigations of thermal properties of spin-crossover nanoparticles based on Ising-like model

Jorge Linares, Catherine Cazelles, Pierre-Richard Dahoo, Nour-El-Islam Belmouri, Kamel Boukheddaden

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

## To cite this version:

Jorge Linares, Catherine Cazelles, Pierre-Richard Dahoo, Nour-El-Islam Belmouri, Kamel Boukheddaden. A comparative study between Monte Carlo entropic sampling method and local mean field investigations of thermal properties of spin-crossover nanoparticles based on Ising-like model. Chemical Synthesis, 2023, 3 (4), 10.20517/cs.2023.23. hal-04447558

HAL Id: hal-04447558

https://uvsq.hal.science/hal-04447558v1

Submitted on 8 Feb 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research Article

Open Access

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

# A comparative study between Monte Carlo entropic sampling method and local mean field investigations of thermal properties of spin-crossover nanoparticles based on Ising-like model

Jorge Linares<fnref n="1" />,<fnref n="2" />,* ID, Catherine Cazelles<fnref n="3" />,*, Pierre Richard Dahoo<fnref n="4" />, Nour-el-islam Belmouri<fnref n="1" />, Kamel Boukheddaden<fnref n="1" />

*Correspondence to: Profs. Jorge Linares, Kamel Boukheddaden, Groupe d'Etudes de la Matière Condensée, Université Paris-Saclay, UVSQ, CNRS, 45 Avenue des Etats Unis, Versailles 78035, France. E-mail: jorge.linares@uvsq.fr; Dr. Catherine Cazelles, IUT de Mantes en Yvelines, Université Paris-Saclay, UVSQ, 7 Rue Jean Hoët, Mantes la Jolie 78200, France. E-mail: cather-ine.cazelles@uvsq.fr

**How to cite this article:** Linares J, Cazelles C, Dahoo PR, Belmouri N, Boukheddaden K. A comparative study between Monte Carlo Entropic sampling Method and Local Mean field investigations of thermal properties of spin-crossover nanoparticles based on Ising-like model. *Chem Synth* 2023;3:45. http://dx.doi.org/10.20517/cs.2023.23

Received: 26 Apr 2023

First Decision: 19 Jul 2023

Revised: 2 Oct 2023

Accepted: 20 Oct 2023

Published: 18 Nov 2023

Academic Editors: Yann Garcia, Bao-Lian Su, Jiang Liu

Copy Editor: Dong-Li Li

Production Editor: Dong-Li Li

## Abstract

The thermally induced transitions between low-spin (LS) and high-spin (HS) configurations of spin-crossover (SCO) nanoparticles are simulated, focusing on the effects of localized surface and bulk interactions on the average magnetization of 2D square lattices. The thermal behaviors and hysteresis cycles are investigated within the framework of the Ising model Hamiltonian and are conducted following two approaches: local mean field approximation (LMFA) and Monte Carlo entropic sampling (MCES) techniques. The results obtained by these two methods are compared for the two square lattice sizes, 6 × 6 and 7 × 7. Thus, when the bulk-surface interaction term is set to zero, the two approaches lead to identical values of the surface and bulk transition temperatures separated by a long intermediate plateau in both cases. Although hysteresis curves exhibit a similar shape, LMFA shows slightly larger widths $\Delta T$ than MCES. On increasing bulk-surface interaction term, the two methods lead to different shifts in equilibrium temperature values for both bulk and surface components, respectively, to lower and higher values by MCES. In general, it is found that LMFA shifts surface equilibrium temperature differently to lower values and enhances the hysteresis effect, particularly for surface molecules. On the other hand, for the 7 × 7 square lattice, the equilibrium temperatures are slightly higher by 1.5% and 3.2% for bulk and surface molecules, respectively, with a narrower hysteresis width in the

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

© The Author(s) 2023. **Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, sharing, adaptation, distribution and reproduction in any medium or format, for any purpose, even commercially, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

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

<footnote n="1">¹Groupe d'Etudes de la Matière Condensée, Université Paris-Saclay, UVSQ, CNRS, Versailles 78035, France.</footnote>

<footnote n="2">²Departamento de Ciencias, Sección Física, Pontificia Universidad Católica del Perú, Lima 15088, Perú.</footnote>

<footnote n="3">³IUT de Mantes en Yvelines, Université Paris-Saclay, UVSQ, Mantes la Jolie 78200, France.</footnote>

<footnote n="4">⁴LATMOS, Université Paris-Saclay, UVSQ, Sorbonne Université, CNRS, Guyancourt 78280, France.</footnote>
Research Article

Open Access

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

# A comparative study between Monte Carlo entropic sampling method and local mean field investigations of thermal properties of spin-crossover nanoparticles based on Ising-like model

Jorge Linares<fnref n="1" />,<fnref n="2" />,* ID, Catherine Cazelles<fnref n="3" />,*, Pierre Richard Dahoo<fnref n="4" />, Nour-el-islam Belmouri<fnref n="1" />, Kamel Boukheddaden<fnref n="1" />

*Correspondence to: Profs. Jorge Linares, Kamel Boukheddaden, Groupe d'Etudes de la Matière Condensée, Université Paris-Saclay, UVSQ, CNRS, 45 Avenue des Etats Unis, Versailles 78035, France. E-mail: jorge.linares@uvsq.fr; Dr. Catherine Cazelles, IUT de Mantes en Yvelines, Université Paris-Saclay, UVSQ, 7 Rue Jean Hoët, Mantes la Jolie 78200, France. E-mail: cather-ine.cazelles@uvsq.fr

**How to cite this article:** Linares J, Cazelles C, Dahoo PR, Belmouri N, Boukheddaden K. A comparative study between Monte Carlo Entropic sampling Method and Local Mean field investigations of thermal properties of spin-crossover nanoparticles based on Ising-like model. *Chem Synth* 2023;3:45. http://dx.doi.org/10.20517/cs.2023.23

Received: 26 Apr 2023

First Decision: 19 Jul 2023

Revised: 2 Oct 2023

Accepted: 20 Oct 2023

Published: 18 Nov 2023

**Academic Editors:** Yann Garcia, Bao-Lian Su, Jiang Liu

Copy Editor: Dong-Li Li

Production Editor: Dong-Li Li

## Abstract

The thermally induced transitions between low-spin (LS) and high-spin (HS) configurations of spin-crossover (SCO) nanoparticles are simulated, focusing on the effects of localized surface and bulk interactions on the average magnetization of 2D square lattices. The thermal behaviors and hysteresis cycles are investigated within the framework of the Ising model Hamiltonian and are conducted following two approaches: local mean field approximation (LMFA) and Monte Carlo entropic sampling (MCES) techniques. The results obtained by these two methods are compared for the two square lattice sizes, 6 × 6 and 7 × 7. Thus, when the bulk-surface interaction term is set to zero, the two approaches lead to identical values of the surface and bulk transition temperatures separated by a long intermediate plateau in both cases. Although hysteresis curves exhibit a similar shape, LMFA shows slightly larger widths $\Delta T$ than MCES. On increasing bulk-surface interaction term, the two methods lead to different shifts in equilibrium temperature values for both bulk and surface components, respectively, to lower and higher values by MCES. In general, it is found that LMFA shifts surface equilibrium temperature differently to lower values and enhances the hysteresis effect, particularly for surface molecules. On the other hand, for the 7 × 7 square lattice, the equilibrium temperatures are slightly higher by 1.5% and 3.2% for bulk and surface molecules, respectively, with a narrower hysteresis width in the

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

© The Author(s) 2023. **Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, sharing, adaptation, distribution and reproduction in any medium or format, for any purpose, even commercially, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

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

<footnote n="1">¹Groupe d'Etudes de la Matière Condensée, Université Paris-Saclay, UVSQ, CNRS, Versailles 78035, France.</footnote>

<footnote n="2">²Departamento de Ciencias, Sección Física, Pontificia Universidad Católica del Perú, Lima 15088, Perú.</footnote>

<footnote n="3">³IUT de Mantes en Yvelines, Université Paris-Saclay, UVSQ, Mantes la Jolie 78200, France.</footnote>

<footnote n="4">⁴LATMOS, Université Paris-Saclay, UVSQ, Sorbonne Université, CNRS, Guyancourt 78280, France.</footnote>
surface. Moreover, with the MCES method, an abrupt transition instead of a hysteresis transition is calculated for surface molecules ($ΔT_{surf}^{MCES} = 0$ K).

**Keywords:** Monte Carlo simulation, local mean-field approach, spin crossover, phase transitions, nanoparticles, bulk-surface interactions

## INTRODUCTION

Spin-crossover (SCO) materials [1–6] are composed of transition metals with an electronic configuration ranging between $3d^4$ and $3d^7$ (chromium, manganese, iron, and cobalt). These metals coordinate with (Sulphur, Oxygen, or Nitrogen) atoms in octahedral symmetry. Iron (II)-based SCO solids constitute interesting materials due to their bistable electronic and vibronic nature that allows them to exhibit, as a function of various stimuli, a spin switching between a diamagnetic low-spin state (LS, $e_g^0 t_{2g}^6$) with total spin $S = 0$ and a paramagnetic high-spin (HS, $e_g^2 t_{2g}^4$) with a spin value, $S = 2$. Common external stimuli are temperature, pressure, electric fields, high-magnetic fields, and light [7–20].

In Fe-based SCO materials, the spin transition between LS and HS is accompanied by an increase of ~ 10% of the Fe-ligand bond lengths due to the weakening of the metal-ligand bonds accompanying the population of the anti-bonding orbitals $e_g^2$. This results in the change of the unit cell volume by ~3%-5%. The SCO molecular solids are excellent candidates as sensors of temperature, pressure [21–24], and gases, notably of hazardous volatile organic compounds, as recently demonstrated in several reports [25,26]. They can also serve as switches of light emission in the electroluminescent devices or photoluminescent [27,28] containing SCO complexes or to enhance plasmonic resonances [29].

To try to explain and reproduce these switching phenomena from the LS to HS states, researchers have proposed different models, some of which are based on Ising-like [30–35], mechano-elastic [36], and electroelastic [37] models, on which switchable core-shell SCO nanoparticles have been investigated in details [38,39].

In the present study, we focus on the modeling of thermal effects in 2D square SCO nanoparticles. The Ising-like model is used here and solved in the framework of local mean field approximation (LMFA) and within the Monte Carlo entropic sampling (MCES) technique in an attempt to reveal the role of interactions with the external environment (matrix). In this work, we discuss and compare the results obtained by these two methods. Experimental results on the effect of the surface or the substrate have been reported [35,40,41].

Although the MCES technique is a good choice to solve the Hamiltonian associated with the extended Ising-like model, it is limited in terms of computing capacity for samples containing a small number of molecules. As we would like to analyze the cases of nanoparticles with large numbers of molecules and even thin films in contact with a substrate, which cannot be studied by MCES, we plan to use an alternative formalism based on LMFA. With this contribution, the purpose of this contribution is to determine the reliability, limitations, and drawbacks of the LMFA method when applied to very small systems by comparing it with the MCES technique.

The manuscript is organized as follows. Section 2 introduces the Model and some details of the parameters used in the model. In section 3, the MCES technique is presented, while the LMFA is introduced in section 4. In section 5, we present the simulated thermal behavior of order parameters and the discussion of the obtained results. Finally, we conclude and outline the possible extensions of this work in section 6.
surface. Moreover, with the MCES method, an abrupt transition instead of a hysteresis transition is calculated for surface molecules ($ΔT_{surf}^{MCES} = 0$ K).

**Keywords:** Monte Carlo simulation, local mean-field approach, spin crossover, phase transitions, nanoparticles, bulk-surface interactions

## INTRODUCTION

Spin-crossover (SCO) materials [1-6] are composed of transition metals with an electronic configuration ranging between 3d⁴ and 3d⁷ (chromium, manganese, iron, and cobalt). These metals coordinate with (Sulphur, Oxygen, or Nitrogen) atoms in octahedral symmetry. Iron (II)-based SCO solids constitute interesting materials due to their bistable electronic and vibronic nature that allows them to exhibit, as a function of various stimuli, a spin switching between a diamagnetic low-spin state (LS, $e_g^0 t_{2g}^6$) with total spin S = 0 and a paramagnetic high-spin (HS, $e_g^2 t_{2g}^4$) with a spin value, S = 2. Common external stimuli are temperature, pressure, electric fields, high-magnetic fields, and light [7-20].

In Fe-based SCO materials, the spin transition between LS and HS is accompanied by an increase of ~ 10% of the Fe-ligand bond lengths due to the weakening of the metal-ligand bonds accompanying the population of the anti-bonding orbitals $e_g^2$. This results in the change of the unit cell volume by ~3%-5%. The SCO molecular solids are excellent candidates as sensors of temperature, pressure [21-24], and gases, notably of hazardous volatile organic compounds, as recently demonstrated in several reports [25,26]. They can also serve as switches of light emission in the electroluminescent devices or photoluminescent [27,28] containing SCO complexes or to enhance plasmonic resonances [29].

To try to explain and reproduce these switching phenomena from the LS to HS states, researchers have proposed different models, some of which are based on Ising-like [30-35], mechano-elastic [36], and electroelastic [37] models, on which switchable core-shell SCO nanoparticles have been investigated in details [38,39].

In the present study, we focus on the modeling of thermal effects in 2D square SCO nanoparticles. The Ising-like model is used here and solved in the framework of local mean field approximation (LMFA) and within the Monte Carlo entropic sampling (MCES) technique in an attempt to reveal the role of interactions with the external environment (matrix). In this work, we discuss and compare the results obtained by these two methods. Experimental results on the effect of the surface or the substrate have been reported [35,40,41].

Although the MCES technique is a good choice to solve the Hamiltonian associated with the extended Ising-like model, it is limited in terms of computing capacity for samples containing a small number of molecules. As we would like to analyze the cases of nanoparticles with large numbers of molecules and even thin films in contact with a substrate, which cannot be studied by MCES, we plan to use an alternative formalism based on LMFA. With this contribution, the purpose of this contribution is to determine the reliability, limitations, and drawbacks of the LMFA method when applied to very small systems by comparing it with the MCES technique.

The manuscript is organized as follows. Section 2 introduces the Model and some details of the parameters used in the model. In section 3, the MCES technique is presented, while the LMFA is introduced in section 4. In section 5, we present the simulated thermal behavior of order parameters and the discussion of the obtained results. Finally, we conclude and outline the possible extensions of this work in section 6.
## THE MODEL

In condensed matter physics, magnetic phenomena observed in a regular lattice are described with the well-known Ising Model based on interactions between quantum spin variables that were first introduced phenomenologically by Pauli [42]. Later on, in 1925, Ising [43] used the concept of spins to describe the thermal behavior of the magnetization of magnetic materials using a 1D model whose statistics were treated in the frame of a transfer matrix method. Pointing quickly the exact resolution of this model (with zero applied field) performed with a remarkable mathematical tour de force of Onsager in 1944 [44], it is also important to notice the considerable efforts of several authors in the realization of approximative methods based on variational treatments of the Ising model, which led to mean-field theory [45], Bethe-Peierls [46], etc., which are today widely used to analytically describe a large panel of cooperative phenomena in solid state physics and beyond.

In the field of molecular switchable materials, the Ising model [30–35] has been adapted to study SCO materials in which two different magnetic states, HS and LS, interact to display at the macroscopic scale various physical properties that give rise to a wide variety of transition behaviors. Among them, one can quote: (i) thermally-induced hysteretic first-order phase transitions [15], gradual or continuous [47] conversions, incomplete [48], and multi-step transformations [49].

To model SCO nanoparticles, a spin fictitious operator with two eigen-values +1 and -1 is associated with each molecule describing the HS and LS states, respectively. Then, the Ising-like Hamiltonian is used to model the energy operator connected to SCO. In this respect, the Hamiltonian $H$ can be written as the sum of three energetic contributions related to the isolated molecule, the interaction between molecules, and finally, between the molecules at the boundary with their immediate external environment. Thus, $H$ is expressed as:

$$H = H_1 + H_2 + H_3 \qquad (1)$$

with

(i) The single molecule contribution:

$$H_1 = \frac{\Delta}{2} \sum_{i=1}^{N} \sigma_i \qquad (2)$$

where $N$ is the number of molecules, and $\Delta$ is the energy gap at zero Kelvin between the fundamental LS state and the excited HS state.

(ii) The interaction term:

$$H_2 = - \sum_{i,j} J_{ij} \sigma_i \sigma_j \qquad (3)$$

where $J_{ij}$ is the coupling constant between the spins of the sites $i$ and $j$.

(iii) The interaction between the surface molecules and their environment:
## THE MODEL

In condensed matter physics, magnetic phenomena observed in a regular lattice are described with the well-known Ising Model based on interactions between quantum spin variables that were first introduced phenomenologically by Pauli [42]. Later on, in 1925, Ising [43] used the concept of spins to describe the thermal behavior of the magnetization of magnetic materials using a 1D model whose statistics were treated in the frame of a transfer matrix method. Pointing quickly the exact resolution of this model (with zero applied field) performed with a remarkable mathematical tour de force of Onsager in 1944 [44], it is also important to notice the considerable efforts of several authors in the realization of approximative methods based on variational treatments of the Ising model, which led to mean-field theory [45], Bethe-Peierls [46], etc., which are today widely used to analytically describe a large panel of cooperative phenomena in solid state physics and beyond.

In the field of molecular switchable materials, the Ising model [30–35] has been adapted to study SCO materials in which two different magnetic states, HS and LS, interact to display at the macroscopic scale various physical properties that give rise to a wide variety of transition behaviors. Among them, one can quote: (i) thermally-induced hysteretic first-order phase transitions [15], gradual or continuous [47] conversions, incomplete [48], and multi-step transformations [49].

To model SCO nanoparticles, a spin fictitious operator with two eigen-values +1 and -1 is associated with each molecule describing the HS and LS states, respectively. Then, the Ising-like Hamiltonian is used to model the energy operator connected to SCO. In this respect, the Hamiltonian $H$ can be written as the sum of three energetic contributions related to the isolated molecule, the interaction between molecules, and finally, between the molecules at the boundary with their immediate external environment. Thus, $H$ is expressed as:

$$H = H_1 + H_2 + H_3 \quad (1)$$

with

(i) The single molecule contribution:

$$H_1 = \frac{\Delta}{2} \sum_{i=1}^{N} \sigma_i \quad (2)$$

where $N$ is the number of molecules, and $\Delta$ is the energy gap at zero Kelvin between the fundamental LS state and the excited HS state.

(ii) The interaction term:

$$H_2 = - \sum_{i,j} J_{ij} \sigma_i \sigma_j \quad (3)$$

where $J_{ij}$ is the coupling constant between the spins of the sites $i$ and $j$.

(iii) The interaction between the surface molecules and their environment:
$$H_3 = -L \sum_{k=1}^{M} \sigma_k \qquad (4)$$

where L is the energy state of each molecule at the surface in addition to the previous ligand $\Delta/2$. The variable, M, is the number of molecules $\sigma_k$ at the surface.

In the following, we give a detailed description of each of these terms.

### Molecular contribution $H_1$

The two states, HS and LS, have different spin and orbital degeneracies, and the spectrum of intramolecular vibrations of the molecule strongly depends on the spin state of the central metal.

And, as it has been reported [1,2], the density of states becomes higher in the HS state than in the LS state, which is then equivalent to considering these two states as having different effective degeneracies, denoted here by $g_{HS}$ and $g_{LS}$. For simplicity, we consider $g_{HS}$ and $g_{LS}$ as independent of temperature. From the energetic point of view, these effective degeneracies enter the Hamiltonian as entropic terms, of which only their ratio plays a relevant role. When these specific characteristics of the two states, HS and LS, are considered, the expression of the molecular term $H_1$ is written as:

$$H_1 = \frac{\Delta - k_B T \ln g}{2} \sum_{i=1}^{N} \sigma_i \qquad (5)$$

Where $g = \frac{g_{HS}}{g_{LS}}$, T is the temperature, and $k_B$ is the usual Boltzmann constant.

At this stage, we see that from the point of view of a single site Hamiltonian, as far as $\Delta - k_B T \ln g > 0$, LS states, for which $\sigma = -1$, are favored. When $\Delta - k_B T \ln g < 0$, the HS states, for which $\sigma = +1$ are favored. At the equilibrium temperature $T_{eq} = \frac{\Delta}{k_B \ln g}$, there is a “crossover” situation between these two configurations. The physical quantity that allows the thermal behavior of the system is the HS fraction, denoted here as $Nhs$, which represents the probability of the population being in the HS state. In terms of spin values, it represents the probability of occupation of the spin value $\sigma = +1$. It is straightforward to establish that $Nhs$ relates to the average value of the fictitious magnetization $<\sigma>$ as follows:

$$Nhs = \frac{1 + <\sigma>}{2} \qquad (6)$$

Simple Boltzmann statistics performed on Hamiltonian (1) leads to the following relations,

$$<\sigma> = \tanh \beta \left( \frac{\Delta - k_B T \ln g}{2} \right) \qquad (7)$$
$$H_3 = -L \sum_{k=1}^{M} \sigma_k \qquad (4)$$

where L is the energy state of each molecule at the surface in addition to the previous ligand $\Delta/2$. The variable, M, is the number of molecules $\sigma_k$ at the surface.

In the following, we give a detailed description of each of these terms.

### Molecular contribution $H_1$

The two states, HS and LS, have different spin and orbital degeneracies, and the spectrum of intramolecular vibrations of the molecule strongly depends on the spin state of the central metal.

And, as it has been reported [1,2], the density of states becomes higher in the HS state than in the LS state, which is then equivalent to considering these two states as having different effective degeneracies, denoted here by $g_{HS}$ and $g_{LS}$. For simplicity, we consider $g_{HS}$ and $g_{LS}$ as independent of temperature. From the energetic point of view, these effective degeneracies enter the Hamiltonian as entropic terms, of which only their ratio plays a relevant role. When these specific characteristics of the two states, HS and LS, are considered, the expression of the molecular term $H_1$ is written as:

$$H_1 = \frac{\Delta - k_B T \ln g}{2} \sum_{i=1}^{N} \sigma_i \qquad (5)$$

Where $g = g_{HS}/g_{LS}$, T is the temperature, and $k_B$ is the usual Boltzmann constant.

At this stage, we see that from the point of view of a single site Hamiltonian, as far as $\Delta - k_B T \ln g > 0$, LS states, for which $\sigma = -1$, are favored. When $\Delta - k_B T \ln g < 0$, the HS states, for which $\sigma = +1$ are favored. At the equilibrium temperature $T_{eq} = \frac{\Delta}{k_B \ln g}$, there is a “crossover” situation between these two configurations. The physical quantity that allows the thermal behavior of the system is the HS fraction, denoted here as $Nhs$, which represents the probability of the population being in the HS state. In terms of spin values, it represents the probability of occupation of the spin value $\sigma = +1$. It is straightforward to establish that $Nhs$ relates to the average value of the fictitious magnetization $<\sigma>$ as follows:

$$Nhs = \frac{1 + <\sigma>}{2} \qquad (6)$$

Simple Boltzmann statistics performed on Hamiltonian (1) leads to the following relations,

$$<\sigma> = \tanh \beta \left( \frac{\Delta - k_B T \ln g}{2} \right) \qquad (7)$$
<figure><img src="image_6.png" /><figcaption>**Figure 1.** Typical thermal evolution (A) of the average magnetization <σ>; (B) of the high spin fraction Nhs. The calculation parameters are: Δ/k<sub>B</sub> = 1, 300 ~ K and ln(g) = 6.01.</figcaption></figure>

and

$$Nhs = \frac{1}{2} \left[ 1 + \tanh \beta \left( \frac{\Delta - k_B T \ln g}{2} \right) \right] \quad (8)$$

from which the curves of Figure 1 showing thermal dependence of the HS magnetization and fraction are easily derived.

### Interaction term H<sub>2</sub>

The interaction term can be split into two contributions: (i) those arising from nearest neighbors, which are parametrized using a short-range interaction term J (assumed here to be independent of the site); and (ii) those including all the other molecules, which are parametrized with a long-range interaction term G, which is written in the frame of the meanfield approach. In this respect, H<sub>2</sub> is written as:

$$H_2 = -J \sum_{<i,j>} \sigma_i \sigma_j - G \langle \sigma \rangle \sum_{i=1}^{N} \sigma_i \quad (9)$$

Considering solely this interaction term H<sub>2</sub>, its resolution leads to the thermal evolution of <σ>, as given in Figure 2, which is typical of second-order phase transition. Indeed, at 0 K, the system is either in the ordered “HS” state with all the spins with σ = 1 or in the “LS” state with σ = -1. When temperature increases, the average value of the spin either decreases or increases, as shown in the upper and lower curves. The temperature at which the average value is zero is called the order-disorder temperature T<sub>O/D</sub> or critical temperature T<sub>C</sub>.

The region on the left-hand side below the upper curve (respectively above the lower curve) with respect to this T<sub>OD</sub> temperature corresponds to an ordered phase, while the right-hand side corresponds to the disordered phase. For an infinite square lattice with G = 0, the critical temperature T<sub>C</sub> = 2.269 ~ J is obtained in line with the “Onsager solution”, where T<sub>C</sub> is given by the equation tanh (2J/k<sub>B</sub>T<sub>C</sub>) = 1/√2.

However, for G ≠ 0, there is no exact analytical solution since the mean-field term acts as an “external field” linearly dependent on the magnetization. On the other hand, for a finite system, T<sub>C</sub> depends on the lattice size. It is worth noting the symmetric character of the thermal dependence of the net magnetization around <σ> = 0 for the ordered phases, obtained by starting from +1 or -1, as shown in Figure 2, calculated using J and G in a 6 × 6 square lattice. It is important to mention that this symmetric character remains valid when
<figure><img src="image_6.png" /><figcaption>**Figure 1.** Typical thermal evolution (A) of the average magnetization <σ>; (B) of the high spin fraction Nhs. The calculation parameters are: Δ/k<sub>B</sub> = 1, 300 ~ K and ln(g) = 6.01.</figcaption></figure>

and

$$Nhs = \frac{1}{2} \left[ 1 + \tanh \beta \left( \frac{\Delta - k_B T \ln g}{2} \right) \right] \quad (8)$$

from which the curves of Figure 1 showing thermal dependence of the HS magnetization and fraction are easily derived.

### Interaction term $H_2$

The interaction term can be split into two contributions: (i) those arising from nearest neighbors, which are parametrized using a short-range interaction term J (assumed here to be independent of the site); and (ii) those including all the other molecules, which are parametrized with a long-range interaction term G, which is written in the frame of the meanfield approach. In this respect, $H_2$ is written as:

$$H_2 = -J \sum_{<i,j>} \sigma_i \sigma_j - G <\sigma> \sum_{i=1}^{N} \sigma_i \quad (9)$$

Considering solely this interaction term $H_2$, its resolution leads to the thermal evolution of <σ>, as given in Figure 2, which is typical of second-order phase transition. Indeed, at 0 K, the system is either in the ordered “HS” state with all the spins with σ = 1 or in the “LS” state with σ = -1. When temperature increases, the average value of the spin either decreases or increases, as shown in the upper and lower curves. The temperature at which the average value is zero is called the order-disorder temperature $T_{O/D}$ or critical temperature $T_C$.

The region on the left-hand side below the upper curve (respectively above the lower curve) with respect to this $T_{OD}$ temperature corresponds to an ordered phase, while the right-hand side corresponds to the disordered phase. For an infinite square lattice with G = 0, the critical temperature $T_C = 2.269 \sim J$ is obtained in line with the “Onsager solution”, where $T_C$ is given by the equation $\tanh \frac{2J}{k_B T_C} = \frac{1}{\sqrt{2}}$.

However, for G ≠ 0, there is no exact analytical solution since the mean-field term acts as an “external field” linearly dependent on the magnetization. On the other hand, for a finite system, $T_C$ depends on the lattice size. It is worth noting the symmetric character of the thermal dependence of the net magnetization around <σ> = 0 for the ordered phases, obtained by starting from +1 or -1, as shown in Figure 2, calculated using J and G in a 6 × 6 square lattice. It is important to mention that this symmetric character remains valid when
<figure><img src="image_7.png" /><figcaption>**Figure 2.** Thermal evolution of the average magnetization $<\sigma>$ obtained in a 6 × 6 square lattice with the following parameters: $J/k_B = 60$ K and $G/k_B = 172$ K.</figcaption></figure>

including the long-range interaction term in the Hamiltonian $H_2$, although the expression of $T_C$ now becomes dependent on $J$ and $G$ [14].

### Short- and long-range interaction contributions $H_1 + H_2$

When the first two terms are combined to model SCO physical behavior with temperature, the Hamiltonian writes as:

$$H = \frac{\Delta - k_B T \ln g - 2G < \sigma >}{2} \sum_{i=1}^{N} \sigma_i - J \sum_{\langle i,j \rangle} \sigma_i \sigma_j \quad (10)$$

The symmetry in the configurations consisting of positive and negative sigma values about $<\sigma> = 0$ of the interaction term $H_2$ is broken when the latter is superposed with the molecular term $H_1$ under the condition that $T_{eq} < T_c$.

In this case, negative solutions of $<\sigma>$ occur when $T < T_{eq}$ ($= \frac{\Delta}{k_B \ln g}$) while positive values are favored when $T > T_{eq}$. Thus, the molecular term $H_1$ drives the configurations under the constraint imposed by the term $H_2$.

Therefore, if $T_{eq} < T_c$, the $<\sigma>$ solution shifts from a $<\sigma> < 0$ at low temperatures to the solution with $<\sigma> > 0$ at higher temperatures, which indicates the existence of a phase transition.

### Surface contribution $H_3$

This $H_3$ term is restricted to the molecules at the surface. The surface Hamiltonian $H_S$ can thus be written as:

$$H_S = \frac{\Delta - k_B T \ln g - 2G < \sigma >}{2} \sum_{k=1}^{M} \sigma_k - J \sum_{\langle i,j \rangle \in \text{surface}} \sigma_i \sigma_j - L \sum_{k=1}^{M} \sigma_k \quad (11)$$

with $M$ the number of molecules at the surface. Here, for simplicity, the short-range interaction term, $J$, between nn spins is considered to be the same at the surface and in the bulk material. Upon factorization, we obtain:
<figure><img src="image_7.png" /><figcaption>**Figure 2.** Thermal evolution of the average magnetization $<\sigma>$ obtained in a 6 × 6 square lattice with the following parameters: $J/k_B = 60$ K and $G/k_B = 172$ K.</figcaption></figure>

including the long-range interaction term in the Hamiltonian $H_2$, although the expression of $T_C$ now becomes dependent on $J$ and $G$ [14].

### Short- and long-range interaction contributions $H_1 + H_2$

When the first two terms are combined to model SCO physical behavior with temperature, the Hamiltonian writes as:

$$H = \frac{\Delta - k_B T \ln g - 2G < \sigma >}{2} \sum_{i=1}^{N} \sigma_i - J \sum_{\langle i,j \rangle} \sigma_i \sigma_j \quad (10)$$

The symmetry in the configurations consisting of positive and negative sigma values about $<\sigma> = 0$ of the interaction term $H_2$ is broken when the latter is superposed with the molecular term $H_1$ under the condition that $T_{eq} < T_c$.

In this case, negative solutions of $<\sigma>$ occur when $T < T_{eq}$ ($= \frac{\Delta}{k_B \ln g}$) while positive values are favored when $T > T_{eq}$. Thus, the molecular term $H_1$ drives the configurations under the constraint imposed by the term $H_2$.

Therefore, if $T_{eq} < T_c$, the $<\sigma>$ solution shifts from a $<\sigma> < 0$ at low temperatures to the solution with $<\sigma> > 0$ at higher temperatures, which indicates the existence of a phase transition.

### Surface contribution $H_3$

This $H_3$ term is restricted to the molecules at the surface. The surface Hamiltonian $H_S$ can thus be written as:

$$H_S = \frac{\Delta - k_B T \ln g - 2G < \sigma >}{2} \sum_{k=1}^{M} \sigma_k - J \sum_{\langle i,j \rangle \in \text{surface}} \sigma_i \sigma_j - L \sum_{k=1}^{M} \sigma_k \quad (11)$$

with $M$ the number of molecules at the surface. Here, for simplicity, the short-range interaction term, $J$, between nn spins is considered to be the same at the surface and in the bulk material. Upon factorization, we obtain:
$$H_S = \frac{\Delta - k_B T \ln g - 2L - 2G < \sigma >}{2} \sum_{k=1}^{M} \sigma_i - J \sum_{\langle i,j \rangle} \sigma_i \sigma_j \quad (12)$$

In order to follow how the equilibrium temperature of the surface molecules ($T_{eq}$)$_S$ changes with the value of the external interaction $L$, it is calculated under the condition that the longrange coupling term is set to zero ($G = 0$), and then it can be expressed as:

$$(T_{eq})_{S,G=0} = \frac{\Delta - 2L}{k_B \ln g} \quad (13)$$

However, it must be noted that the long-range interaction $G$ applies to all molecules, and it may happen that the average at the surface is equal to zero, while at the bulk, this is not the case.

### The total Hamiltonian $H_T$

Now, as one considers the total Hamiltonian, one gets:

$$H = H_1 + H_2 + H_3 = \frac{\Delta - k_B T \ln g}{2} \sum_{i=1}^{N} \sigma_i - J \sum_{\langle i,j \rangle} \sigma_i \sigma_j - G \sum_{i=1}^{N} \sigma_i < \sigma > - L \sum_{k=1}^{M} \sigma_k \quad (14)$$

which can be written as:

$$H_T = \frac{\Delta - k_B T \ln g - 2G < \sigma >}{2} \sum_{i=1}^{N} \sigma_i - J \sum_{\langle i,j \rangle} \sigma_i \sigma_j - L \sum_{k=1}^{M} \sigma_k \quad (15)$$

In the present work, the Hamiltonian $H_T$ is solved following two different approaches, namely the (i) LMFA and (ii) MCES techniques, in order to determine the thermal evolution of the average fictitious magnetization, $< \sigma >$.

### RESOLUTION BY MONTE CARLO ENTROPIC SAMPLING METHOD

In this section, we will distinguish three kinds of short-range interaction terms, namely, $J_{bb}$, the interaction between molecules in the bulk, $J_{ss}$, the interaction between the molecules at the surface, and $J_{bs}$, the interaction between the molecules at the surface and in the bulk. The corresponding Hamiltonian is written as:

$$H = \frac{\Delta - k_B T \ln g - 2G < \sigma >}{2} \sum_{i=1}^{N} \sigma_i - J_{bb} \sum_{bulk} \sigma_i \sigma_j - J_{ss} \sum_{surf} \sigma_i \sigma_j - J_{bs} \sum_{B-5} \sigma_i \sigma_j - L \sum_{k=1}^{M} \sigma_k \quad (16)$$
$$H_S = \frac{\Delta - k_B T \ln g - 2L - 2G < \sigma >}{2} \sum_{k=1}^{M} \sigma_i - J \sum_{\langle i,j \rangle} \sigma_i \sigma_j \quad (12)$$

In order to follow how the equilibrium temperature of the surface molecules ($T_{eq}$)$_S$ changes with the value of the external interaction $L$, it is calculated under the condition that the longrange coupling term is set to zero ($G = 0$), and then it can be expressed as:

$$(T_{eq})_{S,G=0} = \frac{\Delta - 2L}{k_B \ln g} \quad (13)$$

However, it must be noted that the long-range interaction $G$ applies to all molecules, and it may happen that the average at the surface is equal to zero, while at the bulk, this is not the case.

## The total Hamiltonian $H_T$

Now, as one considers the total Hamiltonian, one gets:

$$H = H_1 + H_2 + H_3 = \frac{\Delta - k_B T \ln g}{2} \sum_{i=1}^{N} \sigma_i - J \sum_{\langle i,j \rangle} \sigma_i \sigma_j - G \sum_{i=1}^{N} \sigma_i < \sigma > - L \sum_{k=1}^{M} \sigma_k \quad (14)$$

which can be written as:

$$H_T = \frac{\Delta - k_B T \ln g - 2G < \sigma >}{2} \sum_{i=1}^{N} \sigma_i - J \sum_{\langle i,j \rangle} \sigma_i \sigma_j - L \sum_{k=1}^{M} \sigma_k \quad (15)$$

In the present work, the Hamiltonian $H_T$ is solved following two different approaches, namely the (i) LMFA and (ii) MCES techniques, in order to determine the thermal evolution of the average fictitious magnetization, $< \sigma >$.

## RESOLUTION BY MONTE CARLO ENTROPIC SAMPLING METHOD

In this section, we will distinguish three kinds of short-range interaction terms, namely, $J_{bb}$, the interaction between molecules in the bulk, $J_{ss}$, the interaction between the molecules at the surface, and $J_{bs}$, the interaction between the molecules at the surface and in the bulk. The corresponding Hamiltonian is written as:

$$H = \frac{\Delta - k_B T \ln g - 2G < \sigma >}{2} \sum_{i=1}^{N} \sigma_i - J_{bb} \sum_{bulk} \sigma_i \sigma_j - J_{ss} \sum_{surf} \sigma_i \sigma_j - J_{bs} \sum_{B-5} \sigma_i \sigma_j - L \sum_{k=1}^{M} \sigma_k \quad (16)$$
with

$$h = -\frac{\Delta - k_B T \ln(g) - 2G \frac{m_t}{N}}{2} \qquad (17)$$

where the respective total and surface magnetizations, $m_s$ and $m_t$, respectively, are given by:

$$m_t = \sum_{i=1}^{N} \sigma_i \qquad (18)$$

$$m_s = \sum_{k=1}^{M} \sigma_k \qquad (19)$$

H can be re-written as:

$$H = -hm_t - J_{bb} \times s_b - J_{ss} \times s_s - J_{bs} \times s_{bs} - Lm_s \qquad (20)$$

Where $s_b$, $s_s$, and $s_{bs}$ are the short-range correlation functions related to the bulk, surface, and bulk-surface neighboring sites, respectively. Their expressions are given by the following relations:

$$s_b = \sum_{\text{bulk}} \sigma_i \sigma_j \qquad (21)$$

$$s_s = \sum_{\text{surface}} \sigma_i \sigma_j \qquad (22)$$

$$s_{bs} = \sum_{B-s} \sigma_i \sigma_j \qquad (23)$$

The thermal average value of the fictitious magnetization, $\langle \sigma \rangle$, is calculated by the following expression:

$$\langle \sigma \rangle = \frac{\sum_{i=1}^{N_L} \frac{m_i}{N} d(m_{t_i}, s_{b_i}, s_{s_i}, s_{bs_i}, m_{s_i}) \exp(-\beta E_i)}{Z} \qquad (24)$$
with

$$h = -\frac{\Delta - k_B T \ln(g) - 2G \frac{m_t}{N}}{2} \qquad (17)$$

where the respective total and surface magnetizations, $m_s$ and $m_t$, respectively, are given by:

$$m_t = \sum_{i=1}^{N} \sigma_i \qquad (18)$$

$$m_s = \sum_{k=1}^{M} \sigma_k \qquad (19)$$

$H$ can be re-written as:

$$H = -hm_t - J_{bb} \times s_b - J_{ss} \times s_s - J_{bs} \times s_{bs} - Lm_s \qquad (20)$$

Where $s_b$, $s_s$, and $s_{bs}$ are the short-range correlation functions related to the bulk, surface, and bulk-surface neighboring sites, respectively. Their expressions are given by the following relations:

$$s_b = \sum_{\text{bulk}} \sigma_i \sigma_j \qquad (21)$$

$$s_s = \sum_{\text{surface}} \sigma_i \sigma_j \qquad (22)$$

$$s_{bs} = \sum_{B-s} \sigma_i \sigma_j \qquad (23)$$

The thermal average value of the fictitious magnetization, $\langle \sigma \rangle$, is calculated by the following expression:

$$\langle \sigma \rangle = \frac{\sum_{i=1}^{N_L} \frac{m_i}{N} d(m_{t_i}, s_{b_i}, s_{s_i}, s_{bs_i}, m_{s_i}) \exp(-\beta E_i)}{Z} \qquad (24)$$
Where Z is the partition function $E_i$ is the energy of the macrostate $i$, and $d(m_{t_i}, s_{b_i}, s_{s_i}, s_{bs_i}, m_{s_i})$ is the density of this macrostate. Here, Z is given by:

$$Z = \sum_{i=1}^{N_L} d(m_{t_i}, s_{b_i}, s_{s_i}, s_{bs_i}, m_{s_i}) \exp(-\beta E_i) \quad (25)$$

and

$$E_i = -hm_{t_i} - J_{bb}s_{b_i} - J_{ss}s_{s_i} - J_{bs}s_{bs_i} - Lm_{s_i} \quad (26)$$

In Equations (24) and (25), $N_L$ is the number of different configurations with the same five values $m_t, s_b, s_s, s_{bs},$ and $m_s$. The density of states $d(m_t, s_b, s_s, s_{bs}, m_s)$ is calculated by entropic sampling [48–50], and Equation (24) is solved by numerical techniques such as bisection. From the thermal average values $<\sigma>$, the HS fraction Nhs is calculated using the relation (6).

Let us recall the basic principles of MCES [50–52]. It consists in introducing the detailed balance equation of the Monte Carlo (MC) procedure, as shown in:

$$P_i W(i \to j) = P_j W(j \to i) \quad (27)$$

This suited biased distribution P is designed to favor configurations belonging to weakly degenerated macrostates while dampening those within the highly degenerated macrostates. It can be expressed as:

$$P_i = \frac{1}{d(m_{t_i}, s_{b_i}, s_{s_i}, s_{bs_i}, m_{s_i})} \quad (28)$$

In this case, the balance equation can be written as:

$$\frac{W(i \to j)}{W(j \to i)} = \frac{P_j}{P_i} = \frac{d(m_{t_i}, s_{b_i}, s_{s_i}, s_{bs_i}, m_{s_i})}{d(m_{t_j}, s_{b_j}, s_{s_j}, s_{bs_j}, m_{s_j})} \quad (29)$$

Using $d(m_{t_i}, s_{b_i}, s_{s_i}, s_{bs_i}, m_{s_i})$ as a bias, a MC sampling is run; it is termed a “Monte Carlo step” and yields a histogram of the frequency of the macrostates $H(m_{t_i}, s_{b_i}, s_{s_i}, s_{bs_i}, m_{s_i})$, also written as $H_i(m_{t_i}, s_{b_i}, s_{s_i}, s_{bs_i}, m_{s_i})$.

By construction,
Where Z is the partition function $E_i$ is the energy of the macrostate $i$, and $d(m_{t_i}, s_{b_i}, s_{s_i}, s_{bs_i}, m_{s_i})$ is the density of this macrostate. Here, Z is given by:

$$Z = \sum_{i=1}^{N_L} d(m_{t_i}, s_{b_i}, s_{s_i}, s_{bs_i}, m_{s_i}) \exp(-\beta E_i) \quad (25)$$

and

$$E_i = -hm_{t_i} - J_{bb}s_{b_i} - J_{ss}s_{s_i} - J_{bs}s_{bs_i} - Lm_{s_i} \quad (26)$$

In Equations (24) and (25), $N_L$ is the number of different configurations with the same five values $m_t, s_b, s_s, s_{bs},$ and $m_s$. The density of states $d(m_t, s_b, s_s, s_{bs}, m_s)$ is calculated by entropic sampling [48–50], and Equation (24) is solved by numerical techniques such as bisection. From the thermal average values $<\sigma>$, the HS fraction $Nhs$ is calculated using the relation (6).

Let us recall the basic principles of MCES [50–52]. It consists in introducing the detailed balance equation of the Monte Carlo (MC) procedure, as shown in:

$$P_i W(i \to j) = P_j W(j \to i) \quad (27)$$

This suited biased distribution $P$ is designed to favor configurations belonging to weakly degenerated macrostates while dampening those within the highly degenerated macrostates. It can be expressed as:

$$P_i = \frac{1}{d(m_{t_i}, s_{b_i}, s_{s_i}, s_{bs_i}, m_{s_i})} \quad (28)$$

In this case, the balance equation can be written as:

$$\frac{W(i \to j)}{W(j \to i)} = \frac{P_j}{P_i} = \frac{d(m_{t_i}, s_{b_i}, s_{s_i}, s_{bs_i}, m_{s_i})}{d(m_{t_j}, s_{b_j}, s_{s_j}, s_{bs_j}, m_{s_j})} \quad (29)$$

Using $d(m_{t_i}, s_{b_i}, s_{s_i}, s_{bs_i}, m_{s_i})$ as a bias, a MC sampling is run; it is termed a “Monte Carlo step” and yields a histogram of the frequency of the macrostates $H(m_{t_i}, s_{b_i}, s_{s_i}, s_{bs_i}, m_{s_i})$, also written as $H_i(m_{t_i}, s_{b_i}, s_{s_i}, s_{bs_i}, m_{s_i})$.

By construction,
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Schematic representation of a 2D 6 × 6 square-shaped lattice: filled blue circles represent bulk sites (N_b); filled green and red circles represent edge (N_e) and corner (N_c) sites, respectively, interacting with their immediate environment (matrix).</figcaption></figure>

$$H_i(m_t, s_b, s_s, s_{bs}, m_s) \propto d_{i+1}(m_t, s_b, s_s, s_{bs}, m_s) \frac{1}{d_i(m_t, s_b, s_s, s_{bs}, m_s)} \quad (30)$$

The resulting restricted density of states is calculated after the correction for the bias was applied:

$$d_{i+1}(m_t, s_b, s_s, s_{bs}, m_s) \propto d_i(m_t, s_b, s_s, s_{bs}, m_s) \times H_i(m_t, s_b, s_s, s_{bs}, m_s) \quad (31)$$

The flat character of the histogram $H(m_t, s_b, s_s, s_{bs}, m_s)$ has a convenient convergence criterion to obtain a very good estimation of the density of the states. In our case, seven iterations of the 10⁶ Carlo Steps were necessary. After the table of the values of $(m_t, s_b, s_s, s_{bs}, m_s)$ and $d(m_t, s_b, s_s, s_{bs}, m_s)$ is obtained, the quasi-exact partition function can be calculated using the expression (25).

## LOCAL MEAN-FIELD APPROACH

In the present work, three types of sites related to three lattice regions are considered: the atoms located in the bulk (N_b), which are surrounded by four first-neighbors and those atoms located on the surface and, more precisely, on the edge (N_e) and on the corner (N_c), which interact with three and two first-neighbors, respectively. The molecules located at the surface (edge and corner) interact with an external environment (matrix effect), which gives them specific properties.

The case of a system comprising 36 molecules is represented in Figure 3.

Note that in the framework of a classical mean-field approximation [53] and the case of a 2D square-shaped lattice, the coordination number is $q = 4$ [Figure 4].

The consideration of three types of sites leads us to distinguish three local average order parameters $<\sigma_\alpha>$ and three coordination numbers $q_\alpha$. In the LMFA, which amounts to taking into account a local Hamiltonian for each type of site $\alpha$ with $\alpha = b, e, c$, corresponding to bulk, edge, and corner, respectively.

For each region, $\alpha$, we apply the LMFA. In the case of regions connecting each other only via the long-range interaction, each molecule $\sigma_i$ is surrounded by $q_\alpha < \sigma_\alpha >$, and the local Hamiltonian, given by Equation (14), can be expressed as Equation (32) below:
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Schematic representation of a 2D 6 × 6 square-shaped lattice: filled blue circles represent bulk sites (N_b); filled green and red circles represent edge (N_e) and corner (N_c) sites, respectively, interacting with their immediate environment (matrix).</figcaption></figure>

$$H_i(m_t, s_b, s_s, s_{bs}, m_s) \propto d_{i+1}(m_t, s_b, s_s, s_{bs}, m_s) \frac{1}{d_i(m_t, s_b, s_s, s_{bs}, m_s)} \quad (30)$$

The resulting restricted density of states is calculated after the correction for the bias was applied:

$$d_{i+1}(m_t, s_b, s_s, s_{bs}, m_s) \propto d_i(m_t, s_b, s_s, s_{bs}, m_s) \times H_i(m_t, s_b, s_s, s_{bs}, m_s) \quad (31)$$

The flat character of the histogram $H(m_t, s_b, s_s, s_{bs}, m_s)$ has a convenient convergence criterion to obtain a very good estimation of the density of the states. In our case, seven iterations of the 10⁶ Carlo Steps were necessary. After the table of the values of $(m_t, s_b, s_s, s_{bs}, m_s)$ and $d(m_t, s_b, s_s, s_{bs}, m_s)$ is obtained, the quasi-exact partition function can be calculated using the expression (25).

## LOCAL MEAN-FIELD APPROACH

In the present work, three types of sites related to three lattice regions are considered: the atoms located in the bulk (N_b), which are surrounded by four first-neighbors and those atoms located on the surface and, more precisely, on the edge (N_e) and on the corner (N_c), which interact with three and two first-neighbors, respectively. The molecules located at the surface (edge and corner) interact with an external environment (matrix effect), which gives them specific properties.

The case of a system comprising 36 molecules is represented in Figure 3.

Note that in the framework of a classical mean-field approximation [53] and the case of a 2D square-shaped lattice, the coordination number is $q = 4$ [Figure 4].

The consideration of three types of sites leads us to distinguish three local average order parameters $<\sigma_\alpha>$ and three coordination numbers $q_\alpha$. In the LMFA, which amounts to taking into account a local Hamiltonian for each type of site $\alpha$ with $\alpha = b, e, c$, corresponding to bulk, edge, and corner, respectively.

For each region, $\alpha$, we apply the LMFA. In the case of regions connecting each other only via the long-range interaction, each molecule $\sigma_i$ is surrounded by $q_\alpha < \sigma_\alpha >$, and the local Hamiltonian, given by Equation (14), can be expressed as Equation (32) below:
<figure><img src="image_9.png" /><figcaption>**Figure 4.** Mean-field approximation in a 2D square-shaped lattice. The $\sigma_i$ spin is surrounded by four nearest neighbors with the average magnetization $\langle \sigma \rangle$.</figcaption></figure>

<table><caption>**Table 1.** The ligand field contributions and the number of short-range interactions as a function of various localizations in the square-shaped lattice</caption><thead><tr><th>Site</th><th>Bulk</th><th>Edge</th><th>Corner</th></tr></thead><tbody><tr><td>q<sub>α</sub></td><td>4</td><td>3</td><td>2</td></tr><tr><td>z<sub>α</sub></td><td>0</td><td>1</td><td>1</td></tr><tr><td>Ligand-field</td><td>(Δ - k<sub>B</sub>T *ln(g)*) / 2</td><td>(Δ - k<sub>B</sub>T *ln(g)* - 2L) / 2</td><td>(Δ - k<sub>B</sub>T *ln(g)* - 2L) / 2</td></tr></tbody></table>

$$H_{\alpha} = \sum_{i=1}^{N_{\alpha}} \left( \frac{\Delta - k_{B}T \ln(g)}{2} - J_{\alpha}q_{\alpha} \langle \sigma_{\alpha} \rangle - G \langle \sigma \rangle - Lz_{\alpha} \right) \sigma_{i} = - \sum_{i=1}^{N_{\alpha}} h_{\alpha} \sigma_{i} \quad (32)$$

with

$$h_{\alpha} = - \frac{\Delta - k_{B}T \ln(g) - 2J_{\alpha}q_{\alpha} \langle \sigma_{\alpha} \rangle - 2G \langle \sigma \rangle - 2Lz_{\alpha}}{2} \quad (33)$$

where $q_α$ is the number of interactions between a molecule and its first-neighbors, $z_α$ is the number of interactions $L$ between a molecule and the external environment, and $J_α$ is the interaction between the nearest neighbor (nn) molecules. For molecules located in the edge and in the corner, the average value of the number of interactions with the external environment is set equal to 1. Table 1 summarizes the characteristics of each site.

Here, for simplicity of calculations, we have taken $z = 1$ (the same coordination number) everywhere in the surface (edge and corner), meaning that we are considering the same connectivity for all molecules at the surface.

The short-range interaction term $J_α$ depends on the location of the molecules. Two interaction constants are considered in our calculations:

$J_{\text{bulk}}^{\text{LMFA}}$ for the molecules located in the bulk; $J_{\text{surf}}^{\text{LMFA}}$ for the molecules located at the surface.

These two interaction terms, $J_{\text{bulk}}^{\text{LMFA}}$ and $J_{\text{surf}}^{\text{LMFA}}$, are related to the interaction terms used in the MCES method by weighting the various interactions of each molecule in each region. Let us consider the case of the bulk region. We denote by $N_{bb}$, $N_{bs}$, the respective bulk-bulk and the bulk-surface bonds, leading to a total number of bonds including at least one bulk site ($N_{bb} + N_{bs}$).
<figure><img src="image_9.png" /><figcaption>**Figure 4.** Mean-field approximation in a 2D square-shaped lattice. The $\sigma_i$ spin is surrounded by four nearest neighbors with the average magnetization $\langle \sigma \rangle$.</figcaption></figure>

<table><caption>**Table 1.** The ligand field contributions and the number of short-range interactions as a function of various localizations in the square-shaped lattice</caption><thead><tr><th>Site</th><th>Bulk</th><th>Edge</th><th>Corner</th></tr></thead><tbody><tr><td>q<sub>α</sub></td><td>4</td><td>3</td><td>2</td></tr><tr><td>z<sub>α</sub></td><td>0</td><td>1</td><td>1</td></tr><tr><td>Ligand-field</td><td>(Δ - k<sub>B</sub>T *ln(g)*) / 2</td><td>(Δ - k<sub>B</sub>T *ln(g)* - 2L) / 2</td><td>(Δ - k<sub>B</sub>T *ln(g)* - 2L) / 2</td></tr></tbody></table>

$$H_{\alpha} = \sum_{i=1}^{N_{\alpha}} \left( \frac{\Delta - k_{B}T \ln(g)}{2} - J_{\alpha}q_{\alpha} \langle \sigma_{\alpha} \rangle - G \langle \sigma \rangle - Lz_{\alpha} \right) \sigma_{i} = - \sum_{i=1}^{N_{\alpha}} h_{\alpha} \sigma_{i} \quad (32)$$

with

$$h_{\alpha} = - \frac{\Delta - k_{B}T \ln(g) - 2J_{\alpha}q_{\alpha} \langle \sigma_{\alpha} \rangle - 2G \langle \sigma \rangle - 2Lz_{\alpha}}{2} \quad (33)$$

where $q_\alpha$ is the number of interactions between a molecule and its first-neighbors, $z_\alpha$ is the number of interactions $L$ between a molecule and the external environment, and $J_\alpha$ is the interaction between the nearest neighbor (nn) molecules. For molecules located in the edge and in the corner, the average value of the number of interactions with the external environment is set equal to 1. Table 1 summarizes the characteristics of each site.

Here, for simplicity of calculations, we have taken $z = 1$ (the same coordination number) everywhere in the surface (edge and corner), meaning that we are considering the same connectivity for all molecules at the surface.

The short-range interaction term $J_\alpha$ depends on the location of the molecules. Two interaction constants are considered in our calculations:

$J_{\text{bulk}}^{\text{LMFA}}$ for the molecules located in the bulk; $J_{\text{surf}}^{\text{LMFA}}$ for the molecules located at the surface.

These two interaction terms, $J_{\text{bulk}}^{\text{LMFA}}$ and $J_{\text{surf}}^{\text{LMFA}}$, are related to the interaction terms used in the MCES method by weighting the various interactions of each molecule in each region. Let us consider the case of the bulk region. We denote by $N_{bb}$, $N_{bs}$, the respective bulk-bulk and the bulk-surface bonds, leading to a total number of bonds including at least one bulk site ($N_{bb} + N_{bs}$).
Let $J_{bb}^{MCES}$ and $J_{bs}^{MCES}$ be the bulk-bulk and bulk-surface short-range interactions used in MCES method. The total interaction energy term in the bulk region is given by $(N_{bb} \times J_{bb}^{MCES}) + (N_{bs} \times J_{bs}^{MCES})$. Therefore, the bulk value $J_{bulk}^{LMFA}$ used in Equation (32) is obtained by averaging over $J_{bb}^{MCES}$, $J_{bs}^{MCES}$ weighted by the numbers of their corresponding bonds, as follows:

$$J_{bulk}^{LMFA} = \frac{(N_{bb} \times J_{bb}^{MCES}) + (N_{bs} \times J_{bs}^{MCES})}{(N_{bb} + N_{bs})} \qquad (34)$$

In a similar way, for the surface sites, Equation (35) is written as:

$$J_{surf}^{LMFA} = \frac{(N_{ss} \times J_{ss}^{MCES}) + (N_{bs} \times J_{bs}^{MCES})}{(N_{ss} + N_{bs})} \qquad (35)$$

where $N_{ss}$ and $J_{ss}^{MCES}$ are the number surface-surface bonds and the surface-surface short-range interaction term used in the MCES method, respectively.

The average spin state can be written as follows:

$$<\sigma_{\alpha}> = \tanh(\beta h_{\alpha}) = \tanh\beta\left(-\frac{\Delta - k_{B}T \ln(g) - 2J_{\alpha}q_{\alpha} <\sigma_{\alpha}> - 2G <\sigma> - 2Lz_{\alpha}}{2}\right), \alpha = b, e, c \qquad (36)$$

So, for each region b (bulk), e (edge), and c (corner), the average spin state is given by:

$$<\sigma_b> = \tanh\left(-\frac{\Delta - k_B T \ln(g) - 2 \times J_{bulk}^{LMFA} \times 4 <\sigma_b> - 2G <\sigma>}{2k_B T}\right) \qquad (37)$$

$$<\sigma_e> = \tanh\left(-\frac{\Delta - k_B T \ln(g) - 2 \times J_{surf}^{LMFA} \times 3 <\sigma_e> - 2 \times 1 \times L - 2G <\sigma>}{2k_B T}\right) \qquad (38)$$

$$<\sigma_c> = \tanh\left(-\frac{\Delta - k_B T \ln(g) - 2 \times J_{surf}^{LMFA} \times 2 <\sigma_c> - 2 \times 1 \times L - 2G <\sigma>}{2k_B T}\right) \qquad (39)$$

which is equivalent to:
Let $J_{bb}^{MCES}$ and $J_{bs}^{MCES}$ be the bulk-bulk and bulk-surface short-range interactions used in MCES method. The total interaction energy term in the bulk region is given by $(N_{bb} \times J_{bb}^{MCES}) + (N_{bs} \times J_{bs}^{MCES})$. Therefore, the bulk value $J_{bulk}^{LMFA}$ used in Equation (32) is obtained by averaging over $J_{bb}^{MCES}$, $J_{bs}^{MCES}$ weighted by the numbers of their corresponding bonds, as follows:

$$J_{bulk}^{LMFA} = \frac{(N_{bb} \times J_{bb}^{MCES}) + (N_{bs} \times J_{bs}^{MCES})}{(N_{bb} + N_{bs})} \quad (34)$$

In a similar way, for the surface sites, Equation (35) is written as:

$$J_{surf}^{LMFA} = \frac{(N_{ss} \times J_{ss}^{MCES}) + (N_{bs} \times J_{bs}^{MCES})}{(N_{ss} + N_{bs})} \quad (35)$$

where $N_{ss}$ and $J_{ss}^{MCES}$ are the number surface-surface bonds and the surface-surface short-range interaction term used in the MCES method, respectively.

The average spin state can be written as follows:

$$<\sigma_{\alpha}> = \tanh(\beta h_{\alpha}) = \tanh \beta \left( -\frac{\Delta - k_{B}T \ln(g) - 2J_{\alpha}q_{\alpha} <\sigma_{\alpha}> - 2G <\sigma> - 2Lz_{\alpha}}{2} \right), \quad \alpha = b, e, c \quad (36)$$

So, for each region b (bulk), e (edge), and c (corner), the average spin state is given by:

$$<\sigma_b> = \tanh \left( -\frac{\Delta - k_B T \ln(g) - 2 \times J_{bulk}^{LMFA} \times 4 <\sigma_b> - 2G <\sigma>}{2k_B T} \right) \quad (37)$$

$$<\sigma_e> = \tanh \left( -\frac{\Delta - k_B T \ln(g) - 2 \times J_{surf}^{LMFA} \times 3 <\sigma_e> - 2 \times 1 \times L - 2G <\sigma>}{2k_B T} \right) \quad (38)$$

$$<\sigma_c> = \tanh \left( -\frac{\Delta - k_B T \ln(g) - 2 \times J_{surf}^{LMFA} \times 2 <\sigma_c> - 2 \times 1 \times L - 2G <\sigma>}{2k_B T} \right) \quad (39)$$

which is equivalent to:
$$ \langle \sigma_b \rangle - \tanh \left( - \frac{\Delta - k_B T \ln(g) - 2 \times J_{\text{bulk}}^{\text{LMFA}} \times 4 \langle \sigma_b \rangle - 2G \langle \sigma \rangle}{2k_B T} \right) = \xi_1 (\langle \sigma_c \rangle, \langle \sigma_e \rangle, \langle \sigma_b \rangle) = 0 \quad (40) $$

$$ \langle \sigma_e \rangle - \tanh \left( - \frac{\Delta - k_B T \ln(g) - 2 \times J_{\text{surf}}^{\text{LMFA}} \times 3 \langle \sigma_e \rangle - 2 \times 1 \times L - 2G \langle \sigma \rangle}{2k_B T} \right) = \xi_2 (\langle \sigma_c \rangle, \langle \sigma_e \rangle, \langle \sigma_b \rangle) = 0 \quad (41) $$

$$ \langle \sigma_c \rangle - \tanh \left( - \frac{\Delta - k_B T \ln(g) - 2 \times J_{\text{surf}}^{\text{LMFA}} \times 2 \langle \sigma_c \rangle - 2 \times 1 \times L - 2G \langle \sigma \rangle}{2k_B T} \right) = \xi_3 (\langle \sigma_c \rangle, \langle \sigma_e \rangle, \langle \sigma_b \rangle) = 0 \quad (42) $$

The weighted average of the three local order parameters associated with the three regions, bulk, edge, and corner, leads to the average fictitious spin of the system expressed as:

$$ \langle \sigma \rangle = \frac{N_b \langle \sigma_b \rangle + N_e \langle \sigma_e \rangle + N_c \langle \sigma_c \rangle}{N_T} = \frac{1}{N_T} \sum_{\alpha} N_{\alpha} \langle \sigma_{\alpha} \rangle, \quad \alpha = b, e, c \quad (43) $$

and the high-spin fraction, Nhs, which is the probability that the HS state is occupied and which is written as Equation (6): $Nhs = \frac{1+\langle \sigma \rangle}{2}$.

The purpose is to solve a system of three equations $\xi_1$, $\xi_2$, and $\xi_3$, whose three order parameters $\langle \sigma_b \rangle$, $\langle \sigma_e \rangle$, and $\langle \sigma_c \rangle$, connected by Equation (43), are simultaneous solutions. The Newton method, which is appropriate for rapid convergences, has been used. Calculations are based on three points $x_1$, $x_2$, and $x_3$ close to the solution and which represent $\langle \sigma_b \rangle$, $\langle \sigma_e \rangle$, and $\langle \sigma_c \rangle$, respectively.

The calculation of the three points “close to the solution” is performed in two steps: Let X be the space in which values of ($x_1$, $x_2$, $x_3$) are scanned. Let $\xi$ be the space from which ($\xi_1$, $\xi_1$, $\xi_1$) are calculated as a function of ($x_1$, $x_2$, $x_3$). Then, the first step consists of mapping space X on space $\xi$. Those values of ($x_1$, $x_2$, $x_3$) that map onto ($\xi_1$, $\xi_2$, $\xi_3$) close to zero within $10_{-6}$, are then selected to perform the second step: application of the Newton-Raphson technique.

Solving the system amounts to calculating $h_1$, $h_2$, and $h_3$ such that:

$$ \xi_i(x_1 + h_1, x_2 + h_2, x_3 + h_3) = 0, \quad i = 1, 2, 3 \quad (44) $$

which allows to write in the vicinity of points $x_1$, $x_2$, and $x_3$:
$$<\sigma_b> - \tanh\left(-\frac{\Delta - k_B T \ln(g) - 2 \times J_{\text{bulk}}^{\text{LMFA}} \times 4 <\sigma_b> - 2G <\sigma>}{2k_B T}\right) = \xi_1 (<\sigma_c>, <\sigma_e>, <\sigma_b>) = 0 \qquad (40)$$

$$<\sigma_e> - \tanh\left(-\frac{\Delta - k_B T \ln(g) - 2 \times J_{\text{surf}}^{\text{LMFA}} \times 3 <\sigma_e> - 2 \times 1 \times L - 2G <\sigma>}{2k_B T}\right) = \xi_2 (<\sigma_c>, <\sigma_e>, <\sigma_b>) = 0 \qquad (41)$$

$$<\sigma_c> - \tanh\left(-\frac{\Delta - k_B T \ln(g) - 2 \times J_{\text{surf}}^{\text{LMFA}} \times 2 <\sigma_c> - 2 \times 1 \times L - 2G <\sigma>}{2k_B T}\right) = \xi_3 (<\sigma_c>, <\sigma_e>, <\sigma_b>) = 0 \qquad (42)$$

The weighted average of the three local order parameters associated with the three regions, bulk, edge, and corner, leads to the average fictitious spin of the system expressed as:

$$<\sigma> = \frac{N_b <\sigma_b> + N_e <\sigma_e> + N_c <\sigma_c>}{N_T} = \frac{1}{N_T} \sum_{\alpha} N_{\alpha} <\sigma_{\alpha)}, \alpha = b, e, c \qquad (43)$$

and the high-spin fraction, Nhs, which is the probability that the HS state is occupied and which is written as Equation (6): $Nhs = \frac{1+<\sigma>}{2}$.

The purpose is to solve a system of three equations $\xi_1$, $\xi_2$, and $\xi_3$, whose three order parameters $<\sigma_b>$, $<\sigma_e>$, and $<\sigma_c>$, connected by Equation (43), are simultaneous solutions. The Newton method, which is appropriate for rapid convergences, has been used. Calculations are based on three points $x_1$, $x_2$, and $x_3$ close to the solution and which represent $<\sigma_b>$, $<\sigma_e>$, and $<\sigma_c>$, respectively.

The calculation of the three points “close to the solution” is performed in two steps: Let X be the space in which values of ($x_1$, $x_2$, $x_3$) are scanned. Let $\xi$ be the space from which ($\xi_1$, $\xi_1$, $\xi_1$) are calculated as a function of ($x_1$, $x_2$, $x_3$). Then, the first step consists of mapping space X on space $\xi$. Those values of ($x_1$, $x_2$, $x_3$) that map onto ($\xi_1$, $\xi_2$, $\xi_3$) close to zero within $10_{-6}$, are then selected to perform the second step: application of the Newton-Raphson technique.

Solving the system amounts to calculating $h_1$, $h_2$, and $h_3$ such that:

$$\xi_i(x_1 + h_1, x_2 + h_2, x_3 + h_3) = 0, i = 1, 2, 3 \qquad (44)$$

which allows to write in the vicinity of points $x_1$, $x_2$, and $x_3$:
$$ \xi_i (x_1 + h_1, x_2 + h_2, x_3 + h_3) \cong \xi_i (x_1, x_2, x_3) + h_1 \frac{\partial \xi_i}{\partial x_1} + h_2 \frac{\partial \xi_i}{\partial x_2} + h_3 \frac{\partial \xi_i}{\partial x_3} = 0, i = 1, 2, 3 \quad (45) $$

The following system of three equations with three unknowns is thus obtained as:

$$ h_1 \frac{\partial \xi_1}{\partial x_1} + h_2 \frac{\partial \xi_1}{\partial x_2} + h_3 \frac{\partial \xi_1}{\partial x_3} = -\xi_1 (x_1, x_2, x_3) \quad (46) $$

$$ h_1 \frac{\partial \xi_2}{\partial x_1} + h_2 \frac{\partial \xi_2}{\partial x_2} + h_3 \frac{\partial \xi_2}{\partial x_3} = -\xi_2 (x_1, x_2, x_3) \quad (47) $$

$$ h_1 \frac{\partial \xi_3}{\partial x_1} + h_2 \frac{\partial \xi_3}{\partial x_2} + h_3 \frac{\partial \xi_3}{\partial x_3} = -\xi_3 (x_1, x_2, x_3) \quad (48) $$

And where $h_1$, $h_2$, and $h_3$ are calculated as follows:

$$ h_1 = \left| \begin{array}{ccc} -\xi_1 & \frac{\partial \xi_1}{\partial x_2} & \frac{\partial \xi_1}{\partial x_3} \\ -\xi_2 & \frac{\partial \xi_2}{\partial x_2} & \frac{\partial \xi_2}{\partial x_3} \\ -\xi_3 & \frac{\partial \xi_3}{\partial x_2} & \frac{\partial \xi_3}{\partial x_3} \end{array} \right|, h_2 = \left| \begin{array}{ccc} \frac{\partial \xi_1}{\partial x_1} & -\xi & \frac{\partial \xi_1}{\partial x_3} \\ \frac{\partial \xi_2}{\partial x_1} & -\xi_2 & \frac{\partial \xi_2}{\partial x_3} \\ \frac{\partial \xi_3}{\partial x_1} & -\xi_3 & \frac{\partial \xi_3}{\partial x_3} \end{array} \right|, h_3 = \left| \begin{array}{ccc} \frac{\partial \xi_1}{\partial x_1} & \frac{\partial \xi_1}{\partial x_2} & -\xi_1 \\ \frac{\partial \xi_2}{\partial x_1} & \frac{\partial \xi_2}{\partial x_2} & -\xi_2 \\ \frac{\partial \xi_3}{\partial x_1} & \frac{\partial \xi_3}{\partial x_2} & -\xi_3 \end{array} \right| \quad (49) $$

In a general way, the $\xi_i$ functions with $i = 1$ (bulk), 2 (corner), and 3 (edge) and their derivatives can be expressed as follows:

$$ \xi_i = x_i - \tanh \left( - \frac{\Delta - T \ln(g) - 2 * q_i J_i x_i - 2G < \sigma > - 2 * z_i L}{2k_B T} \right) \quad (50) $$

$$ \frac{\partial \xi_i}{\partial x_i} = 1 - \left[ 1 - \tanh^2 \left( - \frac{\Delta - T \ln(g) - 2 * q_i J_i x_i - 2G < \sigma > - 2 * z_i L}{2k_B T} \right) \right] \left( \frac{q_i J_i}{k_B T} + \frac{G N_i}{N_T k_B T} \right) \quad (51) $$

$$ \frac{\partial \xi_i}{\partial x_j} = - \left[ 1 - \tanh^2 \left( - \frac{\Delta - T \ln(g) - 2 * q_i J_i x_i - 2G < \sigma > - 2 * z_i L}{2k_B T} \right) \right] \left( \frac{G N_i}{N_T k_B T} \right) \quad (52) $$
$$ \xi_i (x_1 + h_1, x_2 + h_2, x_3 + h_3) \cong \xi_i (x_1, x_2, x_3) + h_1 \frac{\partial \xi_i}{\partial x_1} + h_2 \frac{\partial \xi_i}{\partial x_2} + h_3 \frac{\partial \xi_i}{\partial x_3} = 0, i = 1, 2, 3 \quad (45) $$

The following system of three equations with three unknowns is thus obtained as:

$$ h_1 \frac{\partial \xi_1}{\partial x_1} + h_2 \frac{\partial \xi_1}{\partial x_2} + h_3 \frac{\partial \xi_1}{\partial x_3} = -\xi_1 (x_1, x_2, x_3) \quad (46) $$

$$ h_1 \frac{\partial \xi_2}{\partial x_1} + h_2 \frac{\partial \xi_2}{\partial x_2} + h_3 \frac{\partial \xi_2}{\partial x_3} = -\xi_2 (x_1, x_2, x_3) \quad (47) $$

$$ h_1 \frac{\partial \xi_3}{\partial x_1} + h_2 \frac{\partial \xi_3}{\partial x_2} + h_3 \frac{\partial \xi_3}{\partial x_3} = -\xi_3 (x_1, x_2, x_3) \quad (48) $$

And where $h_1$, $h_2$, and $h_3$ are calculated as follows:

$$ h_1 = \left| \begin{array}{ccc} -\xi_1 & \frac{\partial \xi_1}{\partial x_2} & \frac{\partial \xi_1}{\partial x_3} \\ -\xi_2 & \frac{\partial \xi_2}{\partial x_2} & \frac{\partial \xi_2}{\partial x_3} \\ -\xi_3 & \frac{\partial \xi_3}{\partial x_2} & \frac{\partial \xi_3}{\partial x_3} \end{array} \right|, h_2 = \left| \begin{array}{ccc} \frac{\partial \xi_1}{\partial x_1} & -\xi & \frac{\partial \xi_1}{\partial x_3} \\ \frac{\partial \xi_2}{\partial x_1} & -\xi_2 & \frac{\partial \xi_2}{\partial x_3} \\ \frac{\partial \xi_3}{\partial x_1} & -\xi_3 & \frac{\partial \xi_3}{\partial x_3} \end{array} \right|, h_3 = \left| \begin{array}{ccc} \frac{\partial \xi_1}{\partial x_1} & \frac{\partial \xi_1}{\partial x_2} & -\xi_1 \\ \frac{\partial \xi_2}{\partial x_1} & \frac{\partial \xi_2}{\partial x_2} & -\xi_2 \\ \frac{\partial \xi_3}{\partial x_1} & \frac{\partial \xi_3}{\partial x_2} & -\xi_3 \end{array} \right| \quad (49) $$

In a general way, the $\xi_i$ functions with i = 1 (bulk), 2 (corner), and 3 (edge) and their derivatives can be expressed as follows:

$$ \xi_i = x_i - \tanh \left( - \frac{\Delta - T \ln(g) - 2 * q_i J_i x_i - 2G < \sigma > - 2 * z_i L}{2k_B T} \right) \quad (50) $$

$$ \frac{\partial \xi_i}{\partial x_i} = 1 - \left[ 1 - \tanh^2 \left( - \frac{\Delta - T \ln(g) - 2 * q_i J_i x_i - 2G < \sigma > - 2 * z_i L}{2k_B T} \right) \right] \left( \frac{q_i J_i}{k_B T} + \frac{G N_i}{N_T k_B T} \right) \quad (51) $$

$$ \frac{\partial \xi_i}{\partial x_j} = - \left[ 1 - \tanh^2 \left( - \frac{\Delta - T \ln(g) - 2 * q_i J_i x_i - 2G < \sigma > - 2 * z_i L}{2k_B T} \right) \right] \left( \frac{G N_i}{N_T k_B T} \right) \quad (52) $$
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Schematic structure of the SCO compound iron (II) [Fe(btr)₂(NCS)₂], comprising a bis-triazole (btr) whose molecular formula are C₄N₆H₄ and a cyano group NCS.</figcaption></figure>

<table><caption>**Table 2.** Correspondence between the values of short-range interactions *J*/*k*<sub>B</sub> used by LMFA and by MCES method in the case of a 6 × 6 square lattice when *J*<sub>bs</sub><sup>MCES</sup>/*k*<sub>B</sub> = 0 K</caption><thead><tr><th>J<sup>MCES</sup>/k<sub>B</sub> [K]</th><th>J<sup>MCES</sup>/k<sub>B</sub> [K]</th><th>J<sup>MCES</sup>/k<sub>B</sub> [K]</th><th>J<sup>LMFA</sup>/k<sub>B</sub> [K]</th><th>J<sup>LMFA</sup>/k<sub>B</sub> [K]</th></tr></thead><tbody><tr><td>60</td><td>60</td><td>0</td><td>36</td><td>33.33</td></tr></tbody></table>

With the *h*<sub>i</sub> obtained in this way, the solutions (*x*<sub>i</sub>)<sub>sol</sub> are deduced by the relation:

$$ (x_i)_{sol} = x_i + h_i \qquad (53) $$

## RESULTS AND DISCUSSION

We are interested in the effect of the coupling between surface and bulk molecules, and for that, we perform numerical simulations by MCES and by LMFA. To proceed, we choose a set of realistic SCO thermodynamic parameters, and for this purpose, we chose the data derived from experimental data of one of these SCO compounds: [Fe(btr)₂(NCS)₂], btr = 4,4'-bis-1,2,4-triazole [54], whose molecular structure is shown in Figure 5. The energy gap Δ and ln(g) are derived from the enthalpy change (ΔH ≈ 11 kJ/mol) and the entropy change (ΔS ≈ 50 J/mol/K), respectively: Δ = ΔH/R ≈ 1,300 K), ln(g) = ΔS/R ≈ 6.01, and R is the perfect gas constant. The equilibrium temperature of the system T<sub>eq</sub> is deduced by T<sub>eq</sub> = Δ/k<sub>B</sub>/ln(g), which leads to ≈ 216.3 K. In the first step, in the case of a 6 × 6 configuration, the values of the short-range interaction parameters J<sub>bb</sub><sup>MCES</sup> and J<sub>ss</sub><sup>MCES</sup> are set to 60 K, and the intensity of the short-range interaction parameter J<sub>bs</sub><sup>MCES</sup> between bulk and surface is gradually increased from 0 to 60 K. In the second step, the size effect is studied with each method. It should be added that the dimensions of the square lattice used in these simulations do not correspond, for the moment, to synthesized nanoparticles.

### The case J<sub>bs</sub><sup>MCES</sup> = 0

Figure 6 displays the results obtained in the case of a 6 × 6 square-shaped lattice, and the values of the interaction parameters *J* used by MCES and by LMFA, calculated from Equations (34) and (35), are gathered in Table 2. The curves highlight a two-step behavior in the form of two well-separated hysteresis loops. The hysteresis width is defined as ΔT = *T*<sub>up</sub> − *T*<sub>down</sub>, where *T*<sub>up</sub> is the ascending thermal transition temperature and *T*<sub>down</sub> is the descending thermal transition temperature. The average temperature between *T*<sub>up</sub> and *T*<sub>down</sub> is the equilibrium temperature *T*<sub>eq</sub>, which corresponds to a HS fraction equal to 1/2(Nhs = 1/2).

Figure 6 illustrates a qualitative and quantitative agreement between the calculations obtained by MCES and LMFA. The hysteresis loops, observable for each method, at lower temperatures, ~ 144 K, are related to the
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Schematic structure of the SCO compound iron (II) [Fe(btr)₂(NCS)₂], comprising a bis-triazole (btr) whose molecular formula are C₄N₆H₄ and a cyano group NCS.</figcaption></figure>

<table><caption>**Table 2.** Correspondence between the values of short-range interactions *J*/*k*<sub>B</sub> used by LMFA and by MCES method in the case of a 6 × 6 square lattice when *J*<sub>bs</sub><sup>MCES</sup>/*k*<sub>B</sub> = 0 K</caption><thead><tr><th>J<sup>MCES</sup>/<i>k</i><sub>B</sub> [K]</th><th>J<sup>MCES</sup>/<i>k</i><sub>B</sub> [K]</th><th>J<sup>MCES</sup>/<i>k</i><sub>B</sub> [K]</th><th>J<sup>LMFA</sup>/<i>k</i><sub>B</sub> [K]</th><th>J<sup>LMFA</sup>/<i>k</i><sub>B</sub> [K]</th></tr></thead><tbody><tr><td>60</td><td>60</td><td>0</td><td>36</td><td>33.33</td></tr></tbody></table>

With the *h*<sub>i</sub> obtained in this way, the solutions (*x*<sub>i</sub>)<sub>sol</sub> are deduced by the relation:

$$ (x_i)_{sol} = x_i + h_i \qquad (53) $$

## RESULTS AND DISCUSSION

We are interested in the effect of the coupling between surface and bulk molecules, and for that, we perform numerical simulations by MCES and by LMFA. To proceed, we choose a set of realistic SCO thermodynamic parameters, and for this purpose, we chose the data derived from experimental data of one of these SCO compounds: [Fe(btr)₂(NCS)₂], btr = 4,4'-bis-1,2,4-triazole [54], whose molecular structure is shown in Figure 5. The energy gap Δ and ln(g) are derived from the enthalpy change (ΔH ≈ 11 kJ/mol) and the entropy change (ΔS ≈ 50 J/mol/K), respectively: Δ = ΔH/R ≈ 1,300 K), ln(g) = ΔS/R ≈ 6.01, and R is the perfect gas constant. The equilibrium temperature of the system T<sub>eq</sub> is deduced by T<sub>eq</sub> = Δ/k<sub>B</sub>/ln(g), which leads to ≈ 216.3 K. In the first step, in the case of a 6 × 6 configuration, the values of the short-range interaction parameters J<sub>bb</sub><sup>MCES</sup> and J<sub>ss</sub><sup>MCES</sup> are set to 60 K, and the intensity of the short-range interaction parameter J<sub>bs</sub><sup>MCES</sup> between bulk and surface is gradually increased from 0 to 60 K. In the second step, the size effect is studied with each method. It should be added that the dimensions of the square lattice used in these simulations do not correspond, for the moment, to synthesized nanoparticles.

### The case J<sub>bs</sub><sup>MCES</sup> = 0

Figure 6 displays the results obtained in the case of a 6 × 6 square-shaped lattice, and the values of the interaction parameters *J* used by MCES and by LMFA, calculated from Equations (34) and (35), are gathered in Table 2. The curves highlight a two-step behavior in the form of two well-separated hysteresis loops. The hysteresis width is defined as ΔT = *T*<sub>up</sub> − *T*<sub>down</sub>, where *T*<sub>up</sub> is the ascending thermal transition temperature and *T*<sub>down</sub> is the descending thermal transition temperature. The average temperature between *T*<sub>up</sub> and *T*<sub>down</sub> is the equilibrium temperature *T*<sub>eq</sub>, which corresponds to a HS fraction equal to 1/2(Nhs = 1/2).

Figure 6 illustrates a qualitative and quantitative agreement between the calculations obtained by MCES and LMFA. The hysteresis loops, observable for each method, at lower temperatures, ~ 144 K, are related to the
<figure><img src="image_11.png" /><figcaption>**Figure 6.** Comparison of the thermal evolution of the HS fraction in a SCO system for a 6 × 6 square lattice: MCES (black squares), LMFA (blue up triangles). The computational parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 0$ K, $J_{bb}^{LMFA}/k_B = 36$ K, $J_{ss}^{LMFA}/k_B = 33.33$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.</figcaption></figure>

<table><caption>**Table 3.** Values of the thermal hysteresis width for molecules located in the bulk and in the surface obtained by MCES, $\Delta T_{bulk}^{MCES}$ and $\Delta T_{surf}^{MCES}$ and by LMFA, $\Delta T_{bulk}^{LMFA}$ and $\Delta T_{surf}^{LMFA}$.</caption><thead><tr><th>ΔT<sup>MCES</sup><sub>bulk</sub> [K]</th><th>ΔT<sup>LMFA</sup><sub>bulk</sub> [K]</th><th>ΔT<sup>MCES</sup><sub>surf</sub> [K]</th><th>ΔT<sup>LMFA</sup><sub>surf</sub> [K]</th></tr></thead><tbody><tr><td>4.8</td><td>7</td><td>7.6</td><td>10.5</td></tr></tbody></table>

The simulation parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 0$, $J_{bb}^{LMFA}/k_B = 36$ K, $J_{ss}^{LMFA}/k_B = 33.33$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.

behavior of surface molecules (edge and corner). These surface molecules commute from LS to HS state before the bulk molecules, and they drive the thermal transition. The hysteresis loops, positioned for each method at 185 K, are related to the bulk molecules. The two methods of simulation lead to identical values of the transition temperature in the surface ($T_{eqsurf}^{LMFA} = T_{eqsurf}^{MCES} \sim 144$ K) and in the bulk ($T_{eqbulk}^{LMFA} = T_{eqbulk}^{MCES} \sim 185$ K). The bulk and the surface appear as two systems that evolve almost independently of each other. Between 150 and 182 K, the curves reveal, both by MCES and by LMFA, a long intermediate plateau (mixture of LS and HS configurations) around $Nhs \approx 0.55$.

It can be noticed that the hysteresis widths $\Delta T$ are slightly larger by LMFA than by MCES. The values of the thermal hysteresis widths are summarized in Table 3.

### The case $J^{MCES} \neq 0$

The value of the short-range interaction parameter $J_{bs}^{MCES}/k_B$ is gradually increased up to 60 K. Figure 7, Figure 8 and Figure 9, show the results obtained respectively for $J_{bs}^{MCES}/k_B = 20$, $J_{bs}^{MCES}/k_B = 40$ and $J_{bs}^{MCES}/k_B = 60$ K in the case of a 6 × 6 lattice. The values of the interaction parameters $J$ used by MCES and by LMFA, obtained from Equations (34) and (35), are gathered in Table 4.

The phase diagrams of the 6 × 6 system, obtained by LMFA and by MCES, are plotted in Figure 10 and Figure 11, in the space coordinates, temperature vs. $J_{bs}^{MCES}/k_B$, for the bulk and for the surface molecules, respectively.

By MCES, as can be seen in Figure 10 and Figure 11, when the interaction parameter $J_{bs}^{MCES}/k_B$ increases, the equilibrium temperature of the surface molecules, $T_{eqsurf}^{MCES}$, is shifted to higher temperatures, whereas the equilibrium temperature of the bulk molecules, $T_{eqbulk}^{MCES}$, is shifted to lower temperatures. The hysteresis loops,
<figure><img src="image_11.png" /><figcaption>**Figure 6.** Comparison of the thermal evolution of the HS fraction in a SCO system for a 6 × 6 square lattice: MCES (black squares), LMFA (blue up triangles). The computational parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 0$ K, $J_{bb}^{LMFA}/k_B = 36$ K, $J_{ss}^{LMFA}/k_B = 33.33$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.</figcaption></figure>

<table><caption>**Table 3.** Values of the thermal hysteresis width for molecules located in the bulk and in the surface obtained by MCES, $\Delta T_{bulk}^{MCES}$ and $\Delta T_{surf}^{MCES}$ and by LMFA, $\Delta T_{bulk}^{LMFA}$ and $\Delta T_{surf}^{LMFA}$.</caption><thead><tr><th>ΔT<sup>MCES</sup><sub>bulk</sub> [K]</th><th>ΔT<sup>LMFA</sup><sub>bulk</sub> [K]</th><th>ΔT<sup>MCES</sup><sub>surf</sub> [K]</th><th>ΔT<sup>LMFA</sup><sub>surf</sub> [K]</th></tr></thead><tbody><tr><td>4.8</td><td>7</td><td>7.6</td><td>10.5</td></tr></tbody></table>

The simulation parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 0$, $J_{bb}^{LMFA}/k_B = 36$ K, $J_{ss}^{LMFA}/k_B = 33.33$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.

behavior of surface molecules (edge and corner). These surface molecules commute from LS to HS state before the bulk molecules, and they drive the thermal transition. The hysteresis loops, positioned for each method at 185 K, are related to the bulk molecules. The two methods of simulation lead to identical values of the transition temperature in the surface ($T_{eqsurf}^{LMFA} = T_{eqsurf}^{MCES} \sim 144$ K) and in the bulk ($T_{eqbulk}^{LMFA} = T_{eqbulk}^{MCES} \sim 185$ K). The bulk and the surface appear as two systems that evolve almost independently of each other. Between 150 and 182 K, the curves reveal, both by MCES and by LMFA, a long intermediate plateau (mixture of LS and HS configurations) around $Nhs \approx 0.55$.

It can be noticed that the hysteresis widths $\Delta T$ are slightly larger by LMFA than by MCES. The values of the thermal hysteresis widths are summarized in Table 3.

### The case $J^{MCES} \neq 0$

The value of the short-range interaction parameter $J_{bs}^{MCES}/k_B$ is gradually increased up to 60 K. Figure 7, Figure 8 and Figure 9, show the results obtained respectively for $J_{bs}^{MCES}/k_B = 20$, $J_{bs}^{MCES}/k_B = 40$ and $J_{bs}^{MCES}/k_B = 60$ K in the case of a 6 × 6 lattice. The values of the interaction parameters $J$ used by MCES and by LMFA, obtained from Equations (34) and (35), are gathered in Table 4.

The phase diagrams of the 6 × 6 system, obtained by LMFA and by MCES, are plotted in Figure 10 and Figure 11, in the space coordinates, temperature vs. $J_{bs}^{MCES}/k_B$, for the bulk and for the surface molecules, respectively.

By MCES, as can be seen in Figure 10 and Figure 11, when the interaction parameter $J_{bs}^{MCES}/k_B$ increases, the equilibrium temperature of the surface molecules, $T_{eqsurf}^{MCES}$, is shifted to higher temperatures, whereas the equilibrium temperature of the bulk molecules, $T_{eqbulk}^{MCES}$, is shifted to lower temperatures. The hysteresis loops,
<table><caption>**Table 4.** Correspondence between the values of short-range interactions *J*/*k*<sub>B</sub> used by MCES and by LMFA for a 6 × 6 lattice when *J*<sub>bs</sub><sup>MCES</sup>/*k*<sub>B</sub> = 20, 40, and 60 K</caption><thead><tr><th>J<sup>MCES</sup><sub>bb</sub>/k<sub>B</sub> [K]</th><th>J<sup>MCES</sup><sub>ss</sub>/k<sub>B</sub> [K]</th><th>J<sup>MCES</sup><sub>bs</sub>/k<sub>B</sub> [K]</th><th>J<sup>LMFA</sup><sub>bulk</sub>/k<sub>B</sub> [K]</th><th>J<sup>LMFA</sup><sub>surf</sub>/k<sub>B</sub> [K]</th></tr></thead><tbody><tr><td>60</td><td>60</td><td>20</td><td>44</td><td>42.22</td></tr><tr><td>60</td><td>60</td><td>40</td><td>52</td><td>51.11</td></tr><tr><td>60</td><td>60</td><td>60</td><td>60</td><td>60</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 7.** Comparison of the thermal evolution of the HS fraction in a SCO system for a 6 × 6 square lattice: MCES (black squares), LMFA (blue up triangles). The computational parameters are: *J*<sub>bb</sub><sup>MCES</sup>/k<sub>B</sub> = 60 K, *J*<sub>ss</sub><sup>MCES</sup>/k<sub>B</sub> = 60 K, *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 20 K, *J*<sub>bb</sub><sup>LMFA</sup>/k<sub>B</sub> = 44 K, *J*<sub>ss</sub><sup>LMFA</sup>/k<sub>B</sub> = 42.22 K, Δ/k<sub>B</sub> = 1, 300 K, G/k<sub>B</sub> = 172 K, L/k<sub>B</sub> = 290 K, and ln(g) = 6.01.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 8.** Comparison of the thermal evolution of the HS fraction in a SCO system for a 6 × 6 square lattice: MCES (black squares), LMFA (blue up triangles). The computational parameters are: *J*<sub>bb</sub><sup>MCES</sup>/k<sub>B</sub> = 60 K, *J*<sub>ss</sub><sup>MCES</sup>/k<sub>B</sub> = 60 K, *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 40 K, *J*<sub>bb</sub><sup>LMFA</sup>/k<sub>B</sub> = 52 K, *J*<sub>ss</sub><sup>LMFA</sup>/k<sub>B</sub> = 51.11 K, Δ/k<sub>B</sub> = 1, 300 K, G/k<sub>B</sub> = 172 K, L/k<sub>B</sub> = 290 K, and ln(g) = 6.01.</figcaption></figure>

quite distinct for *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 0, therefore, tend to overlap. The thermal hysteresis loop Δ*T*<sub>bulk</sub><sup>MCES</sup> increases (~ 4.8 K wide when *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 0 and ~ 8 K wide when *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 60). Simultaneously, the thermal hysteresis loop Δ*T*<sub>surf</sub><sup>MCES</sup> decreases (~ 7.6 K wide for *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 0 and ~ 4 K wide for *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 60).

In LMFA, as shown in Figure 10 and Figure 11, an opposite behavior is observed for surface molecules. First of all, when the interaction parameter *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> increases, the equilibrium temperature of the surface molecules, *T*<sub>eq surf</sub><sup>LMFA</sup>, is shifted to lower temperatures. Then, the thermal hysteresis loop Δ*T*<sub>surf</sub><sup>LMFA</sup> increases (~ 10.5 K wide for *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 0 and ~ 38.5 K wide for *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 60). These results are summarized in Table 5 and the results obtained when *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 0 are recalled for a better understanding of the results. The LMFA
<table><caption>Table 4. Correspondence between the values of short-range interactions *J*/*k*<sub>B</sub> used by MCES and by LMFA for a 6 × 6 lattice when *J*<sub>bs</sub><sup>MCES</sup>/*k*<sub>B</sub> = 20, 40, and 60 K</caption><thead><tr><th>J<sup>MCES</sup><sub>bb</sub>/k<sub>B</sub> [K]</th><th>J<sup>MCES</sup><sub>ss</sub>/k<sub>B</sub> [K]</th><th>J<sup>MCES</sup><sub>bs</sub>/k<sub>B</sub> [K]</th><th>J<sup>LMFA</sup><sub>bulk</sub>/k<sub>B</sub> [K]</th><th>J<sup>LMFA</sup><sub>surf</sub>/k<sub>B</sub> [K]</th></tr></thead><tbody><tr><td>60</td><td>60</td><td>20</td><td>44</td><td>42.22</td></tr><tr><td>60</td><td>60</td><td>40</td><td>52</td><td>51.11</td></tr><tr><td>60</td><td>60</td><td>60</td><td>60</td><td>60</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 7.** Comparison of the thermal evolution of the HS fraction in a SCO system for a 6 × 6 square lattice: MCES (black squares), LMFA (blue up triangles). The computational parameters are: *J*<sub>bb</sub><sup>MCES</sup>/k<sub>B</sub> = 60 K, *J*<sub>ss</sub><sup>MCES</sup>/k<sub>B</sub> = 60 K, *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 20 K, *J*<sub>bb</sub><sup>LMFA</sup>/k<sub>B</sub> = 44 K, *J*<sub>ss</sub><sup>LMFA</sup>/k<sub>B</sub> = 42.22 K, Δ/k<sub>B</sub> = 1, 300 K, G/k<sub>B</sub> = 172 K, L/k<sub>B</sub> = 290 K, and ln(g) = 6.01.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 8.** Comparison of the thermal evolution of the HS fraction in a SCO system for a 6 × 6 square lattice: MCES (black squares), LMFA (blue up triangles). The computational parameters are: *J*<sub>bb</sub><sup>MCES</sup>/k<sub>B</sub> = 60 K, *J*<sub>ss</sub><sup>MCES</sup>/k<sub>B</sub> = 60 K, *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 40 K, *J*<sub>bb</sub><sup>LMFA</sup>/k<sub>B</sub> = 52 K, *J*<sub>ss</sub><sup>LMFA</sup>/k<sub>B</sub> = 51.11 K, Δ/k<sub>B</sub> = 1, 300 K, G/k<sub>B</sub> = 172 K, L/k<sub>B</sub> = 290 K, and ln(g) = 6.01.</figcaption></figure>

quite distinct for *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 0, therefore, tend to overlap. The thermal hysteresis loop Δ*T*<sub>bulk</sub><sup>MCES</sup> increases (~ 4.8 K wide when *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 0 and ~ 8 K wide when *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 60). Simultaneously, the thermal hysteresis loop Δ*T*<sub>surf</sub><sup>MCES</sup> decreases (~ 7.6 K wide for *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 0 and ~ 4 K wide for *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 60).

In LMFA, as shown in Figure 10 and Figure 11, an opposite behavior is observed for surface molecules. First of all, when the interaction parameter *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> increases, the equilibrium temperature of the surface molecules, *T*<sub>eq surf</sub><sup>LMFA</sup>, is shifted to lower temperatures. Then, the thermal hysteresis loop Δ*T*<sub>surf</sub><sup>LMFA</sup> increases (~ 10.5 K wide for *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 0 and ~ 38.5 K wide for *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 60). These results are summarized in Table 5 and the results obtained when *J*<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> = 0 are recalled for a better understanding of the results. The LMFA
<figure><img src="image_14.png" /><figcaption>**Figure 9.** Comparison of the thermal evolution of the HS fraction in a SCO system for a 6 × 6 square lattice: MCES (black squares), LMFA (blue up triangles). The computational parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 60$ K, $J_{bb}^{LMFA}/k_B = 60$ K, $J_{ss}^{LMFA}/k_B = 60$ K, $Δ/k_B = 1, 300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $ln(g) = 6.01$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 10.** Phase diagram T vs. $J_{bs}^{MCES}/k_B$ for a 6 × 6 system. The red, blue and black squares correspond, respectively, to the upper and lower transitions ($T_{up}$) and ($T_{down}$) temperatures of the thermal HS fraction and to the equilibrium temperature of the bulk by LMFA. The red, blue and black stars correspond, respectively, to the upper and lower transitions ($T_{up}$) and ($T_{down}$) temperatures of the thermal HS fraction and to the equilibrium temperature of the bulk by MCES. The computational parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $Δ/k_B = 1, 300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $ln(g) = 6.01$.</figcaption></figure>

<table><caption>**Table 5.** Values of the thermal hysteresis width for molecules located in the bulk and in the surface in MCES, $ΔT_{bulk}^{MCES}$ and $ΔT_{surf}^{MCES}$, and in LMFA, $ΔT_{bulk}^{LMFA}$ and $ΔT_{surf}^{LMFA}$ as a function of $J_{bs}^{MCES}/k_B$</caption><thead><tr><th>J<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> [K]</th><th>ΔT<sub>bulk</sub><sup>MCES</sup> [K]</th><th>ΔT<sub>bulk</sub><sup>MCES</sup> [K]</th><th>ΔT<sub>bulk</sub><sup>LMFA</sup> [K]</th><th>ΔT<sub>surf</sub><sup>LMFA</sup> [K]</th></tr></thead><tbody><tr><td>0</td><td>4.8</td><td>7.6</td><td>7</td><td>10.5</td></tr><tr><td>20</td><td>6.4</td><td>6.2</td><td>17.5</td><td>20.36</td></tr><tr><td>40</td><td>8.0</td><td>5.0</td><td>30.5</td><td>29</td></tr><tr><td>60</td><td>8.8</td><td>4.0</td><td>45.0</td><td>38.5</td></tr></tbody></table>

The computational parameters are: $Δ/k_B = 1, 300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $ln(g) = 6.01$.

drastically accentuates the hysteresis effect, especially for surface molecules. The system seems constrained, and the tendency for the two hysteresis to overlap is not observed.
<figure><img src="image_14.png" /><figcaption>**Figure 9.** Comparison of the thermal evolution of the HS fraction in a SCO system for a 6 × 6 square lattice: MCES (black squares), LMFA (blue up triangles). The computational parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 60$ K, $J_{bb}^{LMFA}/k_B = 60$ K, $J_{ss}^{LMFA}/k_B = 60$ K, $Δ/k_B = 1, 300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $ln(g) = 6.01$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 10.** Phase diagram T vs. $J_{bs}^{MCES}/k_B$ for a 6 × 6 system. The red, blue and black squares correspond, respectively, to the upper and lower transitions ($T_{up}$) and ($T_{down}$) temperatures of the thermal HS fraction and to the equilibrium temperature of the bulk by LMFA. The red, blue and black stars correspond, respectively, to the upper and lower transitions ($T_{up}$) and ($T_{down}$) temperatures of the thermal HS fraction and to the equilibrium temperature of the bulk by MCES. The computational parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $Δ/k_B = 1, 300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $ln(g) = 6.01$.</figcaption></figure>

<table><caption>**Table 5.** Values of the thermal hysteresis width for molecules located in the bulk and in the surface in MCES, $ΔT_{bulk}^{MCES}$ and $ΔT_{surf}^{MCES}$, and in LMFA, $ΔT_{bulk}^{LMFA}$ and $ΔT_{surf}^{LMFA}$ as a function of $J_{bs}^{MCES}/k_B$</caption><thead><tr><th>J<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub> [K]</th><th>ΔT<sub>bulk</sub><sup>MCES</sup> [K]</th><th>ΔT<sub>bulk</sub><sup>MCES</sup> [K]</th><th>ΔT<sub>bulk</sub><sup>LMFA</sup> [K]</th><th>ΔT<sub>surf</sub><sup>LMFA</sup> [K]</th></tr></thead><tbody><tr><td>0</td><td>4.8</td><td>7.6</td><td>7</td><td>10.5</td></tr><tr><td>20</td><td>6.4</td><td>6.2</td><td>17.5</td><td>20.36</td></tr><tr><td>40</td><td>8.0</td><td>5.0</td><td>30.5</td><td>29</td></tr><tr><td>60</td><td>8.8</td><td>4.0</td><td>45.0</td><td>38.5</td></tr></tbody></table>

The computational parameters are: $Δ/k_B = 1, 300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $ln(g) = 6.01$.

drastically accentuates the hysteresis effect, especially for surface molecules. The system seems constrained, and the tendency for the two hysteresis to overlap is not observed.
<figure><img src="image_16.png" /><figcaption>**Figure 11.** Phase diagram T vs. $J_{bs}^{MCES}/k_B$ for a 6 × 6 system. The red, blue and black squares correspond, respectively, to the upper and lower transitions ($T_{up}$) and ($T_{down}$) temperatures of the thermal HS fraction and to the equilibrium temperature of the surface by LMFA. The red, blue and black stars correspond, respectively, to the upper and lower transitions ($T_{up}$) and ($T_{down}$) temperatures of the thermal HS fraction and to the equilibrium temperature of the surface by MCES. The computational parameters are: $J_{bs}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.</figcaption></figure>

<table><caption>**Table 6.** Correspondence between the values of short-range interactions *J*/*k*<sub>B</sub> used by the MCES method and by LMFA in the case of a 7 × 7 square compound</caption><thead><tr><th>J<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub>[K]</th><th>J<sub>ss</sub><sup>MCES</sup>/k<sub>B</sub>[K]</th><th>J<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub>[K]</th><th>J<sub>bulk</sub><sup>LMFA</sup>/k<sub>B</sub>[K]</th><th>J<sub>surf</sub><sup>LMFA</sup>/k<sub>B</sub>[K]</th></tr></thead><tbody><tr><td>60</td><td>60</td><td>0</td><td>40</td><td>32.72</td></tr><tr><td>60</td><td>60</td><td>20</td><td>46.66</td><td>41.81</td></tr><tr><td>60</td><td>60</td><td>40</td><td>52.33</td><td>50.90</td></tr><tr><td>60</td><td>60</td><td>60</td><td>60</td><td>60</td></tr></tbody></table>

<table><caption>**Table 7.** Values of the equilibrium temperatures of the bulk and surface molecules obtained in MCES and in LMFA, as a function of the system size (6 × 6 and 7 × 7 lattices)</caption><thead><tr><th>Latticesize</th><th>T<sub>eq bulk</sub><sup>MCES</sup>[K]</th><th>T<sub>eq bulk</sub><sup>LMFA</sup>[K]</th><th>T<sub>eq surf</sub><sup>MCES</sup>[K]</th><th>T<sub>eq surf</sub><sup>LMFA</sup>[K]</th></tr></thead><tbody><tr><td>6 × 6</td><td>185.0</td><td>185.0</td><td>144.0</td><td>144.0</td></tr><tr><td>7 × 7</td><td>187.75</td><td>187.75</td><td>148.6</td><td>148.6</td></tr></tbody></table>

The computational parameters are: $J_{bs}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 0$ K, $J_{bb}^{LMFA}/k_B = 40$ K, $J_{ss}^{LMFA}/k_B = 32.72$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.

## Size effect

The purpose of this part is to compare the two simulation techniques in the case of a larger system of size 7 × 7. The results are presented and analyzed for $J_{bs}^{MCES}/k_B$ varying from 0 to 60 K. The corresponding values of $J_{bulk}^{LMFA}/k_B$ and $J_{surf}^{LMFA}/k_B$ calculated from Equations (34) and (35) are given in Table 6.

### The case $J_{bs}^{MCES} = 0$

The results obtained and presented in Figure 12, highlight three important points. First of all, as previously observed in the case of the 6 × 6 lattice, when $J_{bs}^{MCES}/k_B = 0$, the equilibrium temperatures obtained in the bulk and in the surface are strictly identical with the two methods. The values of the equilibrium temperatures are reported in Table 7. Then, the equilibrium temperatures, in the bulk and in the surface, are slightly higher than those observed in the case of the 6 × 6 lattice (See also Table 7 for a comparison). This point may be explained by the fact that for larger sizes of particles, $T_{eq}$ and $T_c$. shift to higher temperatures due to the increase in the number of couplings.
<figure><img src="image_16.png" /><figcaption>**Figure 11.** Phase diagram T vs. $J_{bs}^{MCES}/k_B$ for a 6 × 6 system. The red, blue and black squares correspond, respectively, to the upper and lower transitions ($T_{up}$) and ($T_{down}$) temperatures of the thermal HS fraction and to the equilibrium temperature of the surface by LMFA. The red, blue and black stars correspond, respectively, to the upper and lower transitions ($T_{up}$) and ($T_{down}$) temperatures of the thermal HS fraction and to the equilibrium temperature of the surface by MCES. The computational parameters are: $J_{bs}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.</figcaption></figure>

<table><caption>**Table 6.** Correspondence between the values of short-range interactions $J/k_B$ used by the MCES method and by LMFA in the case of a 7 × 7 square compound</caption><thead><tr><th>J<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub>[K]</th><th>J<sub>ss</sub><sup>MCES</sup>/k<sub>B</sub>[K]</th><th>J<sub>bs</sub><sup>MCES</sup>/k<sub>B</sub>[K]</th><th>J<sub>bulk</sub><sup>LMFA</sup>/k<sub>B</sub>[K]</th><th>J<sub>surf</sub><sup>LMFA</sup>/k<sub>B</sub>[K]</th></tr></thead><tbody><tr><td>60</td><td>60</td><td>0</td><td>40</td><td>32.72</td></tr><tr><td>60</td><td>60</td><td>20</td><td>46.66</td><td>41.81</td></tr><tr><td>60</td><td>60</td><td>40</td><td>52.33</td><td>50.90</td></tr><tr><td>60</td><td>60</td><td>60</td><td>60</td><td>60</td></tr></tbody></table>

<table><caption>**Table 7.** Values of the equilibrium temperatures of the bulk and surface molecules obtained in MCES and in LMFA, as a function of the system size (6 × 6 and 7 × 7 lattices)</caption><thead><tr><th>Latticesize</th><th>T<sub>eq bulk</sub><sup>MCES</sup> [K]</th><th>T<sub>eq bulk</sub><sup>LMFA</sup> [K]</th><th>T<sub>eq surf</sub><sup>MCES</sup> [K]</th><th>T<sub>eq surf</sub><sup>LMFA</sup> [K]</th></tr></thead><tbody><tr><td>6 × 6</td><td>185.0</td><td>185.0</td><td>144.0</td><td>144.0</td></tr><tr><td>7 × 7</td><td>187.75</td><td>187.75</td><td>148.6</td><td>148.6</td></tr></tbody></table>

The computational parameters are: $J_{bs}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 0$ K, $J_{bb}^{LMFA}/k_B = 40$ K, $J_{ss}^{LMFA}/k_B = 32.72$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.

## Size effect

The purpose of this part is to compare the two simulation techniques in the case of a larger system of size 7 × 7. The results are presented and analyzed for $J_{bs}^{MCES}/k_B$ varying from 0 to 60 K. The corresponding values of $J_{bulk}^{LMFA}/k_B$ and $J_{surf}^{LMFA}/k_B$ calculated from Equations (34) and (35) are given in Table 6.

### The case $J_{bs}^{MCES} = 0$

The results obtained and presented in Figure 12, highlight three important points. First of all, as previously observed in the case of the 6 × 6 lattice, when $J_{bs}^{MCES}/k_B = 0$, the equilibrium temperatures obtained in the bulk and in the surface are strictly identical with the two methods. The values of the equilibrium temperatures are reported in Table 7. Then, the equilibrium temperatures, in the bulk and in the surface, are slightly higher than those observed in the case of the 6 × 6 lattice (See also Table 7 for a comparison). This point may be explained by the fact that for larger sizes of particles, $T_{eq}$ and $T_c$. shift to higher temperatures due to the increase in the number of couplings.
<figure><img src="image_17.png" /><figcaption>**Figure 12.** Comparison of the thermal evolution of the HS fraction in a SCO system with a size 7 × 7 showing an excellent agreement between the MCES method (black squares) and the LMFA (blue up triangles). The calculation parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 0$ K, $J_{bb}^{LMFA}/k_B = 40$ K, $J_{ss}^{LMFA}/k_B = 32.72$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.</figcaption></figure>

<table><caption>**Table 8.** Values of the thermal hysteresis width for molecules located in the bulk and in the surface obtained by MCES, $\Delta T_{bulk}^{MCES}$ and $\Delta T_{surf}^{MCES}$, and by LMFA, $\Delta T_{bulk}^{LMFA}$ and $\Delta T_{surf}^{LMFA}$.</caption><thead><tr><th>Latticesize</th><th><i>T</i><sub>bulk</sub><sup>MCES</sup> [K]</th><th><i>T</i><sub>bulk</sub><sup>LMFA</sup> [K]</th><th><i>T</i><sub>surf</sub><sup>MCES</sup> [K]</th><th><i>T</i><sub>surf</sub><sup>LMFA</sup> [K]</th></tr></thead><tbody><tr><td>6 × 6</td><td>4.8</td><td>7</td><td>7.6</td><td>10.5</td></tr><tr><td>7 × 7</td><td>10.4</td><td>14</td><td>2.8</td><td>5</td></tr></tbody></table>

The computational parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 0$ K, $J_{bb}^{LMFA}/k_B = 40$ K, $J_{ss}^{LMFA}/k_B = 32.72$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.

Finally, in a larger particle of size 7 × 7, the hysteresis associated with the bulk are wider than the hysteresis associated with the surface molecules by LMFA and by MCES. The proportion of molecules in the bulk (~ 51%) being greater than that of molecules in the surface (~ 49%), the bulk molecules drive the transition. In the 6 × 6 compound, the hysteresis related to surface molecules was greatest in both MCES and LMFA. These results are summarized in Table 8.

### The case $J_{bs}^{MCES} \neq 0$

The phase diagrams of the 7 × 7 system obtained by LMFA and by MCES are plotted in Figure 13 and Figure 14, in the space coordinates, temperature vs. $J_{bs}^{MCES}/k_B$, respectively, for the bulk and for the surface molecules.

When the value of the interaction parameter $J_{bs}^{MCES}/k_B$ is increased up to 60 K, the three following items are especially significant.

First, for each method, the curves presented in Figure 12 and Figure 13 reveal that the equilibrium temperatures of the surface and of the bulk in the 7 × 7 lattice are higher than those obtained in the 6 × 6 lattice [Figure 10 and Figure 11].

In addition, as already observed in the case of a 6 × 6 lattice and by the LMFA, when the short-range interaction $J_{bs}^{MCES}/k_B$ is increased, the equilibrium temperature of the surface molecules decreases. We further note that the increase of the lattice size leads to a decrease of the width of the hysteresis in the surface. By MCES, the hysteresis transition observed for surface molecules disappears ($\Delta T_{surf}^{MCES} = 0$ K) in favor of an abrupt transition [Table 9].
<figure><img src="image_17.png" /><figcaption>**Figure 12.** Comparison of the thermal evolution of the HS fraction in a SCO system with a size 7 × 7 showing an excellent agreement between the MCES method (black squares) and the LMFA (blue up triangles). The calculation parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 0$ K, $J_{bb}^{LMFA}/k_B = 40$ K, $J_{ss}^{LMFA}/k_B = 32.72$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.</figcaption></figure>

<table><caption>**Table 8.** Values of the thermal hysteresis width for molecules located in the bulk and in the surface obtained by MCES, $\Delta T_{bulk}^{MCES}$ and $\Delta T_{surf}^{MCES}$, and by LMFA, $\Delta T_{bulk}^{LMFA}$ and $\Delta T_{surf}^{LMFA}$.</caption><thead><tr><th>Latticesize</th><th><i>T</i><sub>bulk</sub><sup>MCES</sup> [K]</th><th><i>T</i><sub>bulk</sub><sup>LMFA</sup> [K]</th><th><i>T</i><sub>surf</sub><sup>MCES</sup> [K]</th><th><i>T</i><sub>surf</sub><sup>LMFA</sup> [K]</th></tr></thead><tbody><tr><td>6 × 6</td><td>4.8</td><td>7</td><td>7.6</td><td>10.5</td></tr><tr><td>7 × 7</td><td>10.4</td><td>14</td><td>2.8</td><td>5</td></tr></tbody></table>

The computational parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 0$ K, $J_{bb}^{LMFA}/k_B = 40$ K, $J_{ss}^{LMFA}/k_B = 32.72$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.

Finally, in a larger particle of size 7 × 7, the hysteresis associated with the bulk are wider than the hysteresis associated with the surface molecules by LMFA and by MCES. The proportion of molecules in the bulk (~ 51%) being greater than that of molecules in the surface (~ 49%), the bulk molecules drive the transition. In the 6 × 6 compound, the hysteresis related to surface molecules was greatest in both MCES and LMFA. These results are summarized in Table 8.

### The case $J_{bs}^{MCES} \neq 0$

The phase diagrams of the 7 × 7 system obtained by LMFA and by MCES are plotted in Figure 13 and Figure 14, in the space coordinates, temperature vs. $J_{bs}^{MCES}/k_B$, respectively, for the bulk and for the surface molecules.

When the value of the interaction parameter $J_{bs}^{MCES}/k_B$ is increased up to 60 K, the three following items are especially significant.

First, for each method, the curves presented in Figure 12 and Figure 13 reveal that the equilibrium temperatures of the surface and of the bulk in the 7 × 7 lattice are higher than those obtained in the 6 × 6 lattice [Figure 10 and Figure 11].

In addition, as already observed in the case of a 6 × 6 lattice and by the LMFA, when the short-range interaction $J_{bs}^{MCES}/k_B$ is increased, the equilibrium temperature of the surface molecules decreases. We further note that the increase of the lattice size leads to a decrease of the width of the hysteresis in the surface. By MCES, the hysteresis transition observed for surface molecules disappears ($\Delta T_{surf}^{MCES} = 0$ K) in favor of an abrupt transition [Table 9].
<figure><img src="image_18.png" /><figcaption>**Figure 13.** Phase diagram T vs. ($J_{bs}^{MCES}/k_B$) for a 7 × 7 system. The red, blue and black squares correspond, respectively, to the upper and lower transitions ($T_{up}$) and ($T_{down}$) temperatures of the thermal HS fraction and to the equilibrium temperature of the bulk by LMFA. The red, blue and black stars correspond, respectively, to the upper and lower transitions ($T_{up}$) and ($T_{down}$) temperatures of the thermal HS fraction and to the equilibrium temperature of the bulk by MCES. The computational parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 14.** Phase diagram T versus ($J_{bs}^{MCES}/k_B$) for a 7 × 7 system. The red, blue and black squares correspond, respectively, to the upper and lower transitions ($T_{up}$) and ($T_{down}$) temperatures of the thermal HS fraction and to the equilibrium temperature of the surface by LMFA. The red, blue and black stars correspond, respectively, to the upper and lower transitions ($T_{up}$) and ($T_{down}$) temperatures of the thermal HS fraction and to the equilibrium temperature of the surface by MCES. The computational parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.</figcaption></figure>

<table><caption>**Table 9.** Values of the thermal hysteresis width for molecules located in the bulk and in the surface by MCES, $\Delta T_{bulk}^{MCES}$ and $\Delta T_{surf}^{MCES}$, and by LMFA, $\Delta T_{bulk}^{LMFA}$ and $\Delta T_{surf}^{LMFA}$, for a 6 × 6 and a 7 × 7 SCO particle</caption><thead><tr><th>Latticesize</th><th><i>T</i><sub>bulk</sub><sup>MCES</sup> [K]</th><th><i>T</i><sub>bulk</sub><sup>LMFA</sup> [K]</th><th><i>T</i><sub>surf</sub><sup>MCES</sup> [K]</th><th><i>T</i><sub>surf</sub><sup>LMFA</sup> [K]</th></tr></thead><tbody><tr><td>6 × 6</td><td>8.8</td><td>45</td><td>4</td><td>38.5</td></tr><tr><td>7 × 7</td><td>23.2</td><td>48</td><td>0</td><td>32</td></tr></tbody></table>

The computational parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 60$ K, $J_{ss}^{LMFA}/k_B = 60$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.

Simulations using the LMFA were carried for very large lattices, 20 × 20, 50 × 50 and 200 × 200. Figure 15 illustrates the results obtained for two sets of parameters corresponding to $J_{bs}^{MCES}/k_B = 0$ K and $J_{bs}^{MCES}/k_B = 60$ K and shows that, upon increasing the lattice size, the two-step transition vanishes. The curve closely resembles a single hysteresis loop for the 200 × 200 lattice, and the equilibrium temperature tends towards that
<figure><img src="image_18.png" /><figcaption>**Figure 13.** Phase diagram T vs. ($J_{bs}^{MCES}/k_B$) for a 7 × 7 system. The red, blue and black squares correspond, respectively, to the upper and lower transitions ($T_{up}$) and ($T_{down}$) temperatures of the thermal HS fraction and to the equilibrium temperature of the bulk by LMFA. The red, blue and black stars correspond, respectively, to the upper and lower transitions ($T_{up}$) and ($T_{down}$) temperatures of the thermal HS fraction and to the equilibrium temperature of the bulk by MCES. The computational parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 14.** Phase diagram T versus ($J_{bs}^{MCES}/k_B$) for a 7 × 7 system. The red, blue and black squares correspond, respectively, to the upper and lower transitions ($T_{up}$) and ($T_{down}$) temperatures of the thermal HS fraction and to the equilibrium temperature of the surface by LMFA. The red, blue and black stars correspond, respectively, to the upper and lower transitions ($T_{up}$) and ($T_{down}$) temperatures of the thermal HS fraction and to the equilibrium temperature of the surface by MCES. The computational parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.</figcaption></figure>

<table><caption>**Table 9.** Values of the thermal hysteresis width for molecules located in the bulk and in the surface by MCES, $\Delta T_{bulk}^{MCES}$ and $\Delta T_{surf}^{MCES}$, and by LMFA, $\Delta T_{bulk}^{LMFA}$ and $\Delta T_{surf}^{LMFA}$, for a 6 × 6 and a 7 × 7 SCO particle</caption><thead><tr><th>Latticesize</th><th><i>T</i><sub>bulk</sub><sup>MCES</sup> [K]</th><th><i>T</i><sub>bulk</sub><sup>LMFA</sup> [K]</th><th><i>T</i><sub>surf</sub><sup>MCES</sup> [K]</th><th><i>T</i><sub>surf</sub><sup>LMFA</sup> [K]</th></tr></thead><tbody><tr><td>6 × 6</td><td>8.8</td><td>45</td><td>4</td><td>38.5</td></tr><tr><td>7 × 7</td><td>23.2</td><td>48</td><td>0</td><td>32</td></tr></tbody></table>

The computational parameters are: $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 60$ K, $J_{ss}^{LMFA}/k_B = 60$ K, $\Delta/k_B = 1,300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.

Simulations using the LMFA were carried for very large lattices, 20 × 20, 50 × 50 and 200 × 200. **Figure 15** illustrates the results obtained for two sets of parameters corresponding to $J_{bs}^{MCES}/k_B = 0$ K and $J_{bs}^{MCES}/k_B = 60$ K and shows that, upon increasing the lattice size, the two-step transition vanishes. The curve closely resembles a single hysteresis loop for the 200 × 200 lattice, and the equilibrium temperature tends towards that
<figure><img src="image_20.png" /><figcaption>**Figure 15.** Thermal dependence of the HS fraction studied with the LMFA for square-shaped particles of different sizes: 20 × 20 (blue up triangles), 50 × 50 (orange circles), and 200 × 200 (green stars). The computational parameters are: (A) $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 0$ K, $J_{bb}^{LMFA}/k_B = 36$ K, $J_{ss}^{LMFA}/k_B = 33.33$ K; (B) $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 60$ K, $J_{bb}^{LMFA}/k_B = 60$ K, $J_{ss}^{LMFA}/k_B = 60$ K. For all calculations: $\Delta/k_B = 1, 300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.</figcaption></figure>

of the bulk {$T_{eq} = \Delta/[k_B \ln(g) \approx 216$ K)}. The results obtained in the case of the 7 × 7 lattice suggest that if we
could do simulations on larger lattices with MCES, a single hysteresis loop with an equilibrium temperature
close to 216 K would also be obtained.

CONCLUSION

The present contribution, a simulation work, with the purpose of giving some hints to colleagues who synthe-
size spin-transition nanoparticles, is focused on meticulous investigations of the thermal properties of 2D-SCO
nanoparticles using an adapted two-state Ising-like Hamiltonian. At this end, three types of interaction terms
are considered: long and short-range interactions originating from electronic and elastic molecular couplings
inside the lattice, to which we added a specific energetic contribution accounting for the interactions between
surface molecules with their surrounding environment (inside and outside the lattice). The model is solved in
two distinct ways, which involve different levels of sophistication: (i) MCES techniques, which are exact meth-
ods and (ii) LMFA based on mean-field theory but taking into account the lattice shape. In the former case,
three types of two-body short-range interaction terms are involved between bulk-bulk ($J_{bb}$), surface-surface
($J_{ss}$), and bulk-surface ($J_{bs}$) sites. The thermal dependence of the “order parameter”, $\langle \sigma \rangle$, that is, the av-
erage fictitious magnetization, is calculated by using the Boltzmann distribution obtained from the calculated
density of the macrostates. On the other hand, in the LMFA, two interaction terms corresponding to three
types of sites, corner, edge, and bulk, are included in the calculations, taking into account the dependence of
the coordination number on the site location. The results of the two methods have been compared for different
square lattice sizes, and a good quantitative agreement is obtained between them in the case $J_{bs} = 0$, where both
equilibrium temperatures and hysteresis widths were found in very good agreement. When the bulk-surface
interaction $J_{bs}$ is increased, the hysteresis observed by MCES in the bulk and in the surface tend to be close to
each other and finally overlap for strong $J_{bs}$ values. In contrast, by LMFA, this phenomenon is not observed,
and furthermore, the hysteresis effects are accentuated, particularly for surface molecules. Moreover, with the
MCES method, the increase of $J_{bs}$ results in a shift of the surface (resp. bulk) transition temperatures to higher
 resp. lower) values, while in LMFA, an opposite effect is obtained for surface molecules.

Our comparison of the two techniques (MCES and LMFA) shows that when dealing with small samples, LMFA
yields results that are in perfect agreement with MCES in the case where the interaction between the surface and
the bulk is zero ($J_{bs} = 0$). However, as the interaction $J_{bs}$ increases, causing important interferences between the
two regions of the lattice, the results provided by LMFA increasingly differ from those obtained with the MCES
method, although the general trends remain the same. While considering a very weak interaction between the
<figure><img src="image_20.png" /><figcaption>**Figure 15.** Thermal dependence of the HS fraction studied with the LMFA for square-shaped particles of different sizes: 20 × 20 (blue up triangles), 50 × 50 (orange circles), and 200 × 200 (green stars). The computational parameters are: (A) $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 0$ K, $J_{bb}^{LMFA}/k_B = 36$ K, $J_{ss}^{LMFA}/k_B = 33.33$ K; (B) $J_{bb}^{MCES}/k_B = 60$ K, $J_{ss}^{MCES}/k_B = 60$ K, $J_{bs}^{MCES}/k_B = 60$ K, $J_{bb}^{LMFA}/k_B = 60$ K, $J_{ss}^{LMFA}/k_B = 60$ K. For all calculations: $\Delta/k_B = 1, 300$ K, $G/k_B = 172$ K, $L/k_B = 290$ K, and $\ln(g) = 6.01$.</figcaption></figure>

of the bulk {$T_{eq} = \Delta/[k_B \ln(g) \approx 216$ K)}. The results obtained in the case of the 7 × 7 lattice suggest that if we could do simulations on larger lattices with MCES, a single hysteresis loop with an equilibrium temperature close to 216 K would also be obtained.

## CONCLUSION

The present contribution, a simulation work, with the purpose of giving some hints to colleagues who synthesize spin-transition nanoparticles, is focused on meticulous investigations of the thermal properties of 2D-SCO nanoparticles using an adapted two-state Ising-like Hamiltonian. At this end, three types of interaction terms are considered: long and short-range interactions originating from electronic and elastic molecular couplings inside the lattice, to which we added a specific energetic contribution accounting for the interactions between surface molecules with their surrounding environment (inside and outside the lattice). The model is solved in two distinct ways, which involve different levels of sophistication: (i) MCES techniques, which are exact methods and (ii) LMFA based on mean-field theory but taking into account the lattice shape. In the former case, three types of two-body short-range interaction terms are involved between bulk-bulk ($J_{bb}$), surface-surface ($J_{ss}$), and bulk-surface ($J_{bs}$) sites. The thermal dependence of the “order parameter”, $\langle \sigma \rangle$, that is, the average fictitious magnetization, is calculated by using the Boltzmann distribution obtained from the calculated density of the macrostates. On the other hand, in the LMFA, two interaction terms corresponding to three types of sites, corner, edge, and bulk, are included in the calculations, taking into account the dependence of the coordination number on the site location. The results of the two methods have been compared for different square lattice sizes, and a good quantitative agreement is obtained between them in the case $J_{bs} = 0$, where both equilibrium temperatures and hysteresis widths were found in very good agreement. When the bulk-surface interaction $J_{bs}$ is increased, the hysteresis observed by MCES in the bulk and in the surface tend to be close to each other and finally overlap for strong $J_{bs}$ values. In contrast, by LMFA, this phenomenon is not observed, and furthermore, the hysteresis effects are accentuated, particularly for surface molecules. Moreover, with the MCES method, the increase of $J_{bs}$ results in a shift of the surface (resp. bulk) transition temperatures to higher (resp. lower) values, while in LMFA, an opposite effect is obtained for surface molecules.

Our comparison of the two techniques (MCES and LMFA) shows that when dealing with small samples, LMFA yields results that are in perfect agreement with MCES in the case where the interaction between the surface and the bulk is zero ($J_{bs} = 0$). However, as the interaction $J_{bs}$ increases, causing important interferences between the two regions of the lattice, the results provided by LMFA increasingly differ from those obtained with the MCES method, although the general trends remain the same. While considering a very weak interaction between the
surface and the bulk in the same compound may look unnatural, it may be possible in nanostructured core-shell nanoparticles, where the core and the shell are of different compositions. In addition, it is important to keep in mind that the systems considered are particularly small, and as the size of the system increases, the predominance of the surface, and then $J_{bs}$ interaction, vanishes, which, in turn, means that LMFA will most likely still provide correct results given the large number of molecules.

Finally, it is worth mentioning that the surface properties have been limited here to the outer layer of the lattice, while the concept of the surface may integrate some inner layers close to the interface material/air, which also have specific properties different from those of the bulk. Extensions of the present work to include the above-cited aspects are in progress, as well as other packing clusters such as rhombohedral or faced centered cubic (FCC). The results of this simulation work can also be extended to explain the behavior of polycrystalline SCO films on the $Al_2O_3$ substrate[40] and the influence of the substrate on the functionality of SCO molecular materials[41].

# DECLARATIONS

## Acknowledgments

This project has received financial support from the CNRS through the MITI interdisciplinary programs through its exploratory research program, as well as from the ANR project Mol-CoSM No. ANR-20-CE07-0028-02, the Universities of Ver-sailles and Paris-Saclay-UPSAY. We thank all of them for their strong support.

## Authors' contributions

Proposed this contribution: Linares J, Cazelles C, Boukheddaden K

Wrote all the equations for this research work: Linares J, Boukheddaden K, Cazelles C

Calculated the density of states and made the simulations: Cazelles C, Linares J

Analyzed the results: Linares J, Boukheddaden K, Cazelles C, Belmouri N, Dahoo PR

All the authors participated in writing the paper.

All authors have read and agreed to the published version of the manuscript.

## Availability of data and materials

Not applicable.

## Financial support and sponsorship

This project has received financial support from the CNRS through the MITI interdisciplinary programs through its exploratory research program, as well as from the ANR project Mol-CoSM No. ANR-20-CE07-0028-02, the Universities of Ver-sailles and Paris-Saclay-UPSAY.

## Conflicts of interest

All authors declared that there are no conflicts of interest.

## Ethical approval and consent to participate

Not applicable.

## Consent for publication

Not applicable.

## Copyright

© The Author(s) 2023.
surface and the bulk in the same compound may look unnatural, it may be possible in nanostructured core-shell nanoparticles, where the core and the shell are of different compositions. In addition, it is important to keep in mind that the systems considered are particularly small, and as the size of the system increases, the predominance of the surface, and then $J_{bs}$ interaction, vanishes, which, in turn, means that LMFA will most likely still provide correct results given the large number of molecules.

Finally, it is worth mentioning that the surface properties have been limited here to the outer layer of the lattice, while the concept of the surface may integrate some inner layers close to the interface material/air, which also have specific properties different from those of the bulk. Extensions of the present work to include the above-cited aspects are in progress, as well as other packing clusters such as rhombohedral or faced centered cubic (FCC). The results of this simulation work can also be extended to explain the behavior of polycrystalline SCO films on the $Al_2O_3$ substrate[40] and the influence of the substrate on the functionality of SCO molecular materials[41].

# DECLARATIONS

## Acknowledgments

This project has received financial support from the CNRS through the MITI interdisciplinary programs through its exploratory research program, as well as from the ANR project Mol-CoSM No. ANR-20-CE07-0028-02, the Universities of Ver-sailles and Paris-Saclay-UPSAY. We thank all of them for their strong support.

## Authors' contributions

Proposed this contribution: Linares J, Cazelles C, Boukheddaden K

Wrote all the equations for this research work: Linares J, Boukheddaden K, Cazelles C

Calculated the density of states and made the simulations: Cazelles C, Linares J

Analyzed the results: Linares J, Boukheddaden K, Cazelles C, Belmouri N, Dahoo PR

All the authors participated in writing the paper.

All authors have read and agreed to the published version of the manuscript.

## Availability of data and materials

Not applicable.

## Financial support and sponsorship

This project has received financial support from the CNRS through the MITI interdisciplinary programs through its exploratory research program, as well as from the ANR project Mol-CoSM No. ANR-20-CE07-0028-02, the Universities of Ver-sailles and Paris-Saclay-UPSAY.

## Conflicts of interest

All authors declared that there are no conflicts of interest.

## Ethical approval and consent to participate

Not applicable.

## Consent for publication

Not applicable.

## Copyright

© The Author(s) 2023.
REFERENCES

1. Gütlich P, Gaspar AB, Garcia Y. Spin state switching in iron coordination compounds. *Beilstein J Org Chem* 2013;9:342-91. DOI

2. Gütlich P, Goodwin HA. Spin crossover in transition metal compounds III. In: Topics in current chemistry. 2004. DOI

3. Garcia Y, Moscovici J, Michalowicz A, et al. A spin transition molecular material with a wide bistability domain. *Chem Eur J* 2002;8: 4992-5000. DOI

4. Garcia Y, Gütlich P. Thermal spin crossover in Mn(II), Mn(III), Cr(II) and Co(III) coordination compounds. In: Spin crossover in transition metal compounds II. Berlin, Heidelberg: Springer; 2004. p. 49-62. DOI

5. Gütlich P, Gaspar AB, Garcia Y, Ksenofontov V. Pressure effect studies in molecular magnetism. *C R Chim* 2007;10:21-36. DOI

6. Coronado E, Galán-mascarós J, Monrabal-capilla M, García-martínez J, Pardo-ibáñez P. Bistable spin-crossover nanoparticles showing magnetic thermal hysteresis near room temperature. *Adv Mater* 2007;19:1359-61. DOI

7. Hauser A. Spin-crossover materials. Properties and applications. Edited by Malcolm A. Halcrow. *Angew Chem Int Ed* 2013;52:10419. DOI

8. Boukheddaden K, Linares J, Tanasa R, Chong C. Theoretical investigations on an axial next nearest neighbour Ising-like model for spin crossover solids: one- and two-step spin transitions. *J Phys Condens Matter* 2007;19:106201. DOI

9. Rotaru A, Linares J, Varret F, et al. Pressure effect investigated with first-order reversal-curve method on the spin-transition compounds [FeₓZn₁₋ₓ(btr)₂(NCS)₂]·H₂O (x = 0.6, 1). *Phys Rev B* 2011;83:224107. DOI

10. Faivre C, Bellet D, Dolino G. In situ X-ray diffraction investigation of porous silicon strains induced by the freezing of a confined organic fluid. *Eur Phys J B* 2000;16:447-54. DOI

11. Rotaru A, Dírtu MM, Enachescu C, et al. Calorimetric measurements of diluted spin crossover complexes [FeₓZn₁₋ₓ(btr)₂(NCS)₂]·H₂O with MII = Zn and Ni. *Polyhedron* 2009;28:2531-6. DOI

12. Chastanet G, Gaspar AB, Real JA, Létard JF. Photo-switching spin pairs - synergy between LIESST effect and magnetic interaction in an iron(ii) binuclear spin-crossover compound. *Dedicated to the memory of Olivier Kahn. Chem Commun* 2001;819-20. DOI

13. Gütlich P, Hauser A. Thermal and light-induced spin crossover in iron(II) complexes. *Coord Chem Rev* 1990;97:1-22. DOI

14. Linares J, Jureschi C, Boukheddaden K. Surface effects leading to unusual size dependence of the thermal hysteresis behavior in spin-crossover nanoparticles. *Magnetochemistry* 2016;2:24. DOI

15. Kroeber J, Audiere J, Claude R, et al. Spin transitions and thermal hysteresis in the molecular-based materials [Fe(Htrz)2(trz)](BF4) and [Fe(Htrz)3](BF4)2.cntdot.H2O (Htrz = 1,2,4-4H-triazole; trz = 1,2,4-triazolato). *Chem Mater* 1994;6:1404-12. DOI

16. Constant-Machado H, Stancu A, Linares J, Varret F. Thermal hysteresis loops in spin-crossover compounds analyzed in terms of classical Preisach model. *IEEE Trans Magn* 1998;34:2213-9. DOI

17. Kahn O, Martinez CJ. Spin-Transition polymers: from molecular materials toward memory devices. *Science* 1998;279:44-8. DOI

18. Chiruta D, Jureschi C, Linares J, Garcia Y, Rotaru A. Lattice architecture effect on the cooperativity of spin transition coordination polymers. *J Appl Phys* 2014;115:053523. DOI

19. Enachescu C, Tanasa R, Stancu A, Codjovi E, Linares J, Varret F. FORC method applied to the thermal hysteresis of spin transition solids: first approach of static and kinetic properties. *Physica B Condensed Matter* 2004;343:15-9. DOI

20. Enachescu C, -Machado H, Menendez N, et al. Static and light induced hysteresis in spin-crossover compounds: experimental data and application of Preisach-type models. *Physica B Condensed Matter* 2001;306:155-60. DOI

21. Linares J, Codjovi E, Garcia Y. Pressure and temperature spin crossover sensors with optical detection. *Sensors* 2012;12:4479-92. DOI

22. Boukheddaden K, Ritti MH, Bouchez G, et al. Quantitative contact pressure sensor based on spin crossover mechanism for civil security applications. *J Phys Chem C* 2018;122:7597-604. DOI

23. Linares J, Allal SE, Dahoo PR, Garcia Y. Numerical simulation of a device with two spin crossover complexes: application for temperature and pressure sensors. *J Phys Conf Ser* 2017;936:012048. DOI

24. Jureschi CM, Linares J, Boulmaali A, Dahoo PR, Rotaru A, Garcia Y. Pressure and temperature sensors using two spin crossover materials. *Sensors* 2016;16:187. DOI

25. Sun L, Ndiaye M, El Islam Belmouri N, et al. Spin crossover coordination polymers with pyridine-like modification through selective guest molecules. *Cryst Growth Des* 2022;22:7555-63. DOI

26. Sun L, Belmouri NEI, Ndiaye M, et al. Thermal-driven guest-induced spin crossover behavior in 3D Fe(II)-based porous coordination polymers. *Cryst Growth Des* 2023;23:3402-11. DOI

27. Benaicha B, Van Do K, Yangui A, et al. Interplay between spin-crossover and luminescence in a multifunctional single crystal iron(ii) complex: towards a new generation of molecular sensors. *Chem Sci* 2019;10:6791-8. DOI

28. Titos-Padilla S, Herrera JM, Chen X, Delgado JJ, Colacio E. Bifunctional hybrid SiO2 nanoparticles showing synergy between core spin crossover and shell luminescence properties. *Angew Chem Int Ed Engl* 2011;123:3348-51. DOI

29. Palluel M, Tran NM, Daro N, et al. The interplay between surface plasmon resonance and switching properties in Gold@Spin crossover nanocomposites. *Adv Funct Mater* 2020;30:2000447. DOI

30. Wajnflasz J, Pick R. Transitions "low spin" - "high spin" dans les complexes de Fe²⁺. *J Phys Colloques* 1971;32:C1-91-2. DOI

31. Bousseksou A, Nasser J, Linares J, Boukheddaden K, Varret F. Ising-like model for the two-step spin-crossover. *J Phys I France* 1992;2:1381-403. DOI

32. Linares J, Spiering H, Varret F. Analytical solution of 1D Ising-like systems modified by weak long range interaction: application to spin crossover compounds. *Eur Phys J B* 1999;10:271-5. DOI

33. Varret F, Salunke SA, Boukheddaden K, et al. The Ising-like model applied to switchable inorganic solids: discussion of the static properties. *C R Chim* 2003;6:385-93. DOI
REFERENCES

1. Gütlich P, Gaspar AB, Garcia Y. Spin state switching in iron coordination compounds. *Beilstein J Org Chem* 2013;9:342-91. DOI

2. Gütlich P, Goodwin HA. Spin crossover in transition metal compounds III. In: Topics in current chemistry. 2004. DOI

3. Garcia Y, Moscovici J, Michalowicz A, et al. A spin transition molecular material with a wide bistability domain. *Chem Eur J* 2002;8: 4992-5000. DOI

4. Garcia Y, Gütlich P. Thermal spin crossover in Mn(II), Mn(III), Cr(II) and Co(III) coordination compounds. In: Spin crossover in transition metal compounds II. Berlin, Heidelberg: Springer; 2004. p. 49-62. DOI

5. Gütlich P, Gaspar AB, Garcia Y, Ksenofontov V. Pressure effect studies in molecular magnetism. *C R Chim* 2007;10:21-36. DOI

6. Coronado E, Galán-mascarós J, Monrabal-capilla M, García-martínez J, Pardo-ibáñez P. Bistable spin-crossover nanoparticles showing magnetic thermal hysteresis near room temperature. *Adv Mater* 2007;19:1359-61. DOI

7. Hauser A. Spin-crossover materials. Properties and applications. Edited by Malcolm A. Halcrow. *Angew Chem Int Ed* 2013;52:10419. DOI

8. Boukheddaden K, Linares J, Tanasa R, Chong C. Theoretical investigations on an axial next nearest neighbour Ising-like model for spin crossover solids: one- and two-step spin transitions. *J Phys Condens Matter* 2007;19:106201. DOI

9. Rotaru A, Linares J, Varret F, et al. Pressure effect investigated with first-order reversal-curve method on the spin-transition compounds [FeₓZn₁₋ₓ(btr)₂(NCS)₂]·H₂O (x = 0.6, 1). *Phys Rev B* 2011;83:224107. DOI

10. Faivre C, Bellet D, Dolino G. In situ X-ray diffraction investigation of porous silicon strains induced by the freezing of a confined organic fluid. *Eur Phys J B* 2000;16:447-54. DOI

11. Rotaru A, Dírtu MM, Enachescu C, et al. Calorimetric measurements of diluted spin crossover complexes [FeₓZn₁₋ₓ(btr)₂(NCS)₂]·H₂O with MII = Zn and Ni. *Polyhedron* 2009;28:2531-6. DOI

12. Chastanet G, Gaspar AB, Real JA, Létard JF. Photo-switching spin pairs - synergy between LIESST effect and magnetic interaction in an iron(ii) binuclear spin-crossover compound. *Dedicated to the memory of Olivier Kahn. Chem Commun* 2001;819-20. DOI

13. Gütlich P, Hauser A. Thermal and light-induced spin crossover in iron(II) complexes. *Coord Chem Rev* 1990;97:1-22. DOI

14. Linares J, Jureschi C, Boukheddaden K. Surface effects leading to unusual size dependence of the thermal hysteresis behavior in spin-crossover nanoparticles. *Magnetochemistry* 2016;2:24. DOI

15. Kroeber J, Audiere J, Claude R, et al. Spin transitions and thermal hysteresis in the molecular-based materials [Fe(Htrz)2(trz)](BF4) and [Fe(Htrz)3](BF4)2.cntdot.H2O (Htrz = 1,2,4-4H-triazole; trz = 1,2,4-triazolato). *Chem Mater* 1994;6:1404-12. DOI

16. Constant-Machado H, Stancu A, Linares J, Varret F. Thermal hysteresis loops in spin-crossover compounds analyzed in terms of classical Preisach model. *IEEE Trans Magn* 1998;34:2213-9. DOI

17. Kahn O, Martinez CJ. Spin-Transition polymers: from molecular materials toward memory devices. *Science* 1998;279:44-8. DOI

18. Chiruta D, Jureschi C, Linares J, Garcia Y, Rotaru A. Lattice architecture effect on the cooperativity of spin transition coordination polymers. *J Appl Phys* 2014;115:053523. DOI

19. Enachescu C, Tanasa R, Stancu A, Codjovi E, Linares J, Varret F. FORC method applied to the thermal hysteresis of spin transition solids: first approach of static and kinetic properties. *Physica B Condensed Matter* 2004;343:15-9. DOI

20. Enachescu C, -Machado H, Menendez N, et al. Static and light induced hysteresis in spin-crossover compounds: experimental data and application of Preisach-type models. *Physica B Condensed Matter* 2001;306:155-60. DOI

21. Linares J, Codjovi E, Garcia Y. Pressure and temperature spin crossover sensors with optical detection. *Sensors* 2012;12:4479-92. DOI

22. Boukheddaden K, Ritti MH, Bouchez G, et al. Quantitative contact pressure sensor based on spin crossover mechanism for civil security applications. *J Phys Chem C* 2018;122:7597-604. DOI

23. Linares J, Allal SE, Dahoo PR, Garcia Y. Numerical simulation of a device with two spin crossover complexes: application for temperature and pressure sensors. *J Phys Conf Ser* 2017;936:012048. DOI

24. Jureschi CM, Linares J, Boulmaali A, Dahoo PR, Rotaru A, Garcia Y. Pressure and temperature sensors using two spin crossover materials. *Sensors* 2016;16:187. DOI

25. Sun L, Ndiaye M, El Islam Belmouri N, et al. Spin crossover coordination polymers with pyridine-like modification through selective guest molecules. *Cryst Growth Des* 2022;22:7555-63. DOI

26. Sun L, Belmouri NEI, Ndiaye M, et al. Thermal-driven guest-induced spin crossover behavior in 3D Fe(II)-based porous coordination polymers. *Cryst Growth Des* 2023;23:3402-11. DOI

27. Benaicha B, Van Do K, Yangui A, et al. Interplay between spin-crossover and luminescence in a multifunctional single crystal iron(ii) complex: towards a new generation of molecular sensors. *Chem Sci* 2019;10:6791-8. DOI

28. Titos-Padilla S, Herrera JM, Chen X, Delgado JJ, Colacio E. Bifunctional hybrid SiO2 nanoparticles showing synergy between core spin crossover and shell luminescence properties. *Angew Chem Int Ed Engl* 2011;123:3348-51. DOI

29. Palluel M, Tran NM, Daro N, et al. The interplay between surface plasmon resonance and switching properties in Gold@Spin crossover nanocomposites. *Adv Funct Mater* 2020;30:2000447. DOI

30. Wajnflasz J, Pick R. Transitions "low spin" - "high spin" dans les complexes de Fe²⁺. *J Phys Colloques* 1971;32:C1-91-2. DOI

31. Bousseksou A, Nasser J, Linares J, Boukheddaden K, Varret F. Ising-like model for the two-step spin-crossover. *J Phys I France* 1992;2:1381-403. DOI

32. Linares J, Spiering H, Varret F. Analytical solution of 1D Ising-like systems modified by weak long range interaction: application to spin crossover compounds. *Eur Phys J B* 1999;10:271-5. DOI

33. Varret F, Salunke SA, Boukheddaden K, et al. The Ising-like model applied to switchable inorganic solids: discussion of the static properties. *C R Chim* 2003;6:385-93. DOI
34. Chiruta D, Linares J, Dimian M, Garcia Y. Size effect and role of short- and long-range interactions on 1D spin-crossover systems within the framework of an ising-like model. *Eur J Inorg Chem* 2013;2013:951-7. DOI

35. Ekanayaka TK, Kurz H, Dale AS, et al. Probing the unpaired Fe spins across the spin crossover of a coordination polymer. *Mater Adv* 2021;2:760-8. DOI

36. Nishino M, Enachescu C, Miyashita S, Boukheddaden K, Varret F. Intrinsic effects of the boundary condition on switching processes in effective long-range interactions originating from local structural change. *Phys Rev B* 2010;82:020409. DOI

37. Slimani A, Boukheddaden K, Varret F, Oubouchou H, Nishino M, Miyashita S. Microscopic spin-distortion model for switchable molecular solids: spatiotemporal study of the deformation field and local stress at the thermal spin transition. *Phys Rev B* 2013;87:014111. DOI

38. Nishino M, Singh Y, Boukheddaden K, Miyashita S. Tutorial on elastic interaction models for multistep spin-crossover transitions. *J Appl Phys* 2021;130:141102. DOI

39. Singh Y, Oubouchou H, Nishino M, Miyashita S, Boukheddaden K. Elastic-frustration-driven unusual magnetoelastic properties in a switchable core-shell spin-crossover nanostructure. *Phys Rev B* 2020;101:054105. DOI

40. Jiang X, Hao G, Wang X, et al. Tunable spin-state bistability in a spin crossover molecular complex. *J Phys Condens Matter* 2019;31:315401. DOI

41. Yazdani S, Phillips J, Ekanayaka TK, Cheng R, Dowben PA. The influence of the substrate on the functionality of spin crossover molecular materials. *Molecules* 2023;28:3735. DOI

42. Pauli W. Zur Frage der Zuordnung der Komplexstrukturterme in starken und in schwachen äußeren Feldern. *Z Physik* 1923;20:371-87. DOI

43. Ising E. Beitrag zur theorie des ferromagnetismus. *Z Physik* 1925;31:253-8. DOI

44. Onsager L. Crystal Statistics. I. A two-dimensional model with an order-disorder transition. *Phys Rev* 1944;65:117-49. DOI

45. Néel L. Magnetism and local molecular field. *Science* 1971;174:985-92. DOI

46. Katsura S, Takizawa M. Bethe lattice and the bethe approximation. *Prog Theor Phys* 1974;51:82-98. DOI

47. Constant-machado H, Linares J, Varret F, et al. Dilution effects in a spin crossover system, modelled in terms of direct and indirect intermolecular interactions. *J Phys I France* 1996;6:1203-16. DOI

48. Zhao Q, Xue JP, Liu ZK, Yao ZS, Tao J. Spin-crossover iron(ii) long-chain complex with slow spin equilibrium at low temperatures. *Dalton Trans* 2021;50:11106-12. DOI

49. Fürmeyer F, Carrella LM, Ksenofontov V, Möller A, Rentschler E. Phase trapping in multistep spin crossover compound. *Inorg Chem* 2020;59:2843-52. DOI

50. Shteto I, Linares J, Varret F. Monte Carlo entropic sampling for the study of metastable states and relaxation paths. *Phys Rev E* 1997;56:5128-37. DOI

51. Linares J, Enachescu C, Boukheddaden K, Varret F. Monte Carlo entropic sampling applied to spin crossover solids: the squareness of the thermal hysteresis loop. *Polyhedron* 2003;22:2453-6. DOI

52. Chiruta D, Linares J, Dahoo PR, Dimian M. Analysis of long-range interaction effects on phase transitions in two-step spin-crossover chains by using Ising-type systems and Monte Carlo entropic sampling technique. *J Appl Phys* 2012;112:074906. DOI

53. Cazelles C, Linares J, Ndiaye M, Dahoo P, Boukheddaden K. Hexagonal-shaped spin crossover nanoparticles studied by ising-like model solved by local mean field approximation. *Magnetochemistry* 2021;7:69. DOI

54. Pillet S, Hubsch J, Lecomte C. Single crystal diffraction analysis of the thermal spin conversion in [Fe(btr)₂(NCS)₂]·H₂O: evidence for spin-like domain formation. *Eur Phys J B* 2004;38:541-52. DOI
34. Chiruta D, Linares J, Dimian M, Garcia Y. Size effect and role of short- and long-range interactions on 1D spin-crossover systems within the framework of an ising-like model. *Eur J Inorg Chem* 2013;2013:951-7. DOI

35. Ekanayaka TK, Kurz H, Dale AS, et al. Probing the unpaired Fe spins across the spin crossover of a coordination polymer. *Mater Adv* 2021;2:760-8. DOI

36. Nishino M, Enachescu C, Miyashita S, Boukheddaden K, Varret F. Intrinsic effects of the boundary condition on switching processes in effective long-range interactions originating from local structural change. *Phys Rev B* 2010;82:020409. DOI

37. Slimani A, Boukheddaden K, Varret F, Oubouchou H, Nishino M, Miyashita S. Microscopic spin-distortion model for switchable molecular solids: spatiotemporal study of the deformation field and local stress at the thermal spin transition. *Phys Rev B* 2013;87:014111. DOI

38. Nishino M, Singh Y, Boukheddaden K, Miyashita S. Tutorial on elastic interaction models for multistep spin-crossover transitions. *J Appl Phys* 2021;130:141102. DOI

39. Singh Y, Oubouchou H, Nishino M, Miyashita S, Boukheddaden K. Elastic-frustration-driven unusual magnetoelastic properties in a switchable core-shell spin-crossover nanostructure. *Phys Rev B* 2020;101:054105. DOI

40. Jiang X, Hao G, Wang X, et al. Tunable spin-state bistability in a spin crossover molecular complex. *J Phys Condens Matter* 2019;31:315401. DOI

41. Yazdani S, Phillips J, Ekanayaka TK, Cheng R, Dowben PA. The influence of the substrate on the functionality of spin crossover molecular materials. *Molecules* 2023;28:3735. DOI

42. Pauli W. Zur Frage der Zuordnung der Komplexstrukturterme in starken und in schwachen äußeren Feldern. *Z Physik* 1923;20:371-87. DOI

43. Ising E. Beitrag zur theorie des ferromagnetismus. *Z Physik* 1925;31:253-8. DOI

44. Onsager L. Crystal Statistics. I. A two-dimensional model with an order-disorder transition. *Phys Rev* 1944;65:117-49. DOI

45. Néel L. Magnetism and local molecular field. *Science* 1971;174:985-92. DOI

46. Katsura S, Takizawa M. Bethe lattice and the bethe approximation. *Prog Theor Phys* 1974;51:82-98. DOI

47. Constant-machado H, Linares J, Varret F, et al. Dilution effects in a spin crossover system, modelled in terms of direct and indirect intermolecular interactions. *J Phys I France* 1996;6:1203-16. DOI

48. Zhao Q, Xue JP, Liu ZK, Yao ZS, Tao J. Spin-crossover iron(ii) long-chain complex with slow spin equilibrium at low temperatures. *Dalton Trans* 2021;50:11106-12. DOI

49. Fürmeyer F, Carrella LM, Ksenofontov V, Möller A, Rentschler E. Phase trapping in multistep spin crossover compound. *Inorg Chem* 2020;59:2843-52. DOI

50. Shteto I, Linares J, Varret F. Monte Carlo entropic sampling for the study of metastable states and relaxation paths. *Phys Rev E* 1997;56:5128-37. DOI

51. Linares J, Enachescu C, Boukheddaden K, Varret F. Monte Carlo entropic sampling applied to spin crossover solids: the squareness of the thermal hysteresis loop. *Polyhedron* 2003;22:2453-6. DOI

52. Chiruta D, Linares J, Dahoo PR, Dimian M. Analysis of long-range interaction effects on phase transitions in two-step spin-crossover chains by using Ising-type systems and Monte Carlo entropic sampling technique. *J Appl Phys* 2012;112:074906. DOI

53. Cazelles C, Linares J, Ndiaye M, Dahoo P, Boukheddaden K. Hexagonal-shaped spin crossover nanoparticles studied by ising-like model solved by local mean field approximation. *Magnetochemistry* 2021;7:69. DOI

54. Pillet S, Hubsch J, Lecomte C. Single crystal diffraction analysis of the thermal spin conversion in [Fe(btr)₂(NCS)₂]·H₂O: evidence for spin-like domain formation. *Eur Phys J B* 2004;38:541-52. DOI