Tutorial on the elastic theory of spin crossover materials

There are few systems showing bistability at the molecular scale, quality essential for the miniaturizing of memory devices.¹ Among these, spin-crossover (SC) (or spin transition) systems^2–4 belonging to the larger class of molecular magnets⁵ are one of the most exciting, as their switch between the two states in thermodynamical competition, the diamagnetic low spin (LS) state and the paramagnetic high spin (HS) state, brings a plethora of property changes: magnetic, optical, structural (cell parameters, volume), and vibrational. Spin transition systems are composed of transition metal ions containing four to seven electrons on their last d level and situated in an octahedral ligand field configuration. Actually, the ligand field is responsible for the bistable molecular character by splitting the d orbitals into antibonding e_g and weakly bonding t_2g orbitals.

In Fig. 1, we present the electronic configuration in the case of the most common ion in spin transition compounds—the Fe(II) ion, with six electrons on the last d level. According to the relative values of the ligand field and of the energy of pairing of the electrons: (i) the electrons occupy the orbitals of the lowest energy, by violating the Hundt rule, if the ligand field energy is larger than the pairing energy. In this case, the ion is in the LS state. (ii) The electrons occupy the largest possible number of orbitals, according to the Hundt rule, if the pairing energy is larger than the ligand field energy. Then the spin is maximum and the ion is in the HS state. The most interesting case is when the pairing energy and the ligand field energy are on the same magnitude order: then, the spin conversion can easily take place by variation of external stimuli such as temperature (high temperature favors the HS state due to its larger degeneracy), pressure (high pressure favors the LS state due to its smaller volume), magnetic fields, or by exposing the sample to an electromagnetic radiation (mainly the LIESST Effect—Light Induced Excited Spin State Trapping—which is the trapping at low temperature of the HS metastable state by the light).

In Fig. 2, we present the molecular structure of a typical two-dimensional spin-crossover crystal [Fe(bbtr)₃](ClO₄)₂]⁶ in which the Fe(II) central ion is surrounded by six ligands, together with an image of a hexagonal shaped crystal of the same compound. As we shall see later, this kind of structure was considered an experimental basis for the development of the mechanoelastic model.

When a spin flips, the length of the Fe-ligand bond varies by about 10% [typically ∼0.2 Å on 2 Å for the ion Fe(II)] and the ligand field dramatically changes. The spin transition phenomenon is thus not correctly described by the simplistic condition that the ligand field energy should be on the same order of magnitude as the pairing energy, but by the fact that the variation of ligand field energy is crossing the value of pairing energy.

Because of the difference in metal ligand distances, the frequencies of vibration are lower in the HS state, while the density of the vibration modes and thus the entropy are higher in the HS state. Therefore, the thermal transition always takes place from the LS state at low temperature (which is the fundamental state) to the HS state at high temperature.

The system is, in general, described by the high spin fraction, here denoted as $n H S$ (other symbols are used elsewhere such as $X H S$⁠, $f H S$⁠, or $γ H S$⁠), which represents the proportion of molecules in the HS state. Reversely, $n L S=1− n H S$ corresponds to the fraction of molecules in the LS state.

Below this temperature, the enthalpy factor dominates and the stable state is LS; above $T 1 / 2$⁠, the entropy factor dominates and the HS state is the stable one. Usually, the population of the LS state is negligible at high temperature due to the strong entropy term, which derives from the “effective” degeneration ratio of the system in two states, $g= g H S/ g L S= exp (ΔS/R)$⁠, where $g H S$ and $g L S$ are, respectively, the degeneracies of HS and LS states (it can reach several thousands⁹).

Advances in experimental findings on spin-crossover compounds have been accompanied by the development of appropriate models. Gradual thermal transitions specific to non-interacting (such as solution) compounds have been first reproduced by considering a classical Boltzmann distribution between the vibronic manifolds of the HS and LS states.¹⁰ Further on, thermal transitions showing hysteresis loops have been simulated either in the framework of domain models, considering distributions of switching temperatures,^11,12 or by adding an interaction term to the Gibbs free energy, thus leading to the cooperative mean-field model,¹³ which is useful for simulating hysteresis loops and sigmoidal shaped, relaxation curves after photoexcitation at low temperature. Other phenomena like stepped transition imposed the use of more complex Ising-like models,¹⁴ while the vibrational modes of molecules and the differences of entropy and electronic energy between HS and LS isomers have been obtained by Density Functional Theory (DFT) studies.¹⁵ Using these approaches, all macroscopic spin-crossover phenomena known by the middle of the first decade of this millennium were satisfactory modeled.

These circumstances dramatically changed in 2006 by diffraction studies of single crystals by Pillet et al. in Nancy,¹⁶ which proved the existence of domains during spin transition, followed a couple of years later by optical microscopy experiments performed by spin-crossover research groups Leiden,¹⁷ Versailles,^18,19 and Toulouse.^20,21 All these ground-breaking experiments proved the existence of clusters of the same spin state molecules that nucleate from corners and imposed the development of new models. These approaches considered the elastic interactions between neighbor spin-crossover molecules, which propagate to the whole lattice during the spin transition. In this Tutorial, we shall present the concepts of elastic models, with illustration of results for the mechanoelastic approach. The paper is organized as follows: first, we discuss the mean-field and Ising-like models, showing their strengths and weaknesses and insisting on the points taken over by the elastic model, then we present the steps for implementing the elastic model, explaining the formation of clusters and the changes induced by anisotropic effects on two-dimensional systems. Further sections are devoted to the study of three-dimensional systems, nano or micro-systems embedded in matrices, such as polymers or surfactants, or situated on substrates and to ultrafast dynamics after photoexcitation.

Using phenomenological interaction parameters, acting similarly for all molecules in a given sample, the mean-field models have successfully reproduced several of the characteristic features of spin-crossover (SC) compounds, such as the occurrence of a thermal hysteresis and its dependence on the strength of intermolecular interactions,² the shift of the thermal transition toward higher temperatures^13,22 when applying external pressure, the low temperature photoexcitation curves and the sigmoidal shape of subsequent relaxation curves,²³ or the light induced bistability.^24–26 

This relation allows us to represent the equilibrium value of $n H S$ as a function of temperature (corresponding to the canonical state) and to study the effect of interaction constant $Γ$ on the shape of transition and the width of the hysteresis loop, if any.

Generally, mean-field models give good results for fitting thermal transitions and not very cooperative relaxation curves.²⁹ But, as one can see in Fig. 3, if the cooperativity factor is too high, then there is a discrepancy between experimental and simulation results, which has been explained as the “tail effect” due to the progressive onset of the correlations associated with short-range interactions.³⁰ Because of the premises on which they are based, mean-field models do not distinguish between short- and long-range interactions, neither are able to explain how the effect of the individual switch of molecules spreads through the entire lattice. In addition, mean-field models cannot be applied for the study of size effects in nanoparticulate spin-crossover systems³¹ or for systems prepared as thin films on substrates,³² whose study has become increasingly important in recent years due to potential applications.

The first steps in the theoretical investigation of the behavior of spin transition compounds using the Ising model were made by Wajnflasz and Pick in the 1970s.³³ They were the first to introduce an Ising-like Hamiltonian with the intermolecular interactions and to predict the presence of a first order transition.³³ Twenty years later, in order to avoid the deficiencies of mean-field models, Bousseksou et al. proposed an Ising-like model with nearest-neighbor interactions, solved in mean-field approximation,¹⁴ which was the starting point for more complex Ising-like models, taking into account simultaneously short- and long-range interactions,^14,34–36 with extensive studies devoted to dynamic³⁵ and static properties³⁷ of SC compounds. A recent review of Ising-like models can be found in Ref. 38.

We should notice that the role of the field in the classical Ising model is taken by the term $D− k BTlng$ (the so-called Gibbs energy), which is temperature dependent. Therefore, this approach is not a pure Ising approach, so the results concerning critical temperatures do not have a correspondent in the present model, which is referred to as “Ising-like.” While the short-range interactions can be correlated to the existence of covalent bonds between neighboring molecules, the long-range interactions are considered purely of elastic nature mediated by the lattice.^39,40

In the Arrhenius algorithm, one decides if a molecule switches from a state to the other by generating a positive subunit random number and comparing it with the corresponding switching probability. If the random number is smaller than the switching probability, the transition is accepted; otherwise, the molecule preserves its initial state. A Monte Carlo Step (MCS) is concluded when a number of molecules equal to the number of molecules in the system have been checked. The Ising-like model has been successful in simulating the behavior of spin-crossover systems, including nanoparticles with open boundaries (see Fig. 4, left).⁴³ As we notice in Fig. 4 right, nucleation is present in the systems modeled within the Ising-like approach. However, the nucleation centers are situated both at the edges and inside the bulk, which does not correspond to optical microscopy experimental data that show clusters starting from edges. This is due to the fact that Ising-like models do not take into account distortions created by the different volumes of the two molecular states—actually the physical meaning of the interaction parameters is not clear. As the elastic interactions responsible for the cooperativity phenomenon in spin-crossover compounds originate exactly in these distortions, a new class of models, the elastic models, has been considered.

Essentially, the different volumes of the molecules in the two spin states are at the origin of the elastic stresses in the crystalline network and of the subsequent nonlinear hysteretic behavior of the spin-crossover compounds. The elastic stresses are responsible for the cooperative effects and thus the variation of the switching probability and for the displacements of molecules during the transition from one state to another, as shown in Fig. 5. This effect was described using the concept of “image pressure,”⁴⁴ or, more simply, an “elastic interaction” between the molecules. The long-range character of the elastic interactions has been recently experimentally proven in the case of cooperative single crystals.⁴⁵ 

The elastic models elaborated in the second part of the first decade of this millennium have been preceded by several other theoretical attempts to consider elastic interactions as premises for spin transition. In a previously mentioned paper, when discussing the Ising-like model, Wajnflasz stated that the interactions between spin-crossover units are not from the magnetic origin, but merely are resulting from the bond length change during the transition; he can be therefore considered predecessor of elastic models.³³ The lattice strains and their couplings with the d-electrons responsible for the spin transition have been taken into account by Kambara,⁴⁶ following the phenomenological description made a few years before by Slichter and Drickammer.⁴⁷ Almost simultaneously, Spiering et al. postulated that elastic interactions arise from stresses induced by the change in volume of the lattice during transition,⁴⁴ while Kambara introduced in his model a dynamic variable, stating for intramolecular distortions. Even if, before 1990, it was stated that the elasticity would explain almost of the spin-crossover behavior, these approaches have been abandoned for almost 20 years, due to the lack of relevant experimental data and of enough high-performing computers, needed for the complexity of calculations, due to the large number of parameters.

As we stated in the introductory part, optical microscopy and diffraction experiments that proved the existence of clusters triggered the implementation of elastic models, sometimes called ball and spring models, which considered the molecules as rigid spheres linked by springs and able to move inside a lattice. These models, independently and almost simultaneously proposed by several groups of researchers have different names (mechanoelastic,⁴⁸ electroelastic,^49,50 anharmonic,⁵¹ two-variable elastic,⁵² or simply elastic^53–55), are different at the first sight as they use various Hamiltonian formalisms, diverse shapes, and internal structures of the samples or multiple analysis methods, but, finally, are based on the same hypothesis and lead to the same conclusions. A review on all elastic models and on the differences and similar aspects between them can be found in Ref. 56.

In the following, we shall mainly refer to the results obtained in the framework of the so-called mechanoelastic model.

Let there be a system of n particles, in which every particle has x neighbors (except those on the edges of the system). As discussed above, in the elastic models, one assumes that the molecules are rigid and mobile spheres connected through springs with the same elastic constant k and the interactions between them are only of elastic origin. The key of any elastic model is how to find the new positions of the molecules after every switch. Different studies have proposed different methods for computing the new molecular positions: either relaxing the lattice by taking into account small displacements in all directions for randomly selected molecules,⁵⁴ by considering equation of motion of all molecules in the framework of a Nose–Hoover formalism,⁵³ or by deterministic computation of the motion of molecules,⁴⁸ by solving a system of differential equation, like in the case of the mechanoelastic model.

The following algorithm is used: (i) initially, all molecules are in either the HS or LS state; (ii) we choose a molecule randomly (let it be $i$⁠); (iii) we calculate the local pressure acting on the $i$th molecule using Eq. (15), then the probability of transition (14), depending on the state of the molecule; (iv) we generate a random number $r∈( 0 , 1)$⁠; (v) we check if the molecule switches—this is the case if the generated number is smaller than the probability. (vi) We repeat steps (ii)–(v) for a number of times equal to the number of molecules; (vii) we compute the new positions of all molecules in the system following the procedure discussed above; when we finish this step, we have performed a Monte Carlo Step (MCS).

Then, if studying the relaxation at constant temperature, we simply continue with a subsequent MCS, while if we are interested in the thermal transition, we modify the temperature after the desired number of MCS.

In Table I, we provide a comparison between mechanoelastic and Ising-like models.

The choice of the lattice is another important point to be discussed here. Usually, the bidimensional lattices characteristic to Ising models are rectangular; however, in the elastic model with the closest neighbors’ interactions, it was proved that rectangular (as well as hexagonal) configurations are not stable. Let us follow the proof in Figs. 6 and 7, where we have represented systems in hexagonal and rectangular configurations, in their initial LS states (open squares) and after the lower halves of the systems were switched from the LS state to the HS state. The new positions of the molecules are represented with circles (red corresponds to HS state, while yellow means LS state—this conventional representation will be used throughout all the paper). We notice that the rectangular and the hexagonal systems present significant shape distortions after the transitions. Moreover, some of the molecules are very close to others, which indicates large elastic forces. In order to surpass this difficulty, diagonal springs (linking the central molecule with its second order neighbors) should be considered. However, if one wants to keep a single elastic constant between molecules, then the use of a triangular configuration is needed. As we can notice in Fig. 8, the triangular system has a different behavior from those mentioned before (it does not have significant distortions nor large elastic forces), which makes it a stable system even in the case of using a single elastic constant between the closest molecules. Therefore, for the simplicity of the discussions, as we intend to use the minimum number of parameters, in the following of the present paper we shall refer mainly to systems in triangular configurations.

A last essential parameter in the model is the damping constant. Even if it does not directly intervene in the expression of probabilities, it has an important role in determining the positions of molecules after the transition, which can affect not only the computation time, but also the simulation results. For moderate ratio between the damping coefficient and the elastic constant, the final position of the molecules is the same; however, one can detect differences in the path toward this position. If the damping is high, the system converges to the equilibrium without oscillating (overdamping regime), while for smaller values of the damping the molecule effectuates an osscilatory motion with the amplitude decreasing to zero (underdamping regime). For very small values of the damping coefficient, the system enters into an uncontrolled oscillatory motion. Therefore, in order to save computation time, it is important to use values of the damping coefficient ensuring a situation in-between an underdamping and an overdamping regime.⁵⁹ 

Before presenting the results of the simulations in the framework of the elastic models, let us discuss the values of the parameters needed in the elastic model. We have considered realistic values for the radius in the HS (0.22 × 10⁻⁹ m) and LS state (0.2 × 10⁻⁹ m), thus giving a difference between HS and LS radius of Δr_HL = 0.02 × 10⁻⁹ m. The unit cell length (the distance between molecules in the LS state) has been taken as c₀ = 10⁻⁹ m. In the first stages of applications of the elastic model, the elastic constant took arbitrary values.^48,60 Later studies linked it to the experimentally determinable bulk modulus B by the formula $k=3 c 0B$⁠,^61,62 which results in a value of k of around several N/m for standard spin-crossover solids (with B around 50 kbar⁶³).

The cooperativity is modulated not only by the elastic constant of the springs, but also by the coefficient $κ$⁠, which establishes the relation between effective pressure in the system and all the other quantities. By comparing Arrhenius and Metropolis approaches, in Refs. 61 and 64 the value of $κ$ was found to be equal with $AΔ r H L$⁠, where A is the molecular area and $Δ r H L$ is the radius variation between HS and LS. Let us notice that the molecular area A appears at the denominator in the expression of the pressure (15), so it will be simplified with A from the expression of $κ$ and its value will not be needed in the model. For computational practical reasons, the values of energy quantities are expressed in K. The conversion from K to cm⁻¹ is then simply obtained as 1K = 0.695 cm⁻¹ and from K to J as 1K = 1.3806 × 10⁻²³ J.

In the following, we will discuss the results of the simulation of the relaxation and thermal transition curves in the framework of mechanoelastic model. To study relaxation, we have applied the model to an open boundary hexagonal shaped system, in a triangular configuration, as described above in which every molecule from the bulk has six neighbors (depending on their positions, the molecules at the edge of the system have fewer neighbors than those in the bulk).

We assumed that all the molecules are in the HS state, at a temperature below thermal transition—therefore, the system will relax toward its metastable state. We analyzed configurations in different moments and the dependence of $n H S$ (the percentage of HS molecules) on time (measured in Monte Carlo Steps) for different values of the constant k. The results in Fig. 9 have been obtained at a temperature T = 50 K for a system with 10 981 molecules, with the intrinsic parameters of material D=1000 K, g=1096, giving a thermal transition at around T_1/2=142 K and considering the probability confinement constant $τ=2000$⁠. The relaxation curves are sigmoidal, sign of the cooperativity of the system for lower values of the elastic constant k and they transform to abrupt curves as k increases. The nucleation process itself is modulated by the value of k. The nucleation is absent and the molecules switch randomly from the HS state to the LS state for low values of k, such as in the mean-field model. For higher values, we notice that edge molecules switch faster than those from the bulk ones and nucleation starts from the edges of the system, extending inwards, in an avalanche-like phenomenon.⁶⁵ 

The interface between the LS and HS domains is close to a circular shape, as this shape should store the minimum amount of elastic energy. A more precise computation leading to circular shapes for small clusters has been made using the Kawasaki approach.^64,66 However, it was proved by both mechanoelastic⁶⁷ and electroelastic⁵⁰ simulations that for larger and rectangular systems, the interface shape changes from circular to flat as it propagates throughout the sample and the spin transition proceeds making the simulations closer to experimental observations which show a flat boundary between the HS and LS domains.⁶⁸ We should also notice due to thermal fluctuations which are allowed by the probabilities (14), the boundary is blurred by the diffusion in-between the LS and HS domains, with the fuzziness more accentuated at high temperatures.

We have also analyzed the thermal hysteresis for the same system, by starting in a fully HS configuration at a temperature T = 200 K and decreasing it with a constant temperature sweep rate until all the molecules have switched to the LS state and increasing it back. Both the shape and the width of the hysteresis loops depend on the value of k and become steeper and respectively larger with increasing spring constant. As in the relaxation case, the increase of k leads to the transition from a mean-field regime, with randomly switching molecules, corresponding to smoother and gradual transitions, to one of the appearances and growth of nucleation centers; this feature is more obvious in the case of cooling branches of hysteresis, due to smaller thermal fluctuations. Other studies concerning thermal transition using elastic models concerned the influence of the temperature sweep rate (a fast temperature sweep rate results in a broadening of the thermal hysteresis) and of the pressure effects (Fig. 10).⁶⁰ 

The main results (thermal hysteresis, cluster formation) are similar to those obtained in the framework of the Arrhenius approach. In Fig. 11, we have represented the nucleation for systems with moderate and high interactions and we observe the same trend of producing larger clusters for higher elastic constants.

The existence of the cooperativity in these elastic systems is illustrated in Fig. 12 (left) where we have represented the distributions of intermolecular distances for HS molecules for various $n H S$ values at the beginning of HS–LS relaxation. We notice that the average intermolecular distances decrease as relaxation proceeds, which is equivalent of asserting that the average local pressure over HS molecules is larger and their probability to switch to the LS state increases. A further discussion on the distribution of interaction forces over HS and LS molecules during relaxation can be found in Ref. 61.

The nucleation in elastic systems has been discussed using the concept of “macroscopic nucleation,”⁵⁵ which implies that the size of the critical nucleus is proportional to the system size. A main question in the elastic models is why the nucleation starts always from edges and not from bulk?

In order to understand this feature, let us analyze the data in Fig. 12 (right), where we have represented the elastic energy stored in the whole system for clusters hypothetically grown from a corner, an edge or from inside the bulk, together with two lines representing the gain of Gibbs energy, as a function of cluster size (which is actually $n L S$ fraction). It is obvious that the lowest elastic energy corresponds to clusters nucleating from a corner. An intermediate elastic energy value is obtained for clusters starting from anywhere in the edge, where a hypothetical cluster nucleating from the bulk would determine a huge increase in the elastic energy of the system.

Therefore, the gain of free energy Gibbs $Δ W G$ shows a linear dependence on the $n L S$ and it should be compared with the elastic energy in order to establish if a cluster will grow or not. In the case of a temperature below $T 1 / 2$⁠, this line, which always passes through the origin of axes, has a positive slope, increasing with the decrease in the temperature. This observation allows us to discuss the stability of a cluster in a simple manner: only if the gain in the free Gibbs energy is bigger than the increase in the elastic energy of the system due to distortions, the clusters are stable. Therefore, a hypothetical cluster formed inside the bulk will vanish, as the elastic energy it produces will be much larger than the gain of Gibbs energy, and only cluster from edges (especially at low temperature where the gain in Gibbs energy is large) and from corners will be able to spread. Extensive simulations have shown in addition that the spreading of a single large cluster is expected in small systems with high interactions, while several clusters may develop in the case of larger systems and moderate interactions.⁶⁴ Defects and holes in the cluster formation may change these aspects, as they influence the elastic energy in the system.^70,71

The largest amounts of energy and the largest strains are located at the interface between the two phases as observed in experimental data and predicted in several elastic models (Fig. 13).^49,72,73 This will result in the fact that for a cluster of a given size, the wetting angle between the interface and the lattice side will be close to $π/2$⁠, as proved using the Kawasaki extension of the elastic model.⁶⁴ 

In Secs. III C 1 and III C 2, we have assumed that smaller spherical shaped LS molecules switch to larger HS spherical shaped molecules. However, sometimes the LS to HS transition is accompanied by structural distortions involving larger one-dimensional expansions and different bond-length variations on every axis.^74,75 Moreover, some spin-crossover systems (as 1D systems) exhibit larger elastic interactions on one direction than in the others⁷⁶ (see Fig. 14). Both situations lead to anisotropy effects that we discuss in this section.

Let us assume that the elastic force has a stronger effect on the horizontal axis (⁠ $k x$ has a higher value than $k y$⁠, which means that the pressure on the horizontal axis is higher than the pressure on the vertical axis) and analyze the behavior of the system during relaxation and thermal transition, with the same intrinsic parameters as those used in the isotropic case. If the anisotropy is large, the minimum of the energy in the system is obtained when clusters will form along the vertical axis rather than along the horizontal axis and this orientation is more prominent as the difference between $k x$ and $k y$ is increasing (Fig. 15). Moreover, we notice that the transition from the HS state to the LS state is closer to an exponential for a larger difference between the elastic constants $k x$ and $k y$—in this case, the sample acts like each line will be almost independent of others.

As in the previous case, we also analyzed a thermal hysteresis by starting with a temperature of T = 200 K. The intrinsic parameters are the same used for relaxation and the results are shown in Fig. 16. As in Fig. 15, clusters evolve along the vertical axis rather than along the horizontal axis and the transition from a state to another is more abrupt while increasing the difference between the elastic constants $k x$ and $k y$⁠. During the heating processes, the evolution of clusters is slowed down by the larger thermal fluctuations.

Here we considered a transformation from circular shape molecules in the LS state to elliptical shape molecule in the HS state, with a larger semiaxis on the horizontal direction. As experimental data do not, however, show a too large difference between these two semiaxes, we have considered in our simulations a radius of 0.2 nm in the case of the LS state and two semiaxes of 0.225 and 0.215 nm in the HS state.

In this particular case, the procedure for determining the elastic force is somewhat more complex, as we explain below. In order to determine the force between two molecules, we consider that a molecule has six types of neighbors. Each molecule has six anchors attached, with the following coordinates: $A 1(x,y− r y)$⁠, $A 2 ( x − r x cos ⁡ π 6 , y − r y sin ⁡ π 6 )$⁠, $A 3 ( x + r x cos ⁡ π 6 , y − r y sin ⁡ π 6 )$⁠, $A 4 ( x − r x cos ⁡ π 6 , y + r y sin ⁡ π 6 )$⁠, $A 5 ( x + r x cos ⁡ π 6 , y + r y sin ⁡ π 6 )$⁠, and $A 6(x,y+ r y)$⁠, where x and y are the coordinates of the molecule. These anchors change their position depending on the radii $r x$ and $r y$ and on the molecule's coordinates. To calculate the force between molecules and j, we shall first identify which anchor of the $i$th molecule corresponds to the neighbor j and vice versa. Then, we determine the distance between these molecules and, by multiplying the result with the elastic constant, one finds the elastic force.

The difference between these semiaxes leads to an anisotropic effect, but, keeping it reasonable, the observed anisotropic effect presented in Fig. 17 is less important than if considering a large difference between elastic constants $k x$ and $k y$ (as in Fig. 16 bottom snapshots). However, even in a less extent, we notice in Fig. 17 that spreading of clusters is produced along vertical lines.

As the mechanoelastic model is applied on open systems, without periodic boundary conditions, it is important to establish to which extent the size of the system influences the relaxation and the thermal transition. First, we have computed several relaxation curves for the same system, considering various sequences of random numbers. As the system size is smaller, the distributions of relaxation times spread over a larger range, according to the well-known central limit theorem that states that the re-averaged sum of a sufficiently large number of identically distributed independent random variables each with finite mean and variance will be approximately normally distributed.⁷⁷ As the number of molecules in the system increases, clustering becomes a statistical phenomenon, which leads to smaller fluctuations of the switching time.

By analogy with the Ising model applied for magnetic systems, one could expect that the relaxation time increases with the system size. However, in some situations, in the elastic model, by fine tuning the elastic constants and other system parameters, it is possible to obtain similar relaxation times for different sizes, even if in most cases one will obtain faster relaxation for smaller systems. In Fig. 18, one can notice that the relaxation time for the three systems with different sizes is approximately the same, while the shape is different. In small systems, spatiotemporal fluctuations determine faster build-up of clusters, but at the same time they prevent the cluster to propagate. In larger systems, the first cluster needs a larger time to nucleate; but once formed it will expand very fast and will result in a more pronounced sigmoidal shape.

The decrease in the width of the hysteresis loop with the system size is illustrated in the inset of Fig. 18. As in the Ising-like model, one can find a critical dimension of the system for which the thermal hysteresis vanishes⁶⁰ and which can be considered the superparamagnetic limit of spin-crossover complexes.

Some of the typical experimental spin-crossover solids can be characterized as two-dimensional systems composed of stack of layers connected by weak interactions, but others are purely three-dimensional. As the cooperative phenomenon implies the expansion of clusters toward the inside material, the study of clustering in three-dimensional systems is important for understanding the real nature of the whole process. Moreover, the spin-crossover nanoparticles have different shapes; therefore, it is necessary to study how the cluster behaves depending on the size of samples. The triangular structure previously used in two dimensions is translated in three dimensions in two possible simple regular lattices: face-centered cubic (fcc) and hexagonal close-packed (hcp).⁷⁸ The two structures are equivalent: every bulk molecule has 12 neighbors: six in the same plane and three in upper and lower planes; they differ by the fact that while the fcc system has a triangular pyramid shape (a regular tetrahedron), the hcp system is a hexagonal pyramid. In order to increase the symmetry of the system, we doubled the pyramids in the lower side of the base plan. These two structures correspond to the most stable three-dimensional system, preserving the condition that the distance between two neighboring molecules in pure HS or LS states is the same and requires the use of only a single elastic constant.

Depending on the value of the elastic spring constant, not only the relaxation time varies, but also the shapes of the relaxation curve change. However, for the same spring constant, the relaxation curves look similar for both fcc and hcp systems, with a slightly faster relaxation for hcp systems. In the case of moderate interacting systems, the relaxation curves are sigmoidal as expected for cooperative processes when the energy barrier decreases proportionally to $n H S$⁠. In the case of stronger interactions, abrupt relaxation curves with discontinuous jumps of $n H S$ can be observed. It should be noted that a comparison between the relaxation times for different spring constants is not straightforward as the spin-transition temperature itself changes with the spring constant.⁷⁹ 

In both hcp and fcc systems, there are clusters growing from the corners situated at the base of the pyramids and then advance toward the bulk, as one can see in Fig. 19. As there are several corners in the six faces pyramid, the nucleation probability is somewhat higher than in the case of the tetrahedron; therefore, the relaxation is faster, depending on the number of corners.

A similar study has been made for the thermal hysteresis for both fcc and hcp systems (Fig. 20). In this case, both hcp and fcc structures with similar elastic constants lead to hysteresis loops with the same width, but the nucleation takes place differently than during relaxation. In this situation, both LS and HS clusters start from up and bottom corners. We can understand this fact as follows: in the two-dimensional rectangular or hexagonal samples, we expect that cluster should start from any of the corners, with equal probability. However, in the 3D case, the corners are no more equivalent and the solid angle for the up and down corners is smaller than for the other corners. Therefore, the energy for a cluster to start from up and down corners will be smaller than in the case of a cluster starting from middle corners. The system is especially sensitive to the value of the elastic energy in the vicinity of the thermal transition, whereas one has to compare this elastic energy to the Gibbs energy difference between the states. Consequently, during the thermal transition, the clusters will predominantly start from up and down corners. In the case of a relaxation, at a much lower temperature than the thermal transition, the Gibbs energy difference between LS and HS states is high enough, so the system may allow higher elastic energies as implied by the formation of clusters from middle corners. In addition, there are more middle corners than up and down corners and consequently the cumulative probability for a cluster to start from one of the middle corners during relaxation is higher. A similar situation was observed during the study of cluster behavior in the case of ellipses: even if one expects that the clusters start always from the vertices, sometimes they start form co-vertices, where the number of molecules is larger.⁶⁴ 

In the following, to minimize the influence of the irregular shape of the above 3D lattices, we have performed similar simulations for spheres. The spheres have been obtained by cutting large hcp samples, keeping to a minimum the roughness of the surface. In this case, similar to circular situations, the absence of corners will determine the nucleation development from anywhere on the circular surface (Fig. 21).

Small clusters are of spherical shape (as the minimum interface area is obtained in the case of spheres). The exact shape can be figured out using the Eden algorithm:⁸⁰ (i) we first generate a seed formed by a few connected molecules at the edge of the sphere; (ii) we let the seed nucleate, by only allowing the switch of the molecules neighboring an already switched molecule. We notice that first the shape of the growing cluster is circular with the convexity decreasing while advancing toward the bulk.

From a historical point of view, the introduction in the model of a matrix embedding the spin-crossover sample has been determined by the experimental data showing a different behavior for spin-crossover nanoparticles compared to bulk samples; i.e., a smoother thermal transition, sometimes incomplete, shifted toward lower temperatures accompanied by a narrower or even no hysteresis. A first try to reproduce these features, made in the framework of classical elastic models and taking into account different sizes for open boundaries spin-crossover systems, led to unsatisfactory results: only the less steep transition and the smaller hysteresis that vanishes for a critical size could be reproduced (Fig. 18, inset; similar results have been obtained for the Ising model—Fig. 4, left). Therefore, in order to account for all experimental features, one considered samples embedded in a polymer surrounding environment.^60,81 Similar studies have been performed in the framework of Ising-like models, and the results have been satisfactory, in the limits we already discussed.^43,82,83 For practical reasons, the polymer was represented as a rigid shell of non-switching molecules that interact with the molecules situated at the edge of the lattice by way of springs, with a given elastic constant $k p o l$—which should correspond to hydrogen bonding or van der Waals interactions. Other approaches in the framework of electroelastic models implied the use of an environment with similar properties as the HS phase,⁸⁴ the tuning of the elastic misfit parameters between the spin-crossover lattice and the matrix,^85,86 or an analytical approach using a continuum medium mechanism.⁸⁷ During the HS–LS transition, the volume of the spin-crossover nanoparticle has the tendency to decrease, and, for a frozen or not flexible enough matrix, a negative pressure manifests due to the elongations of the springs linking the edge molecules with neighboring molecules from the polymer network. Assuming these hypotheses, one may obtain, in addition to the decrease in the hysteresis width, a shift of the transition toward lower temperatures and, sometimes, for high enough interactions between edge molecules and matrix molecules, incomplete transition with the corresponding residual HS (Figs. 22 and 23).

We have shown in Sec. III C 2 that the nucleation starts from corners or edges. Figure 24 shows that even if the switching probability of edge molecules is low due to the negative pressures induced by the matrix molecules, the cluster formation is similar to the open boundary systems, due to the propagation of this negative pressure throughout the entire sample (which was not possible in Ising-like systems).

In order to reproduce the behavior of nanoparticles on a surface, the model was adapted as in the following: we consider a 3D sample composed of several layers and the molecules on the first layer are linked with the fixed surface molecules, with a different constant $k s u b$⁠. For open boundary systems composed of several layers of molecules or membrane deposited on deformable substrates, buckling effects have been observed,^58,88 while the interactions with a fixed substrate determine the appearance of a residual $n H S$ during thermal transition of [Fe(pz)Pt(CN)₄] nanoparticles dispersed on a sapphire surface⁸⁹ due to the so-called epitaxial strain.³² 

As we can notice in Fig. 25, the results show the same tendency as for the molecule embedded in matrices and the thermal transition is shifted toward lower temperatures, while clustering starts from edges in some cases form middle of the layer opposite to the surface.³² 

As we have explained in the introductory part, the elastic simulation implies two steps: the spin dynamics, driven using transition probabilities and the lattice relaxation during solving of the differential equations system. For comparing 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 considered in solving the differential equation system after spins have been fixed and $n 2$ corresponds to the number of MCS between two successive solutions of differential equations.⁹⁰ In all the situations discussed so far, we have been interested in finding the equilibrium positions of all molecules in the system and, therefore, r (and $n 1$⁠) were large enough to allow the lattice to relax the excess energy generated by the switch.

However, this is not the case for studies following the ultrafast photoexcitation,⁹¹ where “the elastic step” (the continuous increase of $n H S$ after the irradiation is stopped) has been observed. This phenomenon was explained by the reorganization of the HS molecules toward the edges of the sample.⁹⁰ 

In Fig. 26, we represent simulation of the elastic step for two extreme r values. If r is large, then the lattice relaxes completely until the next switching probability checking, the pressure inside the lattice is released and the HS molecules will keep their state. If r is small, the lattice will do not relax completely and, therefore, the pressure inside the sample will be high and several freshly photoexcited molecules will switch back toward the LS state. In this case, there is a delay between the HS evolution and the increase in the surface of the system. The surface shows a faster variation for larger r values, while, for smaller values of r, an incubation time is observed.

The elastic models were, without a doubt, a step forward in the study of spin transitions. With a more complex and realistic approach than mean-field models, deeper, more flexible, and more accurate than Ising-like models, applicable on a larger scale than density functional theory (DFT) studies, elastic models have led to a new stage in understanding spin transition phenomena. Further developments could include integration with DFT models for a more accurate identification of intrinsic parameters, consideration of geometries identical to those of real samples, and use of more precise elastic potentials. Real structures imply the stabilization of the lattice by taking into account, for example, the torque forces that act on every molecule. Realistic elastic potentials should also prevent the molecules from getting too close (or too far) from one to another.

More than ten years after their initial implementation, elastic models still have much to offer in the study of spin transition compounds and other systems in which long-range interactions predominate.

This work was supported by a grant from the Romanian Ministry of Education and Research, CNCS-UEFISCDI, Project No. PN-III-P4-ID-PCE-2020-1946, within PNCDI III. This paper is dedicated to Professor Dr. François Varret on the occasion of his 80th birthday anniversary as a congratulation to his contribution to the development of the spin-crossover research.

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