Properties of the low-spin high-spin interface during the relaxation of spin-crossover materials, investigated through an electro-elastic model A. Slimani, K. Boukheddaden, F. Varret, M. Nishino, and S. Miyashita Citation: The Journal of Chemical Physics 139, 194706 (2013); doi: 10.1063/1.4829462 View online: http://dx.doi.org/10.1063/1.4829462 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/139/19?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Analysis of long-range interaction effects on phase transitions in two-step spin-crossover chains by using Isingtype systems and Monte Carlo entropic sampling technique J. Appl. Phys. 112, 074906 (2012); 10.1063/1.4756994 Study of the relaxation in diluted spin crossover molecular magnets in the framework of the mechano-elastic model J. Appl. Phys. 109, 07B111 (2011); 10.1063/1.3556702 Spin relaxation due to electron–electron magnetic interaction in high Lande g -factor semiconductors J. Appl. Phys. 108, 054505 (2010); 10.1063/1.3481063 Computer simulations: A window on the static and dynamic properties of simple spin models Am. J. Phys. 76, 445 (2008); 10.1119/1.2839563 Monte Carlo modeling of spin injection through a Schottky barrier and spin transport in a semiconductor quantum well J. Appl. Phys. 96, 4319 (2004); 10.1063/1.1794893

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.255.6.125 On: Mon, 10 Aug 2015 23:18:33

THE JOURNAL OF CHEMICAL PHYSICS 139, 194706 (2013)

Properties of the low-spin high-spin interface during the relaxation of spin-crossover materials, investigated through an electro-elastic model A. Slimani,1 K. Boukheddaden,1,a) F. Varret,1 M. Nishino,2,3 and S. Miyashita3,4 1

Groupe d’Etudes de la Matière Condensée, CNRS-Université de Versailles, 45 Avenue des Etats Unis, F-78035 Versailles Cedex, France 2 Computational Materials Science Center, National Institute for Materials Science, Tsukuba, Ibaraki 305-0047, Japan 3 CREST, JST, 4-1-8 Honcho Kawaguchi, Saitama 332-0012, Japan 4 Department of Physics, Graduate School of Science, The University of Tokyo, Bunkyo-Ku, Tokyo, Japan

(Received 11 August 2013; accepted 25 October 2013; published online 20 November 2013) The present work is devoted to the spatio-temporal investigations of spin-crossover lattices during their thermal relaxation from high- to low-spin state. The analysis is performed using Monte Carlo simulations on a distortable 2D lattice the sites of which are occupied by high-spin (HS) or low-spin (LS) atoms. The lattice is circular in shape and the HS to LS transformation results in single domain nucleation followed by growth and propagation processes. The evolution of the LS:HS interface is monitored during the relaxation process, through the mapping of spin states, displacement fields, local stresses, and elastic energy. The results show a curved interface, the curvature of which is reversed at the mid-transformation. The local stresses and elastic energy peak at the vicinity of the HS:LS interface, with sizeable dependence upon the position along the front line which evidences the edge effects. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4829462] I. INTRODUCTION

The high-spin (HS) ↔ low-spin (LS) transition is still the object of various experimental and theoretical investigations. The interplay between the spin crossover and the structural properties, originating from intra-molecular vibronic coupling, can be enhanced by elastic inter-molecular interactions. The contributions of long-range elastic interactions, particularly at the solid state, lead to rather abrupt thermal spin transition and in many cases to hysteresis behavior denoting a first-order phase transition or two-step spin transition.2–4 However, a gradual transition, corresponding to the simple Boltzmann distribution between two states, is generally obtained in highly diluted crystals (i.e., non-cooperative systems). The molecular spin state transition from LS to HS is accompanied by an expansion of the unit cell of the crystal owing to a stretching of the bond lengths (∼10% for FeII ion). Extensive research efforts have been devoted to the observation and/or the visualization of the spatiotemporal aspects on spin crossover single crystal by optical microscopy.5–12 These studies have allowed to monitor a HS:LS transformation front upon the spin crossover transition and to follow the spatiotemporal features of the HS:LS interface. Such results have evidenced the important role of the mechanical stresses on the front shape and on the nucleation and growth process as well as on the interface propagation. In spin crossover solids, the lattice strain, originating from the molecular volume expansion or decompression overall the crystal, gives rise to long range interactions; a fact which was first pointed out by Ohnishi and Sugano.13 Therefore, a volume change of a few molecules, uniformly distributed in the crystal, leads to 1

a) Electronic mail: [email protected]

0021-9606/2013/139(19)/194706/8/$30.00

a pressure which reflects at the surface of the crystal (due to the finite volume), called the image pressure.14, 15 This pressure effectively acts on all molecules of the crystal (independent of distances) with the same strength. The contribution of the elastic interactions has been studied in detail both theoretically16–26 and experimentally by measuring the elastic constants by Brillouin spectroscopy27 and the strain tensors through temperature dependent x-rays diffraction.28 From the general point of view, the present problem is similar to that of statistical mechanics of miscibility and superstructure formation of solids binary alloys which was studied in the past in the frame of lattice-gas models. In fact, these models are known to be often able to fit quantitatively experimental phase diagrams. However, the model limitations are quite obvious: each lattice site interacts with a fixed neighbor through a fixed interaction constant and the interplay between the elastic distortions and the spin composition cannot be properly modeled. Recent works18–26 pointed out the importance of the spin-lattice interaction which can even change the nature of the phase transition. Here, we present an approach aiming to investigate the physics of the HS:LS interface in spin crossover solids. We consider a circular 2D lattice with fixed topology, where spins and positions of each node are changed through the usual Monte Carlo29 procedure and the positions of nodes are assumed to change only inside the plane. To take account of the lattice contribution to the entropy as well as the change in the rigidity of the lattice between the HS and the LS states, we consider an anharmonic inter-atomic potential26 between the spin-crossover nodes. The paper is organized as follows. In Sec. II, we introduce the electro-elastic Hamiltonian of the model and we describe the simulation procedure and the choice of the interaction parameter values. In Sec. III, we present the results of

139, 194706-1

© 2013 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.255.6.125 On: Mon, 10 Aug 2015 23:18:33

194706-2

Slimani et al.

J. Chem. Phys. 139, 194706 (2013)

the model for the circularly shaped lattice during the low temperature relaxation of the metastable HS state, with a detailed mapping of the stress fields and elastic energy distributions. In Sec. IV, we conclude the paper.

R0 (+1, −1) = R0 (−1, +1) = R0H L , and R0 (−1, −1) = R0LL , where R0H H , R0H L , and R0LL are the respective equilibrium distances between HS-HS, HS-LS, and LS-LS sites. It is straightforward to demonstrate that R0 (Si , Sj ) = ρ0 + ρ1 (Si + Sj ) + ρ2 Si Sj ,

II. THE MODEL: CALCULATION DETAILS

R0H H +2R0H L +R0LL

In the present work, we investigate the elastic properties of a discrete two-dimensional circular lattice with open boundary conditions. The choice of the circular lattice as in Ref. 24 is motivated here by the fact that for large systems all nodes located at the surface are equivalent. Such choice avoids preferential nucleation regions like corners in a square lattice. The considered circular system contains N  π R2 (R = 30) nodes placed respecting a square symmetry, where each node has four nearest neighbors and four next-nearest neighbors. The nodes of lattice are considered as spincrossover molecules and may have two states HS and LS described by the respective eigenvalues +1 and −1 of an associated fictitious spin S. The distance and the interactions between the molecules depend on their spin states. In fact, the interactions between the nearest and next-nearest neighbors may be considered through springs whose stiffness depend on the instantaneous distance. It is important to notice at this stage that real materials have 3D structures in which the atoms can move through the three directions of space. So considering a 2D model allows to describe a class of spin-crossover material having a strong two-dimensional character, such as [Fe(bbtr)3 ] (BF4 )2 30 (bbtr=1,4-di(1,2,3-triazol-1-yl)butane), here abbreviated as Fe(bbtr). This system can be described from the crystallographic point of view, as a stacking of spin crossover layers along one crystallographic axis with covalent interaction in the plane and weak interactions between the layers through hydrogen and van der Waals bonding. On the other hand, the experimental visualization5–7 of the spin transition along two crystallographic directions (parallel and perpendicular to the stacking axis) have shown that the LS domain spread more easily in the plane of the covalent interactions than in the perpendicular direction where many defects and cracks appear due to the weak interaction. The total Hamiltonian of the system accounting for electronic and elastic contributions is written as  ( − kB T ln g)Si  + Aij [rij − R0 (Si , Sj )]2 H = 2 i i,j +



Bik [rik − R0 (Si , Sk )]2 .

(1)

i,k

This Hamiltonian was introduced in a previous work in the case of a square lattice.26 We recall that the first term of Hamiltonian (1) contains the ligand-field energy, , and the entropic contribution, −kB T ln g, arising from the degeneracy g (electronic, orbital, and vibrational) of the spin states. The second and third terms account for the elastic interactions between nearest-(nn) and next-nearest neighbors-(nnn) spin-crossover units, respectively. R0 (Si , Sj ) is the equilibrium bond lengths between two nodes i and j depending on the bond type; R0 (+1, +1) = R0H H ,

(2)

R0H H −R0LL

), where ρ0 = ( ρ1 = ( 4 ) , and 4 HH HL LL R0 −2R0 +R0 ) ( ρ2 = . In the Hamiltonian (1), Aij (Si , 4 Sj ) (resp. Bik (Si , Sk )) denotes the local bond stiffness of nn (resp. nnn) bonds and written under the following form so as to decrease the total elastic constant in (rij ) = A0 + A1 (rij − R0H H )2 and the HS spin state; Aij√ Bik (rik ) = B0 + B1 (rik − 2R0H H )2 . A0 (resp. B0 ) and A1 (resp. B1 ) are respectively the harmonic and the anharmonic contributions to the elastic interaction energy between nn (resp. nnn) neighbors, while rij = ri − rj  (resp. rik ) is the instantaneous distance between two nn (resp. nnn) sites. In the present model, we use as far as possible realistic parameter values derived from experimental data. Thus, we take for ligand-field energy  = 900 K and for the degeneracy ratio g = 15031–33 leading to an entropy change at the transition S = NkB ln g = 82 J.K−1 Mol−1 and a transition temperature Teq = kBln g ≈ 90 K. The equilibrium distances in the LS and HS states are, respectively: R0 (−1, −1) ≈ 1 nm and R0 (+1, +1) ≈ 1.2 nm while for the 0 (−1,−1)) . The HL state we consider R0 (+1, −1) ≈ (R0 (+1,+1)+R 2 lattice parameter between nnn HS (LS) atoms are simply √ √ given by 2R0 (+1, +1) ≈ 1.7 nm ( 2R0 (−1, −1) ≈ 1.4 nm) and the values of the elastic parameters A0 (resp. B0 ) are A0 = 38 000 K.nm−2 ≈ 38 meV Å−2 , A1 = 380 000 Knm−4 ≈ 380 meVÅ−4 , and B0 = 0.3A0 (B1 = 0.3A1 ), which leads to an elastic bulk modulus of ∼10–20 GPa. By the similar way, the equilibrium distances and the elastic parameters corresponding to the mixture HS-LS are introduced as the arithmetic mean value of the pure LS and HS states, which gives A0 ( + 1, −1) = 45.6 meV.Å−2 and B0 = 0.3A( + 1, −1). The Monte Carlo procedure was executed on both spin and position variables. The stochastic algorithm is performed in the following way: for a site i randomly selected, with spin (Si = ±1) and position ri , a new spin value Si ( = −Si ) will be set without position change. This spin change is accepted or rejected by the usual Metropolis criterion.34 Once the new spin value is accepted then the lattice is relaxed mechanically by a slight motion of each nodes selected randomly with a quantity δ ri  ri . Afterwards, a new site will be selected randomly and so on . . . Once all the nodes of the lattice are visited for the spin change, we define such step as the unit of Monte Carlo step and denoted “MCS.” It should be remarked that the relaxation process treat the lattice distortions ri as fast variables and the spin ones as slow variables. III. RESULTS AND ANALYSIS

In this study, we examine the relaxation of metastable HS states at low temperature. Such metastable state may be obtained by photo-excitation (LIESST effect) or by rapid quenching of the high-temperature stable HS state. The aim of this work is to understand how the spin lattice couples to

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.255.6.125 On: Mon, 10 Aug 2015 23:18:33

Slimani et al.

194706-3

J. Chem. Phys. 139, 194706 (2013)

the elastic field and how this microscopic coupling results at the macroscopic scale. A. Spin and lattice relaxation

The HS fraction, nHS , that is the fraction of molecules in the HS state is given as a function of the average fictitious magnetization, S as follows: (1 + S ) . (3) 2 Initially, the lattice is prepared in the HS state from the electronic and elastic point of view at temperature T = 1 K. That means all the spin values are equal to +1 and all the distances between the sites are equal to R0 (+1, +1) = R0H H . The latter configuration is metastable at low temperature and consequently the system will be constrained to relax towards the LS state. Like optical microscopy observation, through this simulation we aim to reproduce a single domain behavior which is strongly dependent of the interaction between molecules. Therefore, different values of the harmonic elastic constant A0 were considered. We summarize in Fig. 1 the time dependence of nHS for different elastic constant, A0 . It is observed an increase of the lifetime of the metastable state as a function of A0 , which is straightforwardly correlated to the increase of the elastic energy barrier. Additionally, we compared the mid-transformation lattice configurations (nHS = 1/2) and found a strong impact of the elastic interaction value on the relaxation behavior as reported in the inset to Fig. 1. Indeed, the larger value A0 = 38 000 K nm−2 generates a single droplet nucleation mechanism, while weaker elastic constant values lead to a multi-droplet nucleation mechanism or ramified structures usually related to gradual thermal transitions. Thus, in the following, we will focus only in the single droplet behavior corresponding to A0 = 38 000 K nm−2 . To access to the elastic properties involved in the renHS =

(1)

LS

A 1,0

A0=30000 K/nm

B C

0,8

HS

2

(2) A0=35000 K/nm

(3)

nHS

0,6

(1)

(2)

2

D

laxation process, the lattice configuration was monitored and summarized in Fig. 2(a), where blue (red) atoms represent LS (HS) states. The nucleation of LS domains takes place at the boundary of the lattice and spreads over the whole system which agrees well with optical microscopy data.6, 7 Such behavior evidences the effect of the elastic interactions which generate long-range strains and hinder the nucleation of the LS state from other parts of the lattice. Further more, it should be noted that the interface shape drastically changes during the relaxation process; initially it is circularly shaped and it changes gradually through the propagation to become a straight line at the middle region of the lattice, at which the interface length goes through a maximum. Beyond the midregion of the lattice, the interface length decreases and goes back to the radial shape, but with a negative curvature. These results support the macroscopic crystal shape effect on the front propagation dynamics, as already reported in Ref. 24. Let us now introduce some energetic considerations related to the present relaxation process. Figure 2(b) shows the evolution of the elastic energy distribution in the system during the relaxation process. The energy density peaks in the region of the LS:HS interface and vanishes far away. By comparing Figures 2(a) and 2(b), we notice a similarity between the elastic energy distribution and the shape of the transformation front. The distribution of the elastic energy could be explained on the basis of long range interaction. If a LS nucleus is formed inside a HS metastable state, there will be an extra energy cost due to the strain mismatch between the HS and LS phases. The stress cannot be released only at the boundary of the two electronic states (electronic HS-LS interface), owing to the long-range nature of the strain field,26 the elastic energy cost will be delocalized among several unit cells surrounding the interface region, i.e., in the volume of the transformed phase. The presence of the long-range strain field thus induces a cooperative nucleation and growth process. It is important to remark that non-cooperative systems show an exponentially shaped relaxation curve from the metastable HS state to the stable LS state. The absence of cooperative interactions hinders the domain formation and the HS to LS transformation process occurs everywhere in the lattice, giving rise to a multi-droplet phenomenon, in which macroscopic LS or HS structures do not appear, as depicted in Fig. 2(c).

(3)

0,4

A0=38000 K/nm

E

0,2

2

B. Geometric aspects of the HS-LS interface F

0,0 0

200

400

600

800

1000

time MC FIG. 1. Relaxation curves from HS → LS as function of Monte Carlo step at T = 1 K for different values of harmonic elastic constant A0 = 30 000, 35 000, and 38 000 K nm−2 and represented by the following colors red, green, and black, respectively. Here we set B0 = 0.3A0 and B1 = 0.3A1 . Inset: corresponding spatial distribution of the spin state for the average value of the HS fraction, nHS = 12 . Blue (resp. red) dots represent LS (resp. HS) atoms. The single droplet nucleation is obtained in the case of a strong elastic interactions.

To investigate the effect of the lattice geometry on the interface shape, we consider two analytical interfaces (planar and radial) for a comparative purpose with those derived from above MC simulations. The planar interface is represented in Fig. 3(a) (left) by a portion of a disk intercepted by the angle θ , while, the radial interface results from the intersection of two circles having different radii and then spin domain resembles to an asymmetric lens (see Fig. 3(a) (right)). For the radial interface, a contact angle of π2 is imposed at the intersection points of the two circles, as reported in Ref. 24. This assumption means that the projection of the gradient of

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.255.6.125 On: Mon, 10 Aug 2015 23:18:33

194706-4

(a)

Slimani et al.

A

J. Chem. Phys. 139, 194706 (2013)

B

C

(b)

B

A

C

1200

1000

800

600

400

200

D

E

F

E

D

F

0

1200

(c)

(d) 1000 800 600 400 200 0

FIG. 2. Lattice configuration (a) and distribution of the elastic energy (b) upon the relaxation process through Monte Carlo simulation time corresponding to the elastic constant A0 = 38 000 K nm−2 . [A → F] correspond respectively to the following instants :1, 200, 400, 500, 600, and 850 MCS. (c) and (d) exhibit the lattice configuration and the distribution of the elastic energy respectively upon the relaxation process through Monte Carlo simulation time for a non-cooperative system with the elastic constant A0 = 20 000 K nm−2 .

the HS fraction on the lattice boundary vanishes, which corresponds to the well known Neumann boundary conditions.35 Therefore, the corresponding interface radius, r, must obey the relation r = R tan θ2 , where R is the lattice radius. MC simulations have been performed on each case (planar and radial) separately as follows: (i) a LS domain with the planar or radial shape is formed at the border of the HS lattice and (ii) the system is relaxed mechanically with locking the spin variables of the nodes. To examine the energy distribution and analyze the shape variation of the interface, we monitor the

(a)

evolution of the relaxed elastic energy, Eelas , of the fully relaxed lattice as function of nHS for different LS domain sizes. The results are reported in Fig. 3(b). The elastic energies of the analytic planar and radial interfaces point out the existence of two remarkable points denoted I1 and I2 at which the two interfaces have the same elastic energy. It appears that for small sized LS domains the radial shape of the interface (red line) is more stable compared to the planar one (blue line) till reaching the point I1 (corresponding to nHS ∼ 0.4), where both interfaces become similar, since the radial interface straights

(b)

MC interface Analytic radial interface Analytic planar interface

I1

Elastic energy (K)

100

I2

80 60 40 20 0 0.0

0.2

0.4

0.6

0.8

1.0

nHS FIG. 3. (a) Schematic views of (left) the planar interface defined as a chord making an angle θ and (right) the radial interface defined as the intersection of two perpendicular circles (contact angle = π2 ), leading to a line (interface) known as the radical line. See text for more explanations. (b) Evolution of the lattice elastic energy as function of nHS . The red (blue) curve corresponds to the case of a radial (planar) analytic interface. While, the black curve is assigned to the numeric elastic energy derived from the Monte Carlo simulation of Figure 2(a).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.255.6.125 On: Mon, 10 Aug 2015 23:18:33

Slimani et al.

194706-5

J. Chem. Phys. 139, 194706 (2013)

in this region. Beyond the point I2 (corresponding to nHS ∼ 0.6), a re-entrant phenomenon occurs since the radial interface becomes again the most stable one. As we can notice, such behavior well reproduces the elastic energy derived from the Monte Carlo simulations, reported by the black curve in Fig. 2(a). It should be noted that the deviation between the numeric and analytic curve shown in the area [0.5,1] of Fig. 3(b) is assigned to the contribution of additional small LS domains located at the border of the lattice upon the Monte Carlo process, and such effect is omitted in the analytical calculations. Let us now discuss the geometric aspects of the planar and radial interfaces of Fig. 3(b) and try to predict some general features of the transformation front shape. If the volume change of the lattice is neglected, we can easily calculate analytically the LS domain area and the interface length for planar and radial analytic interfaces. For a planar interface, intercepting an angle θ (see Fig. 3(a)), the LS domain area, AP , and its interface length, LP , are given by R2 θ (θ − sin θ ) and LP = 2R sin . (4) 2 2 While, for the radial interface the domain area, AC , and interface length, LC , are given by AP =

AC =

R2 [θ (1 − t 2 ) − (1 + t 2 ) sin θ + π t 2 ], and 2 LC = R (π − θ ) t,

(5)

in which θ t = tan . (6) 2 At fixed LS domain area or HS fraction, nHS , the most stable interface shape is the one which has the shortest length, as depicted in Fig. 4(a) where we have summarized the HS fraction dependence of the interface length for the analytical radial (red points) and planar (blue points) fronts as well as for those (black points) of MC simulations of Fig. 1. The results clearly show that the planar interfaces do not reproduce the MC data. In contrast, the HS fraction dependence of the radial

MC interface Analytic radial interface Analytic planar interface

60 50

(b)

Interface Length (nm)

Interface Length (nm)

(a)

interface matches very well the MC simulations. Such behavior confirms that the assumed right contact angle ( π2 ) between the interface and the lattice border constitutes a very good assumption. We have also checked that the elastic energy of the analytical radial interface and that of MC simulations are proportional to the interface’s length, as shown in Fig. 4(b). We should mention at this stage that the elastic energy of the MC data of Fig. 4(b) has been corrected from elastic energy contribution of the LS atoms located at the surface of the lattice. Indeed, during the relaxation process of the HS fraction (see Fig. 2), the HS surface atoms switch very early to the LS state, due to their smaller elastic energy barrier, compared to that of bulky atoms. Thus, the HS fraction and the total elastic energy must be corrected from this residual contribution, which is expected to be negligible for large size lattices. One can also see that the two set of data of Fig. 4(b) are almost superimposed, which proves the consistency of our geometric analytical description of the interface properties. The linear relation between the elastic energy and the interface length means that most of the elastic energy induced by the lattice parameter misfit is stored in a thick stripe around the interface region. Consequently, at fixed HS fraction, the most stable interface is the shortest one which respects the boundary conditions, which are reflected here by the fixed contact angle value π2 . So, finally the present problem is reduced to an interface optimization under the constraints imposed by the presence of a free surface (boundary). The existence of different stable interface shapes is a signature of a non-monotonous propagation velocity. It is interesting to recall that, usually experiments7, 8 report that the propagation of the HS:LS interface takes place at constant velocity. However, the latter result is appropriate to plane waves or to circular waves that have grown sufficiently large to have a low curvature. For the initial development of radial interfaces, there is generally a high curvature. In this case, if the elastic interaction is too strong, the formation of an interface may fail as the perturbation will not reach the threshold value required for the electronic energy to overcome the elastic energy barrier. The relationship between the velocity of a planar interface wave, c0 , and that

40 30 20 10 0

MC interface Analytic radial interface

60 50 40 30 20 10 0

0.0

0.2

0.4

0.6

nHS

0.8

1.0

0

20

40

60

80

100

Elastic energy (K)

FIG. 4. (a) Interfaces lengths versus the HS fraction nHS , for the planar (blue), radial (red), and simulated (black) cases. An excellent agreement is found between the MC results and the analytical prediction of the radial interfaces. (b) Dependence of the interface length versus the elastic energy for the analytic case (red) and the MC simulated case (black). An almost linear behavior is obtained denoting that the length interface should be proportional to the elastic energy.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.255.6.125 On: Mon, 10 Aug 2015 23:18:33

194706-6

Slimani et al. Stochastic regime

(b)

Deterministic regime

Accelerated regime

50

0.5

40

0.4 0.3

30

0.2

20

0.1 10 0.0 0 100

200

300

400

500

600

Instantaneous velocity

Front position

(a)

J. Chem. Phys. 139, 194706 (2013)

700

MCS FIG. 5. (a) Instantaneous elastic profiles corresponding to the distance between two successive sites along the propagation direction of the transformation front. From left to right profiles correspond respectively to the following Monte Carlo instants 300 (filled squares), 400 (filled circles), 500 (filled triangles), and 600 (filled stars) MCS. In all the profiles, the LS (HS) region is represented by blue (red) points, allowing to identify the electronic interface. We can notice the presence of compressive (in the LS phase) and tensile (in the HS phase) stresses on both sides of the interface. These profiles provides also information on the width of the elastic interface estimated here as ∼4 nm. (b) (black line) Successive positions of the interface during the propagation process. (blue line) Instantaneous velocity during propagation process which is the derivative of the position’s curve. Three regimes are identified: (i) stochastic, (ii) flow, and (iii) accelerated regime. Parameter values are those of curve (3) of Fig. 1. The temperature is T = 1 K.

of a curved wave, c, can be approximated by the eikonal equation D , (7) r where, D is a kind of effective diffusion coefficient. If the curvature (resp. the radius) becomes too high (resp. too low), the interface velocity may vanishes indicating a propagation failure. Then a critical radius for the growth initiation exists and is obtained by setting c = 0 in Eq. (7), giving rcrit = cD0 . c = c0 −

C. Spatio-temporal aspect of spin-transition

Let us focus now on the interface velocity. We present in Fig. 5(a) the spatial profile of the elastic HS:LS interface at some chosen instantaneous configurations through the spatial evolution  of the distance dij between two successive atoms ( d = (xj − xi )2 + (yj − yi )2 ) along a line which perpendicularly crosses the transformation front. It is noticed that while the electronic interface is sharp, the elastic inter-

face is clearly broad and extends over several nodes from both sides of the electronic interface. Such result is due to the misfit in lattice parameter between the HS:LS states leading to stress region easily observed on the profiles. Indeed, one can see that close to the interface, the HS (LS) region experiences a tensile (compressive) stress. This fact confirms the “localized” nature of the stress around the interface region, a behavior which is in excellent agreement with experimental optical microscopy data of [Fe(btr)2 (NCS)2 ]H2 O spin-crossover single crystal.6 In addition, we can easily notice from the profiles of Fig. 5(a) that the elastic interfaces spread out in the HS region much more than in the LS, a behaviour which is attributed to the soft elastic nature of the HS state compared to the LS state. From the elastic profiles, the interface position could be easily monitored as depicted in Fig. 5(b). The interface propagation depicts three regimes: (i) an incubation regime (stochastic) from 0 to 300 MCS during which the front is not yet established, followed by (ii) a flow (deterministic) regime in the interval [380−540] MCS in which the front moves at constant velocity,

FIG. 6. (a) Maps of the divergence of the displacement field during the relaxation process described in the snapshots of Fig. 2. The reference state is the HS phase, which leads to a negative value of the divergence in the saturated LS state. (b) Snapshots of the rotational of the displacement field at positions A to F of Fig. 2.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.255.6.125 On: Mon, 10 Aug 2015 23:18:33

194706-7

Slimani et al.

J. Chem. Phys. 139, 194706 (2013)

V = 0.075 nm MCS−1 , and (iii) the process is ended by an accelerated regime when the system “feels” the surface. Now we focus on the stresses generated by the LS domain propagation process of this circular lattice. For that, we determined the elastic stresses at the various positions (A ↔ F) through the displacement field calculations as we did in a previous work26 on a 2D square lattice along the thermal hysteresis loop. The displacement field u(ix , iy ) associated with the lattice site (ix , iy ) is defined as u(ix , iy ) = r(ix , iy ) − r0 (ix , iy ) where r0 (ix , iy ) and r(ix , iy ) are the initial and final atomic positions of the site (ix , iy ). Here, we used the positions of nodes in the perfect HS state as a reference state, i.e., r0 (ix , iy ) = ix R0H H , iy R0H H . Deformations are obtained by space derivation of the displacement expressions. In the continuum approximation, uxx =

ux (ix + 1, iy ) − ux (ix , iy ) ∂ux .  ∂x R0H H

(8)

uyy , uxy , and uyx are also given in the same way. We have calculated the displacement field and its spatial distribution upon the relaxation process at the positions A-F along the relaxation curve of Fig. 1. We derived the spatial distribution of the deformation field for which we depict the divergence (=uxx + uyy ) and the rotational (=uxy − uyx ) in Fig. 6, which respectively bring information about the local relative volume change and the pure shear strains. The divergence of the displacement field, shown in Fig. 6(a), mapping the local dilatation, starts from 0 in the HS state (the reference state) and becomes negative in the LS state, denoting a shrinking of the system during the relaxation. It is also remarked that the significant changes of the dilatation field are mainly located around the HS:LS interfaces. In Fig. 6(b), we present the corresponding snapshot of the rotational of the displacement field, denoted here distortion field, which represents the shear stresses acting in the material. The latter field gives evidence for an enhancement of the shear stresses at the intersects of the HS:LS interface and the boundary of the lattice. We draw in Fig. 7 the distribution of local “kinetic energy” of the lattice sites, given by 12 mvi2 , where the mass m of the molecule is taken equal to 1 and vi2 = vx2i + vy2i + vz2i . We observe that the propagation of the interface is faster at the border of the lattice than in the center, one of the reasons for

FIG. 7. Local distribution of the kinetic energy during relaxation process of Fig. 1 showing that the sites at the lattice’s border move much more rapidly than the bulky atoms.

which the interface straights during the process. This behavior could be also assigned to energetic considerations, since the motion of sites at the border is made easier, because the elastic energy of border atoms is merely smaller than that of bulky atoms. IV. CONCLUSION

In this paper, we have considered the relaxation of the metastable HS state at low-temperature using an electroelastic model for spin-crossover solids running on a circular lattice. We pointed the properties of the macroscopic nucleation, growth, and propagation mechanisms and its dependence with the geometry of the lattice. We demonstrated the intimate relation existing between the lattice shape and the interface propagation process. In particular, we found that the elastic energy of the system is mainly stored in the interface region, except a small part which is also distributed at the surface of the lattice, where the energy barrier for nucleation is smaller due to the lack of neighbors. Based on simple geometric considerations, in which we assumed a right contact angle between the HS:LS interface and the lattice border, we predicted analytically the shape of the interface, its length and the surface of the transformed phase as function of the macroscopic geometric parameters of the lattice. Monte Carlo simulations have very well confirmed our predictions, demonstrating that the geometry controls the main features of the interface dynamics in this electro-elastic system. Further more, in this contribution we clarified the process of formation of the interface, which requires a critical size of LS domain, to overcome the elastic energy barrier. As a consequence, the front velocity depends on the inverse of the interface radius, following the eikonal equation. On the other hand, the fine analysis of the Monte Carlo data upon the propagation process revealed the existence of three regimes that we have identified as corresponding to (i) an incubation (ii) flow, and (iii) acceleration process. The latter is clearly attributed to surface effects, which accelerates the relaxation of the spin state when the system feels the surface. These behaviors are in qualitatively good agreement with experimental optical microscopy observations that we recently reported6, 7 on spin-crossover single crystals. We have also analyzed (although not shown here) the effect of the elastic interaction on the propagation of the HS:LS interface and we found that its speed slows down in rigid lattices. Moreover, we found that the position of the interface depends quadratically on MC time, which means that the propagation velocity is not uniform. Such behavior can be clearly ascribed to the shape of the front transformation, which itself is connected to that of the lattice. Such results suggest new experimental investigations, like creating at low-temperature a radially shaped LS droplet by a photo-thermal process and study the transformation process. Despite of the existence of some experimental attempts36 in this direction, in order to control the front propagation, they unfortunately mainly led to create irreversible structural defects (holes) instead of photoproduced HS domains. The presence of the defects, obviously, perturbs the propagation process of the front transformation. But, we think that much more experimental efforts remains to be done

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.255.6.125 On: Mon, 10 Aug 2015 23:18:33

194706-8

Slimani et al.

before to produce reversibly HS droplets having different radii in order to analyze the relation between the velocity and the shape of the domains. From the theoretical point of view, the study of elliptic systems constitutes also an interesting extension to the present investigations. As a general remark, it is important to point out one limitation of the microscopic elastic models presently used to describe the spin transition, which assumes an isotropic expansion of the unit cell (or lattice parameter) upon the transformation from LS to HS. In fact, the experimental situation is much more complex. In the case of the Fe(bbtr), for example, it is observed that a crystallographic transformation accompanies the spin transition leading to a twinned low-temperature phase, while for [Fe(btr)2 (NCS)2 ].H2 O (btr=4,4 -bis-1,2,4-triazole),37 which has a polymeric structure, x-ray diffraction studies38 have shown the existence of a clear anisotropic expansion upon the spin conversion. Indeed, the LS to HS transition results in an expansion of the unit cell along the b and c axis and a shrinking along the a axis, although the total volume increases. Another example comes from the spin crossover complex, [Fe(bapbpy)(NCS)2] where bapby=N-(6-(6-pyridin2ylamino)pyridin-2-yl)pyridin-2-amine, in which the first example of a straight interface at the thermal spin transition was reported.12 It was also the first example of a robust crystal undergoing several successive thermal hysteresis loops without damage. The analysis of the structural data12 during the transition from LS to HS shows an anisotropic expansion along = 1.2%, the three axis, with the following relative values: a a b c = 2.2%, and = 0.7%. Other examples of anisotropic b c expansion of the unit cell volume through the spin transition exist in literature, which imposes an urgent extension of the current microscopic elastic models, designed to mimic this phenomenon. From the point of view of domain propagation, it is expected that such structural anisotropy will affect the domain propagation velocity and the orientation of the HS/LS interface. We are presently handling this problem in order to explain the stable and reproducible interface orientation observed by optical microscopy in several spin-crossover complexes, such as in [Fe(NCSe)(py)2 ]2 (m-bpypz).8 The second useful extension of the present microscopic elastic models concerns the study of three dimensional objects. In this line, a first step was realized in a recent study39 of spin-crossover membranes and thin films in which a third degree of freedom of the atomic displacements has been included: the atoms can move out of the planes. New behaviors, like buckling and crumpling of the sheets have been evidenced and they seem to be in good agreement with recent experimental observations of the spin transition on nanocrystals.40 However, much more effort is still needed to describe real 3D spin crossover materials. 1 Topics

in Current Chemistry, Spin Crossover in Transition Metal Compounds I-III, edited by P. Gütlich and H. A. Goodwin (Springer, Heidelberg/Berlin, 2004), pp. 233–235. 2 C. P. Köhler, R. Jakobi, E. Meissner, L. Wiehl, H. Spiering, and P. Gütlich, J. Phys. Chem. Solids 51, 239 (1990). 3 R. Jakobi, H. Spiering, and P. Gütlich, J. Phys. Chem. Solids 53, 267 (1992).

J. Chem. Phys. 139, 194706 (2013) 4 K.

Boukheddaden, J. Linares, E. Codjovi, F. Varret, V. Niel, and J. A. Real, J. Appl. Phys. 93, 7103 (2003). 5 F. Varret, A. Slimani, K. Boukheddaden, C. Chong, H. Mishra, E. Collet, J. Haasnoot, and S. Pillet, New J. Chem. 35, 2333–234 (2011). 6 A. Slimani, F. Varret, K. Boukheddaden, C. Chong, H. Mishra, J. Haasnoot, and S. Pillet, Phys. Rev. B. 84, 094442 (2011). 7 C. Chong, A. Slimani, F. Varret, K. Boukheddaden, E. Collet, J. C. Ameline, R. Bronisz, and A. Hauser, Chem. Phys. Lett. 504, 29 (2011). 8 A. Slimani, F. Varret, K. Boukheddaden, D. Garrot, H. Oubouchou, and S. Kaizaki, Phys. Rev. Lett. 110, 087208 (2013). 9 F. Varret, C. Chong, A. Slimani, D. Garrot, Y. Garcia, and D. Naik Anil, “Real-time observation of spin-transitions by optical microscopy,” in SpinCrossover Materials: Properties and Applications, 1st ed., edited by Malcolm A. Halcrow (John Wiley & Sons, Ltd., Oxford, UK, 2013). 10 S. Bedoui, G. Molnár, S. Bonnet, C. Quintero, H. J. Shepherd, W. Nicolazzi, L. Salmon, and A. Bousseksou, Chem. Phys. Lett. 499, 94 (2010). 11 S. Bedoui, M. Lopes, W. Nicolazzi, S. Bonnet, S. Zheng, G. Molnár, and A. Bousseksou, Phys. Rev. Lett. 109, 135702 (2012). 12 S. Bonnet, G. Molnár, J. S. Costa, M. A. Siegler, A. L. Spek, A. Bousseksou, W. Fu, P. Gamez, and J. Reedijk, Chem. Mater. 21, 1123 (2009). 13 S. Ohnishi and S. Sugano, J. Phys. C 14, 39 (1981). 14 H. Spiering, K. Boukheddaden, J. Linares, and F. Varret, Phys. Rev. B 70, 184106 (2004). 15 K. Boukheddaden, Phys. Rev. B 88, 134105 (2013). 16 N. Willenbacher and H. Spiering, J. Phys. C. 21, 1423 (1988). 17 H. Spiering and N. Willenbacher, J. Phys. Condens. Matter 1, 10089 (1989). 18 T. Nakada, P. A. Rikvold, T. Mori, M. Nishino, and S. Miyashita, Phys. Rev. B. 84, 054433 (2011). 19 M. Nishino, K. Boukheddaden, Y. Konishi, and S. Miyashita, Phys. Rev. Lett. 98, 247203 (2007). 20 M. Nishino, C. Enachescu, S. Miyashita, K. Boukheddaden, and F. Varret, Phys. Rev. B 82, 020409(R) (2010). 21 C. Enachescu, L. Stoleriu, and A. Stancu, Phys. Rev. Lett. 102, 257204 (2009). 22 C. Enachescu, M. Nishino, S. Miyashita, L. Stoleriu, A. Stancu, and A. Hauser, Eur. Phys. Lett. 91, 27003 (2010). 23 C. Enachescu, M. Nishino, S. Miyashita, L. Stoleriu, and A. Stancu, Phys. Rev. B. 86, 054114 (2012). 24 M. Nishino, C. Enachescu, S. Miyashita, P. A. Rikvold, K. Boukheddaden, and F. Varret, Sci. Rep. 1, 162 (2011). 25 W. Nicolazzi, S. Pillet, and C. Lecomte, Phys. Rev. B. 78, 174401 (2008). 26 A. Slimani, K. Boukheddaden, F. Varret, and H. Oubouchou, Phys. Rev. B. 87, 014111 (2013). 27 J. Jung, F. Bruchhäuser, R. Feile, H. Spiering, and P. Gütlich, Z. Phys. B. 100, 517 (1996). 28 L. Wiehl, H. Spiering, P. Gütlich, and K. Knorr, J. Appl. Crystallogr. 23, 151 (1990). 29 K. Binder, in The Monte Carlo Method in Condensed Matter Physics (Springer, New York, 1995). 30 R. Bronisz, Inorg. Chem. 44, 4463 (2005). 31 A. Bousseksou, J. Nasser, J. Linares, K. Boukheddaden, and F. Varret, J. Phys. I France 2, 1403 (1992). 32 K. Boukheddaden, I. Shteto, B. Hoo, and F. Varret, Phys. Rev. B. 62, 14796 (2000). 33 W. C. Yu and P. Gielisse, Mater. Res. Bull. 6, 621 (1971). 34 N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953). 35 The Neumann boundary condition is one of the three types of boundary conditions commonly encountered in the solution of partial differential equations. It specifies the normal derivative of the function, U, for example,  = f (r , t), where n is the normal vector to the on a surface. It writes, n.∇U surface at r. 36 S. Bedoui, M. Lopes, S. Zheng, S. Bonnet, G. Molnár, and A. Bousseksou, Adv. Mater. 24, 2475 (2012). 37 W. Vreugdenhil, J. G. Haasnoot, O. Kahn, P. Thuéry, and J. Reedijk, Polyhedron 9, 2971 (1990). 38 S. Pillet, J. Hubsh, and C. Lecomte, Eur. Phys. J. B 38, 541 (2004). 39 K. Boukheddaden and A. Bailly-Reyre, Eur. Phys. Lett. 103, 26005 (2013). 40 R. M. van der Veen, O.-H. Kwon, A. Tissot, A. Hauser, and A. H. Zewail, Nat. Chem. 5, 395–402 (2013).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.255.6.125 On: Mon, 10 Aug 2015 23:18:33

Properties of the low-spin high-spin interface during the relaxation of spin-crossover materials, investigated through an electro-elastic model.

The present work is devoted to the spatio-temporal investigations of spin-crossover lattices during their thermal relaxation from high- to low-spin st...
6MB Sizes 0 Downloads 0 Views