# Monte Carlo Results

It is known that methods for quantum spins, such as the spin wave theory or the GF method presented above, are not satisfactory at high temperatures. The spin wave theory, even with magnon- magnon interactions taken into account, cannot go to temperatures close to *T _{c}.* The GF method on the other hand can go up to

*T*but due to the decoupling scheme, it cannot give correct critical behavior at

_{c }*T*Fortunately, we know that quantum spin systems behave as their classical counterparts at high

_{c}.*T.*So, to see if the phase diagram obtained in the previous section for the quantum model is correct or not, we can consider its classical version and use MC simulations to obtain the phase diagram for comparison. MC simulations are excellent means to overcome approximations used in analytic calculations for the high

*T*region as discussed above.

Figure 10.6 Magnetizations of layer 1 (circles) and layer 2 (diamonds) versus temperature *T* in unit of *] /k _{B}* for

*)*= 0.5 with

_{s}*I = l*0.1,

_{s}=*L=*36.

In this paragraph, we show the results obtained by MC simulations with the Hamiltonian (10.1) but the spins are the classical Heisenberg model of magnitude 5 = 1.

The film sizes are Lx Lx *N _{z}* where

*N*4 is the number of layers (film thickness) taken as in the quantum case presented above. We use here L = 24, 36, 48, 60 to study finite-size effects as shown below. Periodic boundary conditions are used in the

_{z}=*X Y*planes. The equilibrating time is about 10

^{6}MC steps per spin and the averaging time is 2 x 10

^{6}MC steps per spin.

*]*= 1 is taken as unit of energy in the following.

Let us show in Fig. 10.6 the layer magnetization of the first two layers as a function of *T,* in the case *J _{s} = 0.5* with

*N*(the third and fourth layers are symmetric). Although we observe a smaller magnetization for the surface layer, there is clearly no surface transition. This is very similar to the quantum case.

_{z}= 4In Fig. 10.7, we show a frustrated case where *J _{s}* = -0.5. The surface layer in this case becomes disordered at a temperature much lower than that for the second layer. Note that the surface magnetization is not saturated to 1 at Г = 0. This is because the

Figure 10.7 Magnetizations of layer 1 (circles) and layer 2 (diamonds) versus temperature *T* in unit of *J/кв* for *J _{s}* = —0.5 with

*I = —I*= 0.1,

_{s}*L=*36.

Figure 10.8 Susceptibilities of layer 1 (circles) and layer 2 (diamonds) versus temperature *T* in unit of *] /k _{B}* for

*)*= 0.5 with

_{s}*I = l*0.1,

_{s}=*L=*36.

surface spins make an angle with the *z* axis so their *z* component is less than 1 in the GS.

Figure 10.8 shows the susceptibilities of the first and second layers in the case where /* = 0.5 with / = /* = 0.1 where one

Figure 10.9 Susceptibility of layer 1 (circles) and layer 2 (diamonds) versus temperature *T* in unit of *J/k _{B}* for

*j*= —0.5 with / =

_{s}*—l*0.1,

_{s}=*L =*36. Note that for clarity, the susceptibility of the layer 2 has been multiplied by a factor 5.

observes the peaks at the same temperature indicating a single transition in contrast to the frustrated case shown in Fig. 10.9. These results confirm the above results of layer magnetizations.

To establish the phase diagram, the transition temperatures are taken at the change of curvature of the layer magnetizations, i.e., at the maxima of layer susceptibilities shown before. Figure 10.10 shows the phase diagram obtained in the space (*J _{s}, T*). It is interesting to see that this phase diagram resembles remarkably that obtained for the quantum counterpart model shown in Fig. 10.5.

Let us study the finite-size effect of the phase transitions shown in Fig. 10.10. To this end, we use the histogram technique, which has been proved so far to be excellent for the calculation of critical exponents [110-112]. The principle is as follows (see details in Chapter 6). Using the Metropolis algorithm to determine approximately the critical temperature region, then choosing a temperature *T _{0}* as close as possible to the presupposed transition temperature. We then make a very long run at

*T*to establish an energy histogram. From formulas established using the statistical canonical distribution, we can calculate physical quantities in a continuous manner for temperatures around

_{0}*T*[110-112]. We can

_{0}Figure 10.10 Phase diagram in the space *(J _{s}, T*) for the classical Heisenberg model with

*N*= |/

_{z}= 4, I_{s}| = 0.1. Phases I to III have the same meanings as those in Fig. 10.5.

Figure 10.11 Susceptibility versus *T* for *L* = 36, 48, 60 with *J _{s} =* 0.5 and

*I = Is =* 0.1.

then identify with precision the transition temperature as well as the maximal values of fluctuation quantities such as specific heat and susceptibility.

Figure 10.11 shows the susceptibility versus *T for L=* 36, 48, 60 in the case of/_{s} = 0.5. For presentation convenience, the size L= 24 has been removed since the peak for this case is rather flat in the

Figure 10.12 Transition temperature versus *L* with *) _{s}* = 0.5 and

*I*=

*l*= 0.1. The inset shows the enlarged scale.

_{s}scale of the figure. However, it shall be used to calculate the critical exponent *у* for the transition. As seen, the maximum /^{max} of the susceptibility increases with increasing *L.*

For completeness, we show in Fig. 10.12 the transition temperature as a function of *L.* A rough extrapolation to infinity gives *T£°* ~ 1.86 ±0.02.

In the frustrated case, i.e., *) _{s} <* //, we perform the same calculation to estimate the finite-size effect. Note that in this case there are two phase transitions. We show in Fig. 10.13 the layer susceptibilities as functions of

*T*for different

*L.*As seen, both surface and second-layer transitions have a strong size dependence.

We show in Figs. 10.14 and 10.15 the size dependence of the transition temperatures *T* (surface transition) and *T** _{2}* (second-layer transition).

The size dependence of the maxima observed above allows us to estimate the ratio *y/v* (see Chapter 6 for details on the finite- size scaling laws). We show now ln/^{max} as a function of InL for the different cases studied above. Figures 10.16a and 10.16b correspond, respectively, to the transitions of surface and second layer occurring in the frustrated case with *J _{s}=* - 0.5, while Fig. 10.16c corresponds to the unique transition occurring in the non frustrated case with

*J*= 0.5. The slopes of these straight lines give

_{s}*y/v*~ 1.864 ± 0.034 (a), 1.878 ± 0.027 (b), 1.801 ± 0.027 (c). The errors were estimated from the mean least-square fitting

Figure 10.13 Layer susceptibilities versus *T* for *L* = 36, 48, 60 with *] _{s} = *—0.5 and / = —

*I*= 0.1. Left (right) figure corresponds to the first (second) layer susceptibility.

_{s}Figure 10.14 Transition temperature versus *L* for the surface layer in the case *] _{s}* = —0.5 with

*1 = — I*0.1. The inset shows the enlarged scale.

_{s}=Figure 10.15 Transition temperature versus *L* for the second layer in the case *] _{s}* = —0.5 with / = —

*I*= 0.1. The inset shows the enlarged scale.

_{s}Figure 10.16 Maximum of surface-layer susceptibility versus *L* for *L = *24, 36, 48, 60 with *] _{s} =* —0.5 (a,b),

*]*0.5 (c) and

_{s}=*I =*|/

_{s}| = 0.1, in the In — In scale. The slope gives

*y/v*indicated in the figure for each case. See text for comments.

and errors on the peak values obtained with different values of *T** _{0 }*(multiple-histogram technique). Within errors, the first two values, which correspond to the frustrated case, can be considered as identical, while the last one corresponding to the non-frustrated case is different. We will return to this point later.

At this stage, we would like to emphasize the following points. First, we observe that these values of *y/v* are found to be between those ofthe two-dimensional (2D) Ising model (*y/v =* 1.75) and the three-dimensional (3D) one (*y/v* = 1.241/0.63 ~ 1.97). A question which naturally arises on the roles of the Ising-like anisotropy term, of the two-fold chiral symmetry and of the film thickness. The roles of the anisotropy term and the chiral symmetry are obvious: the Ising character should be observed in the result (we return to this point below). It is, however, not clear for the effect of the thickness. Some arguments, such as those from the renormalization group, say that the correlation length in the direction perpendicular to the film is finite; hence, it is irrelevant to the criticality; the 2D character, therefore, should be theoretically preserved. We think that such arguments are not always true because it is difficult to conceive that when the film thickness becomes larger and larger the 2D universality should remain. Instead, we believe that that the finite thickness of the film affects the 2D universality in one way or another, giving rise to "effective” critical exponents with values deviated from the 2D ones. The larger the thickness is, the stronger this deviation becomes. The observed values of *y/v* shown above may contain an effect of a 2D-3D crossover. At this point, we would like to emphasize that, in the case of simple surface conditions, i.e., no significant deviation of the surface parameters with respect to those of the bulk, the bulk behavior is observed when the thickness becomes larger than a dozen of atomic layers [75, 78]: surface effects are insignificant on thermodynamic properties of the film. There are, therefore, reasons to believe that there should be a crossover from 2D to 3D at some film thickness. Of course, this is an important issue which needs to be theoretically clarified in the future. We return now to the effect of Ising anisotropy and chiral symmetry. The deviation from the 2D values may result in part from a complex coupling between the Ising symmetry, due to anisotropy and chirality, and the continuous nature of the classical Heisenberg spins studied here. This deviation may be important if the anisotropy constant *1* is small. For the effect chiral symmetry, it is a complex matter. To show the complexity in determining the critical universality with chiral symmetiy, let us discuss about a simpler case with similar chiral symmetry: the XY model on the fully frustrated Villain’s lattice

(see model described in Problem 7 of Chapter 3 and its solution in Section 18.3). There has been a large number of investigations on the nature of the transition observed in this 2D case in the context of the frustration effects [45, 85, 264]. In this model, though the chirality symmetiy argument says that the transition should be of 2D Ising universality class, many investigators found a deviation of critical exponents from those of the 2D case (see review in Ref. [45]). For example, the following values are found for the critical exponent *w. v* = 0.889 [287] and *v* = 0.813 [203]. These values are close to that obtained for the single transition in a mixed ХУ-Ising model which is 0.85 [204, 257]. It is now believed that the XY character of the spins affects the Ising chiral symmetry giving rise to those deviated critical exponents. Similarly, in the case of thin film studied here, we do not deal with the discrete Ising model but rather an Ising- like Heisenberg model. The Ising character due to chiral symmetry of the transition at the surface is believed to be perturbed by the continuous nature of Heisenberg spins. The transition of interior layer shown in Fig. 10.16b suffers similar but not the same effects because of the absence of chiral symmetry on this layer. So the value *y/v* is a little different. In the non frustrated case shown in Fig. 10.16c, the deviation from the 2D Ising universality class is less important because of the absence of the chiral symmetry. This small deviation is believed to stem mainly from the continuous nature of Heisenberg spins.

To conclude this paragraph, we believe, from physical arguments given above, that the deviation from 2D Ising universality class of the transitions observed here is due to the following causes: the effect of the coupling between the continuous degree of freedom of Heisenberg spin to the chiral symmetry, the small Ising-like anisotropy and the film thickness.

# Concluding Remarks

We have used the GF method and MC simulations to study the Heisenberg spin model with an Ising-like interaction anisotropy in thin films of stacked triangular lattices. The two surfaces of the film are frustrated. We found that surface spin configuration is non-collinear when the surface antiferromagnetic interaction is smaller than a critical value //. In the non-collinear regime, the surface layer is disordered at a temperature lower than that for interior layers ("soft" surface) giving rise to the so-called "magnetically dead surface" observed in some materials [36, 374]. The surface transition disappears for *) _{s}* larger than the critical value _//. Using the GF method we have established a phase diagram in the space (T,

*Js)*which is in a good agreement with the MC result. This is due to the fact that at high temperatures where the transition takes place, the quantum nature of spins used in the GF is lost so that we should find results of classical spins obtained by MC simulations. We have also studied by MC histogram technique the critical behavior of the phase transition using the finite-size effects. The result of the ratio of critical exponents

*y/v*shows that the nature of the transition is complicated due to the influence of several physical mechanisms. The symmetry of the ground state alone cannot explain such a result. We have outlined a number of the most relevant mechanisms. We will return to the criticality of thin films in Chapter 16.