W Nicolazzi

S Pillet

Structural aspects of the relaxation process in spin crossover solids : phase separation, mapping of lattice strain and domain wall structure

Keywords: numbers: 05.50.+q, 64.60.De, 75.30.Wx, 75.60.-d I

We present a non-equilibrium study of the relaxation process in spin crossover solids using numerical simulations of a recently introduced two variable elastic Ising-like model. We analyse the structural lattice distortions accompanying the relaxation from the metastable high spin to the ground low spin state as a function of cooperativity. In the highly cooperative case, a sigmoidal relaxation behaviour of the high spin fraction nHS is described, and occurs jointly with a structural phase separation process. The mean lattice spacing follows a similar sigmoidal trend, owing to the interplay between electronic and lattice variables in the Hamiltonian. Weakly cooperative systems are characterized by single exponential relaxations of the high spin fraction, the corresponding structural transformation proceeds homogeneously with a progressive relaxation of the mean lattice spacing. Long relaxation tail effects are also observed. We highlight the development of lattice strain accompanying the spin transition, and show that structural phase rebuilding proceeds in the late stage of the relaxation by releasing residual strain. Under specific conditions, a temporal decoupling between the electronic and lattice variables is observed, which may have direct applications for interpreting time resolved spectroscopic or diffraction experiments, and for elucidating unusual structural behaviours, such as the development of superstructures, modulated structures or transient phases.

INTRODUCTION

Functional molecular materials have been the subject of intense research activities for the last decades and are currently considered as relevant alternative in the next generation of electronic nanodevices 1 . In this context, spin crossover (SCO) materials have attracted much attention 2 , owing to their bistability and switching properties. SCO complexes exhibit a reversible switching between molecular low spin (LS) and high spin (HS) states, which may be triggered by a change of temperature, pressure [START_REF] Meissner | [END_REF]4 , magnetic field 5,6 or using a pulsed 7,8 or continuous optical excitation 9,10 . The bistability property is closely related to strong interactions between the molecules, the so-called cooperativity, which originates from the large HS to LS molecular volume contraction inducing local lattice distortions, coupled to long-range interactions of elastic origin within the solid 11,12 ; these are mediated by intermolecular contacts between the SCO molecules 13 .

In SCO solids, a metastable HS state may be populated at very low temperature by photo-excitation from the ground LS state through the LIESST effect (Light Induced Excited Spin State Trapping) 9,10 . Detailed information on the out of equilibrium behaviour, light induced excitation and subsequent isothermal * corresponding author: sebastien.pillet@crm2.uhp-nancy.fr relaxations, have been obtained using photomagnetic and optical reflectivity measurements. It has been shown that strong cooperativity in the crystal results in characteristic sigmoidal relaxation curves of the HS fraction n HS . These results have been well interpreted in the framework of the macroscopic phenomenological equation of Hauser 14,15 , which considers in a mean field approach that relaxation is a thermally activated unimolecular process whose rate is modulated by the immediate environment. A slowing down at the end of the relaxation process, the so-called "tail effect", has been further detected in some cases 16 , and attributed to the onset of strong short-range correlations 11 or to inhomogeneities, mainly due to chemical impurities and structural defects, resulting in a spatial distribution of activation energies. 17,18 Various Ising-like models, based on a two level fictitious spin formalism with short range interactions, have been introduced to describe SCO phenomena 19,20 ; these models have been investigated by numerical simulations 11 , exact analytical treatment of the partition function 21 or analytical approaches using mean field approximation 22,23 . These models capture the essential features of spin transitions, and have been further extended to interpret dynamic aspects in the relaxation 23,24 , and the photo-excitation 25 regimes. The dynamic Ising-like model in the mean field approximation 23,25 gives a microscopic origin of the macroscopic equation of Hauser. This approach has allowed to develop analytical schemes beyond mean field approximations, by using local equilibrium method 24 or by taking into account short range correlations 26 , thus reproducing the "tail effect".

It is now well accepted that the formation and dynamics of like spin domains (LSD), i.e. the clustering of adjacent molecules in the same spin state 27 , play a major role in the various spin transition phenomena. Phase separation behaviours have been evidenced by reciprocal space mapping using x-ray and neutron diffraction experiments for the most cooperative SCO systems in the thermal spin transition 28 , photoexcitation 29 or light induced bistability 30 regimes. The corresponding kinetics of phase separation have been further interpreted using the Kolmogorov-Johnson-Mehl-Avrami model of phase transformation from kinetic x-ray diffraction measurements 29,31 . In some cases, the subtle interplay between spin and structural degrees of freedom may lead to intriguing and unusual structural features, such as symmetry breaking transitions resulting in the development of superstructures [32][33][34] or incommensurate modulated structural phases 35 , observed under specific experimental conditions of temperature, light irradiation, or under thermal quenching. More recently, LSD have been directly imaged by optical microscopy [36][37][38] and Raman spectroscopic techniques 39 , highlighting new aspects of the spatio-temporal development of the nucleation, domain growth and propagation of the thermodynamically stable phase. A temporal decoupling of spin and crystallographic phase transition may also occur, as has been reported at the thermal spin transition or during the relaxation process of SCO materials 40,41 . Quite recently, detailed information have been derived on the successive temporal steps of photoinduced spin state switching by a combination of time resolved spectroscopic and x-ray diffraction 42,43 or x-ray absorption techniques 44 . It has been shown that the electronic and structural switching span several temporal orders of magnitude from subpicosecond (electronic processes) to nanosecond (volume expansion) and microsecond (thermal switching) time scales.

As a consequence, the spin transition phenomenon is intrinsically a multi-scale process which requires an understanding of atomic scale lattice distortions, as well as mesoscopic structural organisations to interpret these unusual structural behaviours, and LSD formation and dynamics. Obviously, these effects can not be interpreted within mean field and requires appropriate models accounting for crystallographic aspects. In that direction, several microscopic elastic schemes, which introduce explicitly lattice degrees of freedom, have been recently proposed to provide a clear microscopic origin of cooperativity [45][46][47][48][49][50][51][52][53] . Some of these models are based on the atom-phonon treatment initially introduced by Nasser 50 . The first order character of the spin transition is controlled by the strength of elastic constants which can be dependent on the spin state 45,47,48,50 . Analytical solutions for one dimensional chain have shown that the atom phonon model is isomorph to an Ising like model under effective temperature dependent ligand field 45,50 . At higher dimension, lattice deformations lead to elastic long range interactions, suggesting a different nucleation and domain growth process than the original Ising-like model. Monte Carlo (MC) or molecular dynamics simulations on square 47,53 and hexagonal lattice 54 using mechano-elastic models have displayed a spin conversion which nucleates preferentially from the corner and then domains propagate inside the bulk when open boundary conditions is applied to the lattice. On the contrary, purely elastic schemes with periodic boundary conditions do not display clustering process 46 . We have introduced a microscopic elastic Ising-like model (called hereafter anharmonic model), which considers a spin and distance dependant intermolecular coupling, aiming at providing an efficient description of the structural aspects related to the spin transition 47 . Using this scheme, photo-induced crystallographic phase separation phenomena are well reproduced 55 and the corresponding domain growth kinetics follow the Kolmogorov-Johnson-Mehl-Avrami model in agreement with the experimental x-ray diffraction findings 29,31 . Nonlinear effects of the photo-excitation, such as sigmoidal photo-conversion kinetics, and the presence of an incubation time, are also retrieved.

The aim of the present study is to grasp the essential structural aspects at the microscopic and mesoscopic scale of non-equilibrium relaxation processes in SCO solids to attain a clear picture of the spatiotemporal properties using MC simulations of the anharmonic model. The paper is organized as follows. Section II is devoted to the introduction of the dynamic anharmonic model and the computational details of the MC methods. The transformation mechanisms are discussed in section III, considering the evolution of the thermal relaxation curves for the spin and lattice variables, the corresponding spatial correlation functions, as well as the development and mapping of local lattice strain, and microstructural characterisation of domain walls.

II. KINETIC ANHARMONIC ISING-LIKE MODEL AND COMPUTATIONAL DETAILS

A. Anharmonic Ising-like model

In a previous work 47 , we have introduced an anharmonic elastic Ising-like model to describe the equilibrium properties of SCO solids. The corresponding Hamiltonian, adapted from the standard two level Ising-like model 19,20 , writes

H ({σ}, { r}) = ∆ ef f 2 i σ i + i,j V elast r i,j , r 0 i,j × {J 0 + J 1 (σ i + σ j ) + J 2 σ i σ j } , (1)
The two degenerated molecular states are represented by fictitious Ising spin operators σ, whose eigenvalues σ = +1 and σ = -1 are assigned to the HS and the LS states respectively, with g + and g -the corresponding vibronic degeneracies (g + >> g -). The first term on the right hand side in eq. 1 corresponds to the onesite hamiltonian for a system of N SCO entities, where ∆ ef f = ∆ -k B T ln(g + /g -) and ∆ is the energy difference between the HS and LS ground states. To capture the essential intermolecular interactions in SCO molecular solids, the interaction energy (second term in eq. 1) is developped on pairwise 6-3 Lennard-Jones (LJ) potentials with finite range r max as given in equ. 2.

• if r i,j ≤ r max V elast r i,j , r 0 i,j = r 0 i,j r i,j 6 -2 r 0 i,j r i,j 3 • if r i,j > r max V elast r i,j , r 0 i,j = 0. (2)
r i,j = | r i -r j | and r 0 i,j = r 0 i -r 0 j are the neighboring instantaneous and equilibrium distances respectively between site i and j. They correspond to the Fe...Fe distances in SCO molecular crystals and are assimilated to lattice spacings hereafter. In the following, we set r max to 1.5. To account for the different structural (e.g. lattice spacing) and elastic (e.g. Bulk modulus, thermal expansion) properties between purely HS and LS phases, three distinct V elast potentials are considered, with equilibrium distances r 0 HS , r 0 LS and r 0 HL , for HS-HS, LS-LS and HS-LS neighbouring pairs. In this scheme, the interaction energy depends on the respective spin state and separation distance between neighbouring molecules. We consider that this is of paramount importance to account for local structural relaxations and lattice distortions accompanying the spin state switching and therefore to describe structural phase transition and phase separation phenomena. Within this elastic scheme, the spin and lattice variables are coupled directly in the Hamiltonian. The spin state switching of one molecule generates a local lattice distortion owing to the spin dependence of the V elast potential, which in turns affects the whole lattice through the intermolecular interactions. Effective longrange interactions result from this scheme, even though only short-range interactions are formally present in the Hamiltonian. This is at variance from the classical shortrange Ising-like model applied to SCO solids. The parameters of the anharmonic model can fluctuate with temperature when the system is in contact with a thermal bath. We define two observables whose thermal averages are the usual HS fraction n HS

n HS = 1 + σ 2 , (3)
directly related to the mean "magnetization" σ , and a dimensionless and normalized lattice spacing r norm

r norm = r -r 0 LS r 0 HS -r 0 LS , (4)
where r = r i,j is the mean intersite distance. As already introduced for the standard Ising-like model, first neighbour spin spatial correlations are given by

C σ,σ ′ (t) = 1 N i,j σ i σ j . (5)
The HS-LS nearest neighbour pair fraction, noted n HL , is related to the spin spatial correlation function C σ,σ ′ (t) through the relation 24 n

HL (t) = 1 -C σ,σ ′ (t) 4 (6)
In a similar way, we define a spatial correlation function written as:

C r,r ′ (t) = k,l r k r l -M × r 2 M × r 2 , (7)
where k, l corresponds to nearest neighbour inter-site bonds. M is the realizable number of product between nearest neighbour bonds with the free boundary conditions; it normalizes the correlation function. C r,r ′ probes the extent of spatial structural distortion; a zero value of C r,r ′ corresponds to a non distorted lattice with a uniform thermal mean intersite distance r .

The numerical study of the static and quasistatic thermodynamic properties 47 has revealed that the anharmonic model exhibits an order-disorder thermal transition at a critical temperature T C , different from the second order critical temperature of the Ising model under zero field, due to the existence of both short range Ising coupling and effective long range interactions, mediated by elastic couplings 56 . The anharmonic model can reproduce gradual or abrupt thermal transitions, with or without thermal hysteresis, by fine tuning the values of the model parameters, especially J 2 (Fig. 1). It is noteworthy that the J 2 = 0.5 case exhibits a gradual two step transition, resulting from the competition between short-range and long-range interactions; an ordered pattern of HS and LS species is however not formed under such conditions.

B. Kinetic anharmonic model

For investigating the relaxation process, we consider a MC method as follows. Let P ({σ}, { r}) be the probability for the system to adopt the spin and lattice configuration ({σ}, { r}) at time t. We assume that the temporal evolution of P ({σ}, { r}) is governed by the two variable microscopic master equation:

∂P ({σ}, { r}, t) ∂t = N i=1 d r ′ i {σ ′ } Ω {σ ′ }, { r ′ } → {σ}, { r} P {σ ′ }, { r ′ }, t -Ω {σ}, { r} → {σ ′ }, { r ′ } P ({σ}, { r}, t) , (8)
where Ω {σ ′ }, { r ′ } → {σ}, { r} is the transition rate.

The choice of the transition rate Ω is fundamental when non-equilibrium kinetics are considered. Its rigorous establishment should be based on the knowledge of microscopic processes, resulting from quantum mechanical considerations. It has to retrieve the Boltzmann equilibrium distribution in the stationary state ∂P ({σ}, { r}) /∂t = 0, and the detailed balance condition is imposed:

Ω {σ ′ }, { r ′ } → {σ}, { r} Ω {σ}, { r} → {σ ′ }, { r ′ } = exp (-βH ({σ}, { r})) exp -βH {σ ′ }, { r ′ } . (9)
In SCO compounds, the suitable dynamics would account for the transition probability from a HS vibrational state to a LS one for a single molecule in contact with a thermal bath and an intermolecular phonon bath, corresponding to lattice vibrations 57 . In the following, we decouple the total transition rate Ω in separate spin and lattice contributions, noted W spin ({σ ′ } → {σ}) and W elast { r ′ } → { r} . Spin and lattice degrees of freedom interact individually with the heat bath, inducing stochastically their reversal and incremental modification respectively. This particular choice is justified by the different characteristic time scales at which electronic processes (femtosecond) and nuclear displacements (picosecond to nanosecond) occur. It is well known that for SCO materials, the photoexcitation and thermal relaxation outside the tunnel regime are thermally activated processes 15 , related to the crossing of energy barriers in a purely classical point of view. Accordingly, as suggested for the standard dynamic Ising-like model, a transition rate of Arrhenius type is well suited 23,58 , for which a microscopic origin has recently been provided by expressing the "'exchange"' interactions in the Ising-like model as spin-phonon couplings 59 . There exists several forms of Arrhenius dynamics 60 . We adopt here, a "One Step Dynamic" (OSD) 61 , noted hereafter W OSD , which corresponds to the transition probability from an initial state of energy E i to a final state of energy E f , passing through an intermediate state ("saddle point") of energy E T . W OSD can be written in the general form

W OSD ∼ exp [-β (E T -E i )] . (10)
The expression for E T is simply

E T = E i + E f 2 + E barrier , (11)
where E barrier corresponds to a microscopic energy barrier. In the next, we use the same writting proposed in other previous work 23 for the spin dynamic,

W spin ({σ} → {σ ′ }) = 1 τ 0 spin e -β[Espin-E i 2 ] , (12)
and for the lattice dynamic,

W elast { r} → r ′ = 1 τ 0 elast e -β[E elast -E i 2 ] . (13)
1/τ 0 spin and 1/τ 0 elast correspond to an intrinsic intramolecular frequency associated to spin and lattice state switchings respectively. In contrast with previous results 23,58 , E spin and E elast are assimilated to non constant intramolecular and intermolecular energy barriers, which are sensitive to the local environment since they are function of the energy of the final state:

E spin ({σ ′ }, { r}) = E 0 spin + E f ({σ ′ }, { r}) 2 E elast {σ}, { r ′ } = E 0 elast + E f {σ}, { r ′ } 2 , (14)
where E 0 spin and E 0 elast are phenomenological parameters corresponding to a reference for intra-and intermolecular energy barriers.

C. Computational details

The dynamic properties of the Hamiltonian (1) are studied using MC methods. It is important to distinguish the variable t in the MC approach, which defines the time unit called Monte Carlo step (MCS), from the real time from newtonian classical equations or Schrodinger equation in quantum mechanics. The MC time is defined arbitrarily and depends on the different processes included in a MCS, giving a certain disadvantage compared to other algorithms, integrating directly the dynamical equations, as for example, molecular dynamic 48,52 . Especially, the direct comparison with experimental kinetic data is impossible. MC methods are nevertheless advantageous to calculate thermodynamic quantities. In the next, we consider a system of 32 × 32 (N = 1024) spins on a deformable two-dimensional lattice with the free boundary conditions, allowing energy and "volume" fluctuations. The model parameters are set to J 0 = 1000, J 1 = 0.15, ∆ = 7.2 and ln g = 2. Equilibrium lattice spacings of HS, LS and intermediate HS-LS structural phases are set to r 0 HS = 1.2, r 0 LS = 1 and r 0 HL = 1.1 respectively. Thermal relaxations after a quench are investigated through the n HS and r norm temporal evolutions, calculated numerically by ensemble averages; temporal averages are not valid owing to the absence of ergodicity and to the breaking of translation invariance in time. Practically, thermal means are performed on a set of 20000 identical independent configurations of systems. A MCS corresponds to the three following stages: (i) a site i and a normalized number p are chosen randomly; the fictitious spin σ i flip is updated according to the Arrhenius dynamic for the spin variable W spin . If p < W spin , the new spin configuration is accepted. ii) A second site j at the position r j (x j , y j ) and an other normalized number p ′ are chosen randomly; a new position r ′ j x ′ j , y ′ j is proposed as follows:

x ′ j = x j + d xj , y ′ j = y j + d yj , (15)
where d xj and d yj are continuous displacement drawn on a gaussian distribution with zero mean and adjustable variance. The new position is evaluated, and accepted if p ′ < W elast . iii) The two previous sequences are repeated N times. The transition probabilities should be normalized in the Arrhenius dynamic 58 to satisfy the two following inequalities:

0 < W spin < 1 0 < W elast < 1 (16)
In the simulation, all new configuration propositions inducing a transition rate which does not satisfy these two inequalities, will be systematically rejected. In all calculations, we set E 0 spin = 10 and E 0 elast = 5, and τ 0 spin = 0.8 and τ 0 lattice = 1.0.

D. Local strain mapping

Our simulations are performed on a deformable 2D lattice with equilibrium inter-site distances r 0 HS (resp. r 0 LS ) corresponding to purely HS (resp. LS) structural phases. Along the relaxation path, large structural deformations may occur induced by the HS-LS difference in inter-site distances. Strain tensors are the primary measure of such local deformations in continuum mechanics, and may be computed from the gradient of a continuous displacement field with respect to a reference lattice. Local strain tensors may be derived in a discrete form by considering the local distortion around each lattice site with respect to a reference structural configuration 62,63 . We first consider two reference regular square lattices, namely HS and LS, in which each site is surrounded by 4 neighbours with r 0 HS and r 0 LS inter-site separation distance respectively. For a given structural configuration of the lattice, a transformation matrix J i is assigned to each site i. The local affine transformation matrix J i which best describes (in the least-squares sens) the deformation of the nearest neighbour environment of site i with respect to the corresponding HS or LS reference lattice, is calculated by minimization of the quantity

j r 0 ij J i -r ij 2 , ( 17
)
where r ij is the current position vector between sites i and j and r 0 ij is the corresponding reference (HS or LS) position vector ; the summation runs over all nearest neighbours of site i. The local Lagrangian strain matrix associated to site i is then computed as

ǫ i = 1 2 J i J T i -I . (18)
In the following, to map the local strain on the simulation lattice, and therefore quantify local structural distortions, we use the invariant T r (ǫ i ) of the ǫ i matrix; it corresponds to the relative volume variation of the immediate neighbourhood of the considered molecule i.

III. RESULTS OF THERMAL RELAXATION AFTER A QUENCH

Thermal relaxation after a quench to low temperature is investigated as follows. The system is initially prepared in the totally HS electronic (n HS = 1) and HS structural configuration (r norm = 1). At t = 0, the system is quenched to low temperature (T = 1.8), outside the spinodal regime and at which the HS configuration becomes metastable, and is allowed to evolve under the stochastic dynamics. Various aspects of the relaxation to the LS ground state are probed through the non-equilibrium kinetics of the n HS and r norm variables as a function of time (in MCS); these relaxation curves are discussed in the next section.

A. nHS and r norm relaxation curves

The influence of the J 2 parameter on the relaxation curves is investigated. We have shown previously for the thermal transition regime 47 that J 2 drives the abruptness of the thermal transition, and therefore the first order character, while J 0 plays only a minor role. An increase of J 2 leads to an enhancement of the coupling between spin and lattice degrees of freedom, and therefore of the SCO cooperativity. For the standard Ising-like model, the short-range coupling parameter J drives the sigmoidal character of the relaxation of n HS , as shown by kinetic analytical [22][23][24] or MC simulations 58 from dynamic Ising-like model. Similarly for the present anharmonic model, for strong values of the J 2 parameter (1.45 ≤ J 2 ≤ 2.15), the relaxation curves of the HS fraction n HS follow a sigmoidal trend (Fig. 2), characteristic of the self-accelerated phenomenon. Owing to the interplay between spin and lattice variables in the Hamiltonian, the mean lattice spacing follows a similar sigmoidal trend from r norm = 1.0 to r norm → 0. As the cooperative relaxation proceeds, molecules in the HS electronic configuration progressively switch to the LS electronic state, giving rise to local lattice distortions. As a consequence, the mean lattice spacing gradually contracts from a purely HS structural phase to a purely LS phase. We note on Fig. 2 that the relaxation of the lattice spacing is systematically delayed from the n HS relaxation. This results from the deliberate choice of two independant MC dynamics for the spin and lattice degrees of freedom, with different characteristic times (τ 0 spin <τ 0 lattice ), and from the continuous character of the lattice variables. We have indeed analyzed the relaxation process with different ratios of the characteristic times (not shown here), and found that the corresponding relaxation curves may be significantly influenced. As J 2 decreases, the sigmoidal character of the n HS relaxation curve attenuates, corresponding to a decrease of the cooperativity of the system; this result corroborates recent simulations of thermal relaxation with an elastic model on a hexagonal lattice 49 . The lifetime of the metastable state shortens, as clearly illustrated by the short time behaviours, the escape from the HS metastable state becomes easier with the weakening of the cooperativity, corresponding to a lowering of the energy barrier. In parallel, the delay between r norm and n HS tends to increase with the decrease of the J 2 parameter, indicating that the molecule-lattice coupling weakens. The well-known tail effect, corresponding to a slowing down of the relaxation of n HS at long time, originating from the onset of short range correlations 22,26 is clearly observed and tends to disappear as J 2 decreases. The magneto elastic coupling J 2 plays therefore a major role in the short range correlations within the lattice and drives the local structural distortions. As the J 2 parameter further decreases (0.5 ≤ J 2 ≤ 1.20), the n HS relaxation curves become progressively single-exponential or stretched single-exponential (Fig. 3 (a)). On the contrary, the r norm relaxation curves do not show any noticeable evolution for J 2 below 1.20 (Fig. 3 (b)). The n HS relaxation seems to be much more influenced by the decrease of cooperativity than the mean lattice spacing. In this weakly cooperative case, it is therefore expected that the elastic properties are almost independent on the spin state, and correlatively do not vary along the relaxation process. At variance, for strongly cooperative materials, the elastic properties, such as thermal expansion tensor, Bulk modulus, sound velocity or Debye temperatures, may depend on the spin state. As a matter of fact, experimental evidences have been provided by x-ray diffraction 64 that the thermal expansion tensor in the HS and LS phases differ for the cooperative SCO material [Fe(btr) 2 (NCS) 2 ] .H 2 O, while spectroscopic ellipsometry measurements and diffuse reflectivity measurements under pressure 65 have shown that the Bulk modulus, sound velocity and Debye temperature may also change along the charge transfer transition in a Prussian Blue analog molecular solid. Within the present anharmonic model, since the V elast potential is dependent on the spin states of neighbouring molecules as defined in eq. 2, it has been found that the resulting linear thermal expansion coefficients are different in the HS and LS phases 47 .

For such weak J 2 couplings, the time scale at which the thermal relaxation of n HS and r norm occurs is completely different: a thousand MCS is sufficient for the molecules to adopt the LS equilibrium electronic configuration whereas complete relaxation of the lattice requires more than 30000 MCS steps. We anticipate that this behaviour may found interesting applications in time-resolved experiments which have shown that the electronic (spin state) and lattice degrees of freedom respond to photo-excitation with different time scales 42 . However, it has to be kept in mind that the arbitrary character of the "Monte Carlo time" (or Monte Carlo

Step) prevents a quantitative comparison with the real time of an experiment. The interplay between spin and lattice variables is well illustrated by the mutual inspection of the relaxation time τ spin and τ lattice (Fig 4). These are defined as the half-life (in MCS) of the HS phase from the spin and lattice points of view. As J 2 increases, a clear enhancement of the relaxation time occurs, correlated to an increase of the life-time of the HS metastable state. In parallel, spin and lattice variables becomes much more coupled (see Fig 4 (c)). For weak J 2 values, the relaxation time ratio τ spin /τ lattice is close to zero since the spins relax faster than the lattice variable by several orders of magnitude. On the contrary, for strong J 2 values, the ratio reaches almost 0.7, molecules and lattice relaxation takes place in very similar time scales.

B. Like spin domain formation and structural phase separation

Short range correlations and LSD formation and dynamics play a major role in the properties of cooperative spin transitions both in the quasi-static and in the out of equilibrium regimes, resulting for instance in the nonlinearity of light-induced, thermally-induced and relaxation dynamics, and in the presence of long relaxation tails; this is well documented from photomagnetic and optical reflectivity experiments. From the modeling point of view, the effective interaction range in elastic Isinglike models is of paramount interest 66 . Recent focus has been attached to the simulation of nanoparticle finitesize SCO systems 67,68 , for which the extent of long-range interactions with respect to the system (particle) size in addition to edge effects in open boundary conditions are of major concern. As discussed in the section II A, the present anharmonic model relies on a combination of short-range and effective long-range interactions, which results in the nucleation and growth of LSD in the case of highly cooperative thermal transitions, as we recently reported 47 . Experimentally, the presence of LSDs has been emphasized from x-ray and neutron diffraction experiments, evidenced by so called Bragg peak splitting, and imaged by optical microscopy [36][37][38] and Raman spectroscopic techniques 39 . However, these different types of experiment focus on very different aspects of the spin transition phenomenon. Spectroscopic techniques (optical absorption, optical reflectivity, Mossbauer) or magnetic measurements gives a signal directly related to the HS fraction n HS . On the contrary, diffraction experiments probe structurally ordered regions of a single crystal or polycrystalline sample. Spectroscopic or magnetic measurements do not even require the sample to be crystalline. The present two variable anharmonic model affords the possibility of describing the SCO phenomenon using both aspects. In the following, we propose to distinguish Molecular Like-Spin Domain (MLSD) from Structural Like-Spin Domain (SLSD). MLSD defines the clustering of molecules with the same electronic configuration (HS or LS), it has a clearly defined domain boundary. The size of MLSD can be related to the spatial spin correlation function by the Fourier transform of the two site spin correlation functions, the so-called structure function 69 . SLSD is an extended region of the system, in which the structure (intramolecular, lattice spacing and orientation) is well defined and perfectly ordered and differs from the neighbouring regions. SLSD can be characterized from scattering techniques (diffraction and diffusion) and quantitative analysis of the diffraction pattern. For instance, the domain size and domain boundary thickness has been extracted from high resolution x-ray diffraction experiments or electron microscopy for ferroelastic materials 70 .

The new terminology (MLSD and SLSD) is illustrated below in an analysis of the transformation mechanism considering the highly cooperative (J 2 = 2.05), moderately cooperative (J 2 = 1.25) and weakly cooperative cases (J 2 = 0.5).

Snapshots of the configuration of the system in the highly cooperative case are given in figure 5. In the first stage of the relaxation, two nuclei of molecules with the LS electronic configuration are formed at corners of the system, MLSDs subsequently grow until t = 29850M CS. The spin state change is followed by a contraction of the structural lattice from the HS (r norm = 1) to the LS (r norm → 0) lattice spacing values, leading to a progressive decrease of the overall volume which persists till t = 39860M CS, that is to say well after the HS to LS molecular state switching is complete. This is consistent with the delay between n HS and r norm relaxation curves discussed above.

Although the present elastic model contains both short range and effective long range interactions, nucleation phenomenon starts from the corners of the system, which is the signature of boundary conditions effects and the presence of elastic distortions 53,54,71 . Mechanisms of nuclei formation seems to be different for the elastic models by comparison with short range Ising-like schemes. It has already been shown that in the case of a purely elastic model hamiltonian on system with periodic boundary conditions, clusters are suppressed and molecular state switching occurs uniformly in the lattice 56 . The situation may be different with the anharmonic elastic model, where LSDs are expected even with periodic boundary conditions. The propagation of the domain walls is driven by thermal fluctuations. In such cases, the average dimension R of the domains may follow a power law of the form R ∝ t α . For such a highly cooperative case, SLSD superimposed almost with MLSD as can be seen at t = 18500M CS for instance. Domain coalescence occurs around t = 20000 -23000M CS. It is important to stress that a structural reconstruction occurs in the late stage of the HS to LS relaxation; structural coherence is recovered at completeness of the transition. The cohesive potential J 0 × V elast in the Hamiltonian is most probably the driving force for this late stage structural reconstruction. The temporal evolution of the spin and lattice correlation functions, defined in section (II A), brings new insights on the mechanism of the relaxation process. The evolution of the fraction of HS-LS pairs, n HL , is given in Fig. 6 for strong J 2 values. All curves exhibit a similar trend, characterized by first a rapid raise of the HS-LS pairs to a maximum, and then a slower decrease to zero at completeness of the HS to LS relaxation. The position of the maximum is displaced to shorter time and higher value when J 2 decreases. In this case of highly cooperative systems, the evolution of n HL can be interpreted in the framework of the critical droplet classical theory 72 . First, thermal fluctuations lead to the nucleation of the thermodynamically stable LS phase, leading to a rapid increase of HS-LS first neighbours. At this stage, "germ nuclei" are energetically unfavourable and do not reach the critical size above which domain growth occurs. Close to the maximum, critical nuclei are formed through favourable fluctuations 73 , LS domains develop. The maximum of HS-LS pairs fraction n HL is reached when the total length of interface separating the HS and LS MLSD becomes maximal. Then, these growing domains coalesce, n HL finally decreases. As J 2 increases in comparison with thermal fluctuations, the critical size for the "germ nuclei" is hard to attain, the maximum is displaced to longer time. In parallel, the maximum value decreases, owing to the strong ferroelastic coupling which favours HS-HS and LS-LS pair formation, leading to a higher MSLD size and lower length of the interface. These results are comparable to the numerical or analytical results provided by the standard Ising-like model 22,26 for the spin spatial correlation function. As n HL gives information on the distribution of HS and LS electronic states, and therefore on the MLSD, the inspection of the lattice spatial correlation function C r,r ′ , given in figure 7 for strong values of the J 2 parameter, is complementary. As for the HS-LS pair fraction n HL , the correlation function C r,r ′ exhibits a maximum, whose position depends on the cooperativity. A non zero value of C r,r ′ indicates that elastic distortions occur within the system. These maxima are located at longer time with respect to n HL , reflecting the delay of the lattice response with respect to the spin state change. The progressive formation of LS MSLD induces a built up of internal pressure in the system, due to the structural misfit of LS molecules within the HS structural matrix. Structural distortions set up, C r,r ′ therefore increases. After reaching a maximum, C r,r ′ decreases as the structural relaxations span the entire system. The long tail of the structural correlation function corresponds to the crystal lattice reconstruction in the late stage while the system is totally in the LS electronic configuration. Instantaneous configurations of the system in the case of weak J 2 values are represented on Fig. 8 for different instants of the thermal relaxation. We have seen in section III A that, for such a case, the spin and lattice variables are temporally decorrelated. Switching of HS molecules to the LS state occurs within the first 200 MCS, whereas the lattice has not reacted yet. Indeed, we may distinguish already formed LS clusters at t = 247M CS. A close to homogeneous and uniform compression of the intersite distances, inducing a homogeneous decrease of the whole volume, only begins around t = 2500M CS and is completed around t = 45000M CS. No MLSD or SLSD can be distinguished along the relaxation. The amplitude and position of the maximum of the HS-LS pairs fraction n HL (Fig. 9), follow a similar behaviour than those already analyzed on Fig 6, with nevertheless some important differences. First, the value of the HS-LS pairs fraction at the maximum is much higher, and the maximum position occurs at shorter time. In such a weakly cooperative case, as snapshots of the system do not evidence any MLSD, the transition mechanism is dominated by a homogeneous nucleation process of the thermodynamically stable LS phase, which leads to a rapid formation of numerous HS-LS pairs spread uniformely within the system. MLSD growth is hindered by thermal fluctu- The behaviour of the C r,r ′ lattice spatial correlation function (Fig. 10) is more complicated with respect to the case of strong J 2 values. The maximum of C r,r ′ occurs at very short time (t < 1000M CS) followed by a second maximum around t = 7000M CS for J 2 ≤ 0.95. The first maximum proceeds jointly with the nucleation process discussed for n HL , while the second maximum and subsequent long tail in C r,r ′ corresponds more likely to late structural reconstruction.

C. Mapping of local lattice strain

It is evident from figure 5, that the development of SLSD is followed by structural distortions, notably at the domain boundaries. Figures 11-14 depict a quantitative mapping of the local strain as a function of time in the highly (J 2 = 2.05), weakly (J 2 = 0.5), and moderately cooperative (J 2 = 1.25) situations respectively. High strain values are associated to molecules for which the spatial distribution of the nearest neighbours differs severely from the ideal reference positions. For instance, a LS molecule surrounded by four neighbours at a distance of nearly r 0 HS exhibits a high positive strain value. For the highly cooperative case (Fig. 11), by comparison with the configurations of the system depicted in figure 5, the development of high strain is associated to the nucleation and growth of MLSDs. The highest strain values are located at the MLSD boundaries, which correspond to molecules switching their electronic configuration from HS to LS while the lattice spacings in their immediate neighbourhood is still close to r 0 HS . As these MLSD boundaries propagate in the system as a function of time during the HS to LS relaxation, a front of large structural strain propagates meanwhile. Deep in the growing MLSD, the structure is relaxing with spacing distances close to the equilibrium LS value, the local strain reduces in parallel. Interestingly, significant strain develops also in advance of the propagating MLSD boundary (for instance at t = 18500M CS), as a precursor effect for the domain growth. The relaxation in this highly cooperative case occurs with the nucleation and growth of only two large MLSDs, which coalesce between t = 20270M CS and t = 23000M CS. The coalescence leaves high residual strain in the center of the system, which persists well after all the molecules are converted to the LS electronic state. The late stage of the relaxation corresponds to a release of this residual strain, which parallels a global size reduction of the deformable lattice to the equilibrium structural configuration of the purely LS phase. Obviously, the delay discussed above between the relaxation of the HS fraction and the normalized lattice spacing r norm may be attributed to the release of the accumulated residual stain. It is noteworthy that the release of residual strain starts from the corners of the system, resulting most probably from the open boundary conditions used in our simulations. The accumulation of high strains at the domain walls characterized here has been recently observed by microscopy images using crossed polarizer, detecting birefringence of regions deformed by elastic strain and stress, 39 and has been interpreted as an acoustic wave which propagates in an elastic medium. A very different situation occurs for the weakly cooperative case (Fig. 12). Here, molecules are switching their spin state at almost spatial random in the system (see Fig. 8), no MLSD may be detected. For each LS site, a high strain develops purely locally. Contrary to the highly cooperative case, strain is limited to the LS molecule and does not propagates easily to the neighbouring molecules. This is obviously due to the weak inter-molecular interactions (weak J 2 ), a molecule switching its spin state does not perturb the structural configuration of its immediate neighbourhood. At t = 1200M CS, all the molecules are in the LS electron configuration (see Fig. 8), while the entire system exhibits large strain values; this residual strain is relaxed later on, starting from the corners of the system.

It is interesting to compare the structural characteristics of the relaxation in a moderately cooperative case (J 2 = 1.25) with respect to the two situations discussed above. This moderately cooperative case exhibits an abrupt thermal spin transition, without any noticeable hysteresis (Fig. 1). However, the development of MLSD is fundamentally different from the highly cooperative situation. In the latter case, only few domains nucleate at the boundary of the system and grow, while for J 2 = 1.25, several MLSD nucleate within the system and develop (Fig. 13). Around t = 3400M CS, the system is completely in the LS electronic configuration, but the corresponding lattice spacings are still much closer to the HS value. Starting from t = 3400M CS, the size of the system progressively decreases untill all the inter-site distances are quite close to the equilibrium LS value. Large strain values are associated to the forming MLSDs, strain relaxation occurs after the HS to LS electronic configuration relaxation is completed around t = 3400M CS. Slight strain increase is observed ahead of the propagating MLSD boundaries, as for the highly cooperative case.

D. Structure of the domain wall

The structural distortions associated to domain walls in ferroelastic materials, such as wall thickness and strain distribution, are essential parameters which govern the walls energy and the microscopic domain pattern as well as the interactions of the domain walls with crystal defects. The spontaneous strain at domain walls may be estimated from the Landau free energies; the corresponding profile of the order parameter across the wall is found countinuous of the form tanh(x/W ) with W the wall thickness [START_REF] Salje | Phase Transitions in Ferroelastic and Coelastic Crystals[END_REF][START_REF] Novak | [END_REF][76] , which is indeed observed experimentally for twin walls 77 . In the context of SCO materials, the structure of domain walls and their interaction with defects, such as dislocations and impurities, may have important implication on the dynamics of domain nucleation and growth. The impurity effect can be induced by replacement (dilution) of the Fe(II) SCO active molecules with isostructural Zn(II) or Co(II) analogues which do not present SCO. It has been shown by time and temperature dependent crystallographic analysis that the activation energy to domain growth is indeed perturbed with dilution 31 . We have defined above two concepts of domains, with respect to the two variables of our model : MLSD and SLSD. It is relevant to analyse under details the profiles of the spin and lattice variables in the neighbourhood of the domains boundaries. For that purpose, we give in figure 15 the profiles of the spin variable σ i , lattice variable r i,j , and strain value T r (ǫ i ) along the vertical arrow depicted in figure 5 for the highly cooperative case. The MLSD boundary is associated to an abrupt increase of the σ i value from -1 inside the growing LS domain to +1 in the direction of the HS matrix (Fig. 15 (a)). In parallel, the inter-site distance almost abruptly increase from nearly 1.05 to 1.35 with only one intermediate value in the SLSD boundary; the transition region of the SLSD boundary is thus found quite sharp (Fig. 15 (b)). Interestingly the location of the MLSD and SLSD boundaries match perfectly with each other. Ahead of the SLSD boundary in the not yet converted HS region, the inter-site distances are larger than the equilibrium r 0 HS = 1.2 value and correspond to the build up of strain ahead of the propagating SLSD boundary. The decrease of inter-site distances in the forming LS domain is compensated by this increase of the inter-site distances in the HS phase. This is energetically allowed by the anharmonic shape of the inter-site Lennard-Jones potential V elast . The strain profile (Fig. 15 (c)) exhibits a large transition region, with a maximum at the boundary of the domain, a progressive decrease of strain from the boundary to almost 0.0 in the still unconverted HS phase, and a progressive decrease to a significant residual strain value in the forming LS domain.

IV. CONCLUSION

We have studied the high spin to low spin relaxation phenomenon in spin crossover molecular solids using out of equilibrium Monte Carlo simulations of a recently introduced anharmonic Ising-like Hamiltonian. The model is based on coupled electronic (spin state) and structural degrees of freedom of interacting spin crossover entities on a deformable lattice, accounting for the structural lattice distortions. The concept of like spin domain emerges directly from this two variable scheme. A distinction between electronically defined (MLSD) and structurally defined (SLSD) domains is proposed. By fine tuning the elastic coupling J 2 of the model, SCO materials from high cooperativity to weak cooperativity can be described. In the former case, the sigmoidal relaxation kinetics of the high spin fraction n HS , is retrieved, while the mean lattice spacing r norm follows a similar, albeit delayed, trend. As the cooperativity progressively weakens (decreasing J 2 ), the lifetime of the metastable HS state decreases, the sigmoidal kinetics is maintained but the rate of relaxation increases. A further decrease of J 2 induces a loss of the sigmoidal caracter, while a temporal decoupling between the electronic and structural variables occurs. This decoupling has already been observed experimentally, and may found important applications to interpret time resolved spectroscopic or diffraction experiments. In the highly cooperative case, the spin transition proceeds through nucleation and growth of few domains of molecules with the LS electronic configuration, accompanied by lattice distortions. A front of high lattice strain is associated to the MLSD and propagates in parallel to the domain boundary during the HS to LS relaxation. Lattice strain is also observed ahead of the domain boundary. No MLSD nor SLSD are observed in the weakly cooperative case, while the development of lattice strain is essentially local and does not propagate within the system. In all cases, residual strain is progressively built up in the system during the relaxation, and released in the late stage of the relaxation through a structural phase rebuilding process driven by the cohesive potential J 0 × V elast in the Hamiltonian.

Our work provides the basis for interpreting the dynamics of phase transformations in spin transition materials and investigate the nucleation and growth of structural like spin domain. It gives the first relation between mesoscopic and microscopic processes through the local strain mapping. A quantitative comparison between simulated diffraction pattern and experimental diffraction pattern may be very informative with that respect. This work is under progress and will be published in a forthcoming paper.

FIG. 1 :

 1 FIG.1:(color online) HS fraction nHS (open symb.) and normalized lattice spacing r norm (filled symb.) as a function of temperature, for different values of the J2 coupling parameter. Thermal cycles have been computed as in ref.47 , using the Metropolis MC dynamic for the two variables, thermal averages have been computed over 5000 independant MC simulations. Model parameters: J0 = 1000, J1 = 0.15, ∆ = 7.2 and ln g = 2.

FIG. 2 :

 2 FIG.2:(color online) Relaxation curves for the HS fraction nHS (full lines) and the normalized lattice spacing r norm (dashed lines) for different "strong" values of the J2 parameter (1.45 ≤ J2 ≤ 2.15).

50 FIG. 3 :

 503 FIG. 3: (color online) Relaxation curves for (a) the HS fraction nHS and (b) the normalized lattice spacing r norm for different "weak" values of the J2 parameter.

FIG. 4 :

 4 FIG.4:(color online) Evolution of (a) the lattice relaxation time τ lattice , (b) the spin relaxation time τspin and (c) the τspin/τ lattice ratio as a function of the J2 parameter. Relaxation times are defined as the half-life (in MCS) of the metastable HS phase (nHS = 0.5 and r norm = 0.5 for spin and lattice resp.).

FIG. 5 :

 5 FIG. 5: (color online) Instantaneous configuration of the system (spin and lattice) during thermal relaxation in the highly cooperative case (strong elastic coupling (J2 = 2.05)). LS and HS molecules are denoted as blue and red circles respectively. The vertical arrow at t=18500 MCS indicates the position of the profiles depicted in figure 15. The outmost square gives the size of the system at t = 0.

FIG. 6 :

 6 FIG.6:(color online) Temporal evolution of the fraction of HS-LS pairs nHL for "strong" values of the J2 parameter.

FIG. 7 :

 7 FIG.7:(color online) Temporal evolution of the lattice spatial correlation function C r,r ′ for "strong" values of the J2 parameter.

FIG. 8 : 50 FIG. 9 :

 8509 FIG.8:(color online) Evolution of the lattice and spin system configuration at different times during thermal relaxation in the case of a weakly cooperative system (J2 = 0.5). LS and HS molecules are denoted as blue and red circles respectively. The outmost square gives the size of the system at t = 0.

50 FIG. 10 :

 5010 FIG.10:(color online) Temporal evolution of the lattice spatial correlation functions C r,r ′ for different "weak" values of the J2 parameter. The insert corresponds to a zoom of C r,r ′ at short time (0 < t < 1000).

FIG. 11 :

 11 FIG. 11: (color online) Spatial mapping of the strain value (T r (ǫi)) during thermal relaxation in the highly cooperative case (J2 = 2.05).

FIG. 12 :

 12 FIG.12:(color online) Spatial mapping of the strain value (T r (ǫi)) during thermal relaxation in the weakly cooperative case (J2 = 0.5).

FIG. 13 :FIG. 14 :

 1314 FIG.13:(color online) Evolution of the lattice and spin system configuration at different times during thermal relaxation in the case of a moderately cooperative system (J2 = 1.25). LS and HS molecules are denoted as blue and red circles respectively. The outmost square gives the size of the system at t = 0.

FIG. 15 :

 15 FIG. 15: Profile of the (a) spin variable σi, (b) inter-site distance r i,j and (c) strain value T r (ǫi) along the vertical arrow depicted in figure 5 at t=18500 MCS.

Acknowledgments

This work was supported by the European Network of Excellence MAGMANet (FP6-515767-2), the Université Henri Poincaré and the CNRS. The authors would like to thank K. Boukheddaden for helpful discussion. Part of the numerical calculations was performed at the computing center of the Institut Jean Barriol, Nancy, which is acknowledged.