A comparative study between Monte Carlo entropic sampling method and local mean field investigations of thermal properties of spin-crossover nanoparticles based on Ising-like model

A comparative study between Monte Carlo Entropic sampling Method and Local Mean field investigations of thermal properties of spin-crossover nanoparticles based on Ising-like

In Fe-based SCO materials, the spin transition between LS and HS is accompanied by an increase of ∼ 10% of the Fe-ligand bond lengths due to the weakening of the metal-ligand bonds accompanying the population of the anti-bonding orbitals  2  .This results in the change of the unit cell volume by ∼3%-5%.The SCO molecular solids are excellent candidates as sensors of temperature, pressure [21][22][23][24] , and gases, notably of hazardous volatile organic compounds, as recently demonstrated in several reports [25,26] .They can also serve as switches of light emission in the electroluminescent devices or photoluminescent [27,28] containing SCO complexes or to enhance plasmonic resonances [29] .
In the present study, we focus on the modeling of thermal effects in 2D square SCO nanoparticles.The Isinglike model is used here and solved in the framework of local mean field approximation (LMFA) and within the Monte Carlo entropic sampling (MCES) technique in an attempt to reveal the role of interactions with the external environment (matrix).In this work, we discuss and compare the results obtained by these two methods.Experimental results on the effect of the surface or the substrate have been reported [35,40,41] .
Although the MCES technique is a good choice to solve the Hamiltonian associated with the extended Isinglike model, it is limited in terms of computing capacity for samples containing a small number of molecules.
As we would like to analyze the cases of nanoparticles with large numbers of molecules and even thin films in contact with a substrate, which cannot be studied by MCES, we plan to use an alternative formalism based on LMFA.With this contribution, the purpose of this contribution is to determine the reliability, limitations, and drawbacks of the LMFA method when applied to very small systems by comparing it with the MCES technique.
The manuscript is organized as follows.Section 2 introduces the Model and some details of the parameters used in the model.In section 3, the MCES technique is presented, while the LMFA is introduced in section 4. In section 5, we present the simulated thermal behavior of order parameters and the discussion of the obtained results.Finally, we conclude and outline the possible extensions of this work in section 6.

THE MODEL
In condensed matter physics, magnetic phenomena observed in a regular lattice are described with the wellknown Ising Model based on interactions between quantum spin variables that were first introduced phenomenologically by Pauli [42] .Later on, in 1925, Ising [43] used the concept of spins to describe the thermal behavior of the magnetization of magnetic materials using a 1D model whose statistics were treated in the frame of a transfer matrix method.Pointing quickly the exact resolution of this model (with zero applied field) performed with a remarkable mathematical tour de force of Onsager in 1944 [44] , it is also important to notice the considerable efforts of several authors in the realization of approximative methods based on variational treatments of the Ising model, which led to mean-field theory [45] , Bethe-Peierls [46] , etc., which are today widely used to analytically describe a large panel of cooperative phenomena in solid state physics and beyond.
In the field of molecular switchable materials, the Ising model [30][31][32][33][34][35] has been adapted to study SCO materials in which two different magnetic states, HS and LS, interact to display at the macroscopic scale various physical properties that give rise to a wide variety of transition behaviors.Among them, one can quote: (i) thermallyinduced hysteretic first-order phase transitions [15] , gradual or continuous [47] conversions, incomplete [48] , and multi-step transformations [49] .
To model SCO nanoparticles, a spin fictitious operator with two eigen-values +1 and -1 is associated with each molecule describing the HS and LS states, respectively.Then, the Ising-like Hamiltonian is used to model the energy operator connected to SCO.In this respect, the Hamiltonian  can be written as the sum of three energetic contributions related to the isolated molecule, the interaction between molecules, and finally, between the molecules at the boundary with their immediate external environment.Thus,  is expressed as: with (i) The single molecule contribution: where  is the number of molecules, and Δ is the energy gap at zero Kelvin between the fundamental LS state and the excited HS state.
(ii) The interaction term: where    is the coupling constant between the spins of the sites  and .
(iii) The interaction between the surface molecules and their environment: where  is the energy state of each molecule at the surface in addition to the previous ligand Δ 2 .The variable, , is the number of molecules   at the surface.
In the following, we give a detailed description of each of these terms.

Molecular contribution 𝑯 1
The two states, HS and LS, have different spin and orbital degeneracies, and the spectrum of intramolecular vibrations of the molecule strongly depends on the spin state of the central metal.
And, as it has been reported [1,2] , the density of states becomes higher in the HS state than in the LS state, which is then equivalent to considering these two states as having different effective degeneracies, denoted here by   and   .For simplicity, we consider   and   as independent of temperature.From the energetic point of view, these effective degeneracies enter the Hamiltonian as entropic terms, of which only their ratio plays a relevant role.When these specific characteristics of the two states, HS and LS, are considered, the expression of the molecular term  1 is written as: (5)   Where  =     ,  is the temperature, and   is the usual Boltzmann constant.
At this stage, we see that from the point of view of a single site Hamiltonian, as far as Δ−    ln  > 0, LS states, for which  = −1, are favored.When Δ −    ln  < 0, the HS states, for which  = +1 are favored.At the equilibrium temperature   = Δ   ln  , there is a "crossover" situation between these two configurations.The physical quantity that allows the thermal behavior of the system is the HS fraction, denoted here as  ℎ, which represents the probability of the population being in the HS state.In terms of spin values, it represents the probability of occupation of the spin value  = +1.It is straightforward to establish that  ℎ relates to the average value of the fictitious magnetization <  > as follows: Simple Boltzmann statistics performed on Hamiltonian (1) leads to the following relations, and from which the curves of Figure 1 showing thermal dependence of the HS magnetization and fraction are easily derived.

Interaction term 𝑯 2
The interaction term can be split into two contributions: (i) those arising from nearest neighbors, which are parametrized using a short-range interaction term  (assumed here to be independent of the site); and (ii) those including all the other molecules, which are parametrized with a long-range interaction term , which is written in the frame of the meanfield approach.In this respect,  2 is written as: (9)   Considering solely this interaction term H 2 , its resolution leads to the thermal evolution of <  >, as given in Figure 2, which is typical of second-order phase transition.Indeed, at 0 K, the system is either in the ordered "HS" state with all the spins with  = 1 or in the "LS" state with  = −1.When temperature increases, the average value of the spin either decreases or increases, as shown in the upper and lower curves.The temperature at which the average value is zero is called the order-disorder temperature  / or critical temperature   .
The region on the left-hand side below the upper curve (respectively above the lower curve) with respect to this T OD temperature corresponds to an ordered phase, while the right-hand side corresponds to the disordered phase.For an infinite square lattice with  = 0, the critical temperature   = 2.269 ∼ J is obtained in line with the "Onsager solution", where   is given by the equation tanh 2 However, for  ≠ 0, there is no exact analytical solution since the mean-field term acts as an "external field" linearly dependent on the magnetization.On the other hand, for a finite system,   depends on the lattice size.It is worth noting the symmetric character of the thermal dependence of the net magnetization around <  > = 0 for the ordered phases, obtained by starting from +1 or -1, as shown in Figure 2, calculated using J and G in a 6 × 6 square lattice.It is important to mention that this symmetric character remains valid when including the long-range interaction term in the Hamiltonian  2 , although the expression of   now becomes dependent on  and  [14] .

Short-and long-range interaction contributions
When the first two terms are combined to model SCO physical behavior with temperature, the Hamiltonian writes as: (10)   The symmetry in the configurations consisting of positive and negative sigma values about <  > = 0 of the interaction term  2 is broken when the latter is superposed with the molecular term  1 under the condition that   <   .
In this case, negative solutions of <  > occur when  <   = Δ   ln  while positive values are favored when  >   .Thus, the molecular term  1 drives the configurations under the constraint imposed by the term H 2 .
Therefore, if   <   , the <  > solution shifts from a <  > < 0 at low temperatures to the solution with <  > > 0 at higher temperatures, which indicates the existence of a phase transition.

Surface contribution 𝑯 3
This  3 term is restricted to the molecules at the surface.The surface Hamiltonian   can thus be written as: (11)   with  the number of molecules at the surface.Here, for simplicity, the short-range interaction term, , between nn spins is considered to be the same at the surface and in the bulk material.Upon factorization, we obtain: (12)   In order to follow how the equilibrium temperature of the surface molecules (  )  changes with the value of the external interaction , it is calculated under the condition that the longrange coupling term is set to zero ( = 0), and then it can be expressed as: However, it must be noted that the long-range interaction  applies to all molecules, and it may happen that the average at the surface is equal to zero, while at the bulk, this is not the case.

The total Hamiltonian 𝑯 𝑻
Now, as one considers the total Hamiltonian, one gets: (14)   which can be written as: (15)   In the present work, the Hamiltonian   is solved following two different approaches, namely the (i) LMFA and (ii) MCES techniques, in order to determine the thermal evolution of the average fictitious magnetization,

RESOLUTION BY MONTE CARLO ENTROPIC SAMPLING METHOD
In this section, we will distinguish three kinds of short-range interaction terms, namely,   , the interaction between molecules in the bulk,   , the interaction between the molecules at the surface, and   , the interaction between the molecules at the surface and in the bulk.The corresponding Hamiltonian is written as: (16)  with where the respective total and surface magnetizations,   and   , respectively, are given by: can be re-written as: Where   ,   , and   are the short-range correlation functions related to the bulk, surface, and bulk-surface neighboring sites, respectively.Their expressions are given by the following relations: The thermal average value of the fictitious magnetization, <  >, is calculated by the following expression: Where  is the partition function   is the energy of the macrostate , and     ,    ,    ,    ,    is the density of this macrostate.Here,  is given by: and In Equations ( 24) and ( 25),   is the number of different configurations with the same five values   ,   ,   ,   , and   .The density of states  (  ,   ,   ,   ,   ) is calculated by entropic sampling [48][49][50] , and Equation ( 24) is solved by numerical techniques such as bisection.From the thermal average values <  >, the HS fraction  ℎ is calculated using the relation (6).
Let us recall the basic principles of MCES [50][51][52] .It consists in introducing the detailed balance equation of the Monte Carlo (MC) procedure, as shown in: This suited biased distribution  is designed to favor configurations belonging to weakly degenerated macrostates while dampening those within the highly degenerated macrostates.It can be expressed as: (28)   In this case, the balance equation can be written as: The resulting restricted density of states is calculated after the correction for the bias was applied: The flat character of the histogram  (  ,   ,   ,   ,   ) has a convenient convergence criterion to obtain a very good estimation of the density of the states.In our case, seven iterations of the 10 6 Carlo Steps were necessary.After the table of the values of (  ,   ,   ,   ,   ) and d (  ,   ,   ,   ,   ) is obtained, the quasiexact partition function can be calculated using the expression (25).

LOCAL MEAN-FIELD APPROACH
In the present work, three types of sites related to three lattice regions are considered: the atoms located in the bulk (N b ), which are surrounded by four first-neighbors and those atoms located on the surface and, more precisely, on the edge (N e ) and on the corner (N c ), which interact with three and two first-neighbors, respectively.The molecules located at the surface (edge and corner) interact with an external environment (matrix effect), which gives them specific properties.
The case of a system comprising 36 molecules is represented in Figure 3.
Note that in the framework of a classical mean-field approximation [53] and the case of a 2D square-shaped lattice, the coordination number is q = 4 [Figure 4].
The consideration of three types of sites leads us to distinguish three local average order parameters <   > and three coordination numbers q  .In the LMFA, which amounts to taking into account a local Hamiltonian for each type of site  with  = , , , corresponding to bulk, edge, and corner, respectively.
For each region, , we apply the LMFA.In the case of regions connecting each other only via the long-range interaction, each molecule   is surrounded by   <   >, and the local Hamiltonian, given by Equation ( 14), can be expressed as Equation (32) below: ℎ    (32)   with where   is the number of interactions between a molecule and its first-neighbors,   is the number of interactions  between a molecule and the external environment, and   is the interaction between the nearest neighbor (nn) molecules.For molecules located in the edge and in the corner, the average value of the number of interactions with the external environment is set equal to 1 .Table 1 summarizes the characteristics of each site.
Here, for simplicity of calculations, we have taken z = 1 (the same coordination number) everywhere in the surface (edge and corner), meaning that we are considering the same connectivity for all molecules at the surface.
The short-range interaction term   depends on the location of the molecules.Two interaction constants are considered in our calculations: for the molecules located in the bulk; for the molecules located at the surface.
These two interaction terms,       and       , are related to the interaction terms used in the MCES method by weighting the various interactions of each molecule in each region.Let us consider the case of the bulk region.We denote by   ,   , the respective bulk-bulk and the bulk-surface bonds, leading to a total number of bonds including at least one bulk site (  +   ).
Let     and     be the bulk-bulk and bulk-surface short-range interactions used in MCES method.The total interaction energy term in the bulk region is given by (  ×    ) + (  ×    ).Therefore, the bulk value        used in Equation ( 32) is obtained by averaging over     ,     weighted by the numbers of their corresponding bonds, as follows: In a similar way, for the surface sites, Equation ( 35) is written as: where   and     are the number surface-surface bonds and the surface-surface shortrange interaction term used in the MCES method, respectively.
The average spin state can be written as follows: So, for each region  (bulk),  (edge), and  (corner), the average spin state is given by: (39)   which is equivalent to: The weighted average of the three local order parameters associated with the three regions, bulk, edge, and corner, leads to the average fictitious spin of the system expressed as: and the high-spin fraction, Nhs, which is the probability that the HS state is occupied and which is written as Equation ( 6): The purpose is to solve a system of three equations  1 ,  2 , and  3 , whose three order parameters <   >, <   >, and <   >, connected by Equation (43), are simultaneous solutions.The Newton method, which is appropriate for rapid convergences, has been used.Calculations are based on three points  1 ,  2 , and  3 close to the solution and which represent <   >, <   >, and <   >, respectively.
The calculation of the three points "close to the solution" is performed in two steps: Let X be the space in which values of ( 1 ,  2 ,  3 ) are scanned.Let  be the space from which ( 1 ,  1 ,  1 ) are calculated as a function of ( 1 ,  2 ,  3 ).Then, the first step consists of mapping space X on space .Those values of ( 1 ,  2 ,  3 ) that map onto ( 1 ,  2 ,  3 ) close to zero within 10 −6 , are then selected to perform the second step: application of the Newton-Raphson technique.
Solving the system amounts to calculating ℎ 1 , ℎ 2 , and ℎ 3 such that: which allows to write in the vicinity of points  1 ,  2 , and  3 : The following system of three equations with three unknowns is thus obtained as: And where ℎ 1 , ℎ 2 , and ℎ 3 are calculated as follows: In a general way, the   functions with i = 1 (bulk), 2 (corner), and 3 (edge) and their derivatives can be expressed as follows:   With the ℎ  obtained in this way, the solutions (  )  are deduced by the relation:

RESULTS AND DISCUSSION
We are interested in the effect of the coupling between surface and bulk molecules, and for that, we perform numerical simulations by MCES and by LMFA.To proceed, we choose a set of realistic SCO thermodynamic parameters, and for this purpose, we chose the data derived from experimental data of one of these SCO compounds: [Fe(btr) 2 (NCS) 2 ], btr = 4,4-bis-1,2,4-triazole [54] , whose molecular structure is shown in Figure 5.
The energy gap Δ and ln() are derived from the enthalpy change (Δ ≈ 11 kJ/mol) and the entropy change (Δ ≈ 50 J/mol/K), respectively: Δ = Δ/ ≈ 1, 300 K), ln () = Δ/ ≈ 6.01, and  is the perfect gas constant.The equilibrium temperature of the system   is deduced by   = Δ/  () , which leads to ≈ 216.3 K.In the first step, in the case of a 6 × 6 configuration, the values of the short-range interaction parameters     and     are set to 60 K, and the intensity of the short-range interaction parameter     between bulk and surface is gradually increased from 0 to 60 K.In the second step, the size effect is studied with each method.It should be added that the dimensions of the square lattice used in these simulations do not correspond, for the moment, to synthesized nanoparticles.34) and (35), are gathered in Table 2.The curves highlight a two-step behavior in the form of two well-separated hysteresis loops.The hysteresis width is defined as Δ =    −  , where    is the ascending thermal transition temperature and   is the descending thermal transition temperature.The average temperature between    and   is the equilibrium temperature   , which corresponds to a HS fraction equal to 1/2(ℎ = 1/2).
Figure 6 illustrates a qualitative and quantitative agreement between the calculations obtained by MCES and LMFA.The hysteresis loops, observable for each method, at lower temperatures, ∼ 144 K, are related to the 4.8 The simulation parameters are: behavior of surface molecules (edge and corner).These surface molecules commute from LS to HS state before the bulk molecules, and they drive the thermal transition.The hysteresis loops, positioned for each method at 185 K, are related to the bulk molecules.The two methods of simulation lead to identical values of the transition temperature in the surface (      =      ∼ 144 K) and in the bulk (      =      ∼ 185 K).The bulk and the surface appear as two systems that evolve almost independently of each other.Between 150 and 182 K, the curves reveal, both by MCES and by LMFA, a long intermediate plateau (mixture of LS and HS configurations) around  ℎ ≈ 0.55.
It can be noticed that the hysteresis widths Δ are slightly larger by LMFA than by MCES.The values of the thermal hysteresis widths are summarized in Table 3.
The value of the short-range interaction parameter     /  is gradually increased up to 60 K. Figure 7, Figure 8 and Figure 9, show the results obtained respectively for     /  = 20,     /  = 40 and     /  = 60 K in the case of a 6 × 6 lattice.The values of the interaction parameters  used by MCES and by LMFA, obtained from Equations ( 34) and (35), are gathered in Table 4.
The phase diagrams of the 6 × 6 system, obtained by LMFA and by MCES, are plotted in Figure 10 and Figure 11, in the space coordinates, temperature vs.     /  , for the bulk and for the surface molecules, respectively.By MCES, as can be seen in Figure 10 and Figure 11, when the interaction parameter     /  increases, the equilibrium temperature of the surface molecules,       , is shifted to higher temperatures, whereas the equilibrium temperature of the bulk molecules,       , is shifted to lower temperatures.The hysteresis loops,   The computational parameters are: Δ/  = 1, 300 K, /  = 172 K, /  = 290 K, and ln() = 6.01.
drastically accentuates the hysteresis effect, especially for surface molecules.The system seems constrained, and the tendency for the two hysteresis to overlap is not observed.The computational parameters are: Finally, in a larger particle of size 7 × 7, the hysteresis associated with the bulk are wider than the hysteresis associated with the surface molecules by LMFA and by MCES.The proportion of molecules in the bulk (∼ 51%) being greater than that of molecules in the surface (∼ 49%), the bulk molecules drive the transition.In the 6 × 6 compound, the hysteresis related to surface molecules was greatest in both MCES and LMFA.These results are summarized in Table 8.

The case 𝑱
The phase diagrams of the 7 × 7 system obtained by LMFA and by MCES are plotted in Figure 13 and Figure 14, in the space coordinates, temperature vs.     /  , respectively, for the bulk and for the surface molecules.
When the value of the interaction parameter     /  is increased up to 60 K, the three following items are especially significant.
First, for each method, the curves presented in Figure 12 and Figure 13 reveal that the equilibrium temperatures of the surface and of the bulk in the 7 × 7 lattice are higher than those obtained in the 6 × 6 lattice [Figure 10 and Figure 11].
In addition, as already observed in the case of a 6 × 6 lattice and by the LMFA, when the short-range interaction     /  is increased, the equilibrium temperature of the surface molecules decreases.We further note that the increase of the lattice size leads to a decrease of the width of the hysteresis in the surface.By MCES, the hysteresis transition observed for surface molecules disappears (Δ     = 0 K) in favor of an abrupt transition [Table 9].The computational parameters are:  of the bulk {  = Δ/[  ln() ≈ 216 K}.The results obtained in the case of the 7 × 7 lattice suggest that if we could do simulations on larger lattices with MCES, a single hysteresis loop with an equilibrium temperature close to 216 K would also be obtained.

CONCLUSION
The present contribution, a simulation work, with the purpose of giving some hints to colleagues who synthesize spin-transition nanoparticles, is focused on meticulous investigations of the thermal properties of 2D-SCO nanoparticles using an adapted two-state Ising-like Hamiltonian.At this end, three types of interaction terms are considered: long and short-range interactions originating from electronic and elastic molecular couplings inside the lattice, to which we added a specific energetic contribution accounting for the interactions between surface molecules with their surrounding environment (inside and outside the lattice).The model is solved in two distinct ways, which involve different levels of sophistication: (i) MCES techniques, which are exact methods and (ii) LMFA based on mean-field theory but taking into account the lattice shape.In the former case, three types of two-body short-range interaction terms are involved between bulk-bulk (  ), surface-surface (  ), and bulk-surface (  ) sites.The thermal dependence of the "order parameter", <  >, that is, the average fictitious magnetization, is calculated by using the Boltzmann distribution obtained from the calculated density of the macrostates.On the other hand, in the LMFA, two interaction terms corresponding to three types of sites, corner, edge, and bulk, are included in the calculations, taking into account the dependence of the coordination number on the site location.The results of the two methods have been compared for different square lattice sizes, and a good quantitative agreement is obtained between them in the case   = 0, where both equilibrium temperatures and hysteresis widths were found in very good agreement.When the bulk-surface interaction   is increased, the hysteresis observed by MCES in the bulk and in the surface tend to be close to each other and finally overlap for strong Jbs values.In contrast, by LMFA, this phenomenon is not observed, and furthermore, the hysteresis effects are accentuated, particularly for surface molecules.Moreover, with the MCES method, the increase of   results in a shift of the surface (resp.bulk) transition temperatures to higher (resp.lower) values, while in LMFA, an opposite effect is obtained for surface molecules.
Our comparison of the two techniques (MCES and LMFA) shows that when dealing with small samples, LMFA yields results that are in perfect agreement with MCES in the case where the interaction between the surface and the bulk is zero (  = 0).However, as the interaction   increases, causing important interferences between the two regions of the lattice, the results provided by LMFA increasingly differ from those obtained with the MCES method, although the general trends remain the same.While considering a very weak interaction between the surface and the bulk in the same compound may look unnatural, it may be possible in nanostructured coreshell nanoparticles, where the core and the shell are of different compositions.In addition, it is important to keep in mind that the systems considered are particularly small, and as the size of the system increases, the predominance of the surface, and then   interaction, vanishes, which, in turn, means that LMFA will most likely still provide correct results given the large number of molecules.
Finally, it is worth mentioning that the surface properties have been limited here to the outer layer of the lattice, while the concept of the surface may integrate some inner layers close to the interface material/air, which also have specific properties different from those of the bulk.Extensions of the present work to include the abovecited aspects are in progress, as well as other packing clusters such as rhombohedral or faced centered cubic (FCC).The results of this simulation work can also be extended to explain the behavior of polycrystalline SCO films on the Al 2 O 3 substrate [40] and the influence of the substrate on the functionality of SCO molecular materials [41] .

Figure 2 .
Figure 2. Thermal evolution of the average magnetization <  > obtained in a 6 × 6 square lattice with the following parameters: / = 60 K and / = 172 K.

Figure 3 .
Figure 3. Schematic representation of a 2D 6 × 6 square-shaped lattice: filled blue circles represent bulk sites (Nb); filled green and red circles represent edge (Ne) and corner (Nc) sites, respectively, interacting with their immediate environment (matrix).

Figure 4 .
Figure 4. Mean-field approximation in a 2D square-shaped lattice.The  spin is surrounded by four nearest neighbors with the average magnetization <  >.

Table 2 .
Correspondence between the values of short-range interactions / used by LMFA and by MCES method in the case of a 6 × 6 square lattice when

Figure 6
Figure6displays the results obtained in the case of a 6 × 6 square-shaped lattice, and the values of the interaction parameters  used by MCES and by LMFA, calculated from Equations (34) and(35), are gathered in Table2.The curves highlight a two-step behavior in the form of two well-separated hysteresis loops.The hysteresis width is defined as Δ =    −  , where    is the ascending thermal transition temperature and   is the descending thermal transition temperature.The average temperature between    and   is the equilibrium temperature   , which corresponds to a HS fraction equal to 1/2(ℎ = 1/2).

Figure 10 .
Figure 10.Phase diagram T vs.    / for a 6 × 6 system.The red, blue and black squares correspond, respectively, to the upper and lower transitions () and () temperatures of the thermal HS fraction and to the equilibrium temperature of the bulk by LMFA.The red, blue and black stars correspond, respectively, to the upper and lower transitions () and () temperatures of the thermal HS fraction and to the equilibrium temperature of the bulk by MCES.The computational parameters are:    / = 60 K,    / = 60 K,

Figure 13 .
Figure 13.Phase diagram T vs. (   /) for a 7 × 7 system.The red, blue and black squares correspond, respectively, to the upper and lower transitions () and () temperatures of the thermal HS fraction and to the equilibrium temperature of the bulk by LMFA.The red, blue and black stars correspond, respectively, to the upper and lower transitions () and () temperatures of the thermal HS fraction and to the equilibrium temperature of the bulk by MCES.The computational parameters are:    / = 60 K,    / = 60 K,

Figure 14 .
Figure 14.Phase diagram T versus (   /) for a 7 × 7 system.The red, blue and black squares correspond, respectively, to the upper and lower transitions () and () temperatures of the thermal HS fraction and to the equilibrium temperature of the surface by LMFA.The red, blue and black stars correspond, respectively, to the upper and lower transitions () and () temperatures of the thermal HS fraction and to the equilibrium temperature of the surface by MCES.The computational parameters are:    / = 60 K,    / = 60 K, /  = 60 K,    /  = 60 K,    /  = 60 K,      /  = 60 K, Δ/  = 1, 300 K, /  = 172 K, /  = 290 K, and ln() = 6.01.Simulations using the LMFA were carried for very large lattices, 20 × 20, 50 × 50 and 200 × 200. Figure 15 illustrates the results obtained for two sets of parameters corresponding to     /  = 0 K and     /  = 60 K and shows that, upon increasing the lattice size, the two-step transition vanishes.The curve closely resembles a single hysteresis loop for the 200 × 200 lattice, and the equilibrium temperature tends towards that