## 1 Introduction

Spin-crossover (SCO) compounds are transition metal complexes in thermodynamic competition between the low-spin (LS) state, stable at low temperature, and the high-spin (HS) state, stable at higher temperature. These states have not only different magnetic, optical and vibrational properties but also different molecular volumes, with a slightly larger volume for the HS state. The change in volume during the transition results in the building-up of stress in the crystal and therefore is at the origin of elastic interactions and more generally at the origin of all elastic collective behaviours called cooperativity, which, if strong enough, are responsible at the macroscopic scale for the first-order transition accompanied by bistability and hysteresis phenomena. At the microscopic scale, elastic interactions induce nucleation, cluster formation and growth of like-spin domains, that is, domains of adjacent molecules in the same spin state.

The simplest approaches for the study of the macroscopic behaviour of SCO compounds are mean-field models, which imply random distribution of HS and LS molecules during the transition. Mean-field models have successfully reproduced a number of experimental curves (thermal transition, relaxations, and photoexcitation), mainly in the case of less cooperative diluted SCO compounds. These models, which remain widely used, represent a valuable tool in SCO research. However, these models do not include short-range interactions and therefore do not allow an appropriate analysis of the propagation of molecular switching inside the sample. For a deeper microscopic description of nucleation and domain growth processes, different models should be used. One of the simplest models introducing short-range interactions is the Ising model, adapted for the study of SCO compounds by replacing the classical magnetic field with a temperature-dependent field. The Ising model, reviewed in the article of Pavlik and Linares in this special issue, has been used to successfully treat both the dynamic and static properties [1,2] of SCO compounds in the framework of Metropolis or Arrhenius approaches with only short-range interactions or a combination of short-range and long-range interactions. Although Ising-like models can reproduce to some extent both macroscopic and microscopic behaviours, they describe only partially the relations between elastic interactions and cooperativity due to a lack of a simultaneous treatment of local distortions and global deformations. Indeed, the interactions between neighbouring sites are fixed and the lattice cannot be deformed. More complex models aiming to improve this deficiency have been developed. A first model based on the ligand field theory discussing the effect of coupling between d-electrons and lattice strains has been proposed by Kambara [3], following the phenomenological description of Slichter and Drickamer [4]. In his lattice expansion model using the Eshelby defect theory based on continuum medium mechanics considerations [5], Spiering et al. [6] emphasized that elastic interactions arise from stresses because of changes in the volume, shape and elasticity of the lattice during the transition, but still treated the system in a mean-field approach. In addition, Kambara was the first to introduce intramolecular distortions as a dynamic variable and considered the intermolecular coupling in his model [7]. Meanwhile, Zimmermann and Konig [8] took into account the spin–orbit coupling and the effect of low-symmetry ligand fields and considered the volume expansion of the elastic medium, which would induce an additional contribution to the lattice energy. A decade later, it was stated that the interaction constant can be explained by the elasticity theory in a proportion of around 80% [9]. These “realistic” models, despites their huge interest, have not been extensively used by the SCO community for the complexity of calculations and the large number (mostly unknown) material parameters.

An intermediate model between the Ising-like model and elastic defect theory is the atom–phonon coupling model [10,11], in which the SCO compounds are modelled in a similar way as in the Ising-like model, but neighbouring molecules are linked by springs, whose elastic constant can take three different values according to the spin states of the closest molecules (HS–HS, HS–LS or LS–LS).

New experimental results in the last decade called for the development of more advanced and sophisticated models. After years of theoretical suppositions and partial experimental proofs [12–14], a first direct answer to the fundamental question, “does the heterogeneous phase separation exist and does it play a role in the spin crossover phenomena?” [15], has been provided by Pillet et al. [16] by diffraction analysis of high-quality single crystals, constituting an indubitable proof of the existence of a spin-like domain formation. In a pioneering work, Bonnet et al. [17], followed by several articles of Versailles group [18,19] and Toulouse group [20,21] (for an overview see also Ref. [22]), have visualized the phase separation in a single crystal using optical microscopy. The various optical microscopy results showed that the process always starts from a corner of the crystal, and the dependence of the evolution of subsequent clusters on the size, shape and thermal history of the crystal has been carefully analysed. Further experimental works have investigated the velocity of propagation of the interface between the two phases [23,24]. The role of lattice mismatches (elastic stress and strain) [24] and latent heat in this process was also demonstrated [25]. Other recent experiments showing the building-up of an elastic step after photoexcitation using a femtosecond laser have also been explained in terms of the propagation of elastic distortions [26].

To correctly describe the microscopic mechanism of SCO phenomena and to reproduce the new experimental results, the logical development of models implied the replacement of ambiguous short- and long-range interactions with genuine elastic interactions arising from lattice distortions because of the molecular size difference between the LS and HS states. This led to a new class of models—so-called “ball and spring models”, which consider the molecules as rigid and usually unfixed spheres linked by springs that can be compressed or elongated during the transition resulting in practically an infinite number of possible interaction terms between molecules. In the pioneering years, at the end of the first decade of the new millennium, different ball and spring models have been independently and almost simultaneously proposed by a few groups of researchers working in the SCO field [27–33]. Even if all these models are based on the same idea, formal differences can be detected in the form of Hamiltonians used in the theoretical approaches, in the internal structure and shape of the samples, as well as in the methods used in the simulations. The subsequent period has been devoted to applications and development of models for reproducing various experimental behaviours of SCO compounds, including new emerging phenomena observed in the case of micro- and nanoparticles. The present review is organized as follows: first, we present the most important ball and spring models, emphasizing on the similarities and differences in between. Then we discuss their applications with the accent on size reduction effects.

## 2 A journey to various types of elastic models

### 2.1 Hamiltonians

The starting point for all elastic Hamiltonians is the simple Ising-like Hamiltonian [1,34], which is discussed in more details in the review of Pavlik and Linares [35]. By associating a fictitious spin with each molecule (${S}_{i}=\pm 1$ for HS (LS) state), the short-range interactions $J$ and long-range interactions $G$ can be written as exchange interactions in the Hamiltonian:

$$\text{H}=\frac{1}{2}\sum _{i}\left(D-{k}_{\text{B}}T\phantom{\rule{0.25em}{0ex}}\mathrm{ln}\phantom{\rule{0.25em}{0ex}}g\right){S}_{i}-J\sum _{i,j}{S}_{i}{S}_{j}-G\sum _{i}{S}_{i}\langle S\rangle $$ | (1) |

$${n}_{\text{HS}}=\frac{1+\langle S\rangle}{2}$$ | (2) |

In the elastic models, the interaction terms are replaced in the Hamiltonian with an intermolecular elastic potential [27,31–33,36–48]:

$$\text{H}=\frac{1}{2}\sum _{i}\left(D-{k}_{\text{B}}T\phantom{\rule{0.25em}{0ex}}\mathrm{ln}\phantom{\rule{0.25em}{0ex}}g\right){S}_{i}+\sum _{\langle i,j\rangle}{V}_{ij}^{\text{inter}}$$ | (3) |

The most general expression of the elastic potential can be written as

$$\sum _{\langle i,j\rangle}{V}_{ij}^{\text{inter}}=\sum _{i,j}\frac{k}{2}\mathrm{\delta}{x}_{ij}^{2}$$ | (4) |

In the case of a hexagonal lattice [31], all neighbours are equivalent and the individual spring elongations can be written as

$$\mathrm{\delta}{x}_{ij}={d}_{ij}-({l}_{0ij}+{r}_{i}+{r}_{j})$$ | (5) |

$$\sum _{\langle i,j\rangle}{V}_{ij}^{\text{inter}}=\sum _{i,j}\frac{k}{2}{\left[{d}_{ij}-({l}_{0ij}+{r}_{i}+{r}_{j})\right]}^{2}$$ | (6) |

In the case of a rectangular lattice, to maintain the topology of the bond network and to ensure the mechanical stability with respect to the shear distortions, the presence of second neighbour interactions is needed and consequently the potential will be written as

$$\sum _{\langle i,j\rangle}{V}_{ij}^{\text{inter}}=\sum _{i,j}\frac{{k}_{1}}{2}{\left[{d}_{ij}-({l}_{0ij}+{r}_{i}+{r}_{j})\right]}^{2}+\sum _{i,k}\frac{{k}_{2}}{2}{\left[{d}_{ik}-({l}_{0ik}+{r}_{i}+{r}_{k})\right]}^{2}$$ | (7) |

Simple geometrical considerations give the following relation between ${l}_{0ij}$ and ${l}_{0ik}$ (see Fig. 1):

$${l}_{0ik}=\sqrt{2}{l}_{0ij}+\left(\sqrt{2}-1\right)\left({r}_{i}+{r}_{k}\right)$$ | (8) |

A first approach of the interaction potential assumed the closest molecular spheres in contact at equilibrium and therefore ${l}_{0ij}=0$. In this case, taking into account relation (8), the elastic potential will be written in the following form:

$$\sum _{\langle i,j\rangle}{V}_{ij}^{\text{inter}}=\frac{{k}_{1}}{2}{\sum _{i,j}\left[{d}_{ij}-({r}_{i}+{r}_{j})\right]}^{2}+\frac{{k}_{2}}{2}\sum _{i,k}{\left[{d}_{ik}-\sqrt{2}({r}_{i}+{r}_{k})\right]}^{2}$$ | (9) |

This expression of the elastic potential has been used in the original elastic models in SCO compounds [27,28,37,38] (the one-dimensional [1D] case has been proposed in 2001 by J. Nasser in a study [50] who now can be interpreted as a precursor of the current elastic models).

An alternative approach to keep the rectangular structure during simulations without using the second-order neighbour interactions was more recently proposed in the so-called “stretching and bending model” [49]. It considers an angular harmonic potential corresponding to the angle $\theta $ between two pairs of molecules sharing the same vortex (see Fig. 1). In this case, the second term of potential (7) is replaced with a bending oscillator:

$$\sum _{\langle i,j,k\rangle}{V}_{ijk}=\sum _{i,j,k}\frac{{k}_{2}}{2}{\left[{\theta}_{ijk}-{\theta}_{0}\right]}^{2}$$ |

Another expression for the intermolecular potential ${V}_{ij}^{\text{inter}}$ implies the use of exponential functions for nearest neighbours in the following form:

$$\sum _{\langle i,j\rangle}{V}_{ij}^{\text{inter}}=\sum _{i,j}C\left[\mathrm{exp}\left(\mathrm{\alpha}(\mathrm{\delta}{x}_{ij}-{u}_{0})\right)+\mathrm{exp}(-\mathrm{\beta}(\mathrm{\delta}{x}_{ij}-{u}_{0}))\right]$$ |

Considering the equilibrium distances between the mass centres of $\langle i,j\rangle $ molecules

$${R}_{0ij}={l}_{0ij}+{r}_{i}+{r}_{j}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}{R}_{0ik}={l}_{0ik}+{r}_{i}+{r}_{j}$$ | (10) |

$$\sum _{\langle i,j\rangle}{V}_{ij}^{\text{inter}}=\frac{{k}_{1}}{2}\sum _{i,j}{\left({d}_{ij}-{R}_{0ij}\right)}^{2}+\frac{{k}_{2}}{2}\sum _{i,k}{\left({d}_{ik}-{R}_{0ik}\right)}^{2}$$ | (11) |

Then, the total Hamiltonian can be expressed in terms of an Ising-like Hamiltonian (1) as

$$\text{H}=\sum _{i}\left(\frac{1}{2}\left(D-{k}_{\text{B}}T\phantom{\rule{0.25em}{0ex}}\mathrm{ln}\phantom{\rule{0.25em}{0ex}}g\right)-k{\mathrm{\rho}}_{1}\sum _{j=1}^{z}\left({d}_{ij}-{\mathrm{\rho}}_{0}\right)\right){S}_{i}+k{\mathrm{\rho}}_{1}^{2}\sum _{\langle i,j\rangle}{S}_{i}{S}_{j}+\frac{k}{2}\sum _{\langle i,j\rangle}{\left({d}_{ij}-{\mathrm{\rho}}_{0}\right)}^{2}+2kN{\mathrm{\rho}}_{1}^{2}$$ | (12) |

This Hamiltonian in the above form is the basis of the so-called electroelastic approach for SCO complexes [40–42].

In the case of the so-called anharmonic model, the values of ${k}_{1}$ and ${k}_{2}$ depend on the equilibrium intermolecular distances (as the bulk modulus in experimental data) and can be expressed in the following form:

$${k}_{1}={k}_{01}+{k}_{11}{\left[{d}_{ij}-\left({l}_{0ij}+2{r}_{\text{HS}}\right)\right]}^{2}$$ | (13) |

$${k}_{2}\left({d}_{ik}\right)={k}_{02}+{k}_{12}{\left[{d}_{ik}-\sqrt{2}\left({l}_{0ik}+2{r}_{\text{HS}}\right)\right]}^{2}$$ | (14) |

The anharmonicity has been introduced by considering the Lennard-Jones intermolecular interaction potential whose equilibrium distances and minimum energies at zero temperature can be dependent on the spin state of the molecules forming the bonds [32]. This model allows a better control of thermal expansion phenomena, which is known to be different in the two molecular phases, due to different elastic and vibrational properties. The so-called two-variable anharmonic model has been used for the simulation of X-ray diffraction experiments, reproducing the crystallographic phase separation (Bragg peak splitting) observed during the thermal spin transition or during nonequilibrium phase transformation kinetics like photoexcitations or thermal quenching process [33]. In a similar way as the electroelastic model, the anharmonic model has been rewritten as an Ising-like model [48]:

$$\text{H}=\sum _{\langle \mathit{i},\mathit{j}\rangle}\left[\mathit{J}\left({\mathit{d}}_{\mathit{ij}}\right){\mathit{S}}_{i\phantom{\rule{0.25em}{0ex}}}{\mathit{S}}_{j\phantom{\rule{0.25em}{0ex}}}+\mathit{h}\left({\mathit{d}}_{\mathit{ij}},\mathit{T}\right)\left({\mathit{S}}_{i\phantom{\rule{0.25em}{0ex}}}+{\mathit{S}}_{\mathit{j}}\right)+\mathit{D}\left({\mathit{d}}_{\mathit{ij}}\right)\right]$$ | (15) |

_{ij})S

_{i}S

_{j}and h(d

_{ij})(S

_{i}+ S

_{j}). These three quantities can be expressed as a function of the Lennard-Jones interaction potentials of the three possible bonds ${\mathit{A}}_{ij}{\mathit{V}}_{ij}^{\text{inter}}$ with ij = HH, LL, HL (or LH) (Fig. 2):

$$\mathit{J}\left({\mathit{d}}_{\mathit{ij}}\right)=\frac{{\mathit{A}}_{\text{HH}}{\mathit{V}}_{\text{HH}}^{\text{inter}}+\phantom{\rule{0.25em}{0ex}}{\mathit{A}}_{\text{LL}}{\mathit{V}}_{\text{LL}}^{\text{inter}}-2{\mathit{A}}_{\text{HL}}{\mathit{V}}_{\text{HL}}^{\text{inter}}\phantom{\rule{0.25em}{0ex}}}{4}$$ | (16) |

$$\mathit{h}\left({\mathit{d}}_{\mathit{ij}}\right)=\frac{{\mathit{A}}_{\text{HH}}{\mathit{V}}_{\text{HH}}^{\text{inter}}-\phantom{\rule{0.25em}{0ex}}{\mathit{A}}_{\text{LL}}{\mathit{V}}_{\text{LL}}^{\text{inter}}}{4}+\frac{\mathit{D}-{\mathit{k}}_{\text{B}}\mathit{T}\phantom{\rule{0.25em}{0ex}}\text{ln}\left(\mathit{g}\right)}{4}$$ | (17) |

$$\mathit{D}\left({\mathit{d}}_{\mathit{ij}}\right)=\frac{{\mathit{A}}_{\text{HH}}{\mathit{V}}_{\text{HH}}^{\text{inter}}+\phantom{\rule{0.25em}{0ex}}{\mathit{A}}_{\text{LL}}{\mathit{V}}_{\text{LL}}^{\text{inter}}+2{\mathit{A}}_{\text{HL}}{\mathit{V}}_{\text{HL}}^{\text{inter}}\phantom{\rule{0.25em}{0ex}}}{4}$$ | (18) |

A second-order Taylor development of the three Lennard-Jones potentials around the equilibrium intersite spacings permits to retrieve a Hamiltonian with a similar form to Eq. (12). When the minimum energies of the Lennard-Jones potentials are spin state dependent, the zero-order term corresponds to a lattice-independent Ising coupling. A main consequence is the presence of a clustering process during the spin state change in a system with periodic boundary conditions [48].

As a final note on models, let us mention two other possible terms, which have been considered in some cases. In addition to the intermolecular potential (9), an intramolecular adiabatic potential has also been considered [27], which can be written in the framework of the parabolic approximation for LS and HS molecules as follows:

$${V}^{\text{intra}}\left(x\right)=\frac{A}{2}\left(d+b{(c-x)}^{2}+a{x}^{2}-\sqrt{4{J}^{2}+{(d+b{(c-x)}^{2}-a{x}^{2})}^{2}}\right)$$ | (19) |

When elastic models are investigated by molecular dynamics methods, a kinetic term has also to be added, with the following expression:

$$K=\sum _{i=1}^{N}\frac{{p}_{i}^{2}}{2m}+\sum _{i=1}^{N}\frac{{P}_{i}^{2}}{2M}$$ | (20) |

_{i}represent the conjugate momenta of the coordinate of the centre ${x}_{i}$ and of the radius ${r}_{i}$, respectively, for the $i\text{'}\text{th}$ molecule. M and m correspond to the molecular and ligand masses, respectively. The molecular dynamics is often used together with an exponential-type potential; however, it was stated that the use of harmonic or Lennard-Jones potentials does not change the qualitative results of the elastic models.

### 2.2 Simulation procedures

The simulations using elastic approaches imply two different processes: the switch of spins and the change in molecular positions. These two processes can be separated (first one decides if a molecular switch is accepted or rejected and afterwards, if the spin changes, the new positions of molecules are found) or coupled (the decision of acceptance or rejection is taken after both the switch and the new molecular positions are determined).

During the first process, (1) one attributes to the molecule a spin state, taking into account the degeneracy ratios ${g}_{\text{LS}}/\left({g}_{\text{LS}}+{g}_{\text{HS}}\right)$ or ${g}_{\text{HS}}/\left({g}_{\text{LS}}+{g}_{\text{HS}}\right)/$, (2) the spin change is accepted or rejected using the Metropolis criterion (vide infra) or (3) it is also possible to consider that the spin changes with a probability depending on the local pressure.

The use of the Monte Carlo Metropolis standard procedure has following the steps: (1) one calculates the energy of the system in its actual state, according to the chosen Hamiltonian; (2) after considering one molecular switch (and finding or not the new molecular positions, depending on the model), one computes the new energy of the system; (3) one calculates the Metropolis probability

$$P=\mathrm{min}\phantom{\rule{0.25em}{0ex}}\left(1,\mathrm{exp}\left(-\frac{(D-{k}_{\text{B}}T\phantom{\rule{0.25em}{0ex}}\mathrm{ln}\phantom{\rule{0.25em}{0ex}}g){S}_{i}^{\text{final}}+\left({V}_{\text{elastic}}^{\text{final}}-{V}_{\text{elastic}}^{\text{initial}}\right)}{{k}_{\text{B}}T}\right)\right)$$ | (21) |

As the Metropolis procedure is less suitable for the study of the HS–LS relaxation or for other nonequilibrium phenomena, in Refs. [26,31,51,52], the molecular spin switching is decided through a kinetic Monte Carlo algorithm, taking into account the transition rates that depend on the energy barrier between the states rather than on the energy difference between the states. The following HS → LS and LS → HS switching probabilities have been defined based on the Monte Carlo Arrhenius-type dynamics and considering the local pressure influence on every molecule:

$${P}_{\text{HS}\to \text{LS}}^{i}=\frac{1}{\tau}\phantom{\rule{0.25em}{0ex}}\mathrm{exp}\left(\frac{D-{k}_{\text{B}}T\phantom{\rule{0.25em}{0ex}}\mathrm{ln}\phantom{\rule{0.25em}{0ex}}g}{2{k}_{\text{B}}T}\right)\phantom{\rule{0.25em}{0ex}}\mathrm{exp}\left(-\frac{E-\mathrm{\kappa}{p}_{i}}{{k}_{\text{B}}T}\right)$$ | (22) |

$${P}_{\text{LS}\to \text{HS}}^{i}=\frac{1}{\mathrm{\tau}}\phantom{\rule{0.25em}{0ex}}\mathrm{exp}\left(-\frac{D-{k}_{\text{B}}T\phantom{\rule{0.25em}{0ex}}\mathrm{ln}\phantom{\rule{0.25em}{0ex}}g}{2{k}_{\text{B}}T}\right)\phantom{\rule{0.25em}{0ex}}\mathrm{exp}\left(-\frac{E+\mathrm{\kappa}{p}_{i}}{{k}_{\text{B}}T}\right)$$ | (23) |

The new molecular positions in the system are found by (1) relaxing mechanically the lattice, considering small displacements on all axes for randomly selected molecules, (2) considering the equations of motion in the framework of Nosé–Hoover formalism or (3) deterministically calculating the motion of molecules, by solving exactly a system of differential equations.

The Nosé–Hoover formalism [27] implies the following equations of motion, with the Hamiltonian of the thermal reservoir given by ${\text{H}}_{\text{therm}}=({P}_{S}^{2}/2Q)+3N{k}_{\text{B}}T\phantom{\rule{0.25em}{0ex}}\mathrm{ln}\phantom{\rule{0.25em}{0ex}}s$ ($s$ is a scaling factor, ${P}_{S}$ is the conjugated momentum of $s$ and $Q$ is the effective mass), which are then added to the chosen Hamiltonian of the system:

$$\frac{\text{d}{r}_{i}}{\text{d}t}=\frac{{p}_{i}}{m},\frac{\text{d}{p}_{i}}{\text{d}t}=-\frac{\partial {V}^{\text{intra}}}{\partial {r}_{i}}-\frac{\partial {V}^{\text{inter}}}{\partial {r}_{i}}-\frac{{P}_{S}}{Q}{p}_{i},\frac{\text{d}{x}_{i}}{\text{d}t}=\frac{{p}_{i}}{M},\frac{\text{d}{P}_{i}}{\text{d}t}=-\frac{\partial {V}^{\text{inter}}}{\partial {x}_{i}}-\frac{{P}_{S}}{Q}{P}_{i},\frac{\text{d}s}{\text{d}t}=\frac{{P}_{S}}{Q}s,\frac{\text{d}}{\text{d}t}\left(\frac{{P}_{S}}{Q}\right)=\frac{1}{Q}{\left[\sum _{i}\frac{{p}_{i}^{2}}{m}+\sum _{i}\frac{{P}_{i}^{2}}{M}-3N{k}_{\text{B}}T\right]}_{i}$$ | (24) |

In the case of a deterministic calculation of the molecular positions, the motion of molecules is determined by local elastic forces and is accounted for solving the following differential equations for all molecules in the system.

$$\{\begin{array}{c}m\frac{{\text{d}}^{2}{x}_{i}}{\text{d}{t}^{2}}={F}_{xi}-\mathrm{\mu}\frac{\text{d}{x}_{i}}{\text{d}t}\\ m\frac{{\text{d}}^{2}{y}_{i}}{\text{d}{t}^{2}}={F}_{yi}-\mathrm{\mu}\frac{\text{d}{y}_{i}}{\text{d}t}\end{array}$$ | (25) |

$$\begin{array}{l}{F}_{i,x}=\sum _{\text{neighbour}\phantom{\rule{0.5em}{0ex}}\text{spring}}k\mathrm{\delta}{r}_{ij,x}\\ {F}_{i,y}=\sum _{\text{neighbour}\phantom{\rule{0.5em}{0ex}}\text{spring}}k\mathrm{\delta}{r}_{ij,y}\end{array}$$ | (26) |

A slightly different Monte Carlo numerical process has been proposed to study the nonequilibrium process in the anharmonic model [32]. Indeed, two distinct Arrhenius-type transition rates have been used for the fictitious spins and the continuous lattice variables on the basis of the well-known “one-step dynamics”. Two different characteristic scaling constants have been attributed to mimic the difference between electronic transitions at the molecular level and the lattice vibrations at a mesoscopic/macroscopic level [47].

To find the most probable shape for a cluster and of the interface line between HS and LS domains, Eden [53] and Kawasaki [54] approaches have been used. Following the genuine Eden method [53], an LS cluster is forced to grow by only allowing the flipping of molecules, which have at least one LS neighbour. As the Eden-like method provides clusters, which grow dynamically, to find the most probable shape of a cluster for a given HS fraction, Kawasaki [54,55] dynamics has been applied, keeping constant the number of HS and LS molecules in the system. For this, one chooses pairs of neighbouring HS and LS molecules and then exchanges their states. Then the Monte Carlo Metropolis procedure is applied. In this case, as the number of HS molecules is kept constant, the electronic energy does not play any role in finding the lowest energy states, and only the elastic terms in the Hamiltonian are considered.

For tuning the different relative time scales of spin dynamics and lattice relaxation, one can define the ratio $r={n}_{1}/{n}_{2}$, where ${n}_{1}$ is the number of steps to find the positions of molecules with frozen spins, and ${n}_{2}$ is the number of Monte Carlo steps for spin state switching. It should be noted that a large $r$ (large ${n}_{1}$) favours equilibrium distribution, where after each spin flipping the lattice has enough time to relax the excess energy generated by the switch, as it is the case of the thermal transition [56]. On the contrary, a smaller $r$ (high ${n}_{2}$) favours nonequilibrium distribution, like in the case of fast phenomena subsequent to femtosecond photoexcitation experiments [57].

A few warnings with respect to the Monte Carlo methods are necessary here. First, in Monte Carlo simulation, the time evolution of probability distributions and thermal average quantities are governed by a master equation in which the “real” time (i.e., from second Newton's law or Schrödinger equation) is not considered. The Monte Carlo “time” (Monte Carlo step MCs) is arbitrary defined, and therefore the velocity concept cannot be easily grasped by Monte Carlo methods and the experimental data of nonequilibrium dynamical processes cannot be fitted. Another important point concerns the temperature or pressure concepts in Monte Carlo simulations. The simulated system is in contact with an ideal thermal and pressure bath. It is also difficult to consider thermal diffusion and other transport phenomena by defining local temperature or pressure without breaking the microreversibility properties on which Monte Carlo simulations are based. In addition, the Monte Carlo simulations can be largely affected by kinetic effects: the (meta)stable states are found after long computing time, which increases with the system size. On one hand, this feature can give rise to an anomalous hysteresis for systems with weak interactions, which should produce a gradual and reversible transition for large enough waiting time. At a fast sweeping rate, a small variation in the sweep rate determines large changes in the transition temperatures, although the effect is considerably smaller in the case of lower temperature sweeping rates. In this case, the limit temperatures of the hysteresis loops approach constant values, corresponding to spinodal points in the Husimi–Temperley model with long-range interactions [58]. But, on the other hand, this feature could be useful when studying systems composed of nanoparticles or the ultrafast behaviour of SCO systems. Indeed, although for most bulk SCO compounds the intrinsic kinetics of the thermal SCO phase transition are rapid and the experimental thermal hysteresis does not visibly depend on the temperature sweep rate, this is no longer the case of SCO nanocrystals in an analogy to the superparamagnetic behaviour of nanometric or nanostructured magnetic systems. Experimental approaches by first-order reversal curve method [59–61] have shown that most important kinetic effects are accumulated in the proximity of the major hysteresis loop, and a procedure to disentangle between static and kinetic parts has been proposed based on first-order reversal curve simulations using different temperature sweeping rates [62]. Because the kinetic effects are different as a function of the particle size, the averaging of kinetic curves for particles of different sizes is necessary for a good agreement with experimental data.

Although different in appearance, all the elastic models presented above use the same main hypothesis and lead to similar conclusions, being able to accurately describe the different aspects of the spin transition using simple and realistic assumptions. In the following sections, we shall present how they reproduce some of these aspects.

## 3 Applications and links to experiments

### 3.1 Cluster formation

The long-range nature of the elastic interaction is a key factor for the role played by the system boundary in all elastic models. In periodic boundary systems, the domain growth is suppressed and the configuration is uniform, irrespective of the values of elastic constants [27]. However, the situation is by far more interesting in open boundary systems, where clustering occurs in the presence of strong interactions. Studies have been made for different shapes of bidimensional samples (squares [27,28], rectangles [39,40,63], hexagons [43], circles [44], ellipses [45], etc.) and the similarity with experiments is remarkable: like in the simulation, the cluster formation starts from one corner or edge and then spreads towards the bulk, with a contact (wetting) angle [64] close to $\mathrm{\pi}/2$ (see Fig. 3). The same conclusions have been drawn from the first attempts to study clustering phenomena in 3D thin films, where in addition buckling effects were observed [65], and from very recent approaches of pure 3D systems with different shapes (cubic, tetrahedral, octahedral or spherical) in a face-centred cubic configuration [66]. A simple explanation of the spreading of clusters from corners is provided by Fig. 4, where the total elastic energy of the system in the case of clusters of increasing size is plotted, when grown for corner, edge or bulk. The energy is calculated after the systems have relaxed to their minimum energy state and is represented as a function of the ratio between the cluster and the system size, that is, in fact, the LS fraction. It is obvious that the lowest elastic energy is obtained for a cluster formed in the vicinity of a corner, whereas an almost 50% energy increase is observed for a cluster starting from the edge and by several times energy increase in the case of a cluster starting from bulk. Therefore, a hypothetical cluster formed inside the bulk will not be able to spread, at least as long as the lattice symmetry is maintained and there are no defects inside the crystal. In addition, the clusters spread from any part of the surface only at very low temperatures, when the relatively high local energy increase is compensated by the large value of $D-{k}_{\text{B}}T\phantom{\rule{0.25em}{0ex}}\mathrm{ln}\phantom{\rule{0.25em}{0ex}}g$ term in Hamiltonian (3). Simulations have shown that there is a high probability that a single cluster spreads in the case of small systems with high interactions (as the development of a cluster prevents the spreading of other clusters), whereas several clusters may develop in the case of larger systems and moderate interactions [45]. The role of interstitial defects and vacancies in the cluster formation has also been investigated [67,68].

The cluster spreading in such elastic systems has been discussed in the framework of the “macroscopic nucleation” [38] concept, which can be extended for other systems exhibiting long-range interaction and which implies that the size of the critical nucleus is proportional to the system size. During the propagation process, the existence of three regimes has been identified: the incubation, corresponding to the initial formation of a cluster; the flow, corresponding to a linear velocity of the cluster; and the acceleration, corresponding to the last stages of relaxation when the opposite surface is close. The study of the speed dependence on the elastic constant has shown that the speed is lower for more rigid lattices [40]. This can be correlated to the larger velocity experimentally observed in many cases for the heating process, with respect to the cooling process [23]. In addition, it was stated [23] that the interface velocity is a monotonous function of temperature, passing through zero at the equilibrium temperature and if the temperature is constant its evolution cannot be reversed. One should distinguish between the interface propagation and the propagation of elastic interactions, as it was experimentally shown that the domain propagation is much slower than the propagation of sound waves [23]. This feature is particularly important for ultrafast phenomena, where the volume variation does not exactly follow the variation in the ${n}_{\text{HS}}$ fraction, due to different times needed for cluster spreading and for the propagation of elastic distortions towards the edges of the system [26,52].

Excepting a small part distributed at the surface of the lattice, the elastic energy at equilibrium, as well as the higher local pressures, is located at the interface between the two phases [24,40,47]. However, the value of this elastic energy almost vanishes for specific orientations of the border lines, as it was demonstrated by an anisotropic extension of the electroelastic model, which took into account the atomic structure of the crystal via a simple one-site model [24]. Comparable evidence results from the analysis of the maps of the local pressures, which correspond to local deformation fields in a continuous medium. Similar to higher local pressures, higher strain can be detected at the interface between the two phases (see Fig. 5). Therefore, for a fixed size of a cluster, the most stable interface shape is that with the shortest length, and both electroelastic and mechanoelastic models confirmed that the contact (wetting) angle between the interface and the lattice border must be around $\mathrm{\pi}/2$. For example, in the case of a circular lattice, the nucleation starts at the edge of the lattice with a circular shape, which gradually changes through the propagation, becoming a straight line at the middle of the region, whereas the interface length is maximal and returning to a circular shape towards the end of the process with a decrease in the interface length.

### 3.2 Finite size effects

#### 3.2.1 Introduction

In the last decade, an intense research activity has focused on the investigation of the size reduction effect on the physicochemical properties of all classes of switchable molecular materials or coordination networks, with a special attention to SCO systems [69] and Prussian blue analogue (PBA) coordination polymers [70]. In particular, the understanding of the phase transition mechanism at the nanometric scale and the control of the phase stability is of paramount importance for their future potential applications in memory and display devices as molecular switches or in SCO strain-based systems [71].

It is important to note that strictly speaking, a phase transition only exists in the thermodynamic limit, that is, when the number of degrees of freedom and the volume are infinite, their ratio remaining constant. Phase transitions are because of the existence of latent heat or singularities in the Gibbs energy, which produce divergences or discontinuities of their derivatives (order parameters, susceptibilities, etc.). Finite size effects remove such singularities, and thermodynamic or statistical physics approaches introduced for large systems may not be suitable for smaller systems any more. However, for clarity reasons, we abusively call “phase transition” any phase transformation or spin state changes observed at finite sizes.

Common features exist in the size dependence of physical, chemical and thermodynamic properties in different materials in condensed and soft matter ([72] and references therein). In general, size reduction effects are characterized by the occurrence of two phenomena. First, finite size effects lead to confinement effects, which can affect transport phenomena (electrical, mechanical, optical, etc.) and more generally any collective behaviours. It leads to various quantum physical phenomena and cavity effects. On the other hand, the increase in the surface-to-volume ratio, resulting in the predominance of surface/interface properties, modifies mostly the phase stability and phase transformation kinetics. The typical size below which sizeable modifications appear in transport processes, physical and chemical properties and thermodynamic quantities are largely dependent on the structural characteristics of the material, the observed physical phenomenon, the measured observables (magnetization, thermal or/and electrical conductivity, electric polarization, etc.) and the transport mechanism (delocalized electron, acoustic phonons, dislocations, etc.) related to this phenomenon. Moreover, the different phenomena we usually consider as size effects do not take place at the same size range: although modifications of nucleation and domain growth processes may occur at the micrometre scale, quantum confinement effects are expected to be noticeable in very small systems (few nanometres). In the following, the different theoretical approaches and the fundamental concepts introduced to reproduce the experimental observations in SCO nanomaterials are overviewed.

#### 3.2.2 Surface effects: description with local physical property modifications

In a pioneering work, Kawamoto and Abe [73] have pointed out the importance of the surface in the spin transition properties of small systems. The well-known Ising-like model has been treated by Monte Carlo simulations on cubic and spherical systems with open boundary conditions (OBCs). In comparison with periodic boundary conditions (PBC), they showed that the hysteresis width decreases with size reduction in a different manner according to the nanoparticle shape, whereas the phase stability (in particular the equilibrium temperature) remains size independent (Fig. 6). The different observations were attributed to different nucleation processes existing in cubic and spherical systems and to the decrease in the number of interacting molecules when approaching very small systems.

However, experimental observations in some coordination nanoparticles have brought out a strong size dependence of the phase stability, which results in a downshift of the equilibrium temperature and the presence of incomplete transitions [74,75]. Indeed a number of residual HS molecules in the nanoparticles remain at low temperature and the associated HS fraction increases with the rise in the surface-to-volume ratio suggesting the presence of inactive molecules at the surface. The origin of inactive molecules at the surface with different physical and chemical properties may be attributed to surface chemical reactions, structural disorder and/or the alteration of the ligand by the external environment. This finding has been taken into account by adding specific boundary conditions (surface molecules fixed in the HS state [76] or a weakening of the ligand field at the surface [77,78]). The investigation of local surface effects in squared nanoparticles using the Ising-like model with such specific boundary conditions (surface molecules fixed in the HS state) has demonstrated that the HS surface sheet acts as a “negative pressure” [76] (nb: lattice deformations are not explicitly taken into account in the Ising-like model), but pressure effects can be grasped through an effective modification of the ligand field [76]. This leads to a downshift in the equilibrium temperature, which follows an algebraic law with the size (Fig. 7):

$${\mathit{T}}_{\text{eq}}={\mathit{T}}_{\text{eq}}\left(\infty \right)-\frac{8\mathit{J}}{{\mathit{k}}_{\text{B}}\phantom{\rule{0.25em}{0ex}}\mathrm{ln}\phantom{\rule{0.25em}{0ex}}\left(\mathit{g}\right)}\times \frac{1}{\mathit{L}-2}$$ | (27) |

_{eq}($\infty $) is the equilibrium temperature of the bulk material. This algebraic decay seems to be a universal feature of the size evolution of the phase stability and has been already observed in other first-order transitions such as melting process in bimetallic nanoparticles [79].

The presence of surface physical properties different from those of the nanoparticle core has been formulated in terms of surface energies in a simple nanothermodynamic model [80]. Surface energy is a physical ingredient experimentally accessible and despite the fact that different physical and chemical contributions are involved behind this concept, the thermodynamic approach remains a powerful tool to make direct comparisons with experimental observations [81]. Let us consider a spherical SCO nanoparticle with a radius R constituted of a bulk-like core, which is connected to a surface A, the so-called core–shell approach. This core–shell nanoparticle is in contact with a thermostat and a barostat and the Gibbs free energy G can then be written as an extension of the Slichter and Drickamer model:

$$G={n}_{\text{HS}}{G}_{\text{HS}}+\left(1-{n}_{\text{HS}}\right){G}_{\text{LS}}-T\left({S}_{\text{mix}}^{b}+{S}_{\text{mix}}^{s}\right)+A\left({\mathrm{\sigma}}_{\text{HS}}{n}_{\text{HS}}^{s}+{\mathrm{\sigma}}_{\text{LS}}\left(1-{n}_{\text{HS}}^{s}\right)\right)+\Gamma \left(1-{n}_{\text{HS}}\right)$$ | (28) |

_{HS}and σ

_{LS}correspond to the surface energy density costs for the creation of an HS and an LS surface, respectively. This approach is less restrictive than fixed boundary conditions according to the difference between σ

_{HS}and σ

_{LS}, the size reduction leads either to a downshift (σ

_{HS}< σ

_{LS}) or an upshift (σ

_{HS}> σ

_{HS}) of the transition temperature (Fig. 8b). The latter follows the same algebraic decay as presented in Ref. [76] (Fig. 8a), when the nanoparticle radius is reduced. When the difference $\text{\Delta}\mathrm{\sigma}$ becomes important, incomplete transitions occur at low or at high temperatures. In this latter situation, molecules at the surface become inactive and the fixed boundary conditions proposed in Ref. [77] are retrieved; however, the possibility to fix the surface in the LS state is not excluded.

Two main conclusions emerge from this nanothermodynamic model. First, the surface energy can be seen as the driving force of the spin transition at the nanoscale. Then, the downshift of the transition temperature with the occurrence of a residual HS fraction at low temperature in small nanoparticles seems to be the most common observation in SCO coordination nanoparticles, even if an upshift of the transition temperature and residual LS fraction have also been reported [81]. From a microscopic point of view, less constrained boundary conditions have been applied in the framework of Ising-like [82], anharmonic [72] and electromechanical [83] models, considering that chemical and physical (ligand field, vibronic degeneracies, elastic constant and lattice parameters) properties of the surface are different from those of the nanoparticle core.

Beyond the different theoretical descriptions used with the aim to model surface effects, a more fundamental question to ask is “what can be the mechanisms and the associated relevant chemical and physical ingredients at the origin of the different experimental observations?” In surface physics, it is well known that the first trivial consequence of the surface creation is the emergence of coordination defects: the local environment at the surface is different from the core due to broken bonds. Broken LS bonds cost more energy than HS bonds, which results in a stiffer lattice in the LS phase than in the HS phase at a macroscopic level. In Ref. [74], the surface energy density related to the coordination defect of a squared lattice has been expressed as a function of the bond energy of nearest- and next nearest-neighbour molecules:

$$\mathrm{\sigma}=\frac{-{\mathit{v}}_{1}-2{\mathit{v}}_{2}}{2\mathit{r}}-\frac{{\mathit{v}}_{1}}{\sqrt{\mathit{V}}}$$ | (29) |

_{1}and v

_{2}are the cohesion energies of nearest- and next nearest-neighbour molecules, respectively, r is the side of the squared lattice and V its “volume”. Because the cohesion energy of the LS lattice is stiffer (v

_{HS}< v

_{LS}), from Eq. (29) we can deduce that the coordination defect tends to favour the HS state at the surface (σ

_{HS}< σ

_{LS}). Interestingly, this is not enough to explain the drastic change in the phase stability of coordination nanoparticles.

The simulation of hollow SCO nanoparticles with the anharmonic spin–phonon model has highlighted the importance of the surfaces on the spin transition properties: rather than the quantity of matter, it is the surface-to-volume ratio that counts [77]. Indeed, a filled nanoparticle and a hollow nanoparticle having largely different volumes (i.e., different number of molecules) but the same surface-to-volume ratio exhibit similar spin transition curves.

#### 3.2.3 Surface effects: surface and order parameter relaxations

Surface relaxations or/and surface structural reconstructions, which accompany most of the time coordination defects, cannot be neglected in the modelling of the spin transition at the nanoscale. Therefore, the theoretical investigation of size reduction in SCO materials with models taking explicitly into account lattice degrees of freedom and elastic deformations becomes indispensable. Monte Carlo simulations in the (N,P,T) ensemble using a spring–ball Hamiltonian on a cubic lattice have shown that surface physical properties are delocalized in the core due to the surface–bulk mechanical coupling [84]. In other words, the approximation of a surface without depth may not be correct. Indeed, profiles of the mean squared displacement for different nanoparticle sizes show that below a critical size L_{c} (ca. L_{c} ∼ 7–8 nm in this case) the vibrational and elastic properties of the nanoparticle become completely different in comparison with the corresponding bulk materials.

The calculated Debye temperature and the melting temperature drastically decrease below L_{c} in the case of free surfaces, which can be understood by the fact that molecules at the surface can move more easily than those in the bulk and have the higher mean squared displacement and also a lower Debye temperature. On the contrary, when specific boundary conditions are imposed by ad hoc increasing the elastic constants of surface bonds, the Debye temperature increases with the size reduction.

Free surfaces seem to have only a relatively weak effect on the phase stability: the transition temperature remains nearly size invariant and no residual HS fraction occurs by approaching the nanometre size. The modelling and the prediction of the surface structural relaxation is rather complicated because different contributions, such as the coordination number, the crystallographic properties of the materials and so forth are involved in this phenomenon. However, two kinds of surface relaxations can be distinguished: surface structural relaxations and order parameter relaxations. In this present case, the order parameter of the spin transition is the HS fraction, ${n}_{\text{HS}}$. The electron–lattice coupling in SCO materials does not allow one to independently study these two contributions. Therefore, surface relaxations have been investigated within the framework of the Ising-like model. The decomposition of the free energy of an SCO chain has been proposed in the framework of a 1D Ising-like model, solved by transfer matrix techniques [85], comparing open and periodic boundary conditions. In this work, a detailed analysis of the bulk, the surface and the finite size effects on the thermodynamic properties has highlighted that the contribution of the finite size effects in free energy is always smaller than that of the surface free energy.

A similar 3D Ising-like model has been solved in the inhomogeneous mean-field approximation in which only spin-state surface relaxations are present [86]. The different surface thermodynamic quantities have been extracted as follows:

$$\begin{array}{l}{\mathit{F}}_{\text{\Sigma}}=\mathit{F}-{\mathit{F}}_{\text{B}}={\mathit{F}}_{\text{\Sigma}}^{\text{MB}}+{\mathit{F}}_{\text{\Sigma}}^{\text{relax}}+{\mathit{F}}_{\text{\Sigma}}^{\text{size}}\\ {\mathit{E}}_{\text{\Sigma}}=\mathit{E}-{\mathit{E}}_{\text{B}}={\mathit{E}}_{\text{\Sigma}}^{\text{MB}}+{\mathit{E}}_{\text{\Sigma}}^{\text{relax}}+{\mathit{E}}_{\text{\Sigma}}^{\text{size}}\\ {\mathit{S}}_{\text{\Sigma}}=\mathit{S}-{\mathit{S}}_{\text{B}}={\mathit{S}}_{\text{\Sigma}}^{\text{MB}}+{\mathit{S}}_{\text{\Sigma}}^{\text{relax}}+{\mathit{S}}_{\text{\Sigma}}^{\text{size}}\end{array}$$ | (30) |

Interestingly, when the evolution of the surface thermodynamic quantities is inspected, two findings can be established. Far from the transition temperature(s), the main contribution to the surface thermodynamic quantities comes from the coordination defect. However, around the transition temperatures, spin-state surface relaxations become non-negligible and compete with the coordination defects. However, in system exhibiting first-order transition (and also second-order transition), the surface relaxation effects become predominant mainly because bulk–surface correlation length massively increases (diverges) approaching the HS–LS coexistence region (the critical point).

Surface relaxations have also been investigated by considering a hexagonal lattice in which the unit cell contains two different sites: an active SCO and an inactive site. The aim was to mimic a simplified structure of the coordination network [Fe(pyrazine)(Ni(CN)_{4})] with an elastic anharmonic model [87]. The inactive site tends to favour the LS state by exerting a compressive local stress on the active site in the nanoparticle. This stress is relaxed at the surface where a tensile stress exists, favouring the HS state. The nanoparticle exhibits a two-step spin transition where the active molecules at the surface switch at a lower temperature than those located in the core and can be seen as an incomplete transition if the temperature is not sufficiently decreased (Fig. 10). Remarkably, in this study, no specific boundary conditions are imposed. The particular structure combined with OBCs leads to the emergence of this two-step transition. Interestingly, in the frame of this “elastic” model, contributions from structured and spin-state relaxations to the total surface relaxation could be separated.

#### 3.2.4 Impact of size reduction on cooperativity and collective elastic behaviour

An important question concerns the bistability property at the nanoscale. Can the hysteretic behaviour be observed in ultrasmall nanoparticles? This point is very important for the integration of SCO materials in information storage devices and other applications involving a memory effect. In a first approach, the disappearance of the hysteresis was predicted because of the decrease in the number of interacting molecules and the large fluctuations existing in a system containing a limited number of degrees of freedom. By analogy with magnetic nanoparticles, a phenomenon similar to superparamagnetism was expected in SCO nanoparticles. The first theoretical investigations of size effects using the Ising-like model support this assumption [73]. Whatever the numerical (Monte Carlo) or analytical (mean-field approximation) treatment of the Ising-like model and despite some differences because of the presence of short- or long-range interactions in the nanoparticle, a shrink of the hysteresis loop is generally observed with the size reduction [76,78]. However, experimental observations have shown the presence of hysteretic behaviour in the spin transition curves of some rather small SCO nanoparticles. Among different possible explanations, the size-dependence of the lattice stiffness has been evoked. Indeed, the size dependence of the Debye temperature ${\theta}_{\text{D}}$ extracted from Mössbauer data in PBA and SCO coordination networks have revealed a stiffening of particles below ca. 5 nm (Fig. 11a). Inserting these experimental findings in the phenomenological interaction parameter $\text{\Gamma}$ of the nanothermodynamic model ($\text{\Gamma}\phantom{\rule{0.25em}{0ex}}\sim \phantom{\rule{0.25em}{0ex}}{\theta}_{\text{D}}^{2}$), a nonmonotone size evolution of the hysteresis width is observed: for intermediate nanoparticle sizes the hysteresis loop shrinks and disappears before reappearing in very small nanoparticles (Fig. 11b) [80].

This result is in good agreement with experimental observations. The concept of cooperativity is upset by approaching the nanometer scale. There are less interacting entities in SCO nanoparticles, but it may be compensated by stronger bonds. A word of caution is needed here. Indeed, it is not obvious that the elastic theory remains valid in such small nano-objects and quantities such the Debye temperature and the bulk elastic modulus become meaningless. Another important point is the intrinsic character of the elastic/vibrational properties of a solid. Such quantities are not expected to be size dependent and the apparent or effective stiffening arises likely due to a compressive stress exerted by the external environment. The different effects of the environment are discussed in section 3.2.5. Another assumption to explain the observation of the bistability in nanoparticles is the presence of kinetic effects. Indeed, nanoparticles cannot reach the thermodynamic equilibrium because they are constantly crossed by heat and mechanical flows. Any small perturbations can destabilize the nanoparticle stationary state. The nanoparticle may switch permanently and very quickly from one state to another like in superparamagnetism. In this case, the existence of hysteresis has a kinetic origin and only the average time spent by the nanoparticle may be observed.

#### 3.2.5 Matrix effects (external environment)

Analogue to the Ising-like models mentioned above, the classical elastic models applied for open boundary SCO systems led to same unsatisfactory results: only the less steep transition and the smaller hysteresis, which vanishes for a critical size, could be reproduced. A different approach to account for all experimental features consisted in considering samples embedded in a polymer surrounding environment [88]. The motivation for this approach resides in the fact that some methods of preparation of nanoparticles imply the use of either the inverse micelle method or the polymer protected synthesis to obtain different particle sizes. The polymer itself forms a more or less flexible matrix around the nanoparticles, preventing the agglomeration of nanoparticles to form larger ensembles. The SCO molecules at the surface will therefore interact with the polymer matrix via either van der Waals interactions or hydrogen bonding. In a first approach, the systems used in the Ising-like model have been modified by considering an embedding matrix of nonswitching molecules, interacting with edge SCO molecules [78]. The results have satisfactory reproduced experimental data, but for a more realistic approach, the use of elastic models has been necessary. The polymer is represented in this approach as a rigid or flexible shell of nonswitching molecules, which interact with the molecules situated at the edge of the lattice by way of springs, with a given elastic constant denoted as k_{pol}. When the HS to LS transition occurs, the SCO nanoparticle shrinks, and if the polymer is not flexible enough, an increasing negative pressure appears at the edge of the system, originating in pulling forces from the polymer network. With this premise, the hysteresis not only shifts to lower temperatures, but also for high enough interactions between edge molecules and polymer, a residual HS fraction is produced (Fig. 12).

Similar results have been obtained with an Ising-like model containing a long-range interaction terms. The interaction with the external environment has been taken into account through a short-range coupling term between surface SCO molecules and the matrix containing inactive HS molecules. If the Ising coupling is strong enough, the molecules at nanoparticle surfaces are fixed in the HS state. Very recently, the experimental and theoretical investigations of physical properties of hybrid nano-objects and hollow nanoparticles have known a growing interest for practical and fundamental reasons. Indeed, SCO hybrid nano-objects, especially core–shell nanoparticles, allow a control of the spin transition properties and permit to generate synergies between the physical properties of the core and the shell (optical, magnetic, electric, elastic, etc.), which makes them particularly interesting for technological applications. On the other hand, these objects are very suitable for the fundamental investigations of the role of the external environment and interfacial interactions on the phase stability at the nanoscale. Monte Carlo simulations using the electrovibrational model have been performed by Oubouchou et al. [83] considering a nanoparticle embedded in an external environment whose elastic properties and structure were chosen to be similar to the HS phase (“soft” matrix). The analysis of the strain and local pressure mappings show that increasing the thickness of the surrounding medium leads to a release of the elastic strain at the interface nanoparticle/matrix and a decrease in the transition temperature (Fig. 13).

On the basis of a similar approach, Slimani et al. have demonstrated that the external matrix permits to maintain the bistability property in the SCO nanoparticle if the matrix/nanoparticle couplings are strong enough [41]. A further extension of the electroelastic model has been used to simulate SCO heterostructures [89]. In particular, the thickness and the rigidity influence of the surrounding shell on the spin transition properties of the core have been investigated by Monte Carlo simulation. The shell rigidity strongly affects the thermal hysteresis of the SCO core. In addition, resonance behaviour is observed when the shell rigidity is close to that of the core (Fig. 14).

Félix et al. have simulated core–shell nanoparticles using the spin–phonon model solved by the Monte Carlo simulations [77]. The control of the transition temperature of the hybrid nanoparticles can be achieved by tuning the elastic misfit (lattice parameters) between the core and the shell. When the core and the shell have similar elastic properties, the interfacial energy is zero and no phase stability change is observed in comparison with the corresponding nanoparticle or hollow nanoparticle. When the elastic misfit generates a compressive (respectively tensile) stress on the SCO core, an upshift (respectively downshift) of the transition temperature is observed. Interfacial energy is minimized when the core is in the LS state (respectively HS state). Interestingly, when an active core (e.g., SCO materials) is coupled with a shell sensitive to mechanical stress (e.g., martensitic materials) through a sufficiently coherent interface, a novel kind of hysteresis effect can be generated; the transition temperature is upshifted (respectively downshifted) during the warming (respectively cooling) process.

Despite elastic models give an origin of the core–shell coupling mechanism, it remains difficult to establish a link between the Hamiltonian parameters and experimentally accessible quantities. Therefore, an analytical approach modelling a core–shell nanoparticle with a coherent interface using continuum medium mechanism has been developed [90]. Analytical expressions of the interfacial stress $\left({P}_{\text{in}}\right)$ and the interfacial elastic energy density have been derived as a function of the Young's moduli (Y), the Poisson ratio ($\mathrm{\nu}$) as well as the size (a, b) of the core (subscript 1) and the shell (subscript 2):

$${\mathit{P}}_{\text{in}}=\frac{2{\mathit{b}}^{2}{\mathit{a}}_{2}\frac{1-{\mathrm{\upsilon}}_{2}^{2}\phantom{\rule{0.25em}{0ex}}}{{\mathit{Y}}_{2}\left({\mathit{b}}^{2}-{\mathit{a}}_{2}^{2}\right)}+{\mathit{a}}_{1}-{\mathit{a}}_{2}}{\frac{{\mathit{a}}_{1}}{{\mathit{Y}}_{1}}+\left(1+{\mathrm{\upsilon}}_{1}\right)\left(1-2{\mathrm{\upsilon}}_{1}\right)+\frac{{\mathit{a}}_{2}}{{\mathit{Y}}_{2}}\frac{1+{\mathrm{\upsilon}}_{2}}{\left({\mathit{b}}^{2}-{\mathit{a}}_{2}^{2}\right)}\left[\left(1-2{\mathrm{\upsilon}}_{2}\right){\mathit{a}}_{2}^{2}+{\mathit{b}}^{2}\right]}$$ | (31) |

_{1}, a

_{2}and b.

A quantitative investigation of the impact of the interfacial elastic energy on the phase stability can also be performed to optimize the elastic coupling between the core and the shell. The spin transition of the core–shell nanoparticle can be simulated by injecting the interfacial energy density into the nanothermodynamic model. It is important to note that the role of the elastic misfit on the collective behaviour and the bistability cannot be grasped with this model, because the elastic interaction between molecules is not explicitly taken into account in the Slichter and Drickamer [4] model and also the elastic interfacial stress cannot be coupled with the distortions occurring during the spin state change. This constitutes the main drawback of this approach.

#### 3.2.6 Consequences on thermodynamic quantities and vibrational properties

Size reductions have usually a strong impact on the thermodynamic properties of materials. The increasing role of surface (interface) thermodynamic properties destroys the postulates from which classical thermodynamics are built. For instance, thermodynamic quantities or functions are no more proportional to the quantity of matter and become nonextensive. This has been highlighted through the investigation of size dependence of internal energy and entropy of a spring–ball system calculated with histogram Monte Carlo methods (Fig. 16) [77].

As mentioned above, the thermodynamic equilibrium cannot be reached and even a perfect thermostat and a perfect barostat can no more impose its temperature and pressure, which are no longer intensive quantities. As a consequence, the external and internal pressures are related by a Young-Laplace–like relation, which can explain the apparent change in the nanoparticle stiffness:

$${P}_{\text{in}}={P}_{\text{ext}}-\frac{2}{\sqrt{V}}-4\sqrt{V}\phantom{\rule{0.25em}{0ex}}\frac{\partial \mathrm{\sigma}}{\partial V}$$ | (32) |

## 4 Conclusions

Elastic models constitute a step forward for the understanding of the behaviour of cooperative SCO materials. Largely improving previous models, they are, in all their different but convergent forms, both more precise and realistic than mean-field or Ising-like models used for simulating phenomenological behaviours of bulk SCO compounds and faster than ab initio approaches used for finding the electronic and bonding states on single molecules. They are able to reproduce both macroscopic curves and microscopic evolutions in bulk materials and also to simulate finite size effects for nano- and microparticles, with OBCs or embedded in surfactants. Future challenges for these models should include the simulations on structure compatible to real samples, either 2D or 3D, with the final aim to integrate real parameters found in ab initio models.

## Acknowledgements

We thank Dr. Gabor Molnàr for carefully proofreading the manuscript and for his useful comments. W.N. and C.E. thank, respectively, to CNRS France and UEFISCDI Romania for research grants.