Hexagonal-Shaped Spin Crossover Nanoparticles Studied by Ising-Like Model Solved by Local Mean Field Approximation

Hexagonal-Shaped Spin Crossover Nanoparticles Studied by Ising-Like Model Solved by Local Mean Field Approximation. Abstract: The properties of spin crossover (SCO) nanoparticles were studied for ﬁve 2D hexagonal lattice structures of increasing sizes embedded in a matrix, thus affecting the thermal properties of the SCO region. These effects were modeled using the Ising-like model in the framework of local mean ﬁeld approximation (LMFA). The systematic combined effect of the different types of couplings, consisting of (i) bulk short- and long-range interactions and (ii) edge and corner interactions at the surface mediated by the matrix environment, were investigated by using parameter values typical of SCO complexes. Gradual two and three hysteretic transition curves from the LS to HS states were obtained. The results were interpreted in terms of the competition between the structure-dependent order and disorder temperatures ( T O.D. ) of internal coupling origin and the ligand ﬁeld-dependent equilibrium temperatures ( T eq ) of external origin.


Introduction
Over the last few decades, the storage capacity and downsizing of electronic components have been of concern for potential industrial applications in accordance with Moore's law. The search for new components scaled to dimensions at the nanometer level and materials more efficient for data storage is a major economic issue that has driven research towards alternatives to current silicon technologies. This quest was behind the proposal that emerged to use a molecule as an active element in an electronic device. Molecular electronics, based on the taming and manipulation of electrical, optical, or magnetic signals in devices composed of molecules, is continuously being developed at an ever-increasing rate [1]. In this context, spin transition materials have shown great potential and aroused much interest [2][3][4][5]. Indeed, these compounds have two different stable magnetic states termed high spin (HS) and low spin (LS). They can switch reversibly and while under control from one spin state to another. This transition is generally caused by a change in temperature [6][7][8][9]. It has been shown that it is also sensitive to a wide spectrum of external stimuli such as pressure [3,[10][11][12], visible light [13][14][15][16][17][18], and magnetic field [19].
The change in spin state is accompanied by significant changes in the color, volume, magnetic state, and electrical conductivity of the compound, which can be observed by using different characterization techniques (such as magnetometry, Mössbauer, optical spectroscopy, X-ray diffraction, and calorimetry) [20][21][22][23][24].
The phenomenon of molecular bistability can occur for a certain number of metal ions belonging to the first series of transition metals and, more particularly, for metals with the 3d 4 -3d 8 configurations. For example, Fe (II) [25,26], which is probably the most studied ion to date, has a 3d 6 configuration in its free ion state.
When the Fe (II) ion is coordinated in an octahedral ligand field (O h ) to form a metal complex, the degeneracy of the d orbitals is lifted into the t 2g triply degenerate and e g doubly degenerate states. The difference between these two levels depends on the strength of the ligand field, and the six electrons can then be distributed in two different ways: (i) in the case of a strong field, the six electrons pair up in the t 2g orbitals, which leads to an LS diamagnetic and colored state for which the spin S = 0 and for which the electronic configuration is t 6 2g e 0 g ; (ii) in the case of a weak field, the electrons are distributed over the d orbitals according to Hund's rule, leading to a HS paramagnetic and colorless state for which the spin S = 2 and for which the electronic configuration is t 4 2g e 2 g . This second configuration results in an expansion of the Fe-ligand distance with respect to the LS state, and it is therefore accompanied by an increase of the unit cell volume of the material [27].
According to the type (hydrogen bonding, electrostatic, or pi-pi stacking) and intensity of interactions between molecules, spin transition compounds can exhibit different kinds of transitions with a change in temperature: continuous or gradual transformation while respecting the simple Boltzmann-type statistic, sharp first-order with or without hysteresis cycle, two or multi-steps transitions, or incomplete transitions with a quenched residual HS fraction at a low temperature. These five behaviors have been observed in Fe (II) compounds, and they can occur at low, high, or room temperature. First-order phase transitions accompanied by hysteresis are observed in highly cooperative systems in which intermolecular interactions are strong enough. Historically, König et al. [28] were the first to observe a hysteresis loop in 1976 on the thermally-induced spin transition of the [Fe (4, 7-(CH3) 2 -phen) 2 (NCS) 2 ] compound. Later on, a two-step spin transition was revealed in a Fe (III) complex [29]. It is worth mentioning that during the last ten years, SCO Fe(II) complexes have been synthesized that are characterized by spin transition occuring in three steps [30][31][32].
To try to explain and reproduce these switching phenomena from the LS to HS states, researchers have proposed different models based on Ising-like [33][34][35][36] and atomphonon [37] descriptions, and, more recently mechano-elastic [38] and electroelastic [39] models have been developed.
A crucial issue concerning SCO nanoparticles is the reduction of size. Indeed, when the surface-to-volume ratio increases, it results in an increasingly large effect on the thermal properties of SCO nanoparticles dominated by interactions between the surrounding environment and the molecules located on the edges (or surfaces).
Experimentally, it has been observed that a drastic change in the properties and, particularly, in the cooperativity of the nanoparticle can occur when the size is smaller than a threshold "limit" size. For example, in the context of a thermo-induced transition, this can lead to a change in the transition temperature of several tens of Kelvins.
In the context of Ising-like models and 1D SCO systems embedded in matrixes, a theoretical study based on an Ising-like model [33,34,40] was presented by Chiruta et al. [41] to understand the physical origins of the three-step behavior. Three-state thermal behavior was also found by S.E. Allal et al. [42]. In 2D SCO compounds, S.E. Allal et al. [43,44] studied two-step transitions and analyzed the conditions to obtain three-state behavior.
Recently, we proposed a homogenous local mean-field approach to study both shape and size effects on 3D SCO cubic and cuboid lattices [45]. This study highlighted three steps of switching and three states under the effect of temperature, as well as two transition stages and three states under the effect of pressure.
In the present study, we focused on the modeling of thermal effects in 2D hexagonalshaped SCO nanoparticles. The Ising-like model was used and solved in the framework of local mean field approximation (LMFA) in an attempt to reveal the role of interactions with the external environment (matrix). In a hexagonal-shaped lattice, the molecules have more interactions than in a square lattice. This means that due to the lattice topology, all interactions are locally stronger for a hexagonal symmetry than a square symmetry. In this study, we show that interactions with the environment are important, even in a hexagonal-shaped model. As a result, this symmetry allows one to enhance the effects of the interactions and surfaces on the physical properties of the studied nanoparticles. This paper is organized as follows. Section 2 is devoted to the model and the calculation method, Section 3 presents the results and a discussion, and Section 4 is the conclusion.

Model
The first microscopic model simulating the spin transition was introduced by Wajnflasz and Pick [33], who only considered short-range interactions between the spin-crossover sites. This model was reconsidered by Bousseksou et al. [34] in the mean-field approach involving two "antiferromagnetically" coupled sublattices to reproduce a double-step spin transition. Later, Linares et al. [40] added long-range interactions to the Ising-like model to trigger the phenomenon of hysteresis in 1D compounds. Chiruta et al. [41] included an energetic contribution (L). In SCO nanoparticles, this term allows one to explicitly deal with interactions between the molecules located on the periphery and in close contact with the environment.
As a result of these previous contributions, we considered short-range, long-range, and matrix effects. The Hamiltonian of the global system for SCO compounds can be written as follows: The first term of Equation (1) corresponds to the temperature-dependent field. The ligand field ∆ (> 0) is the energy gap between the HS and the LS states that depends on the symmetry and nature of the ligands involved in the complex. The total number of molecules is N T . T is the absolute temperature of the system, k B is the Boltzmann's constant such as β = 1/(k B T), and g = g HS /g LS is the ratio between the degeneracy of the HS state (g HS ) and the degeneracy of the LS state (g LS ). A fictitious spin operator σ is associated with each molecule; its eigenstates are +1 when the molecule is in the HS state and −1 when the molecule is in the LS state [11,46,47]. The second term accounts for the strength of short-range interactions, J, limited to first-neighbors. The third term, G, pertains to the long-range part of the interaction, and the last term, L, equivalent to an applied pressure on the molecules at the surface, is a coupling term related to the interactions between the environment and the atoms located at the corners and at the edges that vanishes for bulk molecules. The J and G parameters are ferromagnetic-type (>0) interactions and thus favor HS-HS and LS-LS pairs. The ligand field of small spin-crossover nanoparticles may continuously vary between the surface and the core of the nanoparticles. Obviously, the ligand field is smaller at the surface due to the immediate environment of SCO sites, which might be composed of water molecules or event impurities (pollution). In contrast, in the bulk, all SCO metal's centers (iron (II)) are surrounded by nitrogen atoms that locally enhance their ligand field. One may expect a variation of the average ligand field of the nanoparticles as a function of size. To simplify, we only considered the change of the ligand field at the surface in this model and assumed that it remained unchanged inside the lattice or that this variation was negligible at the nanometer scale.
In this study, three types of sites were considered: the atoms located in the bulk (N b ) that are surrounded by six first-neighbors and those located on the edge (N e ) and the corner (N c ) that only interact with four and three first-neighbors, respectively. The molecules located on the surface (edge and corner) therefore have specific properties because they also interact with an external environment (matrix effect).
The case of a system comprising 19 molecules (designated H3) is represented in Figure 1 below. Considering the three regions in the lattice (corner, edge, and bulk) led us to distinguish three local average order parameters σ α associated with each of these regions, which brought us to consider a local Hamiltonian for each region α with α = b, c, and e corresponding to the bulk, corner, and edge, respectively, in the local mean field approximation.
The number of interactions, J, between a molecule and its first-neighbors is denoted as q α , and the number of interactions, L, between a molecule and the external environment is referred by z α .
For each region, the Hamiltonian can be cast as Equation (3) below: with The characteristics of the various local situations considered in the model are listed in Table 1.
The partition function of this inhomogeneous mean-field system comprising N b , N c , and N e atoms belonging to the bulk, corner, and edge regions, respectively, is written as: The partition functions Z b , Z c , and Z e -connected to the bulk, corner, and edge, respectively-have the general form: The average value of the fictitious spin of the system is calculated as a weighted average of the three order parameters of the three regions: bulk, corner, and edge.
The high-spin fraction, Nhs, which is the probability that the HS state is occupied, is given by: Using the Gibbs-Bogoliubov-Feynman inequality, a variational approximation of the upper bound Helmholtz free energy could be derived, such that the total mean-field free energy per site is written as: with Then: Thus, for each region of b, c, and e, the average spin state is: which is equivalent to: These equations are combined with the relation (Equation (1)): The problem therefore consisted of finding the three order parametersσ b , σ c , and σ e -that were simultaneous solutions of the three functions f 1 , f 2 , and f 3 .
This system was solved by the Newton method, which is appropriate for rapid convergences. Calculations were initiated from three points (x 1 , x 2 , and x 3 ) close to the solution and that represented, respectively, σ b , σ c , and σ e . Solving the system amounted to calculating 1 , 2 , and 2 so that: which allowed us to write the following equation in the vicinity of points x 1 , x 2 , and x 3 : Thus: The following system of three equations with three unknowns was thus obtained.
and 1 , 2 , and 3 were calculated as follows: In a general way, the f i functions with i = 1 (bulk), 2 (corner), 3 (edge) and their derivatives can be expressed as follows:

Numerical Results and Discussion
In [48], we showed that the size dependence of the equilibrium temperature T eq of such a system could be described analytically, and, more precisely, we reached the conclusion that the equilibrium temperature was the result of a null total effective ligand-field. As a result, and according to the expressions listed in Table 1, it was deduced that the transition takes place at the equilibrium temperature, T eq , thus obeying the following constraint: which led us to the following expression of the equilibrium temperature: Here, T bulk eq , T edge eq , and T corner eq are the respective equilibrium temperatures of the bulk, edge, and corner whose analytical expressions are: This means that, for an infinite size system (N T N b , since N b >> N c and N b >> N e ) with L/k B and ln(g) being constant, T eq tends towards the equilibrium temperature of the bulk T b eq = ∆/k B ln(g) . Equation (23) shows that the equilibrium temperature of the system has two limiting values, namely T corner eq = ∆−6L k B ln(g) and T bulk eq = ∆ k B ln(g) . It also clearly shows that the equilibrium temperature decreases when the number of interactions with the matrix (L term) increases and, therefore, when the size of the lattice decreases. This point is fundamental in understanding and predicting the occurrence of first order transitions with hysteresis loops. Indeed, in a pure Ising model, the order-disorder (or Curie) temperature, T O.D. , is obtained from ∆/k B = 0, L/k B = 0, and g = 1 (ln(g) = 0) in the Hamiltonian of Equation (1). First order transitions with hysteresis occur when T O.D. > T eq ; otherwise, only a gradual transition takes place.
In the following section results of simulations for 2D hexagonal lattices composed of 19, 37, 61, 91 and127 molecules and designated as H3, H4, H5, H6 and H7 respectively are discussed. The ratio parameter r was defined as r = (N c + N e )/N T , where N T = N b + N c + N e is the total number of molecules (Table 2 below).

Case L = 0
First of all, only the short-and long-range interactions of J and G were considered in the calculations, which consisted of setting L/k B = 0 K. The thermodynamic parameters used in the calculations were values typical of the Fe(btr) 2 (NCS) 2 SCO complex [4]: the enthalpy and entropy changes associated with the spin transition were around ∆H ≈11 kJ/mol and ∆S ≈ 50 J/mol/K, which led to ln(g) = ∆S/R ≈ 6.01 (where Ris the perfect gas constant) and ∆ = ∆H/R ≈ 1300 K. The resulting equilibrium temperature value of the system was T eq = ∆ (k B ln(g)) ≈ 216 K (in Figure 2, T eq = 216.3 K), and according to Equation (24), it also corresponded to the bulk transition temperature. At this temperature, there were as many molecules in the LS state as there are molecules in the HS state, and the HS molar fraction was Nhs = 1 2 . The calculations were carried out for systems H3-H7 (19-127 atoms) and are reported in Figure 2.
As can be seen in Figure 2, the equilibrium temperature T eq was independent of the system's size. This result was directly linked to the absence of interactions with the surrounding matrix. On the other hand, different behaviors were observed with regard to the type of spin transition. Though the transition was gradual for the H3 and H4 lattices, a thermal hysteresis loop appeared for the H5, H6, and H7 systems.
Moreover, the width of the thermal hysteresis ∆T = T up − T down increased as the ratio r = N s /N T decreased. In Figure 3, the plot of T against the size H i clearly reveals a linear change of the thermal hysteresis width as a function of size.
These results are explained by the fact that by increasing the size, the coupling between the molecules was increased and the order-disorder transition temperature T O.D. of the hexagonal system (or Curie temperature) therefore increased. For H5, H6, and H7, T eq < T O.D. , which led to a first-order phase transition. By increasing the ratio r (See Table 2), i.e., by increasing the number of molecules on the surface, T O.D. decreased. The condition T O.D. < T eq was also satisfied, and a gradual transition took place for H3 and H4. The hysteretic behavior was therefore strongly impacted by the decrease in the number of molecules of the SCO compound.

Case L = 0
After considering the strength of interactions between the boundary molecules and the matrix, the 2D SCO went from gradual to three-step behaviors. It is interesting to remark that a two-step behavior with three states was also obtained for the same temperature

Three Steps
The value of the L/k B parameter was gradually increased from 0 to 160 K, and calculations were made for the H6 (91 atoms) system It is worth mentioning that modifying the values of L/kB correspond to a change in the interaction strength between the nanoparticle and its environment (matrix effect) which is achieved with a change in the chemical nature of the matrix. The results are reported in Figure 4. Table 3 provides the values of the transition temperatures of the different regions (corner, edge, and bulk) and the overall system for different coupling strengths with the matrix, as derived from the numerical simulations. Table 3. Equilibrium temperatures of the corner, edge, and bulk calculated for different values of the L parameter for the H6 system (91 atoms). The computational parameters were ∆/k B = 1300 K, J/k B = 40 K, G/k B = 70 K, and ln(g) = 6.01.

L/ k B [K]
T As shown in Figure 4, the value L/k B = 0 K led to a single thermal hysteresis loop at the equilibrium temperature T eq = T bulk eq ≈ 216.3 K. For a slightly higher value of L/k B = 50 K, boundary effects increased and a two-step transition emerged. A complete three-step transition was achieved for strong coupling with the environment, such as for L/k B = 160 K. Figure 5 below highlights the three steps associated with three hysteresis loops HYST1, HYST2, and HYST3 obtained in the heating and cooling modes and that were distinct in the sense of the absence of overlapping hysteresis loops. The stable states were: (i) for T < T , the states with Nhs close to 0; (ii) for temperature values between T and T", the states with Nhs close to 0.06; (iii) for the region between T" and T"', the states with Nhs close to 0.31; and (iv) for T > T"', the states with Nhs > 0.8.  We note that, on the one hand, increasing the value of the L/k B parameter increased the width of the HYST3 and HYST2 hysteresis loops; on the other hand, a concomitant shift towards lower temperatures led to the presence of two well-defined intermediate plateaus.
The first hysteresis loop, HYST3, was in the temperature range of 67-80 K and, as can be seen in Table 3, was close to the transition temperature of the corner, T c eq = ∆−6L k B ln(g) ≈ 57 K, with a plateau at Nhs ≈ 0.06. This means that at this temperature, a very small number of molecules switched to the HS state. The second hysteresis, HYST2, in the temperature range of 114-130 K, was close to the temperature transition of the edge, T e eq = ∆−4L k B ln(g) ≈ 110 K, with a plateau at Nhs ≈ 0.31. As for HYST1, between 197 and 218 K, it corresponded to the bulk transition temperature, T b eq = ∆ k B ln(g) = 216.3 K. Figure 6 exhibits two hysteresis loops. The first one, denoted HYST1, was in the temperature range of 132-236 K, with Nhs = 0.32-1.0, and the second one, HYST2, was in the temperature range of 205-217 K, with Nhs = 0.15-0.32. The two hysteresis loops overlapped over a temperature interval of T 2 − T 1 ≈ 12 K, which led to the presence of three stable states for the same temperature. These states could be achieved through an appropriate variation of temperature or some other physical stimuli such as light. The intermediate plateau obtained between 205 and 236 K with Nhs between 0.27 and 0.40 was a mixture of the LS and HS configurations, with the majority of molecules in the low-spin state. The branches (blue crosses in Figure 6) between Nhs = 0.15-0.26 and Nhs = 0.40-0.91 are the unstable states that were not experimentally observable.

Three States
Note that these results were obtained by explicitly considering the edge and corner effects at the surface. These results were different from those presented in [11,48], for which the latter effects were not explicitly considered and which led to an absence of the overlapping of hysteresis loops.

Phase Diagram
The previous results are summarized in Figure 7 by plotting the phase diagram of the system in the space coordinates, temperature versus L/k B . The transition lines corresponding to the limiting branches of the thermal hysteresis when the transitions were of first-order or to the equilibrium temperature when the transitions were continuous (gradual) are represented. The present phase diagram was made for the lattice size H6 (Figures 4 and 5) and has the advantage of clearly highlighting the existence of four regions of thermal behavior as a function of the surface parameter L. Thus, in region (i) corresponding to L/k B < 30 K, the system presented a single first-order phase transition (equilibrium transition temperature: T eq1 ∼ T 1up +T 1down 2 ), while in region (ii), obtained for 30 < L/k B < 90, the transition from LS to HS involved a gradual transition (equilibrium transition temperature: T eq2 ) followed by a hysteretic first-order transition of region (i). Above L/k B = 90 K and below L/k B = 130 K, i.e., region (iii), the previous continuous transition became first-order with limiting transition temperatures, T 2up and T 2down , and a new bifurcation appeared, leading to the emergence of a third transition line of a gradual phase transition with equilibrium temperature T eq3 . In the fourth region (iv), i.e., for L/k B > 130 K, all previous phase transitions became of the first order, and then the system exhibited three consecutive first-order transitions. Remarkably, the phase diagram of Figure 7 shows that as the L/k B value decreased, the widths of the hysteresis loops (T 3up − T 3down ) and (T 2up − T 2down )-respectively associated with HYST3 and HYST2 and corresponding to the corner and edge contributions-linearly decreased with L/k B , while that of the bulk (T 1up − T 1down ) remained unchanged except for very weak L/k B values. Overall, the different regions of this phase diagram confirm that nano-switchable systems combining competing interactions between surface environment and bulk effects may lead to rich thermodynamic behaviors, going from a single first-order transition (region (i)) to consecutive multi-step first-order phase transitions (region (iv)).

Figure 7.
Phase diagram T = f (L/k B ) for a 2D hexagonal-shaped system comprising 91 atoms (H6). The red and blue squares correspond, respectively, to the upper and lower transitions for the heating (T up ) and cooling (T down ) temperatures of the thermal HS fraction. The black squares correspond to the equilibrium temperatures (T eq ) of the gradual transitions. The computational parameters were ∆/k B = 1300 K, J/k B = 40 K, G/k B = 70 K, and ln(g) = 6.01.

Conclusions
Local mean field approximation was applied to study the properties of SCO nanostructures using the Ising-like model. The results showed the possibility of three types of transition for the five 2D hexagonal systems that were studied: continuous, gradual, or with two or three hysteresis cycles in three or four steps, with overlapping cycles in the case of three steps. These different properties, as functions of temperature, were interpreted in terms of the competition between internal interaction coupling, as modeled by the J and G parameters, and external interaction coupling, as modeled by L. The effects of these terms were dependent on both the size and ratio r, and the overall effect was a different shift of the order-disorder transition temperature T O.D. and the equilibrium temperature T eq . These results are to be compared to the re-entrance effect that impacts the type of transition that occurs as a result of this competition. Finally, a phase diagram in the variables of T v/s L/k B clearly highlights the existence of four regions of thermal behavior as a function of the surface parameter L/k B .