PHYSICAL REVIEW B VOLUME 53, NUMBER 14 1 APRIL 1996-II Surface-induced ordering and disordering in face-centered-cubic alloys: A Monte Carlo study W. Schweika Institut fu¨r Festko¨perforschung, Forschungszentrum Ju¨lich, D-52425 Ju¨lich, Germany D. P. Landau Center for Simulational Physics, The University of Georgia, Athens, Georgia 30602; Institut fu¨r Physik, Johannes Gutenberg-Universita¨t Mainz, Staudinger Weg 7, D-55099 Mainz, Germany; and Ho¨chstleistungsrechenzentrum c/o Forschungszentrum Ju¨lich, D-52425 Ju¨lich, Germany K. Binder Institut fu¨r Physik, Johannes Gutenberg-Universita¨t Mainz, Staudinger Weg 7, D-55099 Mainz, Germany Received 19 September 1995 Using extensive Monte Carlo simulations we have studied phase transitions in a fcc model with antiferro- magnetic nearest-neighbor couplings J in the presence of different free surfaces which lead either to surface- induced order or to surface-induced disorder. Our model is a prototype for CuAu-type ordering alloys and shows a strong first-order bulk transition at a temperature kTcb / J 1.738 005(50). For free 100 surfaces, we find a continuous surface transition at a temperature Tcs Tcb exhibiting critical exponents of the two- dimensional Ising model. Surface-induced ordering occurs as the temperature approaches Tcb and the surface excess order and surface excess energy diverges logarithmically. For a free 111 surface, the surface order vanishes continuously at Tcb accompanied by surface-induced disorder SID . In addition to a logarithmic divergence of the excess quantities of order and energy, we find further critical exponents which confirm the actual theory of SID and critical wetting and which can be understood in terms of rough interfaces. For both cases of free surfaces, the asymptotic behavior of the squared interfacial width shows the expected logarithmic divergence. I. INTRODUCTION low-electron electron-diffraction LEED study of a free 100 surface of a Cu3Au alloy,13,15 and the critical surface Phase transitions in real crystals are affected by the pres- exponent 1 was later determined by spin-polarized LEED.16 ence of the surfaces as well as by possible interfaces between Monte Carlo MC simulations14 of an Ising model of this ordered domains called antiphase boundaries APB's . A alloy also indicated a possible continuous decrease of the number of excellent theoretical studies1­7 have led to de- surface order as the bulk transition was approached. tailed predictions for the properties of semi-infinite systems From experiments on the same alloy using evanescent and have been accompanied by some remarkable experimen- x-ray scattering in grazing incidence evidence was found for tal work in recent years. A review of the related experimental the increasing thickness of the disordered surface layer as the investigations has been given by Dosch.8 There are different transition temperature is approached.17 The related phenom- scenarios for how the phase transitions at the surface and in enon of surface melting was observed in lead by ion the near-surface region depend on the order of the bulk phase scattering18 and confirmed by LEED.19 In the case of SID transition. there is no clear quantitative experimental evidence up to For second-order bulk transitions we now have a rather now for the predicted logarithmic divergence of the thickness good picture of the range of possible surface behavior.1,2 In of a wetting layer of the disordered phase at the surface. some situations the surface undergoes a transition at the bulk However, such behavior has been observed by transmission transition temperature, but the surface critical behavior is electron microscopy for the analogous situation in the bulk described by critical exponents which differ from both the of a Cu-Pd 17% alloy where the width of the antiphase two-dimensional 2D as well as bulk 3D values. In other boundaries diverges logarithmically as the transition tem- cases, when the exchange in the surface layer exceeds that perature is approached from below.20 Although discontinuous within the bulk by a sufficient amount, the surface may order order-disorder transitions are quite common in alloy systems, above the bulk transition, exhibiting 2D exponents. most experimental examples consider simple structures The situation in which the bulk undergoes a first-order showing transitions from a cubic ordered (L12 in the ex- transition has been considered by Lipowsky7,9­12 who amples above to a cubic disordered phase thus avoiding the showed that surface-induced disorder SID and surface- complications due to the possible strain between ordered induced order SIO may then be associated with the bulk variants of lower symmetry. However, as has been shown in transition. One of the most notable results is the continuous a previous Monte Carlo study, the fcc alloys having a Cu- decrease of the order parameter in the surface layer as the Au-type order (L10) and tetragonal symmetry are not only bulk transition temperature is approached. Such a behavior candidates for the observervation of SID but also for SIO had already been found in 1973 by Sundaram et al. in a which may occur because a free 100 surface reduces the 0163-1829/96/53 14 /8937 19 /$10.00 53 8937 © 1996 The American Physical Society 8938 W. SCHWEIKA, D. P. LANDAU, AND K. BINDER 53 frustration effects.21 It is the aim of the present work to con- tribute to the understanding of both SID and SIO by Monte Carlo simulation of a suitable model. In this paper we consider a simple model for an AB bi- nary alloy which can, of course, be reinterpreted as an Ising magnet, or as a lattice-gas model for a fluid or alloy spins Si 1 at lattice site i corresponding to the site being occu- pied or empty, or containing an A atom or B atom, respec- tively . For our Monte Carlo simulation, a thick-film geom- etry is used in which there are two parallel free surfaces, and the nature of the phase transitions in these surfaces forms the subject of this study. By ``thick film'' we mean that the two surfaces can be considered as essentially noninteracting and that our study is meant to address the behavior of semi- infinite bulk systems with free surfaces. Preliminary results have been presented21 for the case where both nearest- neighbor and next-nearest-neighbor coupling are present, but only a single orientation of the surface 100 was considered Fig. 1 a . We now present results with only nearest- neighbor coupling and two different orientations of free sur- faces. In the following section we describe the theoretical back- ground on SID, to provide the framework in which our study can be interpreted, and motivate the choice of quantities that will be studied. In Sec. III we describe the model and meth- ods used, and in Secs. IV and V we present results for 100 and 111 surfaces Fig. 1 b , respectively. We conclude in Sec. VI, while the ground-state behavior and the precise es- timation of the bulk transition temperature are discussed in the appendixes. II. THEORETICAL BACKGROUND A. Surface-induced disordering In this section, we recall the basic theoretical predictions FIG. 1. Section of the ordered lattice of the AB alloy in the on surface-induced disordering SID considering only the L10 CuAuI structure. Frustrated nearest-neighbor interactions are shown by dashed lines. a 100 surface oriented such that there is standard situation which is equivalent to critical wetting3,7,22 (2 2) order in the surface layer; b 111 surface shaded . Note and disregarding other situations which correspond to sur- that in this case a (2 1) order in the surface results, where one- face multicritical points.7,9­12 third of the bonds in the surface layer are still frustrated thick The basic phenomenon of interest is the continuous decay broken lines . Unfrustrated nearest-neighbor interactions are not of the local order parameter 1 at the surface layer as the drawn in order to avoid overcrowding of these pictures. first-order transition in the bulk located at Tcb) is ap- proached from below,9 involving a power law with an expo- Since nent s 0 in d 3 dimensions, the divergence of l¯ as 1 t 0 is only logarithmic.3,7,9­12 As a result, the disordered layer increasingly ``screens'' the effective field acting on the 1 t 1, t 1 T/Tcb 0. 1 surface layer due to the still well-ordered bulk, and thus Eq. 1 becomes plausible. The interfacial thickness or roughness This continuous behavior at a discontinuous bulk transition is expected to diverge as well but more slowly as follows occurs because a layer of the disordered phase gradually in- from22 trudes at the surface. For the three-dimensional systems gov- erned by short-range forces it has been predicted that the thickness l¯ of this layer ( d being the correlation length of / d ln 1/t 1/2, t 0, order-parameter fluctuations in the disordered bulk phase t , should diverge as 0 ln . 3 As for more standard critical phenomena, one can define a l¯ dln 1/t , t 0, 2a number of divergent response functions to ``fields'' H con- jugate to the order parameter in the bulk and H where one may define a critical exponent 1 conjugate to 1). Note, that in this case H and H1 are staggered fields. In the theory of surface effects on bulk critical phenomena1,2 l¯ dt s, s 0 ln , t 0. 2b one distinguishes between the susceptibilities 53 SURFACE-INDUCED ORDERING AND DISORDERING IN . . . 8939 wX 2 w 1/2 1/2 w . 11b 1,1 1 H t 1,1, t 0, 4a 1 t In the last step we have redefined the scaling functions to express the singularities in terms of instead of . Now a comparison of the free-energy functionals of SID 1 1 H t 1, t 0, 4b t and critical wetting shows3,7,22 that these problems are equivalent to each other if one approaches wetting criticality and along a special path where ) in the ( , ) plane. Since it turns out that the singularities along these paths are the same as along the path ( 0, variable , we consider for s s H t s, t 0. 4c t simplicity singularities along the latter path. From Eqs. 8a ­ 11b it is easy to derive the critical be- In Eq. 4c we have used the surface excess order s which havior of the various response functions for this wetting in turn is defined as the derivative of the surface excess free problem. We recall that for critical wetting temperaturelike energy variables and H1 scale in the same way, since one can cross the wetting transition line H s Fs T,H,H1 / H T,H , 5a 1 H1c(T) either by variation of 1 T or of H1 . Therefore, taking a derivative of Eq. 8a with while 1 can be obtained as a derivative with respect to the respect to we get local field H1 , w w / w, 12 1 1 s 1 1/ w 0 1 s 1 Fs T,H,H1 / H1 T,H . 5b where we have indicated how the power law for 1 as func- Finally we mention that the correlation function for order- tion of can be inferred by requiring that the corresponding parameter fluctuations in the surface layer can be used to scaling function 1 behaves as a power law for small argu- define a correlation length , ments such that the dependence on cancels. From Eq. 12 we can read off the value of the exponent 1 1 0 1 1 2 exp / , , 6 2 d 1 2/ w where is a coordinate in the surface layer, and w 1 1 s / w . 13 d 1 dt , t 0. 7 Using Eqs. 5a and 8a we obtain the critical behavior of In order to discuss the scaling relations between these critical the excess quantities, in particular of s : exponents for SID, and the theoretical predictions existing w for them, we recall that SID can be interpreted as a critical w / w s 2 s 14 wetting phenomenon.3,7,22 The scaling relations for critical wetting with short-range forces follow from a scaling analy- determining the critical exponent s sis for Fs and for , in terms of scaling functions F and w s 2 w / w d 3 / d 1 . 15 X . We use here a superscript ``w'' in order to avoid confu- s sion between the exponents for critical wetting and the ex- Because of the vanishing s in d 3 one expects a loga- ponents for SID: rithmic divergence of the surface excess order parameter s , as well as of other excess quantities like the surface F w s 2 s F w , 8a excess energy Es and the thickness of the wetting layer l¯ Eq. 2b . wX w , 8b Similarly one obtains the surface excess susceptibility, the layer susceptibilities and their exponents, respectively, where and are temperaturelike and fieldlike variables for the wetting transition, w w s , w and w are the associated criti- / w , cal exponents. For d d* 3, the upper critical dimension- s 2 2 s ality for wetting transitions, one has the standard hyperscal- 2 w 4 ing relation s s 2 w d 1 , 1 in d 3 16a 2 ws d 1 w 9 w / w for surface excess quantities. Now one can show3 that for 1 1 1 s , critical wetting there is in fact a single independent expo- w nent, since w also can be related to w 1 3 d 2/ w s 1 1 w d 1 , 16b w d 1 w/2. 10 w/ w In d 3 dimensions the resulting equations for F , s and can 1,1 s hence be written as follows: 2 d 1 4/ w w/ w F 1,1 s s 2 wF 2 w f 1/2 w 11a d 1 . 16c 8940 W. SCHWEIKA, D. P. LANDAU, AND K. BINDER 53 From Eqs. 13 ­ 16 one can read off the scaling laws kBTcb / 4 l2 , 23 1 1 1, 17a where Tcb is the transition temperature of the bulk, is the interfacial stiffness of the interface between the ordered and s s 1, 17b disordered phase coexisting at Tcb , and l is its intrinsic thickness Lipowsky7 identifies l 2d 6 4/ w d). The exponent 1 is predicted to be 1,1 s d 1 2 1 , 17c 1 1/2 /2, 0 1/2, 24a 2 d 1 1,1 2 1 d 1 1 in d 3 . 17d 1 2 /2, 1/2 2, 24b Some of the above relations can be obtained very directly, of course, if we simply replace by the temperature distance 1 1, 2. 24c t and by H1 in Eqs. 11a and 11b , to find Note that because of Eq. 20 this prediction is simply equivalent to those of Refs. 22­27 for the exponent w dis- Fs tf t 1/2 wH1 , 18a cussed in the theory of critical wetting. Results for other exponents simply follow from the scaling relations written t 1/2 t 1/2 wH1 , 18b down above. It is, however, of interest to show the relation for the logarithmic relations for the thickness of the surface remember that we have restricted ourselves to d 3 here . From Eqs. 8b , 18b one sees that the exponent near layer l¯ and the width of the rough interface7,22 defined in Eq. 7 for SID is l¯/ d 1/2 1 2 ln 1/t , 0 1/2, 25 w/ w 2/ d 1 1/2 in d 3 . 19 From Eq. 18a we can immediately recognize the scaling / d ln 1/t 1/2, 0 1/2. 26 structure of 1 by taking a derivative with respect to H1 , Here we have restricted ourselves to the regime 1/2, which is of most interest. Note that while a direct determi- 1 1 1 t 1/2 wH1 , nation of from Eq. 23 is difficult, since it is hard to calculate , one could infer its value indirectly if Eqs. 25 1 1 1/2 w, d 3, 20 and 26 are used. However, in computer simulation studies and a further derivative yields of wetting in the Ising model29­31 it was already not possible 1,1 , the singular part of which scales as to verify the predictions22­27 for w that correspond to Eqs. 24a ­ 24c . The interpretation of these findings is uncertain: It has been proposed that one must be extremely close to the 1,1 t 1 1,1 t 1/2 wH1 , transition to see the asymptotic behavior,26 or that the wet- ting transition is weakly of first order.27 Boulter and Parry 1,1 1 1/ w, d 3. 21 have recently proposed another explanation.28 This will be These results contain the mean-field theory of SID,7,9­12 as a further discussed later. It is unclear to us whether these latter special case: for the mean-field theory of wetting w 1, and predictions should apply to the present problems as well. hence in d 3) It should be stressed that Eqs. 23 ­ 26 can apply only if w the interface between the ordered and disordered phase is 1 1/2, 1 1/2, 1,1 0, w 2, s 0, rough. Lipowsky7 supposed that Eq. 22a is valid if one is 22a below the roughening transition temperature of this interface. while the following critical exponents are independent of In view of Monte Carlo studies and of molecular field calcu- w) universal and exact in d 3: lations for semi-infinite Potts model33 on wetting near the roughening temperature of the Ising model,32 we would like s 0, s 1, 0, 1/2, s 1. 22b to mention the possibility that for nonrough interfaces SID Note that Eq. 18a can also be interpreted as might be replaced by a sequence of layering transitions. An F effective exponent eff 12 can then be defined only in the s t2 s f (t 1H1), where the exponents s and 1 have the standard meaning as in the usual theory of surface critical sense that the sequence of steps in a log-log plot is approxi- phenomena, see Ref. 1, and 2 mated by a straight line.7 s (d 1) is satisfied with s 1, 1/2 in d 3. We also have 2 s 1 1 and 2 w B. Surface-induced ordering s (2 s )/ w, of course, and s 2 . Since d 3 is the upper critical dimension for critical wetting, all While SID and SIO are commonly discussed in the litera- scaling laws Eqs. 17a ­ 17d are indeed satisfied with these ture to be just equivalent phenomena, one has to note one Landau theory exponents, Eq. 22 . important difference, namely that for SIO the surface neces- Now one knows that fluctuation corrections for critical sarily orders at a different and higher temperature than the wetting are possibly very important, and the renormalization- bulk in order to wet the bulk phase. An estimate for the group theories22­27 predict that w depends on a nonuniver- surface transition temperature can be obtained from a com- sal parameter , parison to the two-dimensional model. Hence, one expects 53 SURFACE-INDUCED ORDERING AND DISORDERING IN . . . 8941 T 2D cb Tcs Tc . One important consequence is that the pre- 1 dictions for the surface properties, and also the scaling rela- x,n 4 m1 m2 m3 m4 tions discussed above, can no longer be valid. Instead, at the surface transition one expects a behavior which falls in the 1 universality class of the corresponding purely 2D system. y,n 4 m1 m2 m3 m4 The exponent describing the singularity of , (T Tcs) , is simply the bulk 2D exponent (d 2), 1 and the exponent 1 describing the order parameter 1 z,n (T T 4 m1 m2 m3 m4 . 29 cs) 1 is the bulk 2D exponent (d 2). Note that at Tcs the susceptibilities 11 , 1 , and s all have the same Similarly, one can define a total order parameter for the singularity T Tcs (d 2), (d 2) being the 2D bulk whole film and a bulk order parameter for the inner part of susceptibility exponent . The scaling relations for the surface the film. However, in most cases, the anisotropy of the sys- exponents at Tcs thus are just standard bulk scaling laws. tems considered is sufficiently strong to single out only one For T Tcb the logarithmic law for the excess quantities, component of order, while the other two components, as well e.g., those of the order, the energy, and the interface position, as their fluctuations, are negligible. This holds throughout for should still be valid. Surface properties are expected to show the case of SID at the free 111 surface. In other words, in an ``extraordinary'' transition1,2 singularities in the second the ordered phase the symmetry is broken, and only one type derivatives, kink in 1). of domain is considered which in the simulation is created by preparation of the initial state . Since the transition is of III. MODEL AND SIMULATION TECHNIQUE first order in the bulk, the system does not ``jump over'' to the orderings corresponding to the other types of domains, at We study the Ising Hamiltonian on the face-centered- least not for the large lattice sizes studied. In these cases, it is cubic lattice in a L L D geometry, applying periodic justifiable to confine ourselves to a two sublattice order pa- boundary conditions in x and y directions. The two free rameter L L surfaces are oriented to be either 100 or 111 faces. The Hamiltonian used is 1 y,n 2 m1 m2 , 30 H J SiSj , Si 1 , 27 which is simply a single layer order parameter and a scalar i,j quantity. Both types of order parameters, and their related where the sum is over all nearest-neighbor pairs (i, j) with susceptibilities, were calculated in the simulations. The sur- antiferromagnetic nearest-neighbor coupling, i.e., J 0. Note face excess order parameter is obtained by that in the present case, the atoms in the surface layers do not D/2 have any modified nearest-neighbor coupling, but they do, of s n b . 31 course, see fewer neighbors than do those atoms in the bulk. n 1 We used a Metropolis, single spin-flip method with preferen- The bi- layer susceptibilities tial layer sampling which was determined by the nature of n and n,n were obtained from the fluctuation relations34 the order-parameter profile. We implemented an efficient, vectorized single spin-flip algorithm on a CRAY-YMP com- puter. Since our aim was to gain an overview of the system n,n L2 n* n n * n /kBT, 32 behavior for a wide range of temperatures, excessively large values for L were avoided. Typically L varied from 32 to n L2D n* tot n * tot /kBT, 33 128, and thicknesses studied varied from D 40 to D 200 while the surface excess susceptibility to ensure that the two surfaces were independent. s follows as In principle, a complete description of the CuAu-type or- D/2 der requires a three-dimensional order parameter s n b . 34 ( n 1 x , y , z) which refers to the three possible orienta- tions of ordered domains, which are alternate layering of Here, pure Cu and pure Au planes in one of the three 100 direc- b denotes the bi- layer susceptibility in the bulk. The surface layer susceptibilities tions. The components of the order parameter are necessarily 1 and the excess susceptibili- ties based on four sublattices and can be defined for bilayers s were also calculated according to their definitions as derivatives in Eqs. 4b and 4c . As discussed in Appendix only. For example, the first component x,n the order param- C, there are no fluctuation relations analogous to Eqs. 32 eter of the bilayer n is defined as and 33 for the excess quantities such as s . Furthermore, we considered the layer energies En nor- 1 2 malized per spin , the bulk and total energies Eb and Etot , x,n L2 Sjeik­rj, where k the surface excess energy Es , the layer specific heat Cn , and j bilayer n a 1,0,0 . the surface excess specific heat C 28 s : D/2 Using the sublattice magnetizations m1 , m2 , m3 , and m4 E E one obtains s n Eb , 35 n 8942 W. SCHWEIKA, D. P. LANDAU, AND K. BINDER 53 Cn L2D EnEtot En Etot /T2, 36 D/2 Cs Cn Cb . 37 n For locating the transition of the first layer in case of SIO, we used the reduced cumulant U 4 2 1,L 1 1 / 3 1 2 . 38 In order to fit the layer profiles we used the following formula: 1 n max 1 exp 2 n n ... 1, 39 where n is the interface position and the interface width. Note, that the inverse of the gradient of n equals twice that of . This formula is analogous to the tanh function of the interfacial profile between two ordered domains at a second- order bulk transition.39 FIG. 2. Temperature dependence of the surface layer order pa- The typical run length varied from an average of 104 rameter 1 broken curve and the bulk order parameter b full curve for the case of our 100 surface. Two linear dimensions Monte Carlo steps/site to an average of 105 Monte Carlo (L 80, 128 are included to show the onset of size effects at steps/site near Tcb for the larger choices of L. We speak of T an ``average'' number of spin flips because the different lay- cs . ers were not sampled equally. For studies of wetting, layer- ing, and surface critical phenomena,29­32,34 we have found it surface at a second-order transition is primarily discussed: useful to consider sites near the free surfaces more often for then 1(T) has no kink, since the singularity of 1 is the a spin flip than sites in the bulk ``preferential surface site same as that of the bulk free energy density, f b selection'' . Here the situation is different; the surface fields 1 T/Tcb 2 b : 1 then has a singularity only in its cur- suppress fluctuations in the layers very close to the surfaces, vature at Tcb .48 At a first-order transition, fb has a kink and the largest and slowest fluctuations do not occur in the singularity, and so does 1.] bulk. In many cases considered they also do not occur at the Furthermore there are obvious finite-size effects at the surface but rather at, or near, the interface which is typically surface ordering transition. The nature of the surface order- a few layers away from the surface . Therefore, we choose ing transition and the critical behavior of our model is ex- either a ``preferential surface site selection'' which decays pected to be same as for a two-dimensional Ising model, exponentially to a constant with n, or one in which the since obviously the (2 2) order of the ordered surface choice is proportional to the gradient plane has this symmetry. Since the surface transition tem- n / n. We also car- ried out multiple, independent runs with different random perature is above the first-order bulk transition temperature number sequences to obtain estimates for the statistical er- Tcb the nature of the surface transition is, of course, indepen- rors. dent of any wetting properties which occur as T Tcb . Previous calculations40­47 of the bulk transition tempera- The manner in which the surface order propagates into the ture T bulk as the temperature is lowered can be seen in Fig. 3 cb for this model do not have the precision required for the present study. Therefore, the transition temperature of the where we present some profiles of the order parameter and bulk alloy was determined by standard thermodynamic inte- layer energy for a thick film of size L L with L 80 and gration of data for large systems with fully periodic bound- thickness 80 D 200. These large thicknesses ensure that ary conditions. We estimate that k the two free surfaces are independent of each other, as can be BTcb / J 1.738 005 0.000 050. A more detailed description of this study is verified very nicely from the pair-correlation function data presented in the Appendix A. shown in Fig. 3 c : correlations parallel and perpendicular to the free surfaces are identical in the bulk. At the interface the IV. RESULTS FOR 100... SURFACES correlations remain anisotropic and do not obey the symme- try properties of the disordered cubic phase. The perpendicu- In Fig. 2 a we show the temperature variation of the lar correlations across the interface determine the wetting surface layer order parameter 1 and the bulk order param- process. The related correlation length should remain finite eter b . This figure clearly shows the large jump in bulk as T approaches the first-order transition at Tcb ; however, order which occurs at Tcb , but it is also obvious from this this correlation length is not simply the correlation length of plot that surface order begins to develop well above the bulk the disordered bulk phase as it is usually assumed in the transition. The surface order increases as the temperature is literature.7 Furthermore, even in the disordered phase the lowered, and it smoothly alters when the bulk finally under- correlations can be anisotropic, and this must be taken into goes a discontinuous transition. At the first-order bulk tran- regard when choosing a ``typical'' bulk correlation length. In sition temperature Tcb one notices a kink in 1(T), as one particular, this is true for the present model as shown in Fig. may expect for an ``extraordinary transition.'' In the 4. The largest correlation length is found for the 100 direc- literature,1,2,48,49 the extraordinary transition of an ordered tions, 100 4.3 in units of a/2 at Tcb. 53 SURFACE-INDUCED ORDERING AND DISORDERING IN . . . 8943 FIG. 4. Anisotropy of the correlations in the disordered bulk phase at Tcb shown for the three directions 100 , 110 , and 111 simulated with full periodic boundary conditions, size: four sc sublattices, each of L3 693). The largest correlation length is seen in 100 directions, 100 4.3 a/2. Tcs the evaluation of the gradient needs a proper analysis of finite-size effects, we found that finite-size effects are not so important for lower T above Tcb at least for the large sys- tems considered . To a good approximation the order- parameter profiles can be fitted to Eq. 39 as shown for example in Fig. 5. Apart from slight deviations close to the surface, this provides an accurate determination of the posi- tion n and width of the interfaces. The results for the thickness of the ordered surface layer l¯ which equals n 1/2 and for the squared width 2 versus the reduced temperature t reveal the theoretically expected asymptotic logarithmic behavior for a rough interface in the limit T Tcb . However, the more detailed predictions given by Eqs. 25 and 26 are in disagreement with the present re- sults. Identifying 100 as obtained in Fig. 4 with d yields 0.16 which is conceivable in view of the results for the interfacial properties in case of SID at the 111 surface dis- cussed below. However, these values and Eq. 26 would result in a slope 2 / lnt 3 which is inconsistent with the data in Fig. 5 c slope 5.2 . In recent work28 Boulter and Parry pointed out that the theory of wetting phenomena should not be formulated only in terms of the interfacial position l¯, but should include the coupling of l¯ to a second length scale l1; this length scale l1 describes the range over which the surface causes deviations from the simple interfa- cial profile, precisely as we see them in Fig. 5 a . It is con- FIG. 3. Profiles for the 100 surface case: a the layer order ceivable that such effects influence the interpretation of some parameter n ; b the normalized layer energies En / J where J is of our results. Of course, one also has to discuss whether the the nearest-neighbor exchange constant; c the nearest-neighbor roughening transition TR for the 100 interface may be correlation function g(r 1) both for r 1 parallel and perpendicular to above Tcb , in which case a qualitatively different behavior the free surface. results, where the width remains finite and is determined by the finite range of correlations of the homogeneous bulk The width of the ordered layer grows with decreasing phases. A crude estimate for TR may be the transition tem- temperature and the profile develops into an S-shaped curve. perature of the purely 2D system, which is indeed above It can be seen, that the gradient gradually decreases as the Tcb , and even above Tcs . For further comparison, model interface moves into the bulk with decreasing temperature. calculations MC and cluster variation method CVM for While for temperatures slightly below the surface transition the same interaction model 50 show that the width of the 8944 W. SCHWEIKA, D. P. LANDAU, AND K. BINDER 53 FIG. 6. Variation of the fourth-order cumulant for the surface layer with temperature for different values of L for the 100 case. 100 interface for the Cu3Au stoichiometry is large 10 fcc cubes but finite at Tcb . Our data are not sufficient to give a clear and definite answer to this aspect whether the interface eventually becomes rough or not. While in our preliminary communication21 on a related model including a next-nearest-neighbor interaction as well the behavior of the profiles of order parameter and energy was qualitatively very similar to the data shown in Fig. 3, the statistical accuracy of these old data simply was not good enough to locate the surface transition temperature Tcs with sufficient precision to make any meaningful statement about the critical behavior of that model. Thus, in the present work, we take up this problem again: Fig. 6 shows that our present data which span a range of linear dimensions from L 32 to L 200 can locate Tcs with reasonable precision, namely kTcs /J 1.896 0.001. 40 The cumulant intersections see Refs. 35­37 for the finite- size scaling background that justifies this method also are compatible with a value U1* 0.61 0.01, i.e., within our accuracy this number is compatible with the value for the two-dimensional Ising universality class.38 This conclusion is supported by a finite-size scaling33­37 analysis: Fig. 7 a shows that the variation of the order parameter at Tcs is compatible with / 1/8, as expected; furthermore, the op- posite curvature for temperatures T Tcs and T Tcs sup- ports our belief that our estimate of Tcs is correct. From the slope of the cumulants at Tcs see Fig. 7 b we extract an exponent 0.9 0.1 which is consistent with the 2D Ising value. Of course, the accuracy of our exponent estimates still FIG. 5. a Fit dotted line of the order-parameter profile clearly is not yet very impressive, but since the data in Figs. for a 100 surface using Eq. 39 . Data refer to a simula- 6 and 7 have already required a large CPU effort, it will not tion at T 1.739 J /k be easy to improve these estimates further. B for L L D 80 80 180. b Semi- log plot of thickness of the ordered surface layer l¯ ( n 1/2) The second layer remains disordered at Tcs which justifies vs t; c semi-log plot of the squared width of the interface 2 the choice of the scalar two sublattice layer order parameter vs t. l¯ and 2 have the expected asymptotic behavior, however, for the consideration of the surface transition and explains the more detailed predictions of Eqs. 25 and 26 are not why the exponents we found agree with those of the 2D Ising confirmed. universality class. Using the same method of analyzing the 53 SURFACE-INDUCED ORDERING AND DISORDERING IN . . . 8945 FIG. 7. a Log-log plot of the 100 surface order parameter 1 vs L where four different temperatures close to Tcs are included; straight line shown has been drawn with the theoretical slope of 1/8. b log-log plot of the derivative of the cumulant versus the ratio of sizes L/L squares: L 128, dots: L 22). Straight line shown has a slope of 1/ with 0.9. finite-size effects we again found, within the accuracy of our simulations, a continuous ordering transition of the second layer at kTc,n 2 /J 1.791 0.003, with exponents and which are still compatible with those of the 2D Ising model. However, the wetting phase does not further proceed via a sequence of continuous layer transitions, since as mentioned above the gradient at the interface becomes finite and de- creases as T approaches Tcb . More precisely, we found that the transition in the first layer is accompanied by the transi- tion in the third layer and further odd layers , although the amplitude is much smaller. With the ordering of the second FIG. 8. Semi-log plot of a the surface excess order layer there is an analogous transition in the subsequent even s vs t and of b the surface excess energy Es vs t total, parallel, and perpen- layers again with exponentially decaying amplitudes. This dicular parts ; log-log plot of c the surface excess specific heat picture appears to be at least plausible, since the competing Cs for the 100 case. For comparison the slope 1 corresponding L10 ordering variants x and y of the full four sublattice to s 1 is given. order parameter lead the same scalar order parameter n for odd layers n, but to n which are different in sign for even detail. Figure 8 shows that the surface excess quantities are layers n. Furthermore, the additional degeneracy of the compatible with the behavior expected theoretically7,9­12 for ground state, which could lead to APB's and would disturb a SIO transition: s and Es are compatible with logarithmic this type of ordering, is already lifted at finite temperatures divergences as t 0 (T Tcb). Here, for the case of SIO we as it is for the infinite system . define the reduced temperature variable as t( ) 1 Tcb /T. Finally, we consider the surface induced ordering in more The excess specific heat was obtained by use of the appro- 8946 W. SCHWEIKA, D. P. LANDAU, AND K. BINDER 53 FIG. 10. Temperature dependence of the surface layer order parameter 1 broken curve and the bulk order parameter b full curve for the case of an 111 surface. FIG. 9. Log-log plot of a the surface excess susceptibility s vs t and of b the surface susceptibility 1 vs t in case of a 100 surface. Full dots open squares represent data based on the scalar two sublattice order parameter the 3D four sublattice order param- eter . s slowly enter the asymptotic region, where the theoretically expected modulus of the slope is given by s 1. priate fluctuation relation see Eqs. 36 and 37 and is compatible with a Curie-Weiss-like divergence, Cs t 1, i.e., s 1. In Fig. 9 we analyze the asymptotic behavior of the sus- ceptibilities s and 1 . The data for s are essentially the same if one considers the fluctuation of a scalar two sublat- tice layer order parameter or of the full 3D order parameter for bilayers. They are in agreement with the theoretical ex- pectation for the exponent s 1. The statistical accuracy of the layer susceptibilities 1 is not overwhelming. However, both layer and bilayer data seem to approach a constant value in the asymptotic limit T Tcb , which means that 1 0. Any comparison with the theoretical predictions dis- cussed in Sec. II A have to be made with care, since here and generally for SIO Tcs and Tcb do not coincide. Instead, FIG. 11. Profiles for the 111 surface case: a layer order pa- there is only an extraordinary transition at Tcb and hence a rameters vanishing exponent n ; b normalized layer energies En / J plotted versus 1 may be expected. layer number n. 53 SURFACE-INDUCED ORDERING AND DISORDERING IN . . . 8947 FIG. 13. Cumulants U1,L for the order parameter at the 111 surface plotted vs temperature, for different choices of L as indi- cated. scaling law 1 L 1 / which is plausible in view of Eqs. 6 and 7 we find 0.60 0.02. Note, that the theoreti- cal value is 1/2 see Eq. 19 , which assumes that the fluctuations in the surface plane are still governed by the fluctuations along the interface. The comparison with the data for a completely disordered film reveals that, at least for large systems, there remains only a trivial finite-size effect, which may have nothing to do with SID. One has to con- clude that 1 is not affected by and that the scaling law assumed above does not apply here. Remember, that 21 (1/L2) i 0 i 1/L2 and that 1 is of the same FIG. 12. a Log-log plot of the surface layer order parameter order as 21 1/2. Therefore finite-size scaling ideas are appli- for the 111 case vs t; b log-log plot of the surface layer order cable only if 1 / 1. parameter at Tcb vs L. Here open or filled symbols refer to runs We also attempted to locate the transition from the cumu- where the bulk is disordered or ordered, respectively. Straight lines lant crossings for different sizes Fig. 13 ; however, we find have slopes of about -1.04 and -2.06, respectively. no well-defined intersection point for any nonzero value of U V. RESULTS FOR 111... SURFACES 1 . This result can already be expected from the missing finite-size effects above, and further, as will be shown below, The behavior in this case is quite different from that seen also because the layer susceptibilities 1,1 do not diverge as in the previous section. There is no surface order above T Tcb . Tcb Fig. 10 and only finite-size short-range order can be In Fig. 14 we analyze the surface excess order parameter seen in the data. As the transition temperature is approached and surface excess energy. Both quantities are nicely com- from below, order-parameter profiles see Fig. 11 reveal that patible with the expected7,9­12 logarithmic variation as T ap- a somewhat disordered layer appears at the surface and the proaches Tcb . Similarly, the behavior of the surface excess interface between this layer and the bulk order moves into specific heat is compatible with Cs t 1, as expected see the bulk as the transition is approached. The disordering be- Eqs. 36 and 37 . havior is reflected in the disappearance of surface order as The critical behavior of the different surface related sus- T Tcb even though the bulk remains quite well ordered. We ceptibilities is analyzed in Fig. 15. Similar to the slow as- have analyzed 1 for power-law behavior and find, from Fig. ymptotic behavior of 1 and the excess quantities s and 12, that as the bulk transition temperature is approached, the Es , for t 0.003 the surface susceptibility 1 and the surface data slowly enter an asymptotic regime where they are well excess susceptibility s become at least compatible with the described by expected exponents, 1 0.36 as determined from 1 and the scaling relation Eq. 17a , and s 1 universal as theoreti- cally predicted7,9­12 see Eq. 16a . In our simulations 1 1 T/Tcb 1 t0.64 2 . 41 agreement was obtained for estimates of these two suscepti- Three different lattice sizes were used in the simulation. bilities from the appropriate fluctuation relations as well as However, even close at Tcb there are no obvious finite-size from the variation with the conjugate field H shown in Figs. effects for 1 parameter. In order to look more closely to 15 c and 15 d according to their definition as derivatives possible finite-size effects, we studied the size dependence of in Eqs. 4b and 4c . We have not presented as much data 1 at Tcb as shown in Fig. 12 b ; assuming a finite-size for the susceptibilities as we did for the order parameters, 8948 W. SCHWEIKA, D. P. LANDAU, AND K. BINDER 53 of 11 , we have subtracted the susceptibility contribution to 11 in a disordered bulk, estimating this constant simply from a simulation of a disordered system at Tcb which is equivalent to the case when the interface disappears towards infinite distance from the surfaces . All results are in accordance with the scaling relations Eqs. 17a ­ 17d although one can see that the asymptotic region is approached rather slowly by the susceptibility data. One may also note that the verification of critical wetting theory for the standard Ising model has been notoriously difficult.29­31 In view of the exponents which we find, our model sys- tem belongs to the first scaling regime of Eq. 24a for a rough interface. At this point, we recall that Lipowsky7,9­12 has predicted universal mean-field exponents for SID for the case of smooth interfaces only, i.e., 1 1/2, 1 1/2, 1,1 0, etc. That the interface is not smooth may be obvious from the layer profiles see Fig. 12 which reveal finite and modest gradients at the positions of the interface. However, as opposed to the situation of an interface in the bulk APB for T Tcb , the interface is bound to the surface and is not rough in the sense of the Kosterlitz-Thouless theory that it makes arbitrary excursions with increasing size of the sys- tem. Therefore, we examined the interfacial roughness in the limit T Tcb where the interface becomes unbound from the surface. From the layer profiles we obtained a fit to Eq. 39 the interface position n and the interface width . A typical example in case of a free 111 surface is shown in Fig. 16 a . The slight deviations near the surface being unimpor- tant for the present considerations stem from the missing bonds at the free surface. At T Tcb ,51 the procedure is car- ried out for a wide range of lateral dimensions L parallel to the surface, and the logarithmic divergence of 2 is clearly seen Fig. 16 d . Of course, this divergence is only possible in the limit T Tcb , where n also diverges, but for the tem- peratures considered the saturation of at a finite value due to the finiteness of n ) is not yet seen for the range of L displayed here. Finally, our analysis of the interface profiles enables us to check the validity of the theoretical predictions7 of the tem- perature dependence of the interface position and thus of the surface layer thickness of the wetting phase, l¯ n 1/2) and the width of the interface itself, as given in Eqs. 25 and 26 . The results shown in Figs. 16 b and 16 c confirm both a logarithmic law for the position which already followed from the behavior of the excess quantities s and Es) as well as a root logarithmic law for the interfacial thickness. How- ever, while Eqs. 25 and 26 can be fitted to the data in this case, there are some problems with the numbers that follow. One may note, that there is an additional constant term FIG. 14. a Semi-log plot of 111 surface excess order param- which is related to the depinning temperature of the inter- eter vs t 1 T/Tcb . Different symbols show different lattice sizes face . While further the proportionality is even consistent L; b same as a but for the surface excess energy; c Log-log with the value for 0.28 using our value of 1 0.64 and plot of the surface excess specific heat vs t. Straight line with slope Eq. 24a , the bulk correlation length, d 7.3 5 in units of 1 corresponding to s 1 is included for comparison. )/4 a , is obviously larger than the bulk correlation length as determined from independent simulations of the disor- because the derivatives need longer MC runs and more CPU dered bulk phase see Fig. 4 . As already discussed for the time. Our value 1 0.64 would imply also 11 0.28, a case of a free 100 surface, the problems may arise from the result which is nicely compatible with the data which are strong anisotropy of the correlations even in the disordered shown in Fig. 15 e . In order to analyze the critical behavior phase. And we again draw attention to the corrections that 53 SURFACE-INDUCED ORDERING AND DISORDERING IN . . . 8949 FIG. 15. Log-log plot of the 111 surface susceptibilities 1 a and s b vs t. Filled dots are data calculated according to the fluctuation relations using Eqs. 33 and 34 ; open squares represent equivalent data obtained as the derivatives of 1(H) and s(H) which are shown in c and d , respectively; e shows the singular part of 1,1 vs t using Eq. 32 . Straight lines indicate the expected exponents a 1 0.36, b s 1, and e 1,1 0.28, respectively. In parts c and d the broken straight lines indicate the error estimates for the fits, while the dotted curves through the data points are intended to guide the eye. 8950 W. SCHWEIKA, D. P. LANDAU, AND K. BINDER 53 FIG. 16. a Determination of the interfacial width for the case of a 111 surface by fitting the observed profile data points, and broken curve to Eq. 39 , for T 1.734. Here a system of size 120 120 160 is used, and fit parameters are n 8.22, 6.94. b The thickness of the disordered surface layer l¯ vs t, and c the squared thickness of the interface vs t exhibit asymptotically a logarithmic divergence. d Semi-log plot of the squared width vs the lateral linear dimension L at Tcb to demonstrate the logarithmic divergence expected for rough interfaces. Here T 1.738 is shown where for L 200, n was about 30 6. may result from the coupling of the interface position l¯ to the temperature. In this case interactions are frustrated both in length scale l1 over which the surface modifies the tanh pro- the bulk and in the surface plane. All of the various quanti- file, Fig. 16 a , which were discussed by Boulter and Parry.28 ties which we studied confirm the essential predictions of the actual wetting theory that there is only a single independent VI. CONCLUSIONS critical exponent, further the values of the expected universal exponents, and also the scaling relations derived for SID. We We have presented an extensive Monte Carlo study of find clear evidence that the exponent 1 differs from its surface phase transitions in a very simple model for AB bi- mean-field value. This fact implies that the fluctuation cor- nary alloys on a face-centered-cubic lattice. For the case of a rections proposed for critical wetting are indeed relevant 100 surface we have found a continuous surface transition here. Concerning the possible interfaces in the bulk, we ex- at Tcs which is above the bulk transition temperature Tcb and pect that slightly below Tcb the APB's in 111 planes to be surface-induced ordering as Tcb is approached from above. more rough than those in 100 planes, while the latter Qualitatively, the occurrence of this SIO is plausible due to should occur more often because of their lower energies. the existence of ``frustrated'' interactions in the bulk ordered We hope that this work will stimulate corresponding ex- structure, while no frustrated interactions occur in the 100 perimental studies. While there are quite a number of experi- surface. We find that the surface transition belongs to the ments related to SID, in particular there is no experimental universality class of the two-dimensional Ising model. work showing SIO in alloys. We have no doubt that this may For 111 surfaces we find instead that the surface order be a quite typical phenomenon and does not result merely parameter vanishes continuously right at the bulk transition from the special choice of the interaction model. SIO may be 53 SURFACE-INDUCED ORDERING AND DISORDERING IN . . . 8951 TABLE I. Results for Tcb in units of J /kb . Tcb Method Reference 1.893 CVM-T cluster variation method, tetrahedron approximation 40 1.746 5 High- and low-temperature series 41 1.71 MC Monte-Carlo 42 1.766 4 MC 43 1.810 CVM-TO tetrahedron-octahedron 44 1.73 MC 45 1.745 ``Mixed'' CVM TO and 13­14 point clusters 46 1.736 1 MC 47 1.738005 50 MC This work expected, for instance, in particular in those alloys where FIG. 17. Energy E versus N 1/3 slightly above Tcb , at long period ordered structures have been observed indicating T 1.7391, for the disordered phase. Here, full periodic boundary the low energy of antiphase boundaries. SIO may occur conditions are used. Finite-size effects become negligible for data at whenever a free surface favors a particular ordering variant, large N. in particular, if the ordered ground state has a lower than cubic symmetry. Therefore, at first sight one would not ex- however, in all cases one must eventually extrapolate to a pect SIO in alloys like Cu system of infinite linear dimension perpendicular to the sur- 3Au. However, the presence of a surface always lowers the symmetry of the bulk state, and faces. Therefore, it is advantageous to consider the system depending on the interaction model SIO may take place for without free surfaces and with fully periodic boundary con- any alternating order-parameter profile perpendicular to a ditions. A standard and well exploited method uses the inte- free surface. As will be discussed in another publication, the gration with inverse temperature of the specific heats.54 presence of sufficiently strong uniform surface fields can Starting from states of known entropy, i.e., T and also lead to SIO. In Appendix B we present arguments, how- T 0, one determines the free energy for the disordered and ever, that weak uniform surface fields which are expected to ordered branch. In order to obtain a reliable value for Tcb , as occur in real alloys would not change the phenomena quali- well as a reliable estimate for its range of confidence, one tatively. needs to carefully consider not only the apparent statistical errors but also all possible sources of systematic errors. We therefore carefully examined size effects, the quality of ran- ACKNOWLEDGMENTS dom numbers, and effects due to the manner in which the This work was partially supported by the National Sci- MC algorithm sweeps and selects sites in the crystal. Hence, ence Foundation NSF under Grants No. DMR-9100692 i we studied the finite-size effects and different choices of and DMR-9405018, by the Deutsche Forschungsgemein- distances between selected sites sligthly above Tcb , ii in schaft DFG under Grant No. SFB 262/D1, by the Alex- the high-T limit we compared the simulated results to exact ander von Humboldt Stiftung, and by the collaborative high-T series expansion, iii we used different integration NATO Grant No. CRG 921202. We are grateful to R. Lip- schemes. owsky for helpful comments on the manuscript. One of the The simulations were performed using a single-site spin- authors W.S. wishes to acknowledge W. Selke for discus- flip algorithm. The size of the system finally used was sions. 4 693, i.e., it consisted of four simple cubic sublattices each of equal linear dimensions L 69, containing in total APPENDIX A: DETERMINATION OF THE BULK more than 1.3 million Ising spins. For this size the correla- TRANSITION TEMPERATURE tion length always remained much smaller than L and finite- size effects did not introduce any significant systematic error A precise determination of the bulk transition temperature cf. Fig. 17, E(N) at T 1.7391 . The number of MC steps Tcb is needed in order for us to accurately extract the various per site ranged from 1000 far from the transition temperature critical exponents as given in this paper. Since the asymp- to 5000 close to Tcb this choice was determined by the totic region is reached very slowly (t 0.01) the uncertainty energy fluctuations . Each update of the entire lattice was in t should be considerably less than 10 4 for there to be two used in computing averages; however, for a proper estimate decades in t for the determination of the exponents. Existing of the statistical error each MC run was analyzed in parts, results in the literature are either incorrect or of insufficient chosen such that the averages taken for these did not show accuracy for this purpose, as can be seen from Table I. any correlations for subsequent data. We used 100 data Various ways exist to determine Tcb for this first-order points for the energies of the disordered high-T branch and transition see Refs. 52­55 for recent discussions . The MC 134 values for the low-T branch. The density of points was data for the finite system having either a 111 or a 100 free chosen approximately proportional to the curvature of the surface could, in principle, be used for this purpose as well; smooth function E(T) or E( ) , such that all contributions 8952 W. SCHWEIKA, D. P. LANDAU, AND K. BINDER 53 These error bars contain statistical as well as the estimated systematic uncertainties at a 95% confidence level twice the standard deviations . Statistical errors are twice as large as the systematic ones. APPENDIX B: DISCUSSION OF GROUND STATES AND THE EFFECT OF A UNIFORM SURFACE LAYER FIELD As is well known, in an Ising-type description of a binary alloy AB three pairwise interactions are present, EAB , EAA , and EBB , but for the description of bulk behavior only the combination EAB (EAA EBB)/2 matters which trans- lates into the exchange constant J of the equivalent Ising magnet , since the other relevant combination EAA EBB can be absorbed in the scale of the chemical potential difference between the species. At a surface, however, this is no longer true.1 The effect of missing neighbors is to generate a term FIG. 18. Free energy versus temperature as calculated by ther- H(1 1) (1 1) 1 Si in the Ising Hamiltonian, with H1 propor- modynamic integration for a large system (4 693 atoms with tional to EAA EBB with the sum restricted to the spins in the fully periodic boundary conditions. The first-order bulk transition surface plane. Note, that we use the superscript ``(1 1)'' temperature Tcb is determined by the crossing of the branches for to distinguish this uniform field from the staggered field H1 disordered and ordered phase. Eq. 4a , which is conjugated to the (2 2) superstructure to the integral have a similar statistical variance. The large in the surface layer. In addition, the surrounding medium number of data is a computational advantage, since by inte- may act with different forces on A and B species in the gration all data contribute to reduce the statistical error. surface layer, and hence this again results in a contribution of The high-T series for the nearest-neighbor pair-correlation the above form, H(1 1) 1 Si . Here we restrict ourselves for function is simplicity to short-range forces only, ignoring possible ef- fects on more distant layers in this qualitative discussion. (1 1) s Since in a real system the effect of a term H1 Si in is j i j 1 wntanhn J , A1 the Hamiltonian is always present, it is important to make where the first leading coefficients are w sure that it does not completely invalidate the description of 1 1, w2 4, w the phenomena of SIO and SID in our model, where such a 3 28 and w4 68. For the integration we used the simplest Newton-Cotes term was ignored so far. We first discuss the ground-state formulas, i.e., the trapezoidal rule and Simpson rule, as well behavior at T 0, assuming our orientation of the ordered as cubic spline approximations. The comparison of the dif- domains as shown in Fig. 1 stacking of planes ABAB along ferent integration methods allows us to estimate any possible the z axis . Now simple bond counting yields the following systematic errors introduced by the integration procedures. results for the energies per spin: The integration based on cubic spline approximations is more precise than the Newton-Cotes rules; however, since Ebulk 2 J , B1 the energy as a function of the inverse temperature is a rather smooth function, the trapezoidal and Simpson rules E 100 E 010 2 J , B2 already yield quite similar results. For comparison, the trap- ezoidal and Simpson rules yield T 1 1 cb 1.738 16 J /kB and E 001 H1 , B3 Tcb 1.738 01 J /kB , respectively . While the cubic splines have no specific requirements to the data basis, this allows us E to estimate systematic errors resulting from different choices 111 3 J /2. B4 for the location and number of nodes. Figure 18 illustrates Note that for 100 , 010 , and 111 surfaces the energies per the behavior of the two branches of the free energy near the spin in the surface planes do not depend on H(1 1) , because intersection point. 1 for the assumed stacking they all contain an equal number of To summarize our results, we obtained the following: A and B atoms. From the fact that E(111) is enhanced in T comparison to E cb 1.738 005 0.000 05 J /kB , bulk , it is already plausible that surface- induced disordering should occur for 111 surfaces, since F they are energetically less stable than the bulk. On the other cb 2.030 97 0.000 02 J , hand, for 100 and 010 surfaces- which are equivalent by E symmetry, of course-the energy of spins in the bulk and at 1.7729 0.0003 J , the surface is the same. Thus, simple ``bond counting'' argu- E 1.3323 0.0002 J ,S 0.1485 0.0003 kB , ments cannot prove that surface-induced ordering should oc- cur, entropic effects are crucial to stabilize the order in the S 0.4020 0.0002 kB . A2 surface plane at higher temperatures than in the bulk. Since 53 SURFACE-INDUCED ORDERING AND DISORDERING IN . . . 8953 this is an entropic effect, it it not surprising that SIO exists in 1 2 our model only at temperatures up to about 10% higher than 2F film mfilm H L2D H2 b D s T,H , T T T cb . For large enough fields the 001 surface is more stable C4 than the 100 and 010 surfaces, namely for with H 1 1 1 1 1 H1c 2 J . B5 2fb b mb , Alternating A, B planes parallel to the surfaces are more H H2 T T stable than the (2 2) structure. Thus, when one cuts an fcc crystal along an 100 plane such that one has a plate geom- 2fs . C5 etry where two surfaces have much larger area than all other s ms H H2 T T surfaces, the degeneracy among the three types of domains stacking of ABAB . . . planes along the x axis, y axis or z Since Ffilm , mfilm , film are densities of additive variables axis is lifted by surface effects: if Eq. B5 holds, the order- F(T,H) and its derivatives are extensive thermodynamic ing of the crystal will be such that the planes of pure A and variables , we can write them as sums over layer quantities: pure B are parallel to the surfaces, while for D H(1 1) (1 1) F 1 1 H1c , the arrangement shown in Fig. 1 a and f f studied at finite temperatures for H(1 1) film L2D D n T,H , C6 1 0 in Sec. IV of n 1 this paper will be the energetically preferred one. where f Of course, at finite temperature a field H(1 1) n(T,H) denotes the free energy per spin in the nth 1 0 will layer. This also holds for the derivatives have interesting effects on the SIO and SID behavior. E.g., in the case of the (2 2) structure in the 100 plane Fig. 1 a 1 D the A-rich and B-rich sublattice of the surface plane are no mfilm mn T,H , longer equivalent just as a magnetic field breaks the sym- D n 1 metry between the two sublattices of a simple antiferromag- net, enhancing the sublattice magnetization along the field 1 D direction and diminishing the sublattice magnetization with film D n T,H , C7 n 1 opposite orientation . However, we expect that the resulting surface enrichment effect will only lead to a shift of T where cs , and a short-range perturbation of the order parameter and energy profiles should result, while qualitatively the phenomena re- m , main the same. n fn H T APPENDIX C: FLUCTUATION RELATIONS FOR THE SURFACE EXCESS SUSCEPTIBILITIES n mn H 2fn H2 . C8 AND SPECIFIC HEATS T T Now we know that in the considered limit (D ) the pro- We consider the free-energy density of a thin film of files f thickness D and a lateral size of L L in absence of any n ,mn , n ultimately reach their bulk values in the cen- ter of the film, and thus it makes sense to rewrite surface fields H1 0 for simplicity : f n(T,H), mn(T,H), n(T,H) as follows: F 2 f n fn fD/2 fD/2 , etc. C9 L2D f b T,H D f s T,H f film T,H . C1 and hence assuming an even number of layers, but this is An analogous decomposition holds trivially for all deriva- not really essential tives, in particular for the magnetization per spin in the film: 2 D/2 1 2 f film f n T,H C10 m F D n 1 film L2D H mb D ms T,H as D , T C2 2 f f with b T,H D n T,H n 1 f b T,H . C11 m b fb H , Putting instead of D/2 in the last step is allowed, because T f n(T,H) f b 0 near n D/2, so we add only zeros. Com- paring Eqs. C1 and C10 we see that m s fs H , C3 T f f and the susceptibility per spin s T,H n T,H f b T,H . C12 n 1 8954 W. SCHWEIKA, D. P. LANDAU, AND K. BINDER 53 Obviously this reasoning goes through for any derivatives, 2 1 D e.g., efilm eb D es D en , C20 n 1 2 m we have for the specific heats per spin film mb T,H D mn T,H mb T,H , n 1 D C13 2 1 cfilm cb D cs D cn C21a n 1 ms T,H mn T,H mb T,H C14 and for the surface excess specific heat n 1 and cs cn cb . C21b n 1 2 film b T,H D n T,H b T,H , Since n 1 C15 e e 1 e c T T T2 s T,H n T,H b T,H . C16 (k n 1 B 1), from e H or L2Den Hn TrHne H/ Tre H we conclude Of course, this is also compatible with the fluctuation rela- tions, since we can simply start from the susceptibility fluc- T2cn L2D 1 HnH Hn H , tuation relation for the total film where H denotes the total Hamiltonian, while both Hn and 1 H are unnormalized Hamiltonians. The correct expression k 2 2 BT film for c L2D s s s hence is 1 1 T2c H s L2D nH Hn H HD/2H L2D s s s s . C17 n Due to the translational invariance in the directions parallel HD/2 H . C22 to the free surfaces we can decompose the summation over Note that cb has to be calculated as in a sum over the spins in the layer, and the sum n over layers: 1 T2cb L2D HnH Hn H n D/2, C23 1 film D n , with and not as n 1 T2c k n,n bT n s0s s0 s . C18 L2D HnHn Hn Hn n D/2 , C24 0, which is the analog of a layer susceptibility in the bulk. For where denotes any site in the system and 0 an element of calculating the specific heat in the bulk, one must correlate the nth layer. The formula for n can also be written as the energy in the bulk layer n D/2 with the total energy k H of the film, not just with the energy in that layers other- bT n L2D n n C19 wise transverse energy-energy correlations would be miss- derivation exactly as shown before for n 1). The order ing . Apparently, neither the surface excess specific heat cs parameter in the nth layer has to be correlated with the order nor the surface excess susceptibility s would be correctly in the whole film. Of course, the exactly analogous formulas obtained by considering instead of Eqs. C15 and C21b hold for derivatives with respect to temperature, and hence the fluctuation relations of the surface excess energy or the since the energy per spin of the whole film surface excess magnetization, respectively. 1 For reviews of surface critical phenomena, see K. Binder, in J.L. Lebowitz Academic, New York, 1988 , Vol. 12, Chap. 1, Phase Transitions and Critical Phenomena, edited by C. Domb and Ref. 2. and J.L. Lebowitz Academic, New York, 1983 , Vol. 8, Chap. 1. 4 D.E. Sullivan and M.M. Telo da Gama, in Fluid Interfacial Phe- 2 H.W. Diehl, in Phase Transitions and Critical Phenomena, edited nomena, edited by C.A. Croxton Wiley, New York, 1986 , p. by C. Domb and J.L. Lebowitz Academic, New York, 1986 , 45. Vol. 10, p. 75. 5 P.G. de Gennes, Rev. Mod. Phys. 57, 825 1985 . 3 For reviews of wetting phenomena, see S. Dietrich, in Phase 6 M.E. Fisher, J. Stat. Phys. 34, 667 1984 ; J. Chem. Soc. Faraday Transitions and Critical Phenomena, edited by C. Domb and Trans. 282, 1569 1986 . 53 SURFACE-INDUCED ORDERING AND DISORDERING IN . . . 8955 7 R. Lipowsky, J. Appl. Phys. 55, 2485 1984 ; Ferroelectrics 73, 31 K. Binder, D.P. Landau, and. S. Wansleben, Phys. Rev. B 40, 69 1987 . 6971 1989 . 8 H. Dosch, in Critical Phenomena at Surfaces and Interfaces, ed- 32 K. Binder and D.P. Landau, Phys. Rev. B 46, 4844 1992 . ited by G. Ho¨hler, Springer Tracts in Modern Physics, Vol. 126 33 R. Lipowsky, Dissertation, University of Mu¨nchen, 1982. Springer, Berlin, 1992 . 34 K. Binder and D.P. Landau, Phys. Rev. Lett. 52, 318 1984 ; D.P. 9 R. Lipowsky, Phys. Rev. Lett. 49, 1575 1982 . Landau and K. Binder, Phys. Rev. B 41, 4786 1990 . 10 R. Lipowsky, Z. Phys. B 51, 165 1983 . 35 K. Binder and D.W. Heermann, The Monte Carlo Method in Sta- 11 R. Lipowsky and W. Speth, Phys. Rev. B 28, 3983 1983 . tistical Physics. An Introduction Springer, Berlin, 1988 . 12 R. Lipowsky, Z. Phys. B 55, 335 1984 . 36 Finite Size Scaling and the Numerical Simulation of Statistical 13 V.S. Sundaram, B. Farrell, R.S. Alben, and W.D. Robertson, Phys. Systems, edited by V. Privman World Scientific, Singapore, Rev. Lett. 31, 1136 1973 . 14 1990 . V.S. Sundaram, R.S. Alben, and W.D. Robertson, Surf. Sci. 46, 37 K. Binder, in Computational Methods in Field Theory, edited by 653 1974 . 15 E.G. McRae and R.A. Malic, Surf. Sci. 148, 551 1984 . C.B. Lang and A. Gausterer Springer, Berlin, 1992 , p. 59. 38 16 S.F. Alvarado, M. Campagna, A. Fattah, and W. Uelhoff, Z. Phys. T.W. Burkhardt and B. Derrida, Phys. Rev. B 32, 7273 1985 ; B 66, 103 1987 . A.D. Bruce, J. Phys. A 18, L873 1985 ; D.P. Landau and D. 17 H. Dosch, L. Maila¨nder, A. Lied, J. Peisl, F. Grey, R.L. Johnson, Stauffer, J. Phys. Paris 50, 509 1989 . 39 and S. Krumnacher, Phys. Rev. Lett. 60, 2382 1988 ; H. Dosch, J.W. Cahn and J.C. Hilliard, J. Chem. Phys. 28, 258 1958 . 40 L. Maila¨nder, H. Reichert, J. Peisl, and R.L. Johnson, Phys. Rev. R. Kikuchi and H. Sato, Acta Metall. 22, 1099 1974 . B 43, 13 172 1991 . 41 D. F. Styer, Phys. Rev. B 17, 2926 1978 . 18 B. Pluis, A.W. Denier van der Gon, J.F. van der Veen, and A.J. 42 K. Binder, Phys. Rev. Lett. 45, 811 1980 . Riemersma, Surf. Sci. 239, 265 1990 ; B. Pluis, D. Frenkel, and 43 M.K. Phani, J.L. Lebowitz, and M.K. Kalos, Phys. Rev. B 21, J.F. van der Veen, ibid. 239, 282 1990 . 4027 1980 . 19 K.C. Prince, U. Breuer, and H.P. Bonzel, Phys. Rev. Lett. 60, 44 A. Finel and F. Ducastelle, Europhys. Lett. 1, 135 1986 . 1146 1988 . 45 H.T. Diep, A. Ghazali, B. Berge, and P. Lallemand, Europhys. 20 Ch. Ricolleau, A. Loiseau, F. Ducastelle, and R. Caudron, Phys. Lett. 2, 603 1986 . Rev. Lett. 68, 3591 1992 . 46 A. Finel, in Alloy Phase Stability, NATO ASI Series E, edited by 21 W. Schweika, K. Binder, and D.P. Landau, Phys. Rev. Lett. 65, G.M. Stocks and A. Gonis Kluwer, Dordrecht, 1989 , p. 278. 3321 1990 . 47 T.L. Polgreen, Phys. Rev. B 29, 1468 1984 . 22 R. Lipowsky, D.M. Kroll, and R.K.P. Zia, Phys. Rev. B 27, 4499 48 D.P. Landau and K. Binder, Phys. Rev. B 41, 4786 1990 . 1983 ; see also D. M. Kroll and R. Lipowsky, ibid. 28, 6435 49 T.W. Burkhardt and H.W. Diehl, Phys. Rev. B 50, 3894 1994 . 1983 . 50 A. Finel, in Statics and Dynamics of Alloy Phase Transforma- 23 E. Brezin, B.I. Halperin, and S. Leibler, Phys. Rev. Lett. 50, 1387 tions, Vol. 319 of NATO ASI Series B: Physics, edited by P.E.A. 1983 . Turchi and A. Gonis Plenum, New York, 1994 , p. 495. 24 D.S. Fisher and D.A. Huse, Phys. Rev. B 32, 247 1985 . 51 Long MC runs are required. Note that the accurate temperature in 25 R. Lipowsky and M.E. Fisher, Phys. Rev. B 36, 2126 1987 . the simulation is also determined from the distribution of the 26 E. Brezin and T. Halpin-Healey, Phys. Rev. Lett. 58, 1220 1987 ; selected random numbers. J. Phys. Paris 48, 757 1987 . 52 K. Binder, K. Vollmayr, H.-P. Deutsch, J.D. Reger, M. Scheucher, 27 M.E. Fisher and A.J. Jin, Phys. Rev. Lett. 69, 792 1992 ; A.J. Jin and D.P. Landau, in Dynamics of First Order Phase Transitions, and M.E. Fisher, Phys. Rev. B 47, 7365 1993 ; M.E. Fisher, edited by H.J. Herrmann, W. Janke, and F. Karrch World Sci- A.J. Jin, and A.D. Parry, Ber. Bunsenges. Phys. Chem. 98, 357 entific, Singapore, 1992 , p. 253. 1994 . 53 W. Janke, in Computer Simulation Studies in Condensed Matter 28 C.J. Boulter and A.O. Parry, Phys. Rev. Lett. 74, 3403 1995 . Physics VII, edited by D.P. Landau, K.K. Mon, and H.-B. Schu¨t- 29 K. Binder, D.P. Landau, and D.M. Kroll, Phys. Rev. Lett. 56, tler Springer, Berlin, 1994 , p. 29. 2276 1986 . 54 Monte Carlo Methods in Statistical Physics, edited by K. Binder 30 K. Binder and D.P. Landau, Phys. Rev. B 37, 1745 1988 ; J. Springer, Berlin, 1979 . Appl. Phys. 57, 3306 1985 . 55 A. Hu¨ller, Z. Phys. B 93, 401 1994 .