Role of Surface Effects in the Vibrational Density of States and the Vibrational Entropy in Spin Crossover Nanomaterials: A Molecular Dynamics Investigation

Abstract

Size reduction effects on the lattice dynamics of spin crossover (SCO) thin films have been investigated through molecular dynamics (MD) simulations of the density of vibrational states. The proposed simple model structure and reduced force field allows us to obtain good orders of magnitude of the sound velocity in both spin states and takes into account the contribution of free surfaces in the vibrational properties of very thin films (below a thickness of 12 nm). The slab method issue from the field of surface physico-chemistry has been employed to extract surface thermodynamic quantities. In combination with the related slab-adapted method, the slab approach provides a powerful numerical tool to separate surface contributions from finite-size effects. Due to the relatively low stiffness of SCO materials, the lattice dynamics seems to be governed by surface instead of confinement effects. The size evolution of thermodynamic quantities is successfully reproduced, especially the increase of the vibrational entropy with the size reduction, in good agreement with experimental observations.

1. Introduction

Undoubtedly, switchable molecular phenomena are destined to play a key role in various societal applications [1,2,3]. In this context, we focus on spin crossover (SCO) compounds, which are able to change their electronic configuration from a low spin (LS) to a high spin (HS) state by controlled thermal, mechanical, optical, electric or magnetic energy supplies [4,5,6]. The electronic reorganization inherent to the spin state change leads to drastic modifications of the physico-chemical (magnetic, vibrational and electronic) properties. At the solid state, collective elastic behaviors, usually called cooperativity, originated from the molecular structural difference between the two spin states are able to induce first order phase transitions with hysteretic phenomena [7,8]. Owing to these properties, SCO materials are expected to be candidates for future applications in data processing, photonic or spintronic devices as well as actuators [9,10,11].

To this aim, a crucial step is the synthesis and the conception of SCO nano-objects [12] under various forms such as nanoparticles, thin films, nanopatterned structures or nano-composites [13,14,15]. However, size reduction effects can drastically affect the switching properties in an unwanted manner: changes in the phase stability, incomplete phase transitions, loss of the hysteresis, etc. [16,17]… Both experimental and theoretical investigations of size reduction effects have pointed out the strong influence of surface properties, the nature of the interface as well as the impact of the external environment on the switching properties on SCO nanomaterials [18,19,20]. In particular, spin transition phenomena in core-shell nanoparticles suggest a possible control of phase stability and spin transition curves through a modulation of elastic/mechanical stress at the interface [21,22,23,24,25,26,27,28,29]. Therefore, the knowledge of both mechanical and vibrational properties of SCO materials as well as their evolution with the size reductions are of paramount importance in order to understand the spin transition phenomenon at the nanoscale [30]. Investigations of lattice dynamics through nuclear inelastic scattering (NIS) experiments allows us to assess the partial density of vibrational states (vDOS) of the iron ions and to extract elastic moduli as well as thermodynamic quantities [31,32,33]. Notably, NIS experiments on SCO coordination nanoparticles of [Fe(pyrazine)(Ni(CN)${}_4$)] and Prussian blue analogues show that size reduction leads to an increase of vibrational entropy and sound velocity, while the internal energy remains almost unchanged [34,35].

Numerical investigation of the lattice dynamics of SCO materials through molecular dynamics simulations is also highly relevant [36]. Until now, however, vibrational properties of SCO coordination nano-objects have been unsatisfactory modeled since surface effects have not been fully taken into account [36]. In this paper, we use the so-called slab method [37,38] to investigate surface vibrational and associated thermodynamical properties of nanoscale SCO films. As a first attempt, we use here a simple model structure and an ad-hoc force field which was parameterized to obtain the good order of magnitude of the Debye sound velocity and of the frequencies of optical modes.

2. Molecular Dynamics Simulations and Slab Method

We simulate a cuboid crystal of dimensions $L_x$, $L_y$ and $L_z$ with a cubic symmetry lattice constituted with octahedral patterns, which mimic in a simple manner the structure of SCO coordination complexes. Previous works were particularly successful in a quantitative modeling of optical phonon spectrum of SCO complexes through numerical simulations combining MD and ab initio calculations [39,40,41]. However the price to pay is to consider very small systems (several patterns/molecules) so that the acoustic part of the density of vibrational states (i.e., phonon modes with long wavelength) has to be disregarded. In the present work, we are interested in the overall phonon spectrum and therefore, we need sufficiently large systems to observe the long range collective vibrational modes. The octahedra are constituted of a “heavy” atom ($M_{Me}=57u$, u being the unified atomic mass unit), corresponding to the metallic iron cation surrounded by six nitrogen atoms ($M_L=14u$), modeling the six ligands in a simple manner. The different atoms of the octahedron are bonded by two kinds of harmonic pairwise potentials called $V_{Me-L}$ and $V_{LL-intra}$ representing Metal-ligand (Me-L) and intramolecular ligand-ligand (LL-intra) interactions, respectively. The equilibrium distances are chosen in order to keep the perfect octahedral symmetry of the pattern. In comparison with the force field employed in reference [36], the intermolecular interactions has been modified. In order to keep the cubic structure and to limit the rotation of octahedrons in the crystal packing, a harmonic angular potential between three consecutive ligands (see Figure 1a) has been introduced. This avoids the occurrence of anomalous unit cell thermal contraction.

In reference [36], the angular potential between two ligands and the metallic center strongly altered the low frequency part of the phonon spectrum, leading to high values of the Debye sound velocity in disagreement with experimental observations. In the present work, the choice for the angular potential aims to limit this non-physical effect. The octahedra are connected by Lennard-Jones pairwise interaction potentials ($V_{LL-inter}$). An anharmonic interaction potential is essential for taking into account surface effects in the lattice dynamics [42]. Indeed, surface creation can be seen as a bond-dissociation process, which cannot be achieved with harmonic potentials, which are only able to simulate small atom displacements around an equilibrium position. Moreover, the employment of harmonic potentials for intermolecular bonds prevents the system to relax to the more stable structure and residual elastic stress can subsist in the crystal packing. The different numerical parameters of the force fields employed for the MD simulations are summarized in

Table 1. Summary of the numerical parameters used in the molecular dynamics simulations.

Lennard-Jones Potential	Cohesion Energy (kcal.mol${}^{-1}$)	Equilibrium Distance (Å)
L-L (inter-octadra) (LS)	5.73	1.0
L-L (inter-octadra) (HS)	5.55	1.0
Harmonic Potential	Force Constant(kcal.mol${}^{-\mathbf{1}}$Å${}^{-\mathbf{2}}$)	Equilibrium Distance(Å)
Me-L (LS)	90	2.0
Me-L (HS)	60	2.1
L-L (octahedron) (LS)	60	2$\sqrt{2}$
L-L (octahedron) (HS)	40	2.1$\sqrt{2}$
Angular Potential	Force Constant (kcal.mol${}^{-\mathbf{1}}$rad${}^{-\mathbf{2}}$)	Equilibrium Angle(Rad)
L-L-L	2	$\pi$

. We recall that these parameters do not correspond to a specific SCO complex. Force constants and equilibrium distances are orders of magnitudes deduced from Raman spectroscopic and X-ray diffraction data. (See more details about the assessment of the force constants in Supplementary Materials Section S3, except the angular potential, which is not related to experimentally measurable physical quantities, vide supra).

In the following, we consider a crystal in contact with a thermal bath and a pressure reservoir (isobaric isothermal ensemble). To this aim, we perform MD simulations using a (T, P, N) Nosé-Hoover thermostat-barostat [43,44]. This deterministic thermostat-barostat modeling allows the control of the expansion/contraction of the simulation box, while keeping the real time from Newton’s laws of motion. Positions and velocities of each atoms are obtained by integration using the velocity Verlet schema implemented in the optimized molecular dynamics simulator LAMMPS [45] with an elementary time step of $\delta t=1$ fs. In order to extract the vDOS, the velocity auto-correlation function ($\gamma \left(\tau \right)$) is calculated and collected during 3000 steps ($\tau =3$ ps) from the MD simulation after 200,000 steps of thermalization /relaxation (${\tau }_{relax}=0.2$ ns). The vDOS noted $g\left(E\right)$ is obtained by a Fourier transform of $\gamma \left(\tau \right)$. More precisely, the vDOS related to the metallic center, called hereafter partial vDOS, will be considered in this work for a better comparison with experimental NIS results. (See more details on the calculation of $\gamma \left(\tau \right)$, $g\left(E\right)$ and the extraction of the thermodynamic quantities in Supplementary Materials, Section S4). To investigate surface effects, we employ the so-called slab method. Periodic boundary conditions are applied in the three spatial directions ($x-$, $y-$ and $z-$ axis), but along the z direction, the simulation box is increased, thus producing a vacuum of size $L_{vacuum}\approx 10L_z$, which makes it possible to generate two free surfaces (see Figure 1b, right). Therefore the lengths $L_x$ and $L_y$ are fixed to 5 nm whereas $L_z$ is the only dimension, which can be tuned. The slab method has the advantage to simplify the analysis of surface effects, because there is no surface whose normal is not co-linear to another one, avoiding thus complicated surface couplings, typical in nanoparticles. The system looks rather as a thin film. Finally, in order to identify surface effects on the vibrational properties of the thin film, MD simulations are performed with a simulation box size, which matches perfectly the modeled material dimensions, i.e., without the presence of vacuum (see Figure 1b, left). This MD simulation is called slab-adapted method and allows to obtain the spectra of bulk materials.

3. Results

First, the partial vDOS of a bulk material is simulated in the two spin states using the slab-adapted method at $P=1$ atm and at $T=200$ K. We consider a system with a thickness $L_z=$$25$ nm (5000 octahedrons), which is large enough so that size effects are negligible. Indeed, even if the slab adapted method assures the absence of surface effects and simulates a virtually infinite medium, the system needs to contain a sufficient number of lattice degrees of freedom to obtain a large number of vibrational states. This aims to simulate energetic band structures similar to macroscopic materials. Nevertheless, the absence of finite size effects can only be ascertained after an investigation of the size evolution of the vibrational properties (vide infra). The corresponding HS and LS partial vDOS of the bulk material are depicted in Figure 2a. When we go from the LS to the HS state, a redshift of the overall vDOS is observed, corresponding to the softening of both optical and acoustic modes. Unlike experimental observations, the vDOS conserved the same shape in the two spin states because the perfect octahedral symmetry is preserved in the HS and LS states. The high frequency peak observed around 21 meV for the HS spectrum (resp. 26 meV for the LS spectrum) corresponds to the stretching mode Me-L of the HS (resp. LS) octahedron, in fair agreement with typical value of such optical modes measured by NIS or Raman spectroscopy techniques (${\nu }_{HS}\approx 176$ cm ${}^{-1}$, ${\nu }_{LS}\approx 208$ cm${}^{-1}$) [30,46]. The extracted values (see

Table 2. Physical quantities extracted from the partial vDOS spectra of the bulk material in the HS and the LS states using the methods described in reference [32].

	LS	HS
$s_{vib}$ ($k_B$)	$1.67\pm 0.01$	$1.90\pm 0.01$
$u\left(\mathrm{meV}\right)$	$30.2\pm 0.1$	$28.90\pm 0.01$
$<C>$ (N/m)	$63\pm 5$	$43\pm 5$
$\rho$ (kg/m${}^3$)	$1804.8$	$1604.5$
$v_D$ (m/s)	$2070\pm 10$	$2004\pm 10$

) of the mean force constants $\langle C\rangle$ characterizing the high frequency part of the vDOS spectra confirms a softening of the optical modes. While the relative variation of the $\langle C\rangle$ (${\langle C\rangle }_{LS}$ /${\langle C\rangle }_{HS}\approx 1.5$) matches with the expected frequency shift of the stretching mode Me-L upon the spin transition, the absolute values of $\langle C\rangle$ are smaller than the experimental observations. Indeed, a large part of the contribution to the value of $\langle C\rangle$ comes not only from the stretching modes Me-L, but also from the high frequency vibrational modes produced the ligand. Obviously, the “unified atom” description of the ligand does not allow to assess this part of the vibrational spectrum. The Debye sound velocities $v_D$ of both spin states have been calculated from the low frequency part of the vDOS spectra (see Figure 2b) [32]. The decrease of $v_D$ upon a LS to HS transition is retrieved, denoting a softer lattice in the HS state [34]. It is important point to underline the good match of the calculated values of $v_D$, with the experimental data ($v_{LS,exp}=2063\pm 27\mathrm{m}{\mathrm{s}}^{-1}$, $v_{HS,exp}=1933\pm 20\mathrm{m}{\mathrm{s}}^{-1},$ [35]) obtained for different SCO complexes. The lattice stiffness is realistic, which is essential to investigate size reduction effects. Indeed, the bond energies and also the lattice rigidity play a key role in surface properties and phonon spatial confinement. We shall note that the vibrational entropy difference between the HS and LS state $\Delta S=S_N^{HS}-S_N^{LS}$ is positive due to the higher density of vibrational states in the HS state. Here, the calculated entropy variation between the two spin states is $1.97\left(1\right)$ J.K.${}^{-1}$mol${}^{-1}$, close to the value obtained by NIS measurements ($\approx 4$ J.K.${}^{-1}$mol${}^{-1}$,35]). Whether obtained experimentally or numerically, the entropy variation is much lower than the expected value for a spin crossover materials (40–80 J.K.${}^{-1}$mol${}^{-1}$). This is likely due to the fact that we analyze only a part of the vibrational density of states, which are allowed by the selection rules of NIS.

In the next step, the lattice dynamics of SCO materials have been simulated using the slab-adapted and slab methods for different film thicknesses $L_z$. If we compare the vDOS spectrum of the previously obtained bulk material with the thinnest film ($L_z=2.5$ nm) (see Figure 3), we note that there is no drastic modification in appearance of the vDOS with the size reduction. The general shape of the curve is preserved for very thin systems. This point is not surprising since in this simulation the crystalline symmetry and the geometry of a perfect octahedron is not strongly affected by approaching the nanometric scale. In particular the position and the area of the optical peak around 21 meV are identical for the bulk material and the thin film. However, if the vDOS is more carefully inspected, some differences can be distinguished. The vDOS of the thin film seems to be more noisy, a signature of a rarefaction of the number of accessible vibrational modes. Moreover, the occurrence of several peaks in the acoustic part of the spectrum, around 5 meV can be observed (upward black arrow on Figure 3). Size reduction implies the occurrence of two phenomena: finite-size effects (reduction of lattice degrees of freedom), which leads to a discretization of the energy spectrum and surface effects related to a breaking of the invariance by translation in space of the system properties. Finite-size effects produce a discontinuous vDOS spectrum and the occurrence of an acoustic gap at low frequencies [47]. In the present work, these phenomena are not observable due to the fact that the surface $L_x\times L_y$ remains sufficiently large, so that the thin film conserves a large number of lattice degrees of freedom.

On the other hand, free surfaces imply the presence of additional vibrational modes, called Rayleigh modes, which propagate along the surface [37]. These modes are located in the low frequency part of the vDOS, close to the bulk acoustic modes. The amplitude of Rayleigh waves is attenuated within the core of the thin film as one moves away from the surface. These surface modes subsist in the continuum limit (large values of $L_z$) and their existence can be predicted by Navier’s equation, considering an elastic continuum medium with zero elastic stress at the surface. Therefore, these surface modes appear as soon as a free surface is created, whatever the characteristic of the intermolecular potential (harmonic or anharmonic). The anharmonic nature of the intermolecular vibrations increases the amplitude of Rayleigh modes and adds new surface modes since octahedra located at the surface can exhibit larger displacements around their equilibrium positions. The main advantage of the present approach is the possibility to separate surface effects from finite-size effects. The contribution of these surface modes in the extracted thermodynamic quantities $a_N$ (a standing for vibrational entropy S or internal energy u) can be also determined by considering as a reference the state functions corresponding to the infinite system (simulated with the slab-adapted method) $a_N^b$. Surface quantities can be then seen as excess thermodynamic quantities $a_N^s$ [48]:

$$a_N^s=\frac{1}{2}(a_N-a_N^b),$$

(1)

where the factor $1/2$ takes into account the fact that two facing surfaces are created by cleaving the bulk material. The vibrational internal energy $u_N$ and the mean force constant $\langle C\rangle$ remain virtually the same with the size reduction (see Supplementary Materials, Figures S2 and S3). No major difference can be detected between the calculations performed on a system with periodic boundary conditions and a thin film. It means that the contribution of surface vibrational modes due to the presence a free surface is not significant for these two quantities. This result can be understood as follows: the definitions of $\langle C\rangle$ and $u_N$ [32] show that the main contribution of these two quantities comes from the high frequency part of the vDOS, which is not affected by the size reduction. The situation seems to be different for the simulation of nanoparticles for which an increase of both quantities has been observed with the size reduction [36,49]. This is likely due to a spatial confinement, whose discretization effects tends to increase the area under the vDOS spectrum in the optical part to the detriment of the area occupied by the acoustic part. At the opposite, we observe an increase of the vibrational entropy with the size reduction for both spin states (see Figure 4a,b), which is clearly imputed to the presence of surface modes in the low frequency part of the vDOS spectrum. The main contribution to the vibrational entropy is given by the acoustic part of the phononic spectrum. Indeed, from the thickest film to a film of thickness $L_z=12$ nm, the total vibrational entropy of the film $S_N$ and the bulk vibrational entropy $S_N^b$ match very well. In other words, the surface entropy $S_N^s$ is zero since surface effects are negligible. Below $L_z=12$ nm, however, the surface vibrational entropy $S_N^s$ increases, while the bulk vibrational entropy $S_N^b$ remains virtually constant. (The variation observed for the thinnest film in the system without free surfaces can be likely due to a numerical artifact attributed to a problem of self-interactions.) The amplitude of Rayleigh waves exponentially decreases within the core of the nanomaterial as one moves away from the free surfaces. For large systems, the lattice dynamics of the thin films are unaffected by surface effects since Rayleigh’s waves are negligible at several layers (c.a. 5 nm) under the surface. When the thickness is of the same order of magnitude than the penetration depth of surface waves, the vibrational properties of the whole system are impacted by surface effects. Indeed, the lower coordination number produced by the presence of free surfaces allows the surface octahedra to vibrate with a greater amplitude and the existence of surface acoustic modes leads to a propagation of this perturbation in the core of the material. This can be seen as a local increase of the temperature (or a decrease of the Debye temperature) resulting in an increase of the entropy.

Therefore, we can conclude that there are two competing effects in the size evolution of the vibrational entropy. In the case of free surfaces, the increase or the decrease of the entropy is dependent on the strength of the spatial confinement (dimensionality and morphology of the nano-object) and on the surface Rayleigh modes whose vibrational amplitude is governed primarily by the stiffness and the crystallographic properties of the nanomaterial. The stiffer the nanomaterial, the smaller the vibration amplitude of the Rayleigh modes. As a result, the thermodynamic properties of a nano-system with a high sound velocity are governed by spatial confinement and the vibrational entropy is expected to remain constant (thin films) or to decrease (nanoparticles), since the discretization leads to an impoverishment of the number of accessible vibrational states. At the opposite, if the amplitude of surface modes is significant in softer nanomaterials, then additional surface modes are able to counterbalance the effect of spatial confinement and an increase of the entropy can be observed. In a previous work [36], the spatial confinement was dominant since the sound velocity was too high and a decrease of the entropy was observed with the size reduction. The experimental observation of an increase of the vibrational entropy in SCO coordination nanoparticles when approaching the nanometer scale suggest that their thermodynamic properties are mainly governed by surface effects due to a lower sound velocity [35], in opposition with dense structures as metallic or metallic oxyde nano-systems [50,51]. However, the impact of a free surface on the phase stability of thin films appears rather small because entropy difference between the HS and LS phases remains almost unchanged. Indeed, the surface vibrational entropy of both spin states behaves identically with size reduction so that $\Delta S_N^s\approx 0$. In the present work, the HS and LS states are structurally very similar so that no difference is expected to appear with the size reduction. The present results corroborate previous theoretical investigations of size effects on the switching properties of spin crossover materials through statistical physics approaches performed on mass-spring models and solved by Monte Carlo methods [52]. The coordination defects at the surface only produce a slight decrease of transition temperature with the size reduction, far from the variation experimentally obtained. On the other hand, the internal energy variation as well as the entropy differences between the HS and the LS states are likely to be strongly influenced by an external environment through an interface stress and energy as well as the presence of surface disorders. This will the object of future investigations.

4. Conclusions

In conclusion, the lattice dynamics of SCO nanomaterials have been studied through molecular dynamics simulations by adapting the so-called slab method largely employed in surface physics. A reconsideration of the force field permits to obtain a good order of magnitude of the sound velocity for SCO complexes. We show that confinement effects are limited in this class of nanomaterials, while the soft character of the lattice leads to the emergence of surface vibrational modes at low frequencies, which significantly contribute to thermodynamical properties. Therefore an increase of the vibrational entropy with the size reduction is obtained in good agreement with experimental observations. The nature and the symmetry of these Rayleigh-like surface modes have to be specified, which cannot be achieved with a projected vDOS. To this aim, the access to the phonon band structure (the dispersion relation) is essential before any further numerical investigations in this sense. The present work confirms the idea that phenomena as change in phase stability taking place at the nanoscale are governed chiefly by surface effects. The knowledge of surface/interface characteristics, as surface energy and interface stress, according to the nature of the external environment and the coherence of the interface constitutes a next challenge for the understanding of the spin transition mechanism at the nanoscale.