A Generalized Ising-like Model for Spin Crossover Nanoparticles

Cooperative spin crossover (SCO) materials exhibit first-order phase transitions in the solid state, between the high-spin (HS) and low-spin (LS) states. Elastic long-range interactions are the basic mechanism for this particular behavior and are described well by the Ising-like model, which allows the reproduction of most of the experimental results in the literature. Until now, this model has been applied with an interaction parameter between the molecules, which is considered to be independent of the states. In this contribution, we extend the Ising-like model to include interaction energy that depends on the spin states and apply it to study SCO nanoparticles. Our research shows that following this new hypothesis, the equilibrium temperature shifts toward higher values.


Introduction
Spin crossover (SCO) nanoparticles [1][2][3][4][5][6][7][8][9][10][11][12] are composed of transition metals with configurations ranging from 3d 4 to 3d 8 embedded in an organic matrix. They are characterized by two different electronic spin configurations, resulting in two stable, magnetic macrostates, termed high spin (HS) and low spin (LS). LS-HS transitions are triggered by an external constraint. Most of the SCO compounds that have been studied are based on the Fe (II) 3d 6 configuration. For the 3d 8 case (Ni 2+ SCO) and following the work of O.A. Qamar et al. [13], the SCO phenomenon is associated with three different lattice distortions, one of them being from a square planar (LS, S = 0) to a tetrahedral geometry (HS, S = 1).
When a Fe (II) cation coordinates in an octahedral ligand field (O h ) as a metal complex, the degenerated orbitals of the cation split into triply degenerate t 2g and doubly degenerate e g states. The macrostates are related to the micro-configurations permitted by the distribution of the six electrons in the t 2g and e g states, mediated by the competition between the pairing electron affinity and the strength of the ligand field, which defines the energy gap between the states. The LS state corresponds to the limiting strong field ligand case, which favors the pairing of the six electrons in the t 2g orbitals (t 6 2g e 0 g and spin S = 0), and the HS state corresponds to the limiting weak field ligand case, which allows the distribution of the electrons on both the t 2g and e g states within the frame of Hund's rule (t 4 2g e 2 g and total spin S = 2). The LS state is thus diamagnetic when neglecting the second-order Zeeman effect, and the SCO appears in a deeply colored state, while the HS state is a paramagnetic, fainter-colored state. On top of that, the Fe-ligand distance is greater in the HS state compared to the LS state.
From Mössbauer studies, the Debye temperature is observed to be higher for the LS phase than for the HS state. Accordingly, this result suggests that the interactions between the metal complex sites, which are in the LS state, are stronger than those pertaining to interacting HS state sites. This result is in line with a greater HS state Fe-ligand distance. Following these experimental results and conclusions, the theoretical work presented in this contribution is based on an Ising-like model [24][25][26][27][28] that is built on three different coupling parameters to simulate the thermodynamic properties of a SCO solid. Thus, three interaction parameters were introduced: J LL when the two interacting sites are in the LS state, J HH when they are in the HS state and J HL when one site is in the HS state and the other is in the LS state. In the first step, to reduce the number of parameters, the J HL was taken as the average value of the J LL and J HH . Then, the effect of the three different interaction parameters on the equilibrium temperature, as well as on the hysteresis width, were analyzed.
This generalized Ising-like model is presented in Section 2, and the Monte Carlo entropic sampling simulations used to calculate the density of the states are developed in Section 3. The results and analysis are summarized in Section 4.

The Model
In the "classical" Ising-like model [24,25], the spin crossover molecules are described by the fictitious operator σ, which has two eigenvalues: +1 assigned to the HS state (with a degeneracy g HS ) and −1 to the LS state (with a degeneracy g LS ). The Hamiltonian for N-isolated molecules is written as: where ∆ is the gap between the fundamental energies of the HS and LS states. Considering the different degeneracies (g HS and g LS ) between the respective HS and LS states arising from their different spin values and phonon spectra, the Hamiltonian is isomorphic to: , where T is the temperature and k B is the Boltzmann constant. To account for the interaction between the molecules, an interaction potential energy, written under an Ising form, is added to the previous Hamiltonian, such that: and considering the presence of short (J) and long-range (G) interactions, as proposed by Linares et al. [26] and demonstrated later in true elastic models by M. Paez-Espejo et al. [29], the Hamiltonian is finally expressed as which re-writes as where the long-range interaction (G) is treated in the mean-field approximation. It is worth remarking that in this contribution, to reduce the parameters, we do not consider different long-range interactions for HS and LS species. We notice in Equation (5) that the interaction term J ij = J for the nearest molecules (four in the square lattice, six in the cubic lattice) is considered to be independent of the spin states in the "classical" Ising-like model.
In this contribution, we generalize the previous Ising-like model and consider that the short-range interaction term J ij between the sites σ i and σ j depends on their spin states. So, three interaction constants are introduced: J HH if both molecules are HS, σ i = σ j = +1; J LL if both molecules are LS, σ i = σ j = −1; J HL if one molecule is HS and the other LS, σ i = −σ j = ±1. For the sake of simplicity, it is assumed that J HL = J HH +J LL 2 . Using the following parameters: Hamiltonian (5) writes as: where the two-sites' correlations s HH , s LL and s LL write: Here, N ++ , N −− and N +− are the numbers of the HS-HS, LS-LS and HS-LS nearest neighbor configurations, respectively, which fulfill the relation N ++ where N is the total number of sites in the considered 2D lattice.
The partition function is then given as a function of these macroscopic variables by: where β = 1 k B T and d(m, s HH , s HL , s LL ), obtained by the Monte Carlo entropic sampling method (MCES, see Section 3), is the number of configurations with the same m, s HH , s HL and s LL values. N is the total number of macrostates characterized by these m, s HH , s HL and s LL values.
The thermal average value of the operator σ, (T) , denoted σ , is given by the following relation: which is related to the high-spin fraction (Nhs) by the expression: Equation (10) is solved using the bisection technique.

Entropic Sampling
The purpose of the Monte Carlo entropic sampling method (MCES) [30][31][32] is to obtain the density of the states d(m, s HH , s HL , s LL ) in an attempt to scan the entire space of configurations. To achieve this, the idea is to introduce in the detailed balance equation: P i the probability of obtaining the macrostate (m, s HH , s HL , s LL ), shown as: Such a choice of P i favors the macrostates with a small density of states and dampens those with a high density.
The detailed balance equation is then expressed as: Using d(m i , (s HH ) i , (s HL ) i , (s LL ) i ) as a bias, a "Monte Carlo step" run and a histogram of the frequency of the macrostates H(m i , (s HH ) i , (s HL ) i , (s LL ) i ) are obtained. Then, the density of states in the next "Monte Carlo step" is improved by the following relation: This technique is used iteratively until an almost constant (flat) The density of states is then used to calculate the thermal average value of the fictitious spin operator.

Numerical Results and Analysis
In the calculations, a set of realistic thermodynamic parameters derived from experimental data on the prototype SCO complex [Fe(btr) 2 (NCS) 2 ], btr = 4,4 -bis-1,2,4-triazole [16] is used. From the enthalpy H measurements, we derived the value of the energy gap as ∆/k B = 1300 K. The equilibrium temperature T eq is deduced from the evolution of ∆H and is fixed at 216.3 K, from which the entropy ∆S = ∆H/T eq is calculated to be ∆S= 50 J/K/mol. Recall that according to the previous model, the expression of T eq is given by T eq = ∆ k B ln(g) . To study the behavior of SCO molecules, the following definitions are used such that T up is the ascending thermal transition temperature, T down is the descending thermal transition temperature and T eq = T 1/2 is the average temperature between T up and T down in which the HS fraction is equal to 1/2. The hysteresis width is defined as ∆T = T up − T down .
In the first part, we study the effect of the interaction parameters on the occurrence of thermal hysteresis and the behavior of its width when it exists. It is worth mentioning that from the experimental point of view, the control of the interactions is usually attempted through metal dilution by replacing a fraction of Fe (II) active SCO atoms with Co, Ni or Zn ions [21,33,34]. Unfortunately, it is usually observed that the effect of the dopant not only breaks down the interactions as expected, but also modulates the ligand field. The results concerning the size effect are given in the second and last parts.

• Effect of the Interaction Parameters
The effect of the interaction is then simulated by monitoring the value of the effective interaction parameter J/k B .
This "average" interaction parameter is defined here as the usual J parameter of the "classical" Ising-like model, and so it relates to J HH , J HL and J LL by the following relation: where s = s HH + s LL − 2s HL .
o The Case x = J HH /J LL = 1.0 Since J HL = J HH +J LL 2 , the case x = 1 corresponds to J HH = J LL = J HL , which will be denoted for simplicity by J. The evolution of the HS fraction, noted as Nhs, is simulated as a function of the temperature corresponding to different values of the average short-range parameter J. Figure 1 illustrates the thermal behavior of the HS fraction for various J/k B values ranging from 12 to 25 K. The curves obtained for the 6 × 6 square system show that the nature of the transition (hysteretic, abrupt or gradual) is closely correlated to the intensity of the short-range interactions.
When the average short-range interaction J/k B is decreased, the thermal hysteresis width reduces from ∆T = 5.41 K (J/k B = 25 K) to ∆T = 0.41 K (J/k B = 16 K) and vanishes below the critical value J/k B 15K. Below this value, the interactions are not strong enough to produce the collective first-order transition accompanied by the hysteretic behavior, and the system's feature changes into a gradual one, in fair agreement with experimental findings [21,33,34] and previous simulations [26].
Moreover, it can be seen in Figure 1 that when all the interaction parameters J HH /k B , J LL /k B and J HL /k B are equal, the transition temperature T eq is independent of the J/k B values and remains equal to T eq = ∆/(k B ln(g)) = 216.3 K.
The previous results are summarized in Figure 2, where we report the phase diagram of the system in the space coordinates T versus J/k B . The presented phase diagram clearly highlights the existence of two regions of thermal behavior: a region corresponding to J/k B ≥ 15 K in which the system presents a single first-order phase transition, while a gradual transition is observed when J/k B < 15 K. The ratio x = J HH /J LL has been set to the value of 0.4, and therefore, the short-range interaction parameters (J HH , J HL , J LL ) are different. In order to simulate the dilution effect, the intensity of the average short-range interaction J/k B is gradually reduced from J/k B = 19 K to J/k B = 11 K.
The obtained curves are plotted in Figure 3 in the case of a 6 × 6 nanoparticle configuration. The values of the thermal hysteresis widths, as well as the equilibrium temperature T eq , are summarized in Table 1 as a function of the J/k B parameter. The inspection of the thermal behavior of the HS fraction in Figure 3 clearly shows a significant difference compared to Figure 1, which remained symmetric around the invariant value of the transition temperature. Here, the change of the interaction parameter results in a large deformation and shift of Nhs(T), which means that the interaction J affects both the cooperativity and the ligand field energy.  Table 1. Values of the thermal hysteresis width ∆T hyst = T up − T down and the equilibrium temperature T eq as a function of J/k B , J HH /k B and J LL /k B in a SCO system with a size of 6 × 6. The computational parameters are: ∆/k B = 1300 K, G/k B = 172.2 K, x = J HH /J LL = 0.4 and ln(g) = 6.01. As can be seen in Table 1, when the interaction parameter J decreases, the values of J HH /k B and J LL /k B also decrease; the gap (J LL − J HH )/k B decreases as well. An important result is that, at the same time, the thermal hysteresis width decreases, and the equilibrium temperature of the system is shifted to lower temperatures.
In Figure 4, we plot the variations of T up , T down and T eq versus J/k B , which leads to the system's phase diagram in the (T, J/k B ) variables' space, highlighting the existence of two regions of thermal behavior. The thermal hysteresis loop (~1.03 K wide for J/k B = 19 K), associated with a single first-order phase transition, rapidly narrows when decreasing J/k B and collapses at the threshold value (J/k B ) c ≈ 16 K. These results agree qualitatively with Mössbauer spectroscopy and magnetic measurements realized on iron (II)-based SCO materials diluted with high-spin cobalt (II) ions [33,34]. There, the Co has a double role: it breaks down the interactions between the Fe (II) SCO species and also affects the ligand field energy through steric effects. Linear regression in the T eq (J/k B ) variation gives: T(J/k B ) = 0.24072 × J/k B + 216.30. Figure 5 shows the thermal evolution of the HS fraction in the case of a 10 × 10 square lattice. The results are also reported in Table 2. The comparaison of the simulation for 6 × 6 and 10 × 10 square lattices are given in Table 3.   The transition temperatures reported in Figure 6, in the heating mode (T up ) and in the cooling mode (T down ), clearly show that the width of the hysteresis loop ∆T = T up − T down increases when the J/k B interaction is stronger. Linear regression leads for the Teq Branch to: T(J/k B ) = 0.26019 × J/k B + 216.30.  Table 4, when the ratio x = J HH /J LL decreases, the difference between J HH /k B and J LL /k B increases, and the J LL /k B value increases while that of J HH /k B decreases. Table 4. Equilibrium temperatures and hysteresis widths calculated for different values of x = J HH /J LL parameter in a 6 × 6 square lattice. The computational parameters are: ∆/k B = 1300 K, J/k B = 14.8 K, G/k B = 172.2 K and ln(g) = 6.01.

J/k B [K]
x The results reported in Figure 7 for a 6 × 6 system show that the decrease in x progressively shifts the equilibrium temperature T eq toward higher values. The equilibrium temperature, which is equal to 216.  In this section, the size effect is analyzed, and calculations are performed in the case of a 10 × 10 square lattice. In the same way, for the 6 × 6 lattice, the value of the ratio x is gradually decreased from the value of 1.0 to the value of 0.2. The results are summarized in Table 5 and in Figure 8 for parameter J/k B set to 14.8 K.  The results of the 10 × 10 lattice are qualitatively similar to those obtained for the 6 × 6 nanoparticle configuration. When the value of the x ratio decreases, the difference between J HH /k B and J LL /k B increases; the interactions between molecules in the LS states (J LL /k B ) are much stronger than those in the HS state (J HH /k B ), and the system remains in an LS configuration for a longer temperature interval. The equilibrium temperature T eq gradually moves toward the order-disorder temperature T O.D. . The hysteresis width decreases and vanishes for weaker J/k B values when T eq becomes greater than T O.D. .

Conclusions
In this study, the behavior of 2D-SCO nanostructures is analyzed with particular attention to the effect of the short-range interactions involved in the model. Indeed, three coupling terms were considered (i) J HH /k B between two HS molecules, (ii) J LL /k B between two LS molecules and (iii) J HL /k B = J/k B between a molecule in the HS state and another one in the LS state with J/k B = (J HH /k B + J LL /k B )/2.
The J HH /J LL ratio effect is parametrized by x. When x = 1, which corresponds to the case where the three short-range interactions are equal, the equilibrium temperature of the system T eq remains constant. The nature of the transition from LS to HS configuration is governed by the intensity of the interactions and, therefore, the value of the J/k B and G/k B parameters. Increasing the short-range interaction parameter J/k B leads to a hysteretic transition. This behavior is explained by the fact that the Curie (or the order-disorder) temperature designed by T O.D. increases. When x is set to a value other than 1, a "pseudodilution" effect is simulated by gradually reducing J/k B . Figures 4 and 6 highlight that the equilibrium temperature shows a linear decrease towards ∆/(k B ln(g)) = 216.3 K. Moreover, the hysteretic behavior vanishes at a threshold value of J c /k B ≈ 16K and J c /k B ≈ 15K for 6 × 6 and 10 × 10 systems, respectively. This feature is connected to the relative positions of T eq and T O.D. and can be explained by the fact that the order-disorder transition T O.D. decreases faster for small lattices. The condition T eq > T O.D. leads to a gradual transition. In the present study, we mainly focused on the effect of varying x when J/k B is constant. An important result is that when the difference between J HH /k B and J LL /k B is increased, the LS state is stabilized and T eq is shifted toward higher temperatures. This behavior is more significant for larger lattices. The effects of the interactions between surface molecules and their environment will also be explored in future work.  0028-02 and the French-Japan international laboratory (LIA IM-LED). We thank all of them for their financial support.

Conflicts of Interest:
The authors declare no conflict of interest.