- Home
- Documents
- Two-variable anharmonic model for spin-crossover solids: A like-spin domains interpretation

prev

next

out of 12

Published on

04-Dec-2016View

214Download

1

Transcript

Two-variable anharmonic model for spin-crossover solids: A like-spin domains interpretation

W. Nicolazzi, S. Pillet,* and C. LecomteLaboratoire de Cristallographie et Modlisation des Matriaux Minraux et Biologiques, UMR CNRS 7036, Institut Jean Barriol,

Nancy-Universit, BP 239, F-54506 Vandoeuvre-ls-Nancy, FranceReceived 9 June 2008; revised manuscript received 4 September 2008; published 4 November 2008

Spin-crossover SC complexes are one of the most fascinating examples of molecular bistability, whosesolid-state properties are tightly connected to cooperative interactions within the crystal lattice. A variety ofmacroscopic and microscopic models have been developed to explore the cooperative nature of the SC phe-nomenon. We present here a two-variable microscopic Ising-like model for SC solids, accounting for the elasticorigin of the cooperativity using coupled spin and translational lattice degrees of freedom. Within our model,the interaction between a pair of neighboring molecules in the crystal is dependent not only on their spin statesbut also on their separation distance, modeled by spin-dependent Lennard-Jones LJ potentials. This schemeleads explicitly to local variations in the interactions, associated to the local strain induced by the moleculeschanging their spin state. In essence, the LJ potentials provide the anharmonicity of the crystal lattice. Theequilibrium quasistatic properties of the proposed Hamiltonian are analyzed by Monte Carlo simulations ona regular deformable square lattice. We show that the spin dependence of the LJ potentials breaks the spin-statesymmetry in the free energy. The interplay between spin and lattice degrees of freedom shows itself in thetemperature evolution of the fraction of high-spin molecules and the mean lattice spacing, as a function of theintersite coupling. For strong coupling, like-spin domains nucleate and develop, evidenced by a double struc-ture in the distribution of lattice spacings; structural relaxation occurs at the domain walls. In the weaklycooperative situation, the mean lattice constant scales directly with the fraction of high-spin species; structuralrelaxation spans the entire system.

DOI: 10.1103/PhysRevB.78.174401 PACS numbers: 05.50.q, 64.60.De, 75.30.Wx, 75.60.d

I. INTRODUCTION

The thermally induced spin-crossover SC phenomenonbetween the low-spin LS and the high-spin HS states ofFeII molecular complexes has been the subject of intenseresearch activities.1,2 The intrinsic bistability that these ma-terials exhibit originates from intramolecular vibronic cou-pling. The higher electronic and vibrational degeneracy inthe HS state leads to a LS to HS entropy increase, whichdrives the spin conversion. In addition to large variations inthe magnetic and optical properties upon the spin-statechange, drastic structural modifications occur, namely, atypical 0.2 Fe-N bond lengthening and a 2%5% unitcell expansion.3 Besides temperature or pressure changes, thespin conversion may alternatively be triggered by continuouslight excitation, albeit at very low temperature light-inducedexcited-spin-state trapping LIESST process,46 and also athigher temperature within the thermal hysteresis loop bymeans of an intense pulsed laser.7,8 Owing to theirtemperature- and light-induced switching potential, it is wellestablished that these molecular systems may have prospec-tive use as sensors or may be integrated in data storage de-vices.

In the solid state, strong electron-lattice coupling is at theorigin of the rich variety of behaviors SC materials exhibit,such as gradual spin crossover or abrupt spin conversionfirst-order character with often hysteresis in the case ofstrongly cooperative materials. Although the concept israther ill defined, cooperativity is often attributed to the largeHS to LS molecular volume change coupled to elastic inter-actions within the solid. In the framework of the continuumelasticity theory, it has been pictured as the buildup of an

internal pressure, proportional to the concentration of LSspecies.911 The local distortion strain of the crystal latticeaccompanying the spin transition has been traced back topoint defects at the iron site of the spin-changing molecule.Cooperativity shows itself in the abruptness of the thermaltransition and low-temperature sigmoidal HS to LS relax-ation curves, attributed to a self-accelerated process.12,13Polymeric SC materials with one-dimensional 1D chain ortwo-dimensional 2D layer structural architectures are em-blematic examples of highly cooperative systems, withsometimes extremely wide thermal hysteresis.14 In that case,the bridging ligands play the role of efficient spin-statepropagators and contribute to the anisotropic elasticity. Onthe other hand, for SC systems built from purely monomo-lecular entities, it is quite clear that the intermolecular con-tacts, through weak van der Waals, hydrogen bonds, or -stacking, transmit the interactions to long range and are thusresponsible for the elastic properties of the solid.15

Intensive theoretical investigations have been devoted tothe development of more or less sophisticated models of spintransition, treating cooperativity through various schemes.Among others, macroscopic models based on the thermody-namic theory of regular solutions were proposed. In the so-called Slichter and Drickamer SD model, cooperativity isdescribed in mean field through phenomenological enthalpyterms related to the intermolecular interactions.16 In the elas-tic model for SC systems introduced by Spiering andWillenbacher,911 SC molecules are considered as entitieswhose volume and shape depend on their spin state, insertedin an isotropic elastic medium. When a molecule undergoesthe spin transition, it creates a local strain which forces thenearest-neighbor molecules to undergo the spin transition aswell ferroelastic coupling. In its more elaborated form, the

PHYSICAL REVIEW B 78, 174401 2008

1098-0121/2008/7817/17440112 2008 The American Physical Society174401-1

Spiering scheme describes the crystal lattice within the De-bye approximation and allows for anharmonicity through theGrneisen approximation.11,17 A variety of microscopicIsing-like models were also designed, introducing a Hamil-tonian of interacting two-level units with different energiesand degeneracies, as proposed by Wajnflasz and Pick18,19 andBousseksou et al.20 These models can reproduce most of theequilibrium properties of SC materials as well as their dy-namic behavior, including the sigmoidal relaxation, the two-step transitions, and the photoinduced effects.2127 The iso-morphism between the macroscopic SD and microscopicIsing-like model in the mean field has been formulated.24,28A direct relation between the Ising coupling constant J andthe phenomenological cooperative factor in the SD schemehas been derived. The 1D version of the Ising-like model hasbeen further solved analytically.29 An alternative one-dimensional model of SC molecules coupled to a phononfield atom-phonon coupling model was established, the in-termolecular interactions being introduced as harmonicsprings.30,31 This 1D model was recently solvedanalytically32 and further extended to anharmoniccontributions.33 Various elastic models, for which the cou-pling constant depends on the intersite separation distance,have been proposed and investigated numerically throughmolecular dynamics MD methods34,35 and Monte CarloMC simulations.36,37

The mechanism by which the spin conversion proceeds inthe solid state is a fundamental question. In the case of co-operative materials, characterized by an abrupt thermal tran-sition and sigmoidal relaxation kinetics, it has been widelyargued that molecule clustering may play a key role. Thevery first experiments which have indirectly observed so-called like-spin domains LSDs are powder and single-crystal diffraction measurements.3841 These techniques areby no way sensitive to the molecular spin component but onthe contrary to lattice structural observables. In these ex-periments, LSD shows itself at the thermally and light-induced transitions by a clear separation of Bragg peaks cor-responding to ordered HS and LS phases with different cellconstants phase separation. Within this framework, a do-main consists of adjacent molecules, crystallographically or-dered to sufficiently long range i.e., comparable to the x-rayor neutron spatial coherence of the experiment, and withwell-defined intermolecular separation distances. The con-cept of LSD has emerged quite early in the macroscopic andmicroscopic model formulations, in which it is usually pic-tured as a set of correlated like-spin adjacent molecules. De-spite their success, Ising-like models are at present unable toexplain the Avrami kinetics of LSD nucleation and growthderived from diffraction techniques.38 This is well under-standable since under the assumption of pseudospins onrigid-lattice sites, the Ising-like models cannot describechanges in the positions of the molecules accompanying theelectronic configuration ordering and the associated local lat-tice distortion. Hence crystallographic LSD would not bedefined. The importance of going beyond a fixed-lattice Isingmodel has been stressed for investigating ordering process inCu-Au alloys42 or lipid bilayers,43 for instance. In the presentstudy, the strategy for introducing magnetoelastic couplingsin the Ising-like model of SC molecular solids is based on

the development of the intersite intermolecular interactionson pairwise potentials of the Lennard-Jones LJ form, con-sidering spin and lattice degrees of freedom through twocoupled on-site variables, which allows for continuous mo-lecular displacements. As a general matter of fact, these po-sition changes are expected to highly affect the characteris-tics and even the nature of the transition since, even at thelevel of global elastic deformations coupled to a nearest-neighbor Ising model, the behavior is different from that of arigid-lattice Ising model.44 Our goal is to capture the essen-tial features of the elastic contribution to the solid-state co-operativity and reconcile experimental observations derivedthrough various techniques sensible either to the spin spec-troscopic or magnetic techniques or the lattice variablesdiffraction techniques.

The paper is organized as follows. Section II is devoted tothe presentation of the Hamiltonian and the technical detailsof the Monte Carlo simulations. The equilibrium propertiesof the model are described in Sec. III, as well as the depen-dence of the transition mechanism spin crossover or first-order transition on the parameters of the model; the condi-tions for the presence of hysteresis are established. Wefinally conclude in Sec. IV.

II. MODEL AND METHOD

A. Anharmonic Ising-like model

We adopt the two-level scheme of the Ising-like model, inwhich the two molecular states are described by fictitiousscalar spin operators Fig. 1. The HS and the LS states arerepresented by = +1 and =1, respectively, with corre-sponding g+ and g degeneracies of both electronic angularand spin and intramolecular vibrational origins. Consideringthe isomorphism between this degenerate Ising-like modeland the Ising model under a temperature-dependent field,45the on-site Hamiltonian which accounts for the inner degreesof freedom of N SC molecules is written as

H0 =effT

2 i=1N

i, 1

where effT=kBT lng is the temperature-dependentfield, with as the difference in ligand-field energy betweenthe two levels and g=g+ /g1 as the effective degeneracyratio, related to the LS to HS electronic and vibrational en-tropy increase S=kB lng. This two-level scheme is a dras-tic but appropriate simplification of the complete vibronicsystem.20 The S term in the effective field is the key ingre-dient which drives the SC process.

In the standard microscopic Ising-like model, neighboringi and j sites are coupled by a common ferroelastic interac-tion Jij, irrespective of the spin states: JHS-HS=JHS-LS=JLS-LS=J. The interaction term Jiji j is thus phenom-enologically added to Eq. 1 by analogy with ferromag-netism. The thermodynamic properties of the Ising-likemodel are governed by the sign and temperature dependenceof the effective field together with this coupling term.

Starting from the formal expression of the Ising-likemodel, we introduce additional degrees of freedom to ac-

NICOLAZZI, PILLET, AND LECOMTE PHYSICAL REVIEW B 78, 174401 2008

174401-2

count explicitly for the lattice energy and the local distortionstrain of the crystal lattice, induced by a single moleculeswitching its spin state. Let rixi ,yi be the instantaneousposition vector of the ith site and ri

0 the equilibrium positionvector in the undistorted lattice. The configuration of thesystem is now characterized by 3N variables in two dimen-sions 1 ,x1 ,y1 , . . . ,i ,xi ,yi , . . . ,N ,xN ,yN. The interactionenergy can be written in the general form as a lattice sum:

Hlat,r = i,j

Ai, jVelastri,j,ri,j

0 , 2

where the sum runs over all nearest-neighbor pairs i , j.ri,j= r j ri and ri,j

0= r j

0ri

0 are the neighboring intersiteinstantaneous distance and equilibrium distance, respec-tively. The interaction energy depends therefore only on thelocal environment of each independent molecule, such as thedistribution of other molecules and their distance to it, aswell as their spin state Fig. 1b; this is at variance with thestandard Ising-like model. In another context, it has beenformulated in a similar way for modeling magneto-elasticcouplings in three-dimensional 3D Ising schemes.4648 An

important issue is that the intersite distance corresponds tothe FeFe separation, as illustrated in Fig. 1, and shouldnot be confused with the true intermolecular distances usu-ally characterized from structural analyses e.g., HO hy-drogen bonds and CC short contacts.

In some cases, it has been recognized that in addition tothe electronic and intramolecular vibrational entropychanges, the latter being mostly related to the Fe-N stretch-ing and octahedron distortion modes, intermolecular latticevibrations may also contribute significantly.49 Zimmermannand Konig50 applied a simple Debye model with two differ-ent Debye temperatures HS and LS to account for theentire vibrational contribution. From accurate single-crystaldiffraction measurements, we have concluded that the 2Dpolymeric SC compound Febtr2NCS2 H2O exhibits lat-tice modes of higher vibration amplitudes in the HS state.Related to this, the HS and LS crystal lattices exhibit a dif-ferent thermal-expansion behavior.51 The thermoelastic prop-erties for such SC systems is obviously dependent on thespin state, and therefore has to be considered in the model.As a consequence, the independence of the interaction pa-rameter J with the spin state in the standard Ising-like modelis a drastic approximation, which we relax here. The equilib-rium distance ri,j

0 in the undistorted lattice and the elasticcoupling Ai , j between a pair of molecules i and j areconsidered as follows:

Ai, j = AHS ri,j

0

= rHS if i = j = 1,=AHL =rHL if i = j ,=ALS =rLS if i = j = 1.

3

The intersite equilibrium distances correspond to the in-termolecular distances in the HS and LS phases called lat-tice spacing for simplicity hereafter. We recall that they arenot relative to the internal Fe-N distance. Although the inter-molecular distances are known to be different in the HS andLS phases, typical values cannot be given a priori.

Condition 3 can be formally rearranged to

Ai, j = J0 + J1i + j + J2i j , 4

with

J0 =AHS + 2AHL + ALS

4,

J1 =AHS ALS

4,

J2 =AHS 2AHL + ALS

4. 5

A relation similar to Eq. 4 was derived separately byBolvin and Kahn28 and Boukheddaden et al.32 In our case,Velast in Eq. 2 takes the form of an anharmonic intersitepotential of the empirical 6-3 LJ form with finite rangermax:

i if ri,j rmax

ri

0

(a)

(b)

|| - ||r ri i0

LS

HS

i g

i g

ri,j

Velast

FIG. 1. Color online Schematic picture of the two-variablemodel. a HS molecule of spin variable i at position ri, displacedfrom its equilibrium position ri

0. b Molecules on a 2D lattice with

separating distance ri,j = r j ri and interacting through intermo-lecular Lennard-Jones potentials Velast.

TWO-VARIABLE ANHARMONIC MODEL FOR SPIN- PHYSICAL REVIEW B 78, 174401 2008

174401-3

Velastri,j,ri,j

0 = ri,j0

ri,j

6 2 ri,j0

ri,j

3,

ii if ri,j

rmax

Velastri,j,ri,j

0 = 0. 6

The first term in Velast results from strong short-range repul-sive interaction, while the second term represents the long-range attractive interaction. Velast is a dimensionless quantity;the real elastic interaction is given by Ai , jVelast, as speci-fied in Eq. 2. The choice of rmax is crucial to ensure thestability of the crystal lattice and retain the structural topol-ogy each site surrounded by four neighbors on the squarelattice. In other words, there is no diffusion of moleculesacross the whole system. It is to be noted that rmax preventsthe algorithm from accessing the entire phase space spannedby the translational degrees of freedom. Three different po-tentials are considered for HS-HS, HS-LS, and LS-LS neigh-boring pairs Fig. 2. The position of the minimum is relativeto the corresponding equilibrium distance rHS

0, rHL

0, rLS

0, de-

fining three lattice spacings. The potential depth is related tothe strength of the intermolecular coupling AHS, AHL, ALS.For all the simulations, the equilibrium distances of the un-distorted lattices take the values rHS

0=1.2, rHL

0=1.1, and rLS

0

=1, which correspond to a more compact LS phase; this isthe experimentally observed usual trend. Intersite equilib-rium distances between LS and HS species are intermediatebetween purely LS and HS ones. These chosen values for theequilibrium distances are not relative to any specific com-pound, but exaggerated to amplify the effects. The essentialcharacteristic of the selected LJ potential is the asymmetry,

leading to the desired anharmonicity. In a first approxima-tion, harmonic potentials can successfully predict manyproperties of the solid. However, anharmonic terms in theintermolecular potentials are responsible for several phenom-ena, important in the context of SC, such as thermal expan-sion, compressibility, and thermal conductivity. The chosenparameters of the LJ potentials ensure a stiffer lattice in theLS state with respect to the HS state, and correspond to alower bulk modulus B in the HS state. B is related to thesecond derivative of the crystal energy with respect to thevolume. In the case of pairwise additive central intermolecu-lar potentials, B is given by the second derivative of thepotential with respect to the intersite distance. Consideringour scheme, B for the pure HS or LS lattices is roughlyrelated to the curvatures of the LJ potentials at the equilib-rium HS and LS distances. It is less straightforward for theintermediate-spin configurations. It is clear from Fig. 2 thatthe curvature is more pronounced in the LS state. For quan-titative confrontation with experimental data, a more effi-cient intermolecular force field would be more appropriate,possibly of the Buckingham exp-6 form, augmented withelectrostatic interactions through distributed multipoleexpansions.52 The optimum parameters of the LJ potentialscould be determined through an appropriate optimizationprocedure by the knowledge of the cohesive crystal energy,lattice parameters, and thermal-expansion tensors of the HSand LS pure phases.53

The total two-variable Hamiltonian H , r can finallybe written as

H,r = i,j

Velastri,j,ri,j

0 J0 + J1i + j + J2i j

+eff

2 i i, 7

where J0Velast takes only negative values, whereas J1Velastand J2Velast may be either positive or negative, depending onthe respective depth of the LJ potentials. When J2Velast0,the system is ferroelastic, while for J2Velast0, the system issaid to be antiferroelastic. Although ferroelastic materials isthe general case, a few examples of the second kind havebeen discovered so far.54 This expression of the interactionenergy as the sum of three terms is similar to a formulationof Ising-like model based on elasticity theory,32,55 althoughfor the latter, the interaction energy is identical for HS-HSand LS-LS pairs and the corresponding J0 term is a constantenergy shift. We chose here the cohesive energy of the pureHS and LS phases as different, which is rationalized by thefact that both phases usually exhibit modifications of thesupramolecular organization, even possibly related to aspace-group change crystallographic phase transition. Con-trasting with the standard Ising-like model, the LJ pair po-tentials lead to coupling parameters which are defined locallyand depend on the pair i , j. The entire system exhibitstherefore a distribution of interaction constants. A similaridea has also been introduced in a thermodynamic model toexplain unusual spin transition features.56 The introductionof lattice degrees of freedom allows the couplings to fluctu-ate when the system is in contact with a thermal bath, which

FIG. 2. Spin-dependent Lennard-Jones potentials Ai , jVelastdescribing elastic interactions between two nearest-neighbor sites iand j. The position and depth of the potentials define three latticespacings rHS

0,rHL

0,rLS

0 and three intermolecular couplingsAHS,AHL,ALS. rmax is specified as a vertical dotted line.

NICOLAZZI, PILLET, AND LECOMTE PHYSICAL REVIEW B 78, 174401 2008

174401-4

we consider in the following. The system can be character-ized by two quantities: the magnetization per site and themean intersite distance r lattice spacing defined as

=1Z i=1N iNRk=1N drk exp H , 8

r =1Z Rk

N

drki,j

Mri,j

M exp H , 9

where =1N corresponds to a sum over allspin configurations and Rk=1

N drk= dr1

drN to

a continuous sum over all position vector configurations ac-cessible for the system. Z=Ri=1N dri expH is thecanonical partition function and expH is the Boltzmannweight with =1 /kBT corresponding to the inverse tempera-

ture kB=1. M is the number of intersite bonds.

B. Dynamical aspects and Monte Carlo simulations

In order to study the static properties of Hamiltonian 7,we perform MC simulations on a simple square deformablelattice with free boundary conditions. We adopt the NVT-MCmethod, where N is the number of sites molecules, V is thevolume of the system, and T is the thermal bath temperature.The use of periodic conditions would have required atreatment in the isobaric N , P ,T statistical ensemble, whereP is the pressure of the system, allowing the area of thesystem to fluctuate through a length rescaling. LetP1 , . . . ,N ,r1 , . . . ,rN , t be the probability for the systemto have the spin configuration 1 , . . . ,N and the positionvector configuration r1 , . . . ,rN at Monte Carlo step MCSt. We assume that the P , r , t evolution is governed bythe master equation

P,r,t + 1 = i=1

N

driN

,r ,rP,r,t ,r ,rP,r,t , 10

where , r , r is the transition rate corre-sponding to the conditional probability for the system tohave the , r configuration at MCS t knowing it was , r at MCS t1. The transition probability is imposedto follow the Boltzmann equilibrium distribution. The sta-tionary condition P , r , t /t=0 leads to the well-known necessary but nonsufficient detailed balance condi-tion

,r ,r,r ,r =

exp H,rexp H,r . 11

In the following, we choose two transition probabilities,Wspin for the spin variables and Welastr r for the displacement vectors, which is justified by thedifferent characteristic time at which electronic and latticedeformation processes in the adiabatic approximation takeplace. As a matter of fact, electron transitions occur withinfemtoseconds, while atomic displacement and vibrational re-laxation occur in the picosecond to nanosecond time scale. Anonconserved order-parameter dynamic is chosen forWspin and Welastr r. For dealing with thenonequilibrium behavior of the SC phenomenon, the appro-priate choice of the transition rate would be of Arrheniustype, which allows one to correctly reproduce the sigmoidalrelaxation curves of the HS fraction at very lowtemperature.21,22,27 In the present case, we restrict the analy-sis to the equilibrium properties of our model; the nonequi-librium study will be the subject of a future publication. Wechoose the Metropolis single-variable-changing dynamics57for the two variables of the Hamiltonian. The different se-

quences of a MCS are as follows: i A lattice site i is chosenrandomly and the fictitious spin i flip is updated accordingto the Metropolis acceptance criterion. ii Another site j atthe position r jxj ,yj is chosen randomly and a new positionr jxj ,yj is proposed as

xj = xj + dxj ,

yj = yj + dyj , 12

where dxj and dyj are random continuous displacementsdrawn on a Gaussian distribution of zero mean and adjust-able variance. The new position of the site is evaluated andupdated using the Metropolis scheme. iii The two previousstages are repeated N times. The spin and lattice variables arethus technically decoupled in the algorithm; the real inter-play comes from the Hamiltonian. MC simulations havebeen performed on 3232 square lattice N=1024 withopen boundary conditions. In addition, for simplicity, J1 hasbeen set to 0 in all calculations, corresponding to ALS beingequal to AHS. Although this assumption simplifies our Hamil-tonian, this does not bring a real severe limitation since J1mostly affects the transition temperature vide supra.

III. EQUILIBRIUM BEHAVIOR OF THE ANHARMONICMODEL

The properties of the model defined by Hamiltonian 7,called hereafter as anharmonic model, are described accord-ing to the behavior of the intrinsic parameters nHS and ri,j

norm.

nHS can be defined as a function of the fictitious magnetiza-tion:

TWO-VARIABLE ANHARMONIC MODEL FOR SPIN- PHYSICAL REVIEW B 78, 174401 2008

174401-5

nHS =1 +

2. 13

The dimensionless and normalized lattice spacing ri,j

norm is

relative to the distance between two neighboring sites i and j.Its thermodynamic average rnorm is given by

rnorm =r rLS

0

rHS0

rLS0 , 14

where r is the average intersite distance in the system.

A. eff=0 case

We discuss in a first step the behavior of the anharmonicmodel without effective field and setting J1=0, and interpretit with respect to the phase diagram of the well-known zero-field Ising model defined by

HIsing = J2i,j

i j , 15

where J2 is the usual coupling constant. Both models anhar-monic and zero-field Ising model have been treated usingidentical system size to have an estimate of finite-size ef-fects, the true Curie temperature of the 2D Ising model at thethermodynamic limit being indeed well known.58 The phasediagram is derived using the following procedure. At t=0 thesystem is prepared in a totally ordered HS nHS=rnorm=1 orLS nHS=rnorm=0 configuration. For t0 the system is incontact with a thermal bath at a temperature T. The observ-ables nHS and rnorm are estimated as the average over 6

105 MCSs, using only configurations after a waiting timew i.e., disregarding configurations in the nonequilibriumtransient regime. For the different simulations, w rangefrom 104 to 105 MCSs per spin for temperatures far from andnear the critical point, respectively. The results are given inFig. 3.

Just as the Ising case, the anharmonic model presents adisordered nHS=1 /2 phase at high temperature and an or-

dered phase nHS1 /2 at low temperature, separated by acritical temperature TC not estimated in this work.59 It is wellknown that in the ordered phase, the free-energy density ofthe Ising model f IsingnHS has two symmetric and degeneratestable solutions nHS+

0 T and nHS0 T for eff0+ and eff

0

, respectively, surrounding the unstable solution nHS=1 /2. Obviously, it is expected that the application of a non-zero field eff would break the symmetry between nHS+

0 andnHS

0. Surprisingly, even if the temperature-dependent field

eff is absent, the anharmonic model does not have this sym-metry between the two states. This is due to the deliberatechoice of different LJ potentials for HS and LS states Eq.6, which corresponds to an additional entropy differenceSadd between the HS and the LS phases. Sadd acts as anadditional field, noted as h , which breaks the symmetry be-tween the two states in the free-energy density fnHS, asschematically represented in Fig. 4. This additional field is akey feature of this anharmonic model. Under such condi-tions, the overall effective field eff+h cannot be easilycanceled since Velastri,j ,ri,j

0 itself depends on the spins iand j due to the fact that ri,j

0 depends on the spin values ofthe sites i and j. Figure 3 indicates that the more compactnHS

0 T phase is the stable one, while nHS+0 T is metastable.

The additional entropy term can be formally written as addi-tional degeneracies, originating from lattice vibrationsphonons, as will be shown by solving analytically Hamil-tonian 7 in a forthcoming paper. For a purely harmonic 1DSC chain, the corresponding degeneracy is simply related tothe ratio of LS and HS elastic constants.32 A similar depen-dence on ALS and AHS, but a more complicated one, doesexist for the present anharmonic model.

B. eff0 case

We now study the eff0 case, setting lng+ /g=4 forall MC simulations. The two-dimensional Ising-like modelunder temperature-dependent field has been widely studiednumerically24 and exactly in the mean-field approximation.21

FIG. 3. Temperature dependences of nHS for the anharmonicsquares and zero-field finite-size Ising model triangles. Systemsprepared in ordered HS and LS phases are represented by openand filled symbols, respectively. Parameters of the model: J0=999,J1=0, and J2=1.

FIG. 4. Schematic representation of the free-energy densityfnHS in the ordered phase TTC for the zero-field Ising fullline and the anharmonic dotted line model with effT=0.

NICOLAZZI, PILLET, AND LECOMTE PHYSICAL REVIEW B 78, 174401 2008

174401-6

Following previous treatment, the corresponding transitiontemperature Tequ, at which =0 and nHS=1 /2, is definedby the condition

eff = kBTequ lng = 0, 16

which gives

Tequ =

kB lng. 17

This definition of Tequ is an approximation; the true transitiontemperature corresponds to the temperature at which the freeenergies of the HS and LS states are equal. However, MCsimulations do not allow the determination of the free energyof the system, which justifies the use of this approximation.For the Ising-like model, the equilibrium temperature isknown to be independent of the coupling parameter J2 seeEq. 17, while the sharpness of the spin transition increaseswith J2. Whether the change in nHS is continuous gradual ordiscontinuous first-order character depends on whether Tequis above or below the critical temperature of the correspond-ing zero field Ising model TIsing.58 The transition line betweenspin crossover and first-order transition in the J2- phasediagram of the anharmonic model has been determined in thefollowing way. Thermal hysteresis loops have been per-formed by varying J2 parameters for different fixed valuesuntil the total disappearance of bistability is observed, whichcorresponds to the transition line. As seen in Fig. 5, the an-harmonic model displays such a critical behavior but thetransition line differs from the finite-size Ising-like model.This suggests that the critical temperature TC of the anhar-monic model is different from that of the simple Ising case.The introduction of the anharmonic coupling favors the first-order transition behavior. In agreement with our findings, ithas been shown that the introduction of anharmonic contri-bution to the 1D atom-phonon model sharpens also the spintransition and drives the first-order character.33 As ex-pected, first-order character is observed for high J2 valuesand weak HS to LS ligand-field energy difference .

The two extreme cases, namely, those of first-order ther-mal spin transition and simple spin crossover, are representedin Figs. 6 and 7, respectively. Figure 6 shows an abrupt tran-sition with hysteresis. Snapshots of the configuration of thesystem at several temperatures display the different stages ofthe LSHS and HSLS phase transformations. The 1and 2 4 and 5 upper snapshots show germ nucleationof the HS LS state within the LS HS phase. The 3 and6 upper snapshots present LS and HS coarsening spin do-mains. Our two-variable model allows also probing the cor-responding configuration of the crystal lattice. On the lowersnapshots in Fig. 6, equilibrium distances in the undistortedLS and HS lattices are represented in white and black, re-spectively, while intermediate distances are represented ingrayscale. Each thermally induced germ nucleation leads to alocal lattice distortion, which grows to a structurally defineddomain of large area lower snapshots 3 and 6. Thesecoarsening domains coexist in the system near Tequ: the tran-sition is heterogeneous with structural phase separation. Theseparation between HS and LS domains is clearly not abrupt,

FIG. 5. J2- phase diagram of the anharmonic filled squaresand finite-size Ising-like models open circles J0=999 and J1=0. spin

lattice

FIG. 6. Top: Temperature dependences of nHS squares andrnorm circles in the anharmonic model with J0=999, J1=0, and=8. J2 has been chosen above the transition line defined in Fig. 5such as J2=0.9. The two observables are reported during the warm-ing filled symbols and cooling modes open symbols. Bottom:Snapshots of the configuration of the system at different tempera-tures in the hysteresis loop. Spin black: HS state; white: LSstate and lattice black: ri,j

norm=rHS

0 ; white: ri,j

norm

=rLS0 ; shaded gray:

rLS0 ri,j

normrHS0 observables are represented on top and bottom,

respectively. The thermal loops have been performed using 2000Monte Carlo transient steps at a heating rate of dT /dt=2

106 K MCS1 in the warming and then cooling modes, the ini-tial system configuration being the final configuration of the previ-ous temperature. Error bars are smaller than the symbols.

TWO-VARIABLE ANHARMONIC MODEL FOR SPIN- PHYSICAL REVIEW B 78, 174401 2008

174401-7

owing to the continuous character of the lattice variable; re-laxation of the lattice spacing occurs, which corresponds todomain walls. Such a behavior is obviously not accessible ina rigid-lattice model such as the standard Ising-like one. Thisdomain structure is expected for a first-order transition, dueto the short-range nature of the interactions in our model.The usual short-range Ising model itself presents such a char-acteristic clustering. Our results can be related and comparedto the elastic Ising-like model of Konishi et al.36 and Mi-yashita et al.,37 in which the interactions are developed on apurely harmonic intersite potential, furthermore identicalfor HS and LS. This latter restriction is dissimilar to ourscheme and leads to a mean-field behavior no short-rangeinteractions.37 This harmonic elastic model does not exhibitcluster growth, and accordingly cannot explain phase-separation effects, as evidenced by diffractiontechniques.3841

In the same way, the case of thermal spin crossover hasbeen investigated Fig. 7. No metastable state exists; theprocess is dominated by nucleation leading to a gradual con-version without hysteresis. All along the transition, the latticespacing is intermediate between the equilibrium distances inthe HS and LS undistorted phases; the transition is homoge-neous. Structural relaxation spreads the whole system. In ad-dition, the lattice spacing directly follows the fraction of HSspecies; in other words, the Vgard law is satisfied.

rnorm is a simplistic descriptor for how the crystal latticereacts to the spin conversion; the analysis of the distribution

of the ri,j

norm distances is more informative. It is represented in

the form of histograms Hri,j

norm in Figs. 8 and 9, correspond-

ing respectively to the J2=0.6 and J2=0.9 cases discussedabove.

In the case of a spin crossover Fig. 8, the distributionconsists of a large single peak whatever the temperature is,which corresponds to a gradual change in the lattice spacingin the crystal. This is related to a homogeneous transforma-tion mechanism without crystallographic phase separation.One can notice that the distribution of lattice spacing is quitedispersed around the mean value. This is typical for agradual spin conversion in our model and is attributed toweak intermolecular couplings. Starting from the LS stateand increasing the temperature, the width of the distributionfirst increases as the spin conversion proceeds, as a conse-quence of increased lattice distortions. In the late steps of theconversion, the distribution finally becomes narrower. Distri-butions 1, 3, and 4 Fig. 9 present a radically different be-havior for the first-order transition case. At low temperaturehistogram 1, a sharp single peak centered around the latticespacing of the LS undistorted phase rnorm=0 is observed.As temperature increases in the bistability area hysteresisloop, a second peak appears centered around the HS latticeparameter rnorm=1. The double peak evidences domain for-mation with phase coexistence the metastable and the ther-modynamically stable phases separated by domain walls. Atthe same time, both distributions widen, as a consequence of

spin

lattice

FIG. 7. Top: Temperature dependences of nHS squares andrnorm circles in the anharmonic model with J0=999, J1=0, and=8, computed as in Fig. 6. J2 has been chosen below the transi-tion line defined in Fig. 5 such as J2=0.6. Bottom: Snapshots of theconfiguration of the system as in Fig. 6.

FIG. 8. Distributions of intersite distance Hri,j

norm correspond-

ing to a gradual spin conversion, relative to snapshots 1, 2, and 3in Fig. 7. The inset gives the mean value rnorm and width w ofeach distribution.

NICOLAZZI, PILLET, AND LECOMTE PHYSICAL REVIEW B 78, 174401 2008

174401-8

structural relaxation, at the domain walls principally. In thefinal step of the transition, the rnorm=1 peak sharpens andlong-range order is restored in the HS phase. It is noteworthythat distributions 2 and 3 correspond to a quite similar HSfraction but exhibit radically different behaviors. This differ-ence can be related to the observations reported in the litera-ture from x-ray- and neutron-diffraction measurements. Inthe case of highly cooperative materials, the diffraction pat-tern usually consists of a superposition of the diffraction pat-terns of the LS and HS phases, with Bragg-peak splitting.This has been attributed to phase coexistence in the bistabil-ity region.3841 On the contrary, for a homogeneous transi-tion, a continuous peak displacement occurs.60

It is noteworthy that owing to the asymmetric anhar-monic form of the intermolecular LJ potentials, the modelspecifically accounts for the thermal expansion of the crystallattice as illustrated in Fig. 10. For this simulation, the spinconfiguration for all the sites of the system has been fixed toHS or LS bypassing the spin Metropolis update in the algo-rithm. Only the lattice variables were allowed to change.The mean lattice spacing has been evaluated as a function oftemperature. As temperature increases, thermal fluctuationslead to a higher probability for r to be greater than rHS

0 andrLS

0 in the HS and LS states, respectively. This directly resultsfrom the repulsive contribution of the LJ potential beingmuch stronger than the attractive part see Fig. 2. In agree-

ment with the chosen shape of the LJ potentials, the derivedlinear thermal-expansion coefficient is higher in the HS statethan in LS.

C. Condition for the presence of hysteresis

In this section, we analyze the dependence of the equilib-rium temperature Tequ and hysteresis width T as a functionof the parameters J0, J1, and J2 of the anharmonic model. Forthe standard Ising-like model with z neighbors, it has beenshown from a mean-field treatment that the condition forhysteresis bistability is simply given by

zJ2

ln g. 18

When J2 is large compared to , the change in nHS is dis-continuous and the LSHS conversion corresponds to afirst-order transition. Conversely, when J2 is small comparedto , the change in nHS is gradual, as is the case for a spincrossover. The steepness of the change in nHS at Tequ dependson the amplitude of J2, the transition becoming more abruptfor large J2.

In the case of the anharmonic model, we have noted inSec. II that the spin dependence of the LJ potentials behavesas an additional phonon contribution h to the temperature-dependent field eff. This suggests a nontrivial dependenceof Tequ and T on the parameters J0, J1, and J2. In the fol-lowing, we restrict the discussion to the more informativecase of highly cooperative system undergoing a first-ordertransition. Tequ and T are defined as follows:

FIG. 10. Temperature dependences of the mean lattice spacingr for the pure LS and HS phases i.e., in the absence of spinconversion. The linear thermal-expansion coefficients have beenfitted by straight lines. The HS and LS corresponding values aregiven in the inset.

FIG. 9. Distributions of intersite distance Hri,j

norm in the case of

first-order transition, relative to snapshots 1, 3, and 4 in Fig. 6 in thewarming process. The inset gives the mean value rnorm and widthw of each distribution.

TWO-VARIABLE ANHARMONIC MODEL FOR SPIN- PHYSICAL REVIEW B 78, 174401 2008

174401-9

Tequ =T + T

2,

T = T T, 19

where T and T are spin transition temperatures in thewarming and cooling modes, respectively, determined suchas d2nHST /d2T=0. The dependences of the temperatureTequ and thermal hysteresis width T on J0, J1, and J2 arereported in Figs. 11 and 12, respectively. J0, which is themean LJ potential depth, does not seem to play a key role inour model. Tequ and T are barely affected by large changesin J0 except for very small values of J0 not shown here. J2drives the abruptness of the transition and the opening of thehysteresis loop. J1, which corresponds to the difference inHS and LS potential depth and therefore to the stiffness dif-ference between the LS and HS phases, controls mainly theposition of the equilibrium temperature. This latter is there-fore independent of the constant AHL. In the mean-field ther-modynamic model of Slichter and Drickamer,16 there is a ccritical value of the phenomenological interaction parameter beyond which a thermal hysteresis can be observed. Asshown in Fig. 12, such a critical value of the interactionparameter J2 also exists in the anharmonic Ising-like model.

IV. DISCUSSION AND CONCLUSION

We have presented an Ising-like model for spin transitionsystems, taking into consideration the long-range elastic in-teractions together with the vibronic degeneracy of the HSand LS molecules. The model allows for molecular displace-ments governed by intermolecular Lennard-Jones potentialson a deformable lattice. By means of Monte Carlo simula-tions, we have analyzed the static behavior of the model; thedynamic properties will be published in a forthcoming paper.The main advantage of our approach with respect to otherIsing-like schemes is the use of two coupled degrees of free-dom, namely, the spin component and the lattice spacing.

Within our model, the interaction between a pair of neigh-boring molecules in the crystal is dependent not only on theirspin states but also on their separation distance, introducingexplicitly inhomogeneities and a distribution of interactionconstants. One key element of the model is the spin depen-dence of the LJ potentials; a broader potential has been at-tached to the HS state, corresponding to a lower bulk modu-lus. The HS-LS difference in potential depth controls thethermal transition temperature and behaves as an additionalentropy term, attributed phenomenologically to a lattice pho-non contribution. The intermolecular coupling constant J2drives the first-order character of the transition and width ofthe hysteresis loop for strong interactions. The anharmonicityasymmetry of the lattice, provided by the chosen LJ formof the pairwise interaction potentials, brings important as-pects. For instance, thermal contraction effects, as observedand characterized experimentally, emerge directly. From a

FIG. 12. Evolutions of the thermal hysteresis width T as afunction of J0, J1, and J2 for =4.

(a) (b)

FIG. 11. Variations in equilibrium temperature Tequ for =8: a with J1 J0=1000 and J2=0.9 and b with J2 J0=999 and J1=0.Insets in each plot correspond to thermal hysteresis loops when a J1 and b J2 vary. J2=0.65 corresponds approximately to the criticalpoint at which the hysteresis loop opens.

NICOLAZZI, PILLET, AND LECOMTE PHYSICAL REVIEW B 78, 174401 2008

174401-10

1D analytical treatment of our anharmonic model, it will beshown in a forthcoming paper that the anharmonic contribu-tion by comparison with harmonic pairwise potentials low-ers the transition temperature, thus stabilizing the HS phase.It corresponds formally to additional degeneracies of latticevibration origin in the effective field. Similar conclusionshave been drawn in the case of the 1D atom-phonon cou-pling model. It has been shown that the anharmonic correc-tion to the harmonic intersite coupling lowers also the tran-sition temperature and drives the abruptness of thetransition.33

We may anticipate that our model would allow an inter-pretation of the first-order transition mechanism throughlike-spin domain nucleation and growth, as clearly evidencedby several x-ray- and neutron-diffraction measurements. As amatter of fact, for strong intermolecular interactions, the dis-tribution of lattice spacing presents a double sharp structure,centered on the HS and LS equilibrium distances. This isdirect evidence for phase separation in the crystallographicsense. On the contrary, in the spin-crossover case weak in-termolecular coupling, the distribution of lattice spacingscales with the fraction of HS molecules. In this formalism,the notion of like-spin domains emerges from the lattice vari-able. This behavior contrasts with other elastic models with

only long-range interactions and mean-field behavior,37which do not exhibit cluster growth and phase separation. Itis well known that the usual Ising model with short-rangeinteractions also exhibits cluster growth, but in this case,only spin variables are concerned: it does not correspond to acrystallographic phase separation such as the one evidencedin the present study.

Our model may be easily adapted to treat dilution ordoping effects, that is, the case for which the crystal latticecontains both SC and non-SC molecules with different latticespacings. Especially the influence of the strain induced bythe doping elements on the dynamic behavior would behighly relevant. This possibility will be explored in a futurework.

ACKNOWLEDGMENTS

This work was supported by the European Network ofExcellence MAGMANet Grant No. FP6-515767-2, theUniversit Henri Poincar, and the CNRS. The authorswould like to thank K. Boukheddaden for helpful discussion.Part of the numerical calculations was performed at the com-puting center of the Institut Jean Barriol, Nancy, which isacknowledged.

*Corresponding author; sebastien.pillet@lcm3b.uhp-nancy.fr1 P. Gtlich, A. Hauser, and H. Spiering, Angew. Chem., Int. Ed.

Engl. 33, 2024 1994.2 Topics in Current Chemistry, edited by P. Gtlich and H. A.

Goodwin Springer-Verlag, Berlin, 2004, Vols. 233235.3 P. Guionneau, M. Marchivie, G. Bravic, J.-F. Ltard, and D.

Chasseau, Topics in Current Chemistry Springer-Verlag, Berlin,2004, Vol. 233, p. 97.

4 S. Decurtins, P. Gtlich, C. P. Khler, and H. Spiering, Chem.Phys. Lett. 105, 1 1984.

5 A. Hauser, Chem. Phys. Lett. 124, 543 1986.6 A. Hauser, Comments Inorg. Chem. 17, 17 1995.7 E. Freysz, S. Montant, S. Ltard, and J.-F. Ltard, Chem. Phys.

Lett. 394, 318 2004.8 S. Bonhommeau, G. Molnar, A. Galet, A. Zwick, J. A. Real, J. J.

McGarvey, and A. Bousseksou, Angew. Chem., Int. Ed. 44,4069 2005.

9 N. Willenbacher and H. Spiering, J. Phys. C 21, 1423 1988.10 H. Spiering and N. Willenbacher, J. Phys.: Condens. Matter 1,

10089 1989.11 H. Spiering, Topics in Current Chemistry Springer-Verlag, Ber-

lin, 2004, Vol. 235, p. 171.12 A. Hauser, Chem. Phys. Lett. 192, 65 1992.13 A. Hauser, J. Jeftic, H. Romstedt, R. Hinek, and H. Spiering,

Coord. Chem. Rev. 190-192, 471 1999.14 J. Kroeber, J. P. Audiere, R. Claude, E. Codjovi, O. Kahn, J. G.

Haasnoot, F. Groliere, C. Jay, A. Bousseksou, J. Linares, F.Varret, and A. G. Vassal, Chem. Mater. 6, 1404 1994.

15 J. A. Real, A. B. Gaspar, V. Niel, and M. C. Munoz, Coord.Chem. Rev. 236, 121 2003.

16 C. T. Slichter and H. G. Drickamer, J. Chem. Phys. 56, 2142

1972.17 H. Spiering, K. Boukheddaden, J. Linares, and F. Varret, Phys.

Rev. B 70, 184106 2004.18 J. Wajnflasz, Phys. Status Solidi 40, 537 1970.19 J. Wajnflasz and R. Pick, J. Phys. Colloq. 32, C1-91 1971.20 A. Bousseksou, H. Constant-Machado, and F. Varret, J. Phys. I

5, 747 1995.21 K. Boukheddaden, I. Shteto, B. Ho, and F. Varret, Phys. Rev. B

62, 14796 2000.22 K. Boukheddaden, I. Shteto, B. Ho, and F. Varret, Phys. Rev. B

62, 14806 2000.23 F. Varret, S. A. Salunke, K. Boukheddaden, A. Bousseksou, E.

Codjovi, C. Enachescu, and J. Linares, C. R. Chim. 6, 3852003.

24 M. Nishino, S. Miyashita, and K. Boukheddaden, J. Chem. Phys.118, 4594 2003.

25 K. Boukheddaden, J. Linares, E. Codjovi, F. Varret, V. Niel, andJ. A. Real, J. Appl. Phys. 93, 7103 2003.

26 K. Boukheddaden, J. Linares, R. Tanasa, and C. Chong, J. Phys.:Condens. Matter 19, 106201 2007.

27 M. Nishino, K. Boukheddaden, S. Miyashita, and F. Varret, Phys.Rev. B 68, 224402 2003.

28 H. Bolvin and O. Kahn, Chem. Phys. 192, 295 1995.29 K. Boukheddaden, J. Linares, H. Spiering, and F. Varret, Eur.

Phys. J. B 15, 317 2000.30 J. A. Nasser, Eur. Phys. J. B 21, 3 2001.31 J. A. Nasser, K. Boukheddaden, and J. Linares, Eur. Phys. J. B

39, 219 2004.32 K. Boukheddaden, S. Miyashita, and M. Nishino, Phys. Rev. B

75, 094112 2007.33 K. Boukheddaden, Prog. Theor. Phys. 112, 205 2004.

TWO-VARIABLE ANHARMONIC MODEL FOR SPIN- PHYSICAL REVIEW B 78, 174401 2008

174401-11

34 M. Nishino, K. Boukheddaden, Y. Konishi, and S. Miyashita,Phys. Rev. Lett. 98, 247203 2007.

35 K. Boukheddaden, M. Nishino, and S. Miyashita., Phys. Rev.Lett. 100, 177206 2008.

36 Y. Konishi, H. Tokoro, M. Nishino, and S. Miyashita, Phys. Rev.Lett. 100, 067206 2008.

37 S. Miyashita, Y. Konishi, M. Nishino, H. Tokoro, and P. A.Rikvold, Phys. Rev. B 77, 014105 2008.

38 S. Pillet, V. Legrand, M. Souhassou, and C. Lecomte, Phys. Rev.B 74, 140101R 2006.

39 N. Huby, L. Gurin, E. Collet, L. Toupet, J. C. Ameline, H.Cailleau, T. Roisnel, T. Tayagaki, and K. Tanaka, Phys. Rev. B69, 020101R 2004.

40 S. Pillet, J. Hubsch, and C. Lecomte, Eur. Phys. J. B 38, 5412004.

41 K. Ichiyanagi, J. Hebert, L. Toupet, H. Cailleau, P. Guionneau,J.-F. Ltard, and E. Collet, Phys. Rev. B 73, 060408R 2006.

42 Z. Xi, B. Chakraborty, K. W. Jacobsen, and J. K. Norskov, J.Phys.: Condens. Matter 4, 7191 1992.

43 M. Nielsen, L. Miao, J. H. Ipsen, O. G. Mouritsen, and M. J.Zuckermann, Phys. Rev. E 54, 6889 1996.

44 D. J. Bergman and B. I. Herpin, Phys. Rev. B 13, 2145 1976.45 S. Doniach, J. Chem. Phys. 68, 4912 1978.46 E. H. Boubcheur and H. T. Diep, J. Appl. Phys. 85, 6085 1999.47 P. Massimino and H. T. Diep, J. Appl. Phys. 87, 7043 2000.48 E. H. Boubcheur, P. Massimino, and H. T. Diep, J. Magn. Magn.

Mater. 223, 163 2001.49 G. Molnar, V. Niel, A. B. Gaspar, J. A. Real, A. Zwick, A.

Bousseksou, and J. J. McGarvey, J. Phys. Chem. B 106, 97012002.

50 R. Zimmermann and E. Konig, J. Phys. Chem. Solids 38, 7791977.

51 V. Legrand, S. Pillet, C. Carbonera, M. Souhassou, J.-F. Ltard,P. Guionneau, and C. Lecomte, Eur. J. Inorg. Chem. 2007, 5693.

52 A. J. Stone, The Theory of Intermolecular Forces Clarendon,Oxford, 1996.

53 T. F. Middleton, J. Hernndez-Rojas, P. N. Mortenson, and D. J.Wales, Phys. Rev. B 64, 184201 2001.

54 K. Nakano, S. Kawata, K. Yoneda, A. Fuyuhiro, T. Yagi, S.Nasu, S. Moritomo, and S. Kaizaki, Chem. Commun. 2004,2892.

55 J. Linares, H. Spiering, and F. Varret, Eur. Phys. J. B 10, 2711999.

56 R. Boca, M. Boca, L. Dlhan, K. Falk, H. Fuess, W. Haase, R.Jarosciak, B. Papankova, F. Renz, M. Vrbova, and R. Werner,Inorg. Chem. 40, 3025 2001.

57 N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H.Teller, and E. Teller, J. Chem. Phys. 21, 1087 1953.

58 L. Onsager, Phys. Rev. 65, 117 1944.59 Finite-size and boundary effects remove the singularity parts of

the free-energy density f Ising. This in turns move away the criti-cal temperature from the critical point which exists rigorously inthe thermodynamic limit, when N.

60 A. Goujon, B. Gillon, A. Debede, A. Cousson, A. Gukasov, J.Jeftic, G. J. McIntyre, and F. Varret, Phys. Rev. B 73, 1044132006.

NICOLAZZI, PILLET, AND LECOMTE PHYSICAL REVIEW B 78, 174401 2008

174401-12