A First Order Phase Transition Studied by an Ising-Like Model Solved by Entropic Sampling Monte Carlo Method

: Two-dimensional (2D) square, rectangular and hexagonal lattices and 3D parallelepipedic lattices of spin crossover (SCO) compounds which represent typical examples of ﬁrst order phase transitions compounds are studied in terms of their size, shape and model through an Ising-like Hamiltonian in which the ﬁctitious spin states are coupled via the respective short and long-range interaction parameters J , and G . Furthermore, an environmental L parameter accounting for surface effects is also introduced. The wealth of SCO transition properties between its bi-stable low spin ( LS ) and high spin ( HS ) states are simulated using Monte Carlo Entropic Sampling (MCES) method which favors the scanning of macro states of weak probability occurrences. For given J and G , the focus is on surface effects through parameter L . It is shown that the combined ﬁrst-order phase transition effects of the parameters of the Hamiltonian can be highlighted through two typical temperatures, T O.D . , the critical order-disorder temperature and T eq the equilibrium temperature that is ﬁxed at zero effective ligand ﬁeld. The relative positions of T O.D. and T eq control the nature of the transition and mediate the width and position of the thermal hysteresis curves with size and shape. When surface effects are negligible ( L = 0), the equilibrium transition temperature, T eq . becomes constant, while the thermal hysteresis’ width increases with size. When surface effects are considered, L (cid:54) = 0, T eq . increases with size and the ﬁrst order transition vanishes in favor of a gradual transition until reaching a threshold size, below which a reentrance phenomenon occurs and the thermal hysteresis reappears again, as shown for hexagonal conﬁguration. 4 × 4 6. These results show that the MCES method, is well adapted to investigate the thermal properties of spin transitions in the 2D and 3D SCO compounds although it is limited to small sizes. Future work will concern cubic centered and faced cubic centered structures and iron (II) spin-crossover binuclear compounds.


Introduction
The quest to improve our knowledge on phenomenon and processes driving the properties of complex materials is challenging for researchers involved not only in experimental but also in theoretical research fields. The term complex is used here to denote molecular systems subjected at the nanometer scale to the influence of a surrounding media. Modeling and simulating physical processes can be achieved either from analytical and or from numerical methods. Whereas the former can solve exactly both finite and infinite model systems, the latter is more appropriate to study finite complex ones. Monte Carlo methods are appropriate to study stochastic systems and as such have been applied to spin crossover (SCO) materials.
Concerning SCO materials [1][2][3][4][5][6][7][8][9][10][11][12], a typical example of a first-order phase transition, the Monte Carlo methods are used in combination with certain models such as: Ising-like model, atom-phonon coupling (APC) model or mechanical-elastic model. All these models are used to simulate the response of SCO materials when they are subjected to external perturbations such as: thermal or/and pressure variation, light irradiation, applied external magnetic or electrical fields. The importance granted to these materials are based on their potential application in displays, nanoelectronics, data storage or as temperature or/and pressure or light sensors. These materials which are characterized by two stable states, a diamagnetic low spin (LS) with a degeneracy g LS state and a paramagnetic high spin (HS) state with a degeneracy g HS (> g LS ). At the solid state, the interactions between the molecules are at the origin of a thermal hysteresis, that is the fingerprint of a first order phase transition.
In this contribution, we apply the Monte Carlo entropic sampling (MCES) method, described in Section 3, to solve the Hamiltonian associated to the Ising-like model, that will be described in Section 2. The results and the discussion are presented in Section 4, and we conclude in Section 5.

Spin Crossover Phenomenon with a Typical First-Order Phase Transition
The spin transition phenomenon was discovered in 1931 by Cambi et al. [13,14]. They revealed an "abnormal" magnetic behavior during the study of a series of Fe(III) compounds based on dithiocarbamate ligand. Later, in 1963, Madeja et al. [15] presented results of measurements of magnetic susceptibilities of the series of compounds [Fe(II)phen 2×2 ] (phen = 1,10-phenanthroline; X − = Cl − , Br − , I − , N 3 − , SCN − , SeCN − ). In 1964, Baker et al. [16] reported the first spin transition of a solid state Fe(II) complex with the [Fe(phen) 2 (NCS) 2 ] complex. The same year, Ewald et al. [17] introduced the notion of spin-crossover.
Spin transition (ST) or spin crossover (SCO) complex magnetic systems form a class of compounds that are thus characterized by two spin states, low spin (LS) and high spin (HS). They are bistable molecular materials having the properties to switch between the LS state and the HS state when they are subjected to an external physical constraint such as pressure, temperature, electric or magnetic field or light-matter interaction. The corresponding physical magnetic, optical, thermal, electrical and mechanical properties depend on the change of the spin state and consequently can be used as pressure, temperature, light or electromagnetic field sensors. Computational physics to simulate these properties has thus emerged as a field of research based on a Hamiltonian operator that accounts for the different coupling scheme present in the complex SCO compounds and a pertinent numerical algorithm.
Wajnflasz and Pick [18] model, developed in the 70s, the so called "Ising-like" model, by the introduction of a fictitious spin (σ = +1 (HS) and σ = −1 (LS) to describe the spin states and a parameter J that accounts for the interactions between the molecules, may be regarded as the initial step towards a comprehensive modeling of the effects studied in this work. The total Hamiltonian was solved in the mean field approximation.
The Hamiltonian associated to the Wajnflasz and Pick model writes as follows: ∆ (>0) is the difference between the fundamental energy of the HS and LS states. N is the total number of molecules.
In the mean field approximation and taking account the different degeneracies g HS and g LS of the two states, the Hamiltonian (1) can be written as: where k B is the Boltzmann constant, g = g HS g LS , z is the number of nearest neighbors and σ is the thermal average value of σ.
This genuine model (2) has a very simple phase diagram. The transition temperature, T eq , of the system, defined as the temperature at which a half of the (HS) spins converts Symmetry 2021, 13, 587 3 of 14 to (LS), is obtained by cancelling the effective ligand field, ∆ − k B T ln(g), which gives T eq = ∆ k B ln (g) , and the transition is of first-order when T C = J k B > T eq and continuous (gradual) otherwise.
It is worth mentioning that the present description of the SCO phenomenon in terms of Ising-like model is only phenomenological. More realistic descriptions, based on elastic models [19][20][21][22] taking into account the volume change at the transition, demonstrated that this type of models is isomorphic with an Ising-like model combining short-and longrange interactions [23].
These interactions, represented by the coupling parameters, J and G in the present extended Wajnflashz model of Equation (3), are due to the elastic strain field propagating along the lattice through acoustic phonons. The long-range interaction contribution is assumed here, for simplicity, as uniform and of infinite-range. This type of interaction, is known to be equivalent of mean-field model, which then justifies the mean-field treatment in Equation (1). Therefore, to go beyond model (1), and in order to reproduce the experimental results of thermal hysteresis shapes, it is important to combine in the same Hamiltonian both long-and short-range interactions.
Among these extensions, in 1999 Linares et al. [24] analyzed the shape of thermal hysteresis in 1D SCO compounds, and proposed significant improvement of model (1) by integrating a long-range interaction term (G) along with a short-range one (exchangelike, J), in the extended 1D Ising-like Hamiltonian, and studied analytically the conditions to obtain first and gradual spin transitions. One year later, the dynamical version of this model was investigated by Boukheddaden et al. [25].
Here, we are interested in the equilibrium thermodynamic properties of the 2D version of this model, whose Hamiltonian writes: where, i, j means that only the term, J (short-range interaction), connects the nearest neighbors sites. As in the previous case, the infinite long-range interaction, expressed through σ in the Equation (3) allows, by self-consistent numerical method, to reproduce in specific conditions, the thermal hysteresis observed in spin-crossover compounds.
To take into account the interaction between the molecules at the surface and the external medium, an additional energetic contribution to the ligand field (L term) has been introduced [20][21][22].
where σ k stands for the spin state of the M molecules at the surface. As previously shown [26][27][28], the parametric form of the total Hamiltonian is expressed as: where: is the effective energy gap between (HS) and (LS) states, and the macroscopic variables m, s and c are defined as follows: From the total number of configurations, 2 N , there are d(m, s, c) configurations with the same m, s and c values. The purpose of the Monte Carlo-entropic sampling method (MCES) is to scan these "density of states" d (m, s, c). Let's introduce N (< 2 N ) the number of different macro-states characterized by the same values of (m, s, c).
The partition function is given by: and β = 1 k B T . The thermal average value of the operator σ is given by the following relation, which is related to the high-spin fraction (Nhs) by the expression: Moreover, by assuming that the transition temperature T eq of the system corresponds to a null total effective ligand-field, then it can be shown that: where T bulk eq and T sur f eq are the bulk and surface transition temperatures respectively. Note that M depends on the type of configuration of the nanoparticles.

Monte-Carlo Entropic Sampling Method
Monte Carlo methods are frequently used in solid state physics to study the equilibrium properties of matter at the macroscopic scale and, to a lesser extent the nonequilibrium properties close to equilibrium. The properties of metastable states remain a challenge to the theoretical physicist. For example, by construction, the Monte Carlo Metropolis algorithm scans the most probable configurations, which are the configurations that shape the Boltzmann equilibrium thermodynamic trajectory in the phase space.
To obtain the density of states, the entire space of configurations must be visited, which is achieved by the entropic sampling method (MCES) [29][30][31]. This method is based on the idea that any desired distribution P can be obtained from a Monte Carlo method, if this distribution appears in the detailed balance equation that writes as: where W i→j is the transition probability from state (i) to state (j). In this contribution, the "states" in Equation (14) are the macro-states characterized by (m i , s i , c i ) and the detailed balance equation is given by the relation, In the framework of the entropic sampling algorithm, the following probability is introduced such that: Such a choice of P i favors the macrostates with small density of states and dampens those with high density.
Using d(m i , s i , c i ) as a bias, a "Monte-Carlo step" run and a histogram of the frequency of the macrostates H(m i , s i , c 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) H(m i , s i , c i ) histogram is obtained. The density of states d(m i , s i , c i ) is then used to calculate the partition function and the thermal average value of the fictitious spin operator.

Results and Discussion
In a first stage, the thermal dependence of the (HS) fraction obtained for 2D configurational topology of nanoparticles is presented and the size, shape and surface (L) effects under temperature are analyzed for square, rectangular and hexagonal shaped lattices. Next, the shape effects and the influence of a ligand-field (L) are studied for 3D compounds.  (6), the set of thermodynamic parameters ∆/k B = 1300 K (enthalpy change) and ln(g) = 6.01 (entropy change) leads to (∆/k B )/ ln(g) ≈ 216.3 K which therefore corresponds to the equilibrium temperature T eq for the 3 × 3 to 12 × 12 lattices, as depicted in Figure 1. Additionally denoted by T 1/2 , it corresponds to an equal distribution of the configurations in the state (BS) and in the state (HS). As can also be seen in Figure 1, the change of spin state is not of the same type depending on the size of the lattice. The 4 × 4 system exhibits a sharp spin transition which is the imprint of a more cooperative system, in which each metal center interacts more As can also be seen in Figure 1, the change of spin state is not of the same type depending on the size of the lattice. The 4 × 4 system exhibits a sharp spin transition which is the imprint of a more cooperative system, in which each metal center interacts more strongly with the others and is influenced by the change in the spin state of its neighbors. A hysteresis loop, typical of a first-order transition, is associated with larger systems 5 × 5, 8 × 8, 10 × 10 and 12 × 12, which are highly cooperative. By increasing the size of the nanoparticle configuration, the coupling between the molecules strengthens and the orderdisorder transition temperature T O.D. of the system (or Curie temperature) increases. At the onset of the following condition T eq < T O.D. , a hysteretic behavior appears. A zoom around the transition temperature, in Figure 2, highlights this phenomenon.     By adding the interaction term L between surface molecules and the environment, it is wise for analysis to distinguish different types of sites in the simulations. In fact, molecules located in the bulk (N b ) and on the surface (N s ) do not have the same number of first-neighbors (q) and the same number of interactions with the environment (z). Then, the expressions of surface and bulk transition temperatures (N = N tot = N b + N s and M = N s ), can be cast as: , T sur f eq = ∆ − 2 L k B ln(g) and T eq = N b T bulk eq + N s T sur f eq N tot (18) Taking into account the interaction term L, a competition appears between T eq and T O.D. . When the lattice size is reduced, the total number of interactions between molecules decreases and the order-disorder temperature decreases. At the same time, the equilibrium temperature decreases. As can be seen in Figure 4, for the 4 × 4 compounds, the equilibrium temperature T eq decreases faster than the order-disorder temperature T O.D. , which leads to the appearance of a first-order transition with a hysteresis loop. This behavior is similar to that reported by Peng Taking into account the interaction term L, a competition appears between Teq and TO.D.. When the lattice size is reduced, the total number of interactions between molecules decreases and the order-disorder temperature decreases. At the same time, the equilibrium temperature decreases. As can be seen in Figure 4, for the 4 × 4 compounds, the equilibrium temperature Teq decreases faster than the order-disorder temperature TO.D., which leads to the appearance of a first-order transition with a hysteresis loop. This behavior is similar to that reported by Peng et al. [32] for a [Fe(pz) {Ni (CN)4 }] SCO nanoparticle.

2D Rectangular Lattice: Shape Effects under Temperature
In this section, the thermal-dependence of the (HS) fraction is analyzed for different lattice shapes comprising 144 molecules. We start with a 12 × 12 square lattice and the length of one side is gradually increased at the expense of the other side, until a 2 × 72 rectangular-shaped lattice is obtained. The ratio t between the number of molecules on the surface and the total number of molecules is calculated for each lattice configuration. The values of the t parameter are gathered in Table 1 and Figure 5 shows the curves obtained for the different shapes 12 × 12, 8 × 18, 6 × 24, 4 × 36, 3 × 48 and 2 × 72. It is worth mentioning that in the case of the 72 × 2 lattice, all the molecules are located at the surface which is equivalent to a value of the parameter t equal to 1 and the width of the thermal hysteresis ∆T = T up − T down is then maximum. It gradually decreases and disappears for the case 6 × 24 and 8 × 18. The thermal hysteresis reappears for the square lattice 12 × 12, indicating the presence of a re-entrance phenomenon. Table 1. Evolution of the width of the thermal hysteresis ∆T = T up − T down as a function of the shape of the 2D system and of the ratio parameter t corresponding to Figure 5. The computational parameters are: ∆/k B = 1300 K, J/k B =13.0 K, G/k B = 179.7 K, L/k B = 120 K and ln(g) = 6.01.

Shape
T down (K) mentioning that in the case of the 72 × 2 lattice, all the molecules are located at the surface which is equivalent to a value of the parameter t equal to 1 and the width of the thermal hysteresis ∆ = − is then maximum. It gradually decreases and disappears for the case 6 × 24 and 8 × 18. The thermal hysteresis reappears for the square lattice 12 × 12, indicating the presence of a re-entrance phenomenon. Table 1. Evolution of the width of the thermal hysteresis ∆ = − as a function of the shape of the 2D system and of the ratio parameter t corresponding to Figure 5. The computational parameters are: Δ/kB = 1300 K, J/kB =13.0 K, G/kB = 179.7 K, L/kB = 120 K and ln(g) = 6.01.

Case of Triangular Interactions
In this section, the case of nanoparticles with a triangular interaction leading to hexagonal-shaped systems is explained. Each molecule located inside the bulk has, as it can be seen in Figure 6, six nearest-neighbours (nn) while those situated at the external surface have three or four nn.

Case of Triangular Interactions
In this section, the case of nanoparticles with a triangular interaction leading to hexagonal-shaped systems is explained. Each molecule located inside the bulk has, as it can be seen in Figure 6, six nearest-neighbours (nn) while those situated at the external surface have three or four nn. Five different sizes of these systems with 19, 37, 61 and 91 molecules have been studied. They are characterized, as it is shown in Table 2, by quite different values for the ratio t. The latter ratio is important to highlight the influence of the environment, correlated to Five different sizes of these systems with 19, 37, 61 and 91 molecules have been studied. They are characterized, as it is shown in Table 2, by quite different values for the ratio t. The latter ratio is important to highlight the influence of the environment, correlated to the L interaction term, for the molecules at the surface. Table 2. Relationship between the size of the system and the ratio parameter t between the surface atoms N s and the total number of atoms N tot in the case of a hexagonal-shaped lattice. Typical values from experimental results of the parameter ∆ (the energy difference between the fundamental energy level of the (HS) and (LS) states) and of g = g HS /g LS have been used: ∆/k B = 1200 K and ln(g) = 6.5.

Size of the System
Firstly, the interactions between the molecules at the surface and the environment are discarded, which amounts to consider L/k B = 0 in the calculations.
As it is shown in Figure 7, the transition from (BS) to (HS) configuration as a function of temperature for the H3 system shows a small thermal hysteresis. The width of the hysteresis loop increases as the size of the crystal lattice increases. This behavior is explained by the fact that the Curie (or the order-disorder) temperature designed by T O.D. decreases when the size of the molecules decreases whereas the transition temperature designed by T eq remains constant (in our case T eq = 184.6 K for the parameters chosen for the calculations). The value of the transition temperature is deduced from the following equation: T eq = ∆ k B ln(g) , as given in Equation (18). The first-order phase transition takes place when the condition T O.D. > T eq is satisfied. This condition is fulfilled for systems H3 to H7. The case corresponding to the interaction between the molecules at the surface and their immediate environment is next investigated with L/kB = 125 K. The results are shown in Figure 8 below.  The case corresponding to the interaction between the molecules at the surface and their immediate environment is next investigated with L/k B = 125 K. The results are shown in Figure 8 below. Figure 7. The simulated thermal behavior of the total high spin (HS) fraction, Nhs, for different system sizes: H3 (green square), H4 (blue square), H5 (cyan square), H6 (magenta square) and H7 (dark yellow square). The computational parameters are: Δ/kB = 1200 K, J/kB = 14.8 K, G/kB = 136 K, L/kB = 0 K and ln(g) = 6.5.
The case corresponding to the interaction between the molecules at the surface and their immediate environment is next investigated with L/kB = 125 K. The results are shown in Figure 8 below.   Figure 9 shows that the material undergoes a spin transition with hysteresis. The thermal hysteresis decreases in width when size decreases (127, 91, 61 atoms) until its disappearance for the case H4 (37 atoms) and then reappears (for 19 atoms). This is a reentrance phase case. To highlight this behavior, the variation of the width of the hysteresis as a function of the total number of molecules has been plotted in Figure 9 below. Symmetry 2021, 13, x FOR PEER REVIEW 12 of 16 Fig. 9 shows that the material undergoes a spin transition with hysteresis. The thermal hysteresis decreases in width when size decreases (127, 91, 61 atoms) until its disappearance for the case H4 (37 atoms) and then reappears (for 19 atoms). This is a re-entrance phase case. To highlight this behavior, the variation of the width of the hysteresis as a function of the total number of molecules has been plotted in Figure 9 below. As observed in the 2D case, in the absence of surface matrix effects (L/k B = 0 K), the equilibrium temperature T eq is a constant whatever the size and shape of the studied system. With the computational parameters ∆/k B = 1300 and ln(g) = 6.01, T eq = (∆/k B )/lng ≈ 216.3 K. The 4 × 4 × 6 system displays a thermal hysteresis which is the sign of a first-order transition. As can be seen in Figure 11, the width of this hysteresis decreases for the 4 × 4 × 5 system and disappears for the cases 4 × 4 × 4 and 4 × 4 × 3. The values of the parameter t defined as the ratio between the number of molecules on the surface (N s ) and the total (N T ) number of molecules, depend on the shape of the system, are given in Table 3 and show that the thermal hysteresis decreases when the ratio t increases. This phenomenon can be explained by the fact that bulk sites favor the first order transition, while surface atoms, which increase in number with decreasing size, promote a gradual transition due to the lack of neighboring sites.   The value of parameter L/kB has been gradually increased from the value 0 K to the value 200 K. The results are reported in Figure 12 for the 4 × 4 × 6 shaping system and show that the increase in L/kB progressively shifts the equilibrium temperature of the system towards lower temperatures and is accompanied by a first-order transition. Indeed, the influence of the increases. It results in a decreased value of and then the appearance of a thermal hysteresis when the condition Teq < TO.D. is satisfied. The value of parameter L/k B has been gradually increased from the value 0 K to the value 200 K. The results are reported in Figure 12 for the 4 × 4 × 6 shaping system and show that the increase in L/k B progressively shifts the equilibrium temperature of the system towards lower temperatures and is accompanied by a first-order transition. Indeed, the influence of the T sur f eq increases. It results in a decreased value of T eq and then the appearance of a thermal hysteresis when the condition T eq < T O.D. is satisfied. The value of parameter L/kB has been gradually increased from the value 0 K to the value 200 K. The results are reported in Figure 12 for the 4 × 4 × 6 shaping system and show that the increase in L/kB progressively shifts the equilibrium temperature of the system towards lower temperatures and is accompanied by a first-order transition. Indeed, the influence of the increases. It results in a decreased value of and then the appearance of a thermal hysteresis when the condition Teq < TO.D. is satisfied.

Conclusions
Monte Carlo entropic sampling (MCES) applied to an Ising-like model of SCO compounds shows that surface atoms have no effect on the system's equilibrium temperature whatever the size or shape, when the coupling parameter of edge molecules L/k B is set to zero. However, a first order transition with thermal hysteresis occurs with an increase in the width of the hysteresis loop when the size increases. When L/k B = 0, opposite effects are observed. The thermal hysteresis increases when the size of the system decreases and when the shape gradually evolves towards an elongated rectangular 2 × 72. This feature is connected to the relative positions of the equilibrium temperature T eq and the order-disorder transition T O.D. (or Curie) temperature. When the ratio t of the number of molecules at the surface to the total number of molecules increases, the energetic contribution L acts as a negative ligand-field, lowering the equilibrium temperature. Hence, for smaller lattice sizes, the condition T eq < T O.D. is fulfilled and leads to a first-order transition. The "competition" between the effects of T eq and T O.D. is particularly striking for the hexagonal system, for which a re-entrance phenomenon is clearly depicted for 3D SCO configuration of size 4 × 4 × 6. These results show that the MCES method, is well adapted to investigate the thermal properties of spin transitions in the 2D and 3D SCO compounds although it is limited to small sizes. Future work will concern cubic centered and faced cubic centered structures and iron (II) spin-crossover binuclear compounds.