Disentangling Surface Energy and Surface/Interface Stress Effects in Spin Crossover Nanomaterials

Understanding phase transformation processes at the nanoscale is an important step for the integration of phase change materials into functional nanodevices. Here, a numerical method to assess surface energy and surface stress contributions to finite size effects in spin crossover nanomaterials is proposed. Starting from their formal definitions, molecular dynamic simulations combined with continuum mechanics and thermodynamic considerations on model thin films are performed. The surface energy values extracted from the simulations are 71 and 79 mJ m−2 for the high spin (HS) and low spin (LS) states, respectively. In the limit of a weak lattice parameter misfit, the calculated isotropic stresses of a HS/LS interface are 94 and 45 mJ m−2 for the LS and HS states, respectively. From these results, a ≈10 K downshift of the spin transition temperature is predicted with the size reduction. Whereas surface energy favors systematically the HS state, it is not the case for the interface stress. Indeed, it is demonstrated that the anisotropy of mechanical stresses can lead to the stabilization of the LS state. These results confirm the possible control of the phase stability in these smart nanomaterials through interface engineering.


Introduction
[3] To this aim, phase transitions in materials consisting of switchable elementary bricks like spin crossover (SCO) compounds are promising candidates. [4]Indeed, SCO complexes are constituted of a transition metal cation of electronic configuration 3d n , with n = 4, …, 7, surrounded by ligands in a pseudo-octahedral symmetry. [5]hen the ligand field strength is of the same order of magnitude DOI: 10.1002/apxr.202200055[8] In the solid state, the thermodynamic competition between the two spin states and the elastic interactions between SCO entities lead to non-linear effects, which are likely to produce solid-solid first order phase transition and abrupt spin state changes with hysteretic behaviors (bistability phenomenon) can be observed. [9,10]oncomitantly with the electronic reorganizations upon a LS to HS transition, a plethora of physical properties changes, including magnetic, optical, structural, vibrational, and elastic properties.Owing to these assets, SCO materials are expected to be candidates for future applications in data processing, photonic, spintronic, and actuator devices. [4,11,12]efore any development and integration into devices, the size of SCO materials has to be reduced.[15] An important issue is the evolution of the switching properties when the size of the SCO system is reduced down to the nanoscale.As a matter of fact, nanosized materials behave differently from their bulk counterpart.In particular, phase stability can be strongly affected by size reduction effects.In first order phase transitions, it results in a change of the transition temperature T eq and the occurrence of incomplete transitions.In SCO materials, it results most of the time in a decrease of the transition temperature, [16][17][18] a common feature with other first order transitions such as the melting process. [19]However, a slight increase of T eq has been observed in some thin films and sometimes more complicated behaviors have been reported, [20] such as non-monotonous behaviors of T eq with the size reduction.From a fundamental viewpoint, one of the key features of nanosized materials is their high surface to volume ratio.[24] They reveal that surface phenomena appears to favor almost systematically the HS phase, leading to an ineluctable decrease of the transition temperature. [25,26]However, another relevant contribution has to be considered: the external environment.38][39][40][41] However, no quantitative investigation of surface/interface effects has been reported yet and it is difficult to predict in which extent it will be possible to modulate surface effects through the control of interface properties.In principle, two quantities are needed to fully describe mechanical properties of surfaces/interfaces: the surface energy and the surface/interface stress, which are experimentally not easily accessible quantities. [42,43]In this work, we propose to combine molecular dynamics and continuum mechanics to obtain some primary estimations of surface energy in both spin states and to extract interface stress properties from a simplified model structure of thin films, which was recently introduced. [44]In Section 2, we introduce the methodology employed to extract thermodynamic and elastic properties of the bulk materials.Then, surface energy and surface/interface stress are defined and the numerical methods to assess them as well as their integration into the Gibbs free energy of thin films are introduced.Section 3 will be devoted to the presentation of the different numerical results, which are compared with the literature.In Section 4, we extend these numerical methods to anisotropic structures.Finally, the different results are summarized in Section 5 and some perspectives to this work are given.

Model Structure and Force Field
We simulate a cuboid SCO crystal of dimensions L x , L y , and L z , constituted of "heavy" metal centers (M Me = 57 u, u being the unified atomic mass unit) occupying the cubic lattice nodes.Each metal is surrounded by six nitrogen atoms (M L = 14 u), the whole forming an octahedral pattern, which models the structure of a SCO compound in a very simple manner (see Figure 1a).The molecular mechanics treatment of this system involves connecting this collection of atoms by means of springs which define the force field (FF).Here, it is important to note that even if electrons are not considered explicitly, the electronic reorganization related to the spin state change is taken into account in the different structural parameters and force constants included in the FF.The total potential energy is expressed as a simple sum of terms describing respectively bond stretching, angle bending, and nonbonding interactions.This simplified structure does not need to introduce the contribution of proper and improper dihedral torsions.More precisely, the intramolecular ligand-ligand (LL-intra) and metal-ligand (Me-L) atoms are bonded by different harmonic pairwise potentials (bond stretching), called respectively V LL − intra and V Me − L , which take the following form: where k is the force constant, r ij (t) is the instantaneous distance between two interacting atoms, and r 0 is the equilibrium distance.The octahedrons are connected by non-bonded interactions, which take the form of a Lennard-Jones pairwise potential (LL-inter) where ϵ inter is the intermolecular cohesion energy and  corresponds to the distance at which V LL − inter is zero.This interaction is limited to first neighbor octahedra in order to limit the computational time.We shall mention that the anharmonic character of the intermolecular potential is essential since it allows large displacements of atoms, far from their respective equilibrium positions.This has been shown to enhance surface effects. [46]More importantly, unlike a harmonic potential, anharmonic interactions permit the bond dissociation and surface energies can thus be calculated.In order to avoid thermal contraction and limit the rotation of molecules, harmonic angular potentials (angle bending energy) between three consecutive ligands (see Figure 1b) have been added where θ0 is the equilibrium angle, C  is the angular bending force constant, and θijk (t) is the instantaneous angle between three interacting ligands.It is legitimate to ask the question as to the relevance of this simplified model.It is important to emphasize that experimental or theoretical information concerning the surface physical properties of SCO nanomaterials are scarce and the effect of spin state change on these properties remains elusive.Before entering into complexity (e.g., low symmetry, large unit cells), it is essential to start with simple, general models, which are able to catch the main features of the structural properties of SCO compounds (e.g., octahedral symmetry).The different numerical parameters of the FF employed for the molecular dynamics (MD) simulations are summarized in Table 1.
These parameters are order of magnitudes deduced from DFT calculations, [47] X-ray diffraction, and Raman spectroscopy, [48,49] except the angular potential, which cannot be easily related to experimentally measurable physical quantities, vide supra.(For more details about the assessment of the FF, see the supplementary materials of ref. [50]).It is interesting to compare the strength of intermolecular interactions of this present FF with those employed to simulate the epitaxial strain at a SCO complexmetallic interface. [51,52]Indeed, an effective spring constant k can be roughly deduced in the harmonic approximation.We obtain In this case, the LS system is made up of LS multi-layers and a single-layer HS residual fraction on the surface.The number of LS layers (N LS ) varies with size, while on the contrary, the number of HS molecules on the surface remains constant:N HS = 1.At the HS/LS interface: the misfit of lattice parameters between the two subsystems generates biaxial residual stresses along the x-and y-directions.At the LS surface, a surface energy exists due to the broken LS bonds, which is calculated by the model II.Then, the stability of the LS phase in model III is not only affected by the surface energy as in model II, but rather by the sum of surface energy and interface stresses at the LS surface and HS/LS interface, respectively.A similar description for the HS system made up of HS multi-layers and a single-layer LS residual fraction on the surface can be achieved.
Table 1.Summary of the parameters employed in the molecular dynamics simulations.
75 eV Å −2 and k HS = 1.21 eV Å −2 for the LS and HS states, respectively.The averaged value 〈k〉 ≈ 1.48 eV Å −2 is very close to the calculated spring constant extracted from a fit of DFT calculations in ref. [51].This indicates that the parametrization of this FF allows a quantitative extraction of mechanical properties of SCO materials, a crucial point for the modeling of surface/interface properties.Note also that electrostatics is not modeled in our FF, because its impact on the vibrational properties is considered negligible. [47]n the following, three model systems are built in order to separate bulk, solid-vacuum surface and solid-solid interface contributions in physical quantities characterizing SCO nano-objects, namely, in the internal energy.First, a SCO system where closed periodic boundary conditions are applied, is considered (see Figure 1b, model I).This situation constitutes a reference and is often called slab-adapted method. [50]Model I allows to extract thermodynamic and mechanical properties of a system where surface and interface effects are absent.Then, a solid-vacuum surface is created by increasing the length h = L z + L vacuum of the simulation box in the z-direction (L vacuum ≈ 20 nm).Numerical results become independent of L vacuum if this latter is sufficiently large to avoid self-interactions due to periodic boundary conditions. [50]his corresponds to the model II (see Figure 1c), from which surface energy can be calculated for the two spin states.Finally, a multi-layer system surrounded by vacuum and constituted of LS and HS phases, separated by a semi-coherent solid-solid interface is built.An elastic misfit m is generated at the interface due to the difference in terms of lattice parameters and mechanical properties between the two phases (see Figure 1d, model III).An interesting limit is to establish a contact between a solid constituted in the LS phase, covered by a HS monolayer, in order to mimic the presence of a residual HS fraction located at the surface.The simulation of a semi-coherent solid-solid interface and the extraction of the corresponding elastic energy stored due to the misfit can thus be performed.Models II and III are the application of the so-called slab method. [50]

Thermodynamic Aspects and Phase Stability of Finite-Size Systems
Let us consider a slab simulating a thin film in contact with a thermal bath and a barostat where the pressure P is set to zero, corresponding to the isothermal isoabaric ensemble (T, P = 0, N).In the case of the model I and for a sufficiently large system (L z → +∞), the variation of Gibbs free energy per molecule ΔG(T) between the HS and the LS phases can be assimilated to that of the bulk material ΔG b = ΔH − TΔS, where ΔH and ΔS stand for the enthalpy and the entropy variations per molecule, respectively.In the case of models II and III, the contributions of surface/interface have to be taken into account in the thermodynamic properties of the thin film.A first approach uses the concept of excess quantity.Indeed, bulk values of extensive quantities such as ΔG, are modified by an excess amount owing to the presence of surfaces/interfaces: where ) corresponds to surface (respectively, in-terface elastic) Gibbs free energy per molecule, which are size dependent quantities.We note here that the interfacial energy is considered negligible in comparison with the surface energy for a semi-coherent interface and this term will be thus omitted.In a recent work, [50] we have shown that the surface entropy difference (per molecule) ΔS s was negligible in the case of this modeled structure.The contribution of solid-solid interface to the entropy is more difficult to estimate, but we can reasonably assume that ΔS int (L z ) ≈ 0 in this work since a semi-coherent interface between two phases having similar mechanical and structural properties are considered, so that the production of configurational entropy is limited.This does not mean that there is no contribution from surface and interface in the entropy difference between the two phases.Indeed, through nuclear inelastic scattering (NIS) experiments, a change in the entropy variation between the two spin states has been observed with the size reduction and a dependence of thermodynamical quantities with the nature of the matrix in which nanoparticles are embedded has been highlighted. [27]However, no experimental information concerning the size evolution of surface entropy variation in thin films are available.In nanoparticles, the entropy increase is mainly attributed to the existence of vibrations of the nanoparticle as a whole [44] (breathing, torsional modes), which either does not exist or are strongly damped in films.Yet, we must be careful with this assertion in absence of clear experimental proofs and the possible role of surface/interface entropy contributions remains an open question.Therefore, this present study focuses on the energetic and mechanical properties of the surface/interface and the entropy variation ΔS is reduced to that of the bulk material in a crude approximation.By definition, the spin equilibrium occurs when the HS and LS phases coexist in the same proportions, that is, when ΔG = 0.The equilibrium temperature T eq is where T b eq , is the equilibrium temperature in the case of the bulk material.T s eq and T int eq represent the change in the transition temperature due to the presence of surfaces and interfaces, respectively.
In the next sections, we will show the different methods that allow the assessment of each quantity in the total Gibbs free energy and the extraction of the equilibrium temperature T eq .

Molecular Dynamics and Lattice Dynamics Investigations
The calculation of bulk thermodynamic properties requires the estimation of both vibrational and electronic properties at a finite temperature, as ΔH and ΔS are the sum of vibrational and electronic contributions.In order to access to vibrational properties, the total vibrational density of states (vDOS) has to be calculated.To this aim, we perform molecular dynamics (MD) in the isothermal isobaric ensemble (T, P, N) using the Nosé-Hoover thermostat-barostat, [53,54] which allows to control the pressure and the temperature (set to T = 200 K) of the external environment.The lattice dynamics properties of the modeled structure have been investigated in both spin states.Initially, the velocities of each atom constituting the slab have been randomly generated according to the Maxwell-Boltzmann distribution, so that the initial temperature of the system is identical to that of the thermal bath.Then, Newton's equations of motion have been integrated using the velocity Verlet algorithm, implemented in the lammps software. [55]This numerical method allows to access to the time evolution of the position and the velocity of each atom.The simulations have been run with an elementary time step t = 1 fs.The first 20 000 steps ( relax = 20 ps) correspond to a transient regime devoted to the relaxation/thermalization. The velocity data have been collected during the 1000 subsequent steps for the determination of the velocity auto-correlation function (), defined as where v i (t) is the velocity at time t of the atom i and ∑ N i is the summation over the N considered atoms.The vDOS corresponds to the real part of the Fourier's transform of (): From the (vDOS), the thermodynamic quantities such as the vibrational entropy per molecule S vib can be extracted.Considering Bose-Einstein statistics applied to a perfect phonon gas (anharmonic effects are not taken into account in the expression of S vib ), S vib can be written as where  = 1/k B T, k B is the Boltzmann constant and ℏ is the reduced Planck constant.The electronic part of the entropy is primarily attributed to the total spin momentum S, which can be written as In the case of the ion Fe(II), the spin momentum is S HS = 2 (respectively S LS = 0) for the paramagnetic HS (respectively diamagnetic LS) state.The enthalpy difference per molecule is calculated as follows where the bulk transition temperature T b eq is fixed at 200 K.The Debye sound velocity v D is extracted using the Debye model, valid in the low-energy region of the phonon spectrum [56]  (11)   where  = ∑ i N i M i /V u.c. is the material density, N i M i is the product of the number of atoms of type i by their respective mass, the sum running over the different type of atoms, and V u.c. is the volume of the unit cell.m = 6M Me M L ∕(M Me + 6M L ) is the effective mass of ligands and metal cations and is the normalized vDOS.The bulk modulus B (or equivalently the Young's modulus Y) can be estimated as follows [56] where 3 are constants, which depend on the Poisson ratio  = 1/3. [57]

Computation of Surface/Interface Physical Quantities
The surface energy  0 is defined as the work required for bond dissociation (here LL-inter), divided by the created surface area.In other words, it is related to the excess of potential energy per unit area in a system having free surfaces in comparison with its bulk counterpart: [58] where E II (respectively E I ) is the potential energy calculated using model II (respectively model I) and A = L x × L y is the area of the surface.It is important to note that, during the cleaving process, two equivalent and identical surfaces are created.This explains the presence of a factor 1/2 in 13.The surface energy is dependent of the number of broken bonds and thus on the number of atoms at the surface.It is convenient to introduce an effective surface energy: with  s = N s /N ≈ 2/L z is the form factor of the thin film, N s is the number of surface/interface molecules at the two faces, and N = N b + N s is the total number of molecules.Hence, the Gibbs free energy difference per molecule related to the creation of surfaces in a thin film is where a s = l 2 is the elementary surface unit cell area (see the red filled square in Figure 1c).
The other important quantity is the mechanical stress and the related energy due to the presence of surfaces/interfaces.Let us consider a homogeneous thin film characterized by homogeneous bulk stress and strain tensors  b ij and  b ij , without the presence of surface/interface.We can thus define the excess of elastic energy whose variation is called interfacial elastic energy change W int : where  ij (z) and ϵ ij (z) stands for the real stress and strain tensors in the presence of the interface.V = A × h is the volume of the simulation box.These tensors can be decomposed into a perpendicular and a parallel contribution and considering that the thin film remains in a mechanical equilibrium in the direction normal to the surface, that is ∂ iz /∂z = 0 for i = x, y, z, the expression of W int can be simplified as follows where Here, s ij is the interfacial excess quantity (per unit area) of the parallel components of the bulk stress tensor and e ij is the interfacial excess quantity of the perpendicular components of the bulk strain tensor.In this work, the SCO slab is surrounded by vacuum.There is neither stress nor strain tensor in this phase.Therefore, the mechanical equilibrium condition leads to  ⟂ ij =  b iz = 0 (i = x, y) and the surface stress tensor is reduced to the following expression The expression of the interfacial elastic energy change becomes It is important to mention that the tensor s ij corresponds to the stored in-plane stress of the surface and should not be confused with the Laplace over-pressure existing in a curved crystal.In this latter case, the normal component of the stress tensor  zz is no longer zero and follows a Young-Laplace like equation called the Gibbs-Thompson equation in the case of solid interfaces.As long as we deal with a free planar surface, only a tangential surface stress is expected.This condition is fulfilled for small elastic misfits m.Indeed, for large m, the thin film bends and a deflection can be observed.The surface stress must then be compensated by bulk shearing stresses, the mechanical equilibrium is not necessary reached and Hooke's law is no more valid.This situation is out of the scope of the present work.
The interfacial elastic Gibbs free energy difference per molecule between the HS and LS phase is where a I = (a HS s + a LS s )∕2 is a mean elementary interface area, assuming a homogeneous distribution of strain.In practice, the interfacial stress is numerically calculated using the Gradient conjugated minimization method [59] which allows the system to reach the mechanical equilibrium [60] sij (m) = h σij (m) (23)   where σij (m) is the spatially averaged residual interface stress, which depends on the misfit m.Then, an effective interface stress can be introduced as follows Finally, the interfacial elastic contribution to the Gibbs free energy difference per molecule can be deduced 3. Results

Calculation of Bulk Thermoelastic Properties Using Model I
The investigation of lattice dynamical properties has been carried out in both spin states through the calculation of the total density of vibrational states (vDOS) on a slab having a thickness of ≈13 nm (2500 octahedrons) in the framework of the model I.This thickness has been shown to be sufficiently large to model the vibrational properties of a bulk material, [50] avoiding the discretization of the vDOS spectrum related to confinement effects.
The slab adapted method assures the absence of surface effects.
The vDOS reveals a red-shift of the HS vibrational modes relative to the LS modes in accordance with the parameters employed in the FF for the two spin states.Above E = 20 meV, several optical modes can be observed.The Me-L stretching mode is located at 26 meV (respectively 21 meV) in the LS (respectively HS) state.The two following peaks, situated at E = 40 meV and E = 60 meV in the LS state (respectively E = 36 meV and E = 48 meV in the HS state), denote the existence of two ligand-ligand bending vibration modes within octahedrons.
As expected, no drastic change of shape could be mentioned when we compare the LS and HS density of vibrational states due to the force field properties.In real systems, upon the spin state switching, the molecule does not change only its molecular volume, but also its shape, leading to more important differences of the vDOS between the two spin states.In the present work, only the volumetric expansion of the octahedron upon a LS to HS transition is taken into account.As a consequence, the variation of the thermodynamic quantities between the two spin states is certainly underestimated.
As mentioned above, the calculation of the vDOS allows us to extract thermodynamic quantities such as the vibrational entropy (Equation ( 8)), which is (S LS = 108.5 k B , S HS = 117.2k B ).We have also calculated the electronic part of the entropy (Equation (9)), which must be summed with the vibrational contribution to be integrated into the thermodynamic model.The results are summarized in Table 2.[63][64][65] From the low frequency part of the vDOS spectrum, the Debye sound velocity, characterizing the mechanical properties of the thin film, has been extracted through the interpolation of g(E)∕E 2 between 0 − 1 meV (see insert of Figure 2). [56]Not only the relative variation of v D between the HS and LS states, but also the absolute values for both spin states are of the good order of magnitude when compared with the data obtained from NIS measurements. [44,66,67] corresponds to the range where the Debye approximation is valid. [56]erefore, the stiffness of the simulated thin film represented by the bulk and the Young moduli ( 12) is realistic for SCO molecular complexes. [57]Such simulation of the mechanical properties is essential before any quantitative theoretical investigations of surface/interface effects on phase stability and switching properties of SCO nanomaterials.We can note that very recently, an empirical molecular mechanics force field for a specific system, the [Fe(pyrazine)] [Ni(CN) 4 ] spin complex has been built. [68]The parametrization of the force field has been performed from Raman spectroscopic and nuclear inelastic scattering (NIS) data and turns out to be rather precise.In particular, the density of vibrational states is well described in both spin states and, as a result, the mechanical (elastic moduli), thermodynamical (entropy), and vibrational (optical and acoustic) properties could be assessed with reasonable accuracy.The simplified model employed in this work keeps in substance the empirical parametrization of the force field developed in ref. [68].Importantly, the calculated me-chanical and thermodynamic quantities are very close in the two cases, providing thus support to the present approach.

Surface Energy Effects on the Phase Stability
The calculation of surface energy related to the existence of dangling bonds requires the assessment of the potential energy excess for a system with free surfaces (model II) (Equation ( 13)).The molecular dynamics simulations show a linear increase of potential energies in both systems (bulk and with free surfaces) while increasing the size in both spin states (see Figure 3a), with a size-independent potential energy excess ΔE = E II − E I = 573.39(respectively 554.75) kcal mol −1 in the LS (respectively HS) state, leading to size-independent values of the surface energy  0 = 79 (respectively 71) mJ m −2 .For comparison, surface energies in both spin states have been estimated using an analytical approach based on the broken bonds methods (see Section SI, Supporting Information), which can be seen as a physico-chemical interpretation of the surface energy at the atomic scale.This approach consisted of estimating the energy of dissociation at T = 0 K of bonds, which have to be broken to create a surface having a certain crystallographic orientation.The values obtained from this analytical expression are very similar to those extracted from molecular dynamics simulations.This result is not surprising since the thermal expansion is negligible due to the presence of harmonic angular potentials in the FF and the crystalline symmetry of the bulk materials is preserved at the surface.This result demonstrates the accuracy of the slab method for the determination of excess quantities as surface energy.Obviously, if the structure and the topology of the surface is the result of an overly complicated chemical and structural reactional path, only the numerical methods will be relevant.
Let us discuss now the order of magnitude of  0 obtained in this work.To our best knowledge, these two values of  0 constitute primary estimations of surface energy in SCO materials.To date, no surface energy measurement has been performed on spin transition nanomaterials and experimental data concerning surface properties of molecular materials and organometallic complexes are in general scarce.It is nevertheless interesting to compare these numerical values with surface energies of organic and polymeric materials.For example, the surface energies of polymeric materials such as polytetrafluoroethylene (PTFE), poly(methyl methacrylate) (PMMA), and poly(vinyl chloride) (PVC), are respectively 20, 35, and 40 mJ m −2 at room temperature. [69]The molecular material Boron subphthalocyanine chloride (SubPc) has a surface energy of ≈48 mJ m −2 . [70]he surface tension or energy is a manifestation of intermolecular forces and is related to other properties derived from intermolecular forces such as the Young modulus or sound velocity.This is highlighted by the broken bonds methods: the surface energy has a linear dependence with the bond energy and also with the Young's modulus of these materials.The Young's moduli of PTFE, PMMA, PVC, and SubPc are respectively 0.5, 2, 3, and 5 GPa.This correspondence of Y with  0 is intuitive: the more rigid the material, the higher the surface energy.The Young's modulus of the simulated structure is around 6-8 GPa, that is to say approximately twice the Young's modulus of the polymeric materials previously mentioned.The surface energy is expected to be nearly twice than that of PVC that is 75-80 mJ m −2 in good agreement with the numerical results coming from MD simulations (the comparison is possible under the condition that the density of these different materials is close to each other).Overall, the different investigations of structural and mechanical properties of SCO materials have shown that the bulk and the Young's moduli of these bistable materials fall in the range of ≈4-13 GPa.According to the above discussion, surface energies in SCO materials are expected thus to be between 30 and 120 mJ m −2 .
Figure 3b shows the size evolution of the effective surface energy related to the broken bonds for both spin states.From the bulk where  eff → 0 (surface effects are negligible) to a size of L z = 1 nm, the effective surface energy increases to reach  eff ≈ 79 (respectively  eff ≈ 71 mJ m −2 ) in the LS (respectively HS) phases.As shown by the Equation ( 15), the relevant parameter for the study of the size dependence of the phase stability change is not the surface energy density itself, but rather the difference in the effective surface energy Δ of the two spin states.This difference is always negative whatever the material size (see Figure 3b, insert).As a consequence, the surface energy favors the HS phase and also the creation of HS surfaces.The difference of surface energy density appears to be relatively weak (Δ 0 ∕ 0 ≈ 10%) and according to Equation ( 15), the increase of unit cell area upon an LS to HS transition may contribute to the reduction of this gap.In the present case, the surface energy difference contributes to decrease the Gibbs free energy variation by ΔG = −4 meV for the thinnest film (see Figure 3c).It results in a downshift of the transition temperature of ≈4.5 K.A neg-ative value for the surface energy difference Δ should be thus the general trend, but it is fair to say that there has been no experimental observation reported.Indeed, we cannot exclude that the great complexity of organometallic compounds, generating complicated crystallographic structures for SCO crystals might induce different trends.For example, the shape of equilibrium crystals (or nanocrystals) results in a minimization of the total energy (together with a maximization of the total entropy) considering the physical constraints (temperature, pressure, chemical potential, size of the systems, crystallinity…).According to all these physical ingredients, the appearance of certain facets can be more favorable than others and this might not be the same in the two spin states.Clearly, intricate couplings can arise between the electronic spin state, crystal structure, particle morphology, and size, leading to deviations from the general trend.In such cases, numerical calculations of surface energy become obviously essential.
Previously, a nanothermodynamic approach has been developed to simulate the switching properties of spherical SCO nanoparticles. [71]This model takes into account the difference in terms of surface properties between the HS and LS phases by introducing a spin state dependent effective surface energy, which includes in a phenomenological manner the surface energy, interface stress, and surface entropy variations.The nanothermodynamic model allows to obtain a quantitative analysis concerning the size evolution of the phase stability since the different quantities can be deduced from calorimetric measurements, except the effective surface energy difference, which was chosen as an adjustable parameter.The value obtained in MD simulations (Δ 0 = −8 mJ m −2 ) is very close to the smallest value used in the nanothermodynamic model (Δ 0 = −5 mJ m −2 ).For this latter case, a slight downshift of the transition temperature with the size reduction was simulated without any presence a residual HS fractions at low temperature.It is important to keep in mind, however, that the contribution of surface energy phase stability is expected to be more important in nanoparticles than in thin films, due to a higher surface to volume ratio.
This work makes it possible to give a physical origin to this phenomenological parameter while waiting to have an experimental access to these quantities.Importantly, the nanothermodynamic model was able to simulate the experimentally observed transition curves of coordination nanoparticles only if this phenomenological parameter was increased significantly (≈−100 mJ −2 ).It means that surface/interface entropy and interface/surface stress appear to be the dominant contributions to surface effects involved in the phase stability change measured in some class of spin transition nanomaterials.

The Interface Stress Effects on the Phase Stability
Although the surface energy favors the HS phase and consequently decreases the transition temperature by 4.5 K, it seems, however, that the downshift of the transition temperature predicted by MD simulation is lower than that expected experimentally. [18]It is therefore necessary to examine other contributions to surface effects, such as the interface stresses related to the existence of an HS or LS residual fraction.It is important to mention, that strictly speaking, there should be a surface stress.The consideration of surface stress is possible if electronic contributions are explicitly included in the force field, which is not the case in this work.This can be achieved using many-body potentials like embedded atom methods (EAM) or modified embedded atomic methods (MEAM), which is out of the scope of this present work.However, the presence of a monolayer on the thin film can generate a solid-solid interface and a mechanical stress appears, due to the lattice parameter "misfit" between the HS and LS phases.The elastic energy associated to this interface could have a significant impact on the phase stability of SCO nanomaterials.The LS (respectively HS) system now consists of a LS (respectively HS) volume, covered with a size-independent single layer of HS (respectively LS) residual fraction (see Figure 1d, model III).
In order to quantitatively study the effects of the surface/interface stresses on the phase stability, the thin film has to reach a partial mechanical equilibrium (stationary state).This step corresponds to the elastic accommodation phenomenon, leading to the emergence of elastic deformations within the thin film.This allows to establish a mechanical equilibrium in the direction normal to the surface along the z-axis.The components of the stress tensor in this direction are close to zero, due to the presence of the open boundary conditions (thin film/monolayer/vacuum).On the other hand, a residual stress remains in the plane (x, y) parallel to the surface.Model III produces this situation of partial mechanical equilibrium for two reasons: 1) In model III, the simulation box is closed along the x-and yaxes.This can be assimilated to a virtual pressure "sensor".In order to "detect" (calculate) the stress in the xOy plane, the "sensor" (the simulation box) must be in contact with the system.At constant volume, the simulation box is indeformable, consequently, the stresses exerted in the plane parallel to the surface do not relax completely by transforming into elastic deformations.The unrelaxed stresses in the system are called residual stresses arising from the existence of the HS/LS interface.By using model III, we therefore have access to these residual HS/LS interface stresses, denoted by the following σxx and σyy , having non-zero values.
2) Model III allows to obtain the equilibrium condition at the surface ( σiz ≈ 0; with i = x, y, z).This condition is necessary to relax the stresses which are not the consequence of the presence of the HS/LS interface.Indeed, the void generated in the z-direction creates free (HS and LS) surfaces.The change in the volume can only be achieved through a variation of the slab thickness.
Figure 4a,b shows the relaxation of the stress tensor components, for an LS system (respectively HS system), constituted with nine LS layers (respectively nine HS layers) and a residual HS (respectively LS) fraction located on the top of the LS (respectively HS) surface, as a function of the number of iterations of the gradient conjugated minimization algorithm (GC). [59]ue to the absence of migration (diffusion) of atoms at the interface with the size reduction, the local properties and stresses of the HS/LS interface are expected to be size independent.Before reaching a partial mechanical equilibrium (<2000 GC steps for the HS state and <700 GC steps for the LS state), the stresses have a noisy and random profile, and consequently, the related interface properties cannot be extracted.At the mechanical equilibrium, the stress components reach a plateau in all directions of space.In particular, we obtain the residual interface stresses ( σxx and σyy ), which are positive in both spin states, due to the isotropic structural properties of the slab.As expected, the stress along the z-axis is zero ( σiz ≈ 0).According to the thermodynamic conventions, a positive component is assigned to a compressive stress, while a negative one corresponds to a tensile stress.An unexpected result is the negligible shear stress ( σxy = 0), meaning that there is no relative displacement of the surface atoms with respect to the atoms of the volume along the interface.This observation is mainly attributed to the choice of a weak value for the lattice parameter misfit between the HS and LS structures (m = 0.8%) for the numerical calculations (see Section 2.4).This numerical misfit m is much lower than the value generally observed experimentally (≈3-5%).The choice of a low elastic misfit makes it possible to work with thin films which remain flat, thus having an infinite curvature.The Laplace like overpressure is then removed and the only constraints existing in the thin film are tangential stresses.(We show some consequences of a large elastic misfit in the Section SII, Supporting Information.) In order to discuss the effect of elastic stresses on the phase stability and to compare with the contribution of surface energy, these interface stresses have to be expressed in (J m −2 ) rather than in (MPa).This avoids an impractical dependence of these quantities with the geometry of the simulation box.To this aim, the interface stress sij is calculated according to 20, where the length of the simulation box h has been set to ≈47 nm.The interface stress value in the LS (respectively HS) case is sxx = syy = 94 mJ m −2 (respectively 45 mJ m −2 ), close to the values obtained for the Evolution of isotropic mean residual stresses σij with the number of steps of the gradient conjugated method (GC) [59]   ) This point can be qualitatively understood by the Shuttleworth relation: s =  0 + ∕, [72] which relates the work to create a surface  0 with work to deform a surface/interface s, whose difference is written as follows where the second term expresses the difference in variation of the surface energy with respect to a strain ϵ exerted in the plane of the surface between the two spin states.This term is negli-gible for liquids because they cannot be deformed at a constant number of atoms.In this case, the work to deform a surface is equal to the work to create a surface.Hence, conceptually, there is no difference between surface tension and surface energy in liquids.However, in SCO materials, this derivative is not negligible and a surface/interface stress or tension is no more equivalent to the surface energy.It now appears obvious that the relatively large difference in the interface stresses between the two spin states (−49 mJ m −2 ), comes from the second term of the Shuttleworth equation, which is certainly negative.Indeed, (∂/∂ϵ) LS is positive in the LS phase, because this latter undergoes a positive strain or an elongation of the LS lattice parameter, due to a tensile stress coming from the HS residual layer, whereas in the HS phase, (∂/∂ϵ) HS is negative, because the HS slab undergoes a compression and then a negative strain is created parallel to the interface.Despite the small elastic misfit employed in our numerical simulations, the difference in the HS/LS interface stresses is larger than the surface energy difference.This confirms that the phase stability of SCO nano-systems is more affected by the effects of solid-solid interfaces than solid-vacuum surface effects.
As mentioned earlier, the stresses in the HS and LS systems are negligible, except for the residual stresses along the x-and y-directions.Then, the stress matrix of the HS/LS interface sij contains only the two diagonal terms (s xx and syy ).Consequently, the computation of the effective interface stress is reduced to the trace of the matrix elements: seff = sxx εxx + syy εyy (27)   We recall that the interface stresses and the corresponding strains are isotropic.Then we obtain The value of the residual stress is obtained at constant volume, which is not the case for the strains.One solution is to use the continuum mechanics in order to transform the previous equation into an expression, which only depends on the stress.To this aim, we consider that the whole system as a continuous medium exerts a permanent residual stress on the boundaries of the simulation box with a low strain ratio.In this case, the stresses and strains can be linked by the following relation : εxx = 2(1 − )s xx ∕L z Y. [73] Then, the effective interface stress becomes This expression results in an increase of the effective interface stress by almost 3.5 (respectively 1.0) mJ m −2 in the LS (respectively HS) state at the smallest sizes, which corresponds to a maximum difference of Δs eff = −2.5 mJ m −2 between the two spin states (see Figure 4c).A decrease in the free enthalpy difference ΔG = −4.1 meV is then attributed to the effects of interface stress.If we consider interface stress effects only in the LS system (G HS int = 0), which is a common experimental observation, the free enthalpy difference drops to reach ΔG = −5.5 meV.As a consequence, the transition temperature decreases from 200 to 193.9 K in this latter case (Figure 4d).The contribution of interface effects to the free enthalpy value and phase stability appears thus greater than the contribution due to surface energy.Overall, the surface energy due to the dangling bonds induces a decrease of 4.5 K in the transition temperature, while the interface stress due to the presence of a HS residual fraction induces a decrease of 6.1 K in the transition temperature.Then, the sum of the two contributions gives a decrease of 10.6 K in the transition temperature of the simulated thin film, which remains lower than the experimental observations (between 15 and 20 K). [18] This is certainly due to the weak elastic misfit introduced in numerical simulations, which is three times lower than the lattice parameter misfits generally observed between HS and LS crystal structures.The current theoretical and numerical tools do not allow to extract quantities of interfaces with high elastic misfits.
A theory of non-linear elasticity must be developed and new numerical tools must be set up to allow future studies of the influence of highly deformed interfaces on the phase stability of SCO nanomaterials.

Structural Anisotropy Effects
In Section 3.3, the numerical calculations revealed the impact of isotropic interface stresses, which lead to a decrease of the transition temperature, giving an origin to the change in phase stability observed in certain classes of spin transition nanomaterials. [14]owever, the experimental results on thin films of the complex [Fe II (HB(tz) 3 ) 2 ](tz = 1,2,4-triazol-1-yl), show an opposite behavior of the transition temperature with size reduction.In particular, a 3 K increase of the transition temperature accompanies the reduction of the thin film thickness. [20]In this particular case, the strong structural anisotropy was invoked to explain the singular behavior.By virtue of the Curie principle, this anisotropy must have repercussions on the physico-chemical properties of this thin film, notably on the mechanical properties.The size evolution of spin transition curves of this thin film was simulated using a nanothermodynamic model, in which the phenomenological surface parameter was ad hoc adjusted to favor the LS phase and therefore induce an increase in the transition temperature.This nanothermodynamic model takes into account the different surface contributions, including the effective elastic stresses.The numerical approach, developed in this work, makes it possible to go further and to quantify the role of these different contributions in the observed phenomena.We propose to simulate the anisotropic stress/strain effects of the HS/LS interface on the phase stability using the slab method.The HS/LS interface is assumed to play a similar role than a substrate/thin film interface, even if the mechanical properties of the substrate are very different from those of the thin film.In particular, the stiffness of the substrate is generally rather high and therefore it remains flat.
It has been reported that during the LS→HS phase transition, the above mentioned thin film undergoes an anisotropic volume change.The lattice parameter a (abitrary oriented in the x-direction) contracts with a relative variation of m x = Δa∕a = −2.3%, while the lattice parameter b (respectively c) oriented in the y-direction (respectively z-direction) exhibits an expansion whose associated strain is m y = Δb∕b = +1% (respectively m z = Δc∕c = +5.6%).A contraction of m x = −0.08% in the x-direction is imposed in the numerical simulations, twice as large as the expansion along the y-axis (m y = 0.04%), so that the experimentally observed ratio m x /m y is preserved.(The misfit values are lower for the reasons already mentioned previously.)We have chosen negligible expansion along the z-axis, because it does not contribute to the interface properties parallel to the plane (x, y).
Self-evidently, Figure 5a,b shows that the interface stresses are now anisotropic in the two spin states, with the tension stress σyy being smaller than the compression stress ( σxx ) due to the relatively smaller expansion of the system in the y-direction.In addition, the tensile and compressive stresses in the HS system are greater in magnitude than the stresses in the LS system.This is because in the HS system there are nine HS layers, which contract/expand, largely influencing the HS/LS interface properties, while in the LS system there is only one HS layer on the surface which undergoes this phenomenon.Using Equation (23), we Evolution of mean anisotropic residual stresses σij with the number of steps of the gradient conjugated method (GC) [59] Figure 5c shows that the stress anisotropy favors the LS phase by increasing the Gibbs free energy difference by almost 7.9 meV.As a result, the transition temperature increases by 8.2 K for the thinnest film (see Figure 5d).The surface energy is by definition always positive, while the interface stress can be positive or negative.Taking into account the strong anisotropy of the structural properties of SCO materials, this work suggests the possible control of the phase stability through the modulation of the interfacial mechanical properties or, the other way around, the possible prediction of the size evolution of the phase stability.

Conclusions
A molecular dynamics approach has been employed to investigate size reduction effects on the phase stability of SCO nanomaterials.Previous works carried on size reduction effects in SCO materials introduced a quantity allowing the phenomenological analysis of the predominance of surface phenomena at the nanometer scale: the effective surface energy.By considering this quantity dependent on the spin state and by injecting it into a nanothermodynamic model, it was then possible to reproduce the evolution of the switching properties with the size reduction.However, this quantity actually combines all the contributions associated with various physico-chemical phenomena likely to occur at a surface/interface, including surface energy, surface entropy, interface stress, and even enthalpies of chemical reactions.In this work, we were able to separate the surface and interface contributions, by giving definitions coming from the solid-state physics and continuum mechanics: the surface energy is the work necessary to create a surface, while the surface/interface stress is related to the work required to deform a surface/interface.The use of these definitions in molecular dynamics simulations allows us to have preliminary estimations of surface energy and interface stress quantities in SCO materials.From a quantitative approach, which consists in extracting the excess of the potential energy of a cubic structure with octahedral patterns, we deduce the energy of {100} surfaces, which is 79 mJ m −2 for the LS state and 71 mJ m −2 for the HS state, in good agreement with the first estimations coming from the scarce available experimental data.These numerical values have been injected into a thermodynamic model in order to study the thickness dependence of the phase stability of a SCO slab.
The surface energy difference between the two spin states (≈− 8 mJ m −2 ) due to the dangling bonds induces a decrease in the transition temperature of 4.5 K, which is lower than that observed experimentally (between 15 and 20 K).We then considered, in addition to surface effects, the effects of mechanical stresses at the HS/LS interface.In the limit of weak misfits, we calculated the isotropic stresses of the HS/LS interface in the two spin states, which is 94 (respectively 45) mJ m −2 in the LS state (respectively, HS state).Despite a weak misfit used in the numerical calculations (to preserve the validity of the elastic theory employed in this work), the interface stresses are of the same order of magnitude as the surface energy values.This suggests that with greater elastic misfits (≈3-4%), the interface contributions to the change of the transition temperature should become dominant.The interface stresses decrease the transition temperature of the simulated thin film by 6.1 K.By summing both surface energy and interface stress contributions, we show that the transition temperature of the SCO slab decreases by almost 10.6 K.
In summary, we can deduce that the mechanical interface stresses exerted by solid environments are responsible for the drastic change in phase stability and transition temperature with the size reduction.Moreover, we have demonstrated that the anisotropy of mechanical stresses can produce an increase of the transition temperature.This result confirms the possible control of the phase stability of spin transition nanomaterials through interface engineering.We shall underline that the present work represents only an initial approach, using a simple model based on coordination octahedra, to compute surface thermodynamic quantities in SCO nanomaterials and it is important to restate its main limitations as follows.Notably, using the present approach, the surface stress cannot be assessed and the interfacial stress is certainly underestimated due to the imposed small misfit.In certain systems, the surface/interface entropy and interface energy might also represent relevant contributions, contrary to our model system where the role of these quantities on the phase stability can be neglected.The analysis of these issues using more realistic models and more sophisticated force fields represents thus an important scope for future works.Finally, let us note also that there are certainly some aspects of the surface/interface phenomena encountered in these materials, such as the development of interfacial dislocations for large misfits, which should be addressed by techniques other than MD simulations.

Figure 1 .
Figure1.a) 3D captured image of the cuboid SCO structure extracted from the ovito visualization tool,[45] illustrating the metal (Me) coordinated by six ligands (L).b) 2D schematic representation of the model system I (model I) with closed periodic boundary conditions in the three directions of space.The transparent atoms are the reflection of atoms inside the simulation box.Within the octahedra: the dashed-dotted and dotted straight lines represent the harmonic interactions between LL-intra and Me-L, respectively.Between the octahedra: the solid straight line represents the Lennard-Jones (LJ) potential connecting neighboring octahedra and the dashed curved lines symbolize the harmonic angular interaction between three consecutive ligands.c) Model II : the slab method is applied by increasing the size of the simulation box by a length h in the z-direction, creating an LS or HS slab with two equivalent solid-vacuum (100) surfaces.In this model, only coordination defects and the corresponding surface energy, produced by the creation of two solid-vacuum surfaces, affect the phase stability of SCO nanomaterials.The red filled square whose corners are constituted by four metal atoms, represents the elementary unit cell of area a s = l 2 , where l = (2a MeL + a LL-inter ) is the lattice parameter.d) In model III, one of the two solid-vacuum surfaces shown in model II is replaced by a solid-solid HS/LS interface.In this case, the LS system is made up of LS multi-layers and a single-layer HS residual fraction on the surface.The number of LS layers (N LS ) varies with size, while on the contrary, the number of HS molecules on the surface remains constant:N HS = 1.At the HS/LS interface: the misfit of lattice parameters between the two subsystems generates biaxial residual stresses along the x-and y-directions.At the LS surface, a surface energy exists due to the broken LS bonds, which is calculated by the model II.Then, the stability of the LS phase in model III is not only affected by the surface energy as in model II, but rather by the sum of surface energy and interface stresses at the LS surface and HS/LS interface, respectively.A similar description for the HS system made up of HS multi-layers and a single-layer LS residual fraction on the surface can be achieved.

Figure 2 .
Figure 2. Normalized vDOS of the simulated bulk material in the HS (red circles) and the LS (blue stars) states.Insert: g(E)∕E 2 as a function of E for the two spin states.The lowest energy part of the curves (black lines) corresponds to the range where the Debye approximation is valid.[56]

Figure 3 .
Figure 3. Size evolution of : a) the potential energy for the two systems, bulk (model I) and with free surfaces (model II), in the two spin states; b) the effective surface energy related to broken bonds  eff (Equation (14)), in the HS (red circles) and LS (blue stars) states.Insert: The negative difference of the effective surface energy Δ eff =  HS eff −  LS eff favors the creation of HS surfaces; c) the Gibbs free energy difference related to the surface energy ΔG s = HS s − G LS s < 0 (Equation (15)); and d) the transition temperature in presence of surface energy related to broken bonds (interface elastic effects are absent, ΔG int = 0 meV).

Figure 4 .
Figure 4. Evolution of isotropic mean residual stresses σij with the number of steps of the gradient conjugated method (GC)[59] in a) the LS state and b) the HS state.c) Size evolution of the effective interface stress seff , in the HS (red circles) and LS (blue stars) states.Insert: The difference of the effective interface stress Δs eff = sHS eff − sLS eff < 0 favors the HS surface.d) Size evolution of the Gibbs free difference (ΔG int = G HS int − G LS int < 0, 25) related to the interface stress in the case of HS and LS residual monolayers situated, respectively, on the top of an LS and an HS slab.The red dots represent the case of the HS residual monolayer effect on the top of a LS system (G HS int = 0).e) The interface stress effect on the transition temperature (Equation (5)) with the size reduction in the case of LS and HS residual monolayers (black stars) or only HS residual monolayer (red dots).
Figure 4. Evolution of isotropic mean residual stresses σij with the number of steps of the gradient conjugated method (GC)[59] in a) the LS state and b) the HS state.c) Size evolution of the effective interface stress seff , in the HS (red circles) and LS (blue stars) states.Insert: The difference of the effective interface stress Δs eff = sHS eff − sLS eff < 0 favors the HS surface.d) Size evolution of the Gibbs free difference (ΔG int = G HS int − G LS int < 0, 25) related to the interface stress in the case of HS and LS residual monolayers situated, respectively, on the top of an LS and an HS slab.The red dots represent the case of the HS residual monolayer effect on the top of a LS system (G HS int = 0).e) The interface stress effect on the transition temperature (Equation (5)) with the size reduction in the case of LS and HS residual monolayers (black stars) or only HS residual monolayer (red dots).

Figure 5 .
Figure5.Evolution of mean anisotropic residual stresses σij with the number of steps of the gradient conjugated method (GC)[59] in a) the LS and b) the HS phase.c) Size evolution of the Gibbs free energy difference (ΔG int = G HS int − G LS int > 0, 30) related to the interface anisotropic stress.d) Size evolution of the transition temperature with the effect of anisotropic interfacial stresses.
Figure5.Evolution of mean anisotropic residual stresses σij with the number of steps of the gradient conjugated method (GC)[59] in a) the LS and b) the HS phase.c) Size evolution of the Gibbs free energy difference (ΔG int = G HS int − G LS int > 0, 30) related to the interface anisotropic stress.d) Size evolution of the transition temperature with the effect of anisotropic interfacial stresses.

Table 2 .
Summary of structural, thermoelastic, and mechanical properties of the bulk material in the two spin states.