Theoretical methods for understanding advanced magnetic materials: the case of frustrated thin films

Materials science has been intensively developed during the last 30 years. This is due, on the one hand, to an increasing demand of new materials for new applications and, on the other hand, to technological progress which allows for the synthesis of materials of desired characteristics and to investigate their properties with sophisticated experimental apparatus. Among these advanced materials, magnetic materials at nanometric scale such as ultra thin films or ultra fine aggregates are no doubt among the most important for electronic devices.In this review, we show advanced theoretical methods and solved examples that help understand microscopic mechanisms leading to experimental observations in magnetic thin films. Attention is paid to the case of magnetically frustrated systems in which two or more magnetic interactions are present and competing. The interplay between spin frustration and surface effects is the origin of spectacular phenomena which often occur at boundaries of phases with different symmetries: reentrance, disorder lines, coexistence of order and disorder at equilibrium. These phenomena are shown and explained using of some exact methods, the Green's function and Monte Carlo simulation. We show in particular how to calculate surface spin-wave modes, surface magnetization, surface reorientation transition and spin transport.


I. INTRODUCTION
Material science has made a rapid and spectacular progress during the last 30 years, thanks to the advance of experimental investigation methods and a strong desire of scientific community to search for new and high-performance materials for new applications. In parallel to this intensive development, many efforts have been devoted to understanding theoretically microscopic mechanisms at the origin of the properties of new materials. Each kind of material needs specific theoretical methods in spite of the fact that there is a large number of common basic principles that govern main properties of each material family.
In this paper, we confine our attention to the case of magnetic thin films. We would like to show basic physical principles that help us understand their general properties. The main purpose of the paper is not to present technical details of each of them, but rather to show what can be understood using each of them. For technical details of a particular method, the reader is referred to numerous references given in the paper. For demonstration purpose, we shall use magnetically frustrated thin films throughout the paper. These systems combine two difficult subjects: frustrated spin systems and surface physics. Frustrated spin systems have been subject of intensive studies during the last 30 years [1]. Thanks to these efforts many points have been well understood in spite of the fact that there remains a large number of issues which are in debate. As seen below, frustrated spin systems contain many exotic properties such as high ground-state degeneracy, new symmetries, successive phase transitions, reentrant phase and disorder lines. Frustrated spin systems serve as ideal testing grounds for theories and approximations. On the other hand, during the same period surface physics has also been widely investigated both experimentally and theoretically. Thanks to technological progress, films and surfaces with desired properties could be fabricated and characterized with precision. As a consequence, one has seen over the years numerous technological applications of thin films, coupled thin films and super-lattices, in various domains such as magnetic sensors, magnetic recording and data storage. One of the spectacular effects is the colossal magnetoresistance [2,3] which yields very interesting transport properties. The search for new effects with new mechanisms in other kinds of materials continues intensively nowadays as never before.
Section II is devoted to the presentation of the main theoretical background and concepts to understand frustrated spin systems and surface effects in magnetic materials. Needless to say, one cannot cover all recent developments in magnetic materials but an effort is made to outline the most important ones in our point of view. Section III is devoted to a few examples to illustrate striking effects due to the frustration and to the presence of a surface. Concluding remarks are given in Section IV.

A. Theory of Phase Transition
Many materials exhibit a phase transition. There are several kinds of transition, each transition is driven by the change of a physical parameter such as pressure, applied field, temperature (T ), ... The most popular and most studied transition is no doubt the one corresponding to the passage from a disordered phase to an ordered phase at the so-called magnetic ordering temperature or Curie temperature T c . The transition is accompanied by a symmetry breaking. In general when the symmetry of one phase is a subgroup of the other phase the transition is continuous, namely the first derivatives of the free energy such as internal energy and magnetization are continuous functions of T . The second derivatives such as specific heat and susceptibility, on the other hand, diverge at T c . The correlation length is infinite at T c . When the symmetry of one phase is not a symmetry subgroup of the other, the transition is in general of first order: the first derivatives of the free energy are discontinuous at T c . At the transition, the correlation length is finite and often there is a coexistence of the two phases. For continuous transitions, also called second-order transition, the nature of the transition is characterized by a set of critical exponents which defines its "universality class". Transitions in different systems may belong to the same universality class.
Why is the study of a phase transition interesting? As the theory shows it, the characteristics of a transition are intimately connected to microscopic interactions between particles in the system.
The theory of phase transitions and critical phenomena has been intensively developed by Landau and co-workers since the 50's in the framework of the mean-field theory. Microscopic concepts have been introduced only in the early 70's with the renormalization group [4][5][6]. We have since then a clear picture of the transition mechanism and a clear identification of principal ingredients which determine the nature of the phase transition. In fact, there is a small number of them such as the space dimension, the nature of interaction and the symmetry of the order parameter.

B. Frustrated Spin Systems
A spin is said "frustrated" when it cannot fully satisfy all the interactions with its neighbors. Let us take a triangle with an antiferromagnetic interaction J(< 0) between two sites: we see that we cannot arrange three Ising spins (±1) to satisfy all three bonds. Among them, one spin satisfies one neighbor but not the other. It is frustrated. Note that any of the three spins can be in this situation. There are thus three equivalent configurations and three reverse configurations, making 6 the number of "degenerate states". If we put XY spins on such a triangle, the configuration with a minimum energy is the so-called "120-degree structure" where the two neighboring spins make a 120 0 angle. In this case, each interaction bond has an energy equal to |J| cos(2π/3) = −|J|/2, namely half of the full interaction: the frustration is equally shared by three spins, unlike the Ising case. Note that if we go from one spin to the neighboring spin in the trigonometric sense we can choose cos(2π/3) or − cos(2π/3) for the turn angle: there is thus a two-fold degeneracy in the XY spin case. The left and right turn angles are called left and right chiralities. In an antiferromagnetic triangular lattice, one can construct the spin configuration from triangle to triangle. The frustration in lattices with triangular plaquettes as unit such as in face-centered cubic and hexagonal-close-packed lattices is called "geometry frustration". Another category of frustration is when there is a competition between different kinds of incompatible interactions which results in a situation where no interaction is fully satisfied. We take for example a square with three ferromagnetic bonds J(> 0) and one antiferromagnetic bond -J, we see that we cannot "fully" satisfy all bonds with Ising or XY spins put on the corners.
Frustrated spin systems are therefore very unstable systems with often very high ground-state degeneracy. In addition, novel symmetries can be induced such as the two-fold chirality seen above. Breaking this symmetry results in an Ising-like transition in a system of XY spins [7,8]. As will be seen in some examples below, the frustration is the origin of many spectacular effects such as non collinear ground-state configurations, coexistence of order and disorder, reentrance, disorder lines, multiple phase transitions, etc.

C. Surface Magnetism
In thin films the lateral sizes are supposed to be infinite while the thickness is composed of a few dozens of atomic layers. Spins at the two surfaces of a film lack a number of neighbors and as a consequence surfaces have physical properties different from the bulk. Of course, the difference is more pronounced if, in addition to the lack of neighbors, there are deviations of bulk parameters such as exchange interaction, spin-orbit coupling and magnetic anisotropy, and the presence of surface defects and impurities. Such changes at the surface can lead to surface phase transition separated from the bulk transition, and surface reconstruction, namely change in lattice structure, lattice constant [9], magnetic ordering, ... at the surface [10][11][12].
Thin films of different materials, different geometries, different lattice structures, different thicknesses ... when coupled give surprising results such as colossal magnetoresistance [2,3]. Microscopic mechanisms leading to these striking effects are multiple. Investigations on new artificial architectures for new applications are more and more intensive today. In the following section, we will give some basic microscopic mechanisms based on elementary excitations due to the film surface which allows for understanding macroscopic behaviors of physical quantities such as surface magnetization, surface phase transition and transition temperature.

D. Methods
To study properties of materials one uses various theories in condensed matter physics constructed from quantum mechanics and statistical physics [13,14]. Depending on the purpose of the investigation, we can choose many standard methods at hand (see details in Refs. [15,16]): (i) For a quick obtention of a phase diagram in the space of physical parameters such as temperature, interaction strengths, ... one can use a mean-field theory if the system is simple with no frustration, no disorder, ... Results are reasonable in three dimensions, though critical properties cannot be correctly obtained (ii) For the nature of phase transitions and their criticality, the renormalization group [4][5][6] is no doubt the best tool. However for complicated systems such as frustrated systems, films and dots, this method is not easy to use (iii) For low-dimensional systems with discrete spin models, exact methods can be used (iv) For elementary excitations such as spin waves, one can use the classical or quantum spin-wave theory to get the spin-wave spectrum. The advantage of the spin-wave theory is that one can keep track of the microscopic effect of a given parameter on macroscopic properties of a magnetic system at low temperatures with a correct precision (v) For quantum magnetic systems, the Green's function method allows one to calculate at ease the spin-wave spectrum, quantum fluctuations and thermodynamic properties up to rather high temperatures in magnetically ordered materials. This method can be used for collinear spin states and non collinear (or canted) spin configurations as seen below (vi) For all systems, in particular for complicated systems where analytical methods cannot be easily applied, Monte Carlo simulations can be used to calculate numerous physical properties, specially in the domain of phase transitions and critical phenomena as well as in the spin transport as seen below.
In the next section, we will show some of these methods and how they are practically applied to study various properties of thin films.

A. Exactly Solved Two-dimensional Models
Why are exactly solved models interesting? There are several reasons to study such models: • Many hidden properties of a model cannot be revealed without exact mathematical demonstration • We do not know of any real material which corresponds to an exactly solved model, but we know that real materials should bear physical features which are not far from properties described in some exactly solved models if similar interactions are thought to exist in these materials.
• Macroscopic effects observed in experiments cannot always allow us to find their origins if we do not already have some theoretical pictures provided by exact solutions in mind.
To date, only systems of discrete spins in one and two dimensions (2D) with short-range interactions can be exactly solved. Discrete spin models include Ising spin, q-state Potts models and some Potts clock-models. The reader is referred to the book by Baxter [17] for principal exactly solved models. In general, one-dimensional (1D) models with short-range interaction do not have a phase transition at a finite temperature. If infinite-range interactions are taken into account, then they have, though not exactly solved, a transition of second or first order depending on the decaying power of the interaction [18][19][20][21]. In 2D, most systems of discrete spins have a transition at a finite temperature. The most famous model is the 2D Ising model with the Onsager's solution [22].
In this paper, we are also interested in frustrated 2D systems because thin films in a sense are quasi two-dimensional. We have exactly solved a number of frustrated Ising models such as the Kagomé lattice [23], the generalized Kagomé lattice [24], the generalized honeycomb lattice [25] and various dilute centered square lattices [26][27][28].
For illustration, let us show the case of a Kagomé lattice with nearest-neighbor (NN) and next-nearest neighbor (NNN) interactions. As seen below this Kagomé model possesses all interesting properties of the other frustrated models mentioned above.
In general, 2D Ising models without crossing interactions can be mapped onto the 16-vertex model or the 32-vertex model which satisfy the free-fermion condition automatically as shown below with an Ising model defined on a Kagomé lattice with interactions between NN and between NNN, J 1 and J 2 , respectively, as shown in  We consider the following Hamiltonian where σ i = ±1 and the first and second sums run over the spin pairs connected by single and double bonds, respectively. Note that the Kagomé original model, with antiferromagnetic J 1 and without J 2 interaction, has been exactly solved a long time ago showing no phase transition at finite temperatures [29]. The ground state (GS) of this model can be easily determined by an energy minimization. It is shown in Fig. 2 where one sees that only in zone I the GS is ferromagnetic. In other zones the central spin is undetermined because it has two up and two down neighbors, making its interaction energy zero: it is therefore free to flip. The GS spin configurations in these zones are thus "partially disordered". Around the line J 2 /J 1 = −1 separating zone I and zone IV we will show below that many interesting effects occur when T increases from zero.
The partition function is written as where K 1,2 = J 1,2 /k B T and where the sum is performed over all spin configurations and the product is taken over all elementary cells of the lattice. To solve this model, we first decimate the central spin of each elementary cell of the lattice and obtain a checkerboard Ising model with multi-spin interactions (see Fig. 3).
The Boltzmann weight of each shaded square is given by  The partition function of this checkerboard Ising model is thus where the sum is performed over all spin configurations and the product is taken over all the shaded squares of the lattice. To map this model onto the 16-vertex model, we need to introduce another square lattice where each site is placed at the center of each shaded square of the checkerboard lattice, as shown in Fig. 4. At each bond of this lattice we associate an arrow pointing out of the site if the Ising spin that is traversed by this bond is equal to +1, and pointing into the site if the Ising spin is equal to -1, as it is shown in Fig. 5. In this way, we have a 16-vertex model on the associated square lattice [17]. The Boltzmann weights of this vertex model are expressed in terms of the Boltzmann  00  00  11  11  00  00  11  11  000  000  000  111  111  111   00  00  00  11  11  11 00 00 00 11 11 11 FIG. 4: The checkerboard lattice and the associated square lattice with their bonds indicated by dashed lines.
weights of the checkerboard Ising model, as follows Generally, a vertex model is soluble if the vertex weights satisfy the free-fermion conditions so that the partition function is reducible to the S matrix of a many-fermion system [30]. In the present problem the free-fermion conditions are the following As can be easily verified, Eqs. (7) are identically satisfied by the Boltzmann weights Eqs. (6), for arbitrary values of K 1 and K 2 . Using Eqs. (6) for the 16-vertex model and calculating the free energy of the model [17,23] we obtain the critical condition for this system This equation has up to four critical lines depending on the values of J 1 and J 2 . For the whole phase diagram, the reader is referred to Ref. [23]. We show in Fig. 6 only the small region of J 2 /J 1 in the phase diagram which has two striking phenomena: the reentrant paramagnetic phase and a disorder line.
The reentrant phase is defined as a paramagnetic phase which is located between two ordered phases on the temperature axis as seen in the region −1 < J 2 /J 1 < −0.91: if we take for instance J 2 /J 1 = −0.94 and we go up on the temperature axis, we will pass through the ferromagnetic phase F, enter the "reentrant" paramagnetic phase, cross the disorder line, enter the partial disordered phase X where the central spins are free, and finally enter the paramagnetic phase P [23]. The reentrant paramagnetic phase takes place thus between a low-T ferromagnetic phase and a partially disordered phase. Note that in phase X, all central spins denoted by the number 5 in Fig. 1 are free to flip at all T while other spins are ordered up to the transition at T c . This result shows an example where order and disorder coexists in an equilibrium state.
It is important to note that though we get the exact solution for the critical surface, namely the exact location of the phase transition temperature in the space of parameters as shown in Fig. 6, we do not have the exact expression of the magnetization as a function of temperature. To verify the coexistence of order and disorder mentioned above we have to recourse to Monte Carlo simulations. This is easily done and the results for the order parameters and the susceptibility of one of them are shown in Fig. 7 for phases F and X at J 2 /J 1 = −0.94. As seen, the F phase disappears at T 1 ≃ 0.47 and phase IV (defined in Fig. 2) sets in at T 2 ≃ 0.50 and disappears for T > 1.14. T is measured in the unit of J 1 /k B . The paramagnetic region between T 1 and T 2 is the reentrant phase. Note that the disorder line discussed below cannot be seen by Monte Carlo simulations. Let us now give the equation of the disorder line shown in Fig. 6: Usually, one defines each point on the disorder line as the temperature where there is an effective reduction of dimensionality in such a way that physical quantities become simplified spectacularly. Along the disorder line, the partition function is zero-dimensional and the correlation functions behave as in one dimension (dimension reduction). The disorder line is very important in understanding the reentrance phenomenon. This type of line is necessary for the change of ordering from the high-T ordered phase to the low-T one. In the narrow reentrant paramagnetic region, preordering fluctuations with different symmetries exist near each critical line. Therefore the correlation functions change their behavior when crossing the "dividing line" as the temperature is varied in the reentrant paramagnetic region. On this dividing line, or disorder line, the system "forgets" one dimension in order to adjust itself to the symmetry of the other side. As a consequence of the change of symmetries there exist spins for which the two-point correlation function (between NN spins) has different signs, near the two critical lines , in the reentrant paramagnetic region. Hence it is reasonable to expect that it has to vanish at a disorder temperature T D . This point can be considered as a non-critical transition point which separates two different paramagnetic phases. The two-point correlation function defined above may be thought of as a non-local "disorder parameter". This particular point is just the one which has been called a disorder point by Stephenson [31] in analyzing the behavior of correlation functions for systems with competing interactions. Other models we solved have several disorder lines with dimension reduction [25,26] except the case of the centered square lattice where there is a disorder line without dimension reduction [27]. We believe that results of the exactly solved model in 2D shown above should also exist in three dimensions (3D), though we cannot exactly solve 3D models. To see this, we have studied a 3D version of the 2D Kagomé lattice which is a kind of body-centered lattice where the central spin in the lattice cell is free if the corner spins are in an antiferromagnetic order: the central spin has four up and four down neighbors making its energy zero as in the Kagomé lattice. We have shown that the partial disorder exists [32,33] and the reentrant zone between phase F and phase X in Fig. 6 closes up giving rise to a line of first-order transition [34].
To close this paragraph, we note that for other exactly solved frustrated models, the reader is referred to the review by Diep and Giacomini [35].

B. Elementary Excitations: Surface Magnons
We consider a thin film of N T layers with the Heisenberg quantum spin model. The Hamiltonian is written as where J ij is positive (ferromagnetic) and D ij > 0 denotes an exchange anisotropy. When D ij is very large with respect to J ij , the spins have an Ising-like behavior. For simplicity, let us suppose for the moment that all surface parameters are the same as the bulk ones with no defects and impurities. One of the microscopic mechanisms which govern thermodynamic properties of magnetic materials at low temperatures is the spin waves. The presence of a surface often causes spin-wave modes localized at and near the surface. These modes cause in turn a diminution of the surface magnetization and the magnetic transition temperature. The methods to calculate the spin-wave spectrum from simple to more complicated are (see examples given in Ref. [15]): (i) the equation of motion written for spin operators S ± i of spin S i occupying the lattice site i of a given layer. These operators are coupled to those of neighboring layers. Writing an equation of motion for each layer, one obtains a system of coupled equations. Performing the Fourier transform in the xy plane, one obtains the solution for the spin-wave spectrum.
(ii) the spin-wave theory using for example the Holstein-Primakoff spin operators for an expansion of the Hamiltonian. This is the second-quantization method. The harmonic spin-wave spectrum and nonlinear corrections can be obtained by diagonalizing the matrix written for operators of all layers.
(iii) the Green's function method using a correlation function between two spin operators. From this function one can deduce various thermodynamic quantities such as layer magnetizations and susceptibilities. The advantage of this method is one can calculate properties up to rather high temperatures. However, with increasing temperature one looses the precision.
We summarize briefly here the principle of the Green's function method for illustration (see details in Ref. [36,37]). We define one Green's function for each layer, numbering the surface as the first layer. We write next the equation of motion for each of the Green's functions. We obtain a system of coupled equations. We linearize these equations to reduce higher-order Green's functions by using the Tyablikov decoupling scheme [38]. We are then ready to make the Fourier transforms for all Green's functions in the xy planes. We obtain a system of equations in the space ( k xy , ω) where k xy is the wave vector parallel to the xy plane and ω is the spin-wave frequency (pulsation). Solving this system we obtain the Green's functions and ω as functions of k xy . Using the spectral theorem, we calculate the layer magnetization. Concretely, we define the following Green's function for two spins S i and S j as The equation of motion of G i,j (t, t ′ ) is written as where [...] is the boson commutator and ... the thermal average in the canonical ensemble defined as with β = 1/k B T . The commutator of the right-hand side of Eq. (12) generates functions of higher orders. In a first approximation, we can reduce these functions with the help of the Tyablikov decoupling [38] as follows We obtain then the same kind of Green's function defined in Eq. (11). As the system is translation invariant in the xy plane, we use the following Fourier transforms where ω is the magnon pulsation (frequency), k xy the wave vector parallel to the surface, R i the position of the spin at the site i, n and n ′ are respectively the indices of the planes to which i and j belong (n = 1 is the index of the surface). The integration on k xy is performed within the first Brillouin zone in the xy plane. Let ∆ be the surface of that zone. Equation (12) becomes ( ω − A n )g n,n ′ + B n (1 − δ n,1 )g n−1,n ′ + C n (1 − δ n,NT )g n+1,n ′ = 2δ n,n ′ < S z n > where the factors (1 − δ n,1 ) and (1 − δ n,NT ) are added to ensure that there are no C n and B n terms for the first and the last layer. The coefficients A n , B n and C n depend on the crystalline lattice of the film. We give here some examples: • Film of simple cubic lattice A n = −2J n < S z n > Cγ k + 2C(J n + D n ) < S z n > +2(J n,n+1 + D n,n+1 ) < S z n+1 > +2(J n,n−1 + D n,n−1 ) < S z n−1 > (17) B n = 2J n,n−1 < S z n > (18) where C = 4 and γ k = 1 2 [cos(k x a) + cos(k y a)]. • Film of body-centered cubic lattice where γ k = cos(k x a/2) cos(k y a/2) Writing Eq. (16) for n = 1, 2, ..., N T , we obtain a system of N T equations which can be put in a matrix form where u is a column matrix whose n-th element is 2δ n,n ′ < S z n >. For a given k xy the magnon dispersion relation ω( k xy ) can be obtained by solving the secular equation det|M| = 0. There are N T eigenvalues ω i (i = 1, ..., N T ) for each k xy . It is obvious that ω i depends on all S z n contained in the coefficients A n , B n and C n .
To calculate the thermal average of the magnetization of the layer n in the case where S = 1 2 , we use the following relation (see chapter 6 of Ref. [15]): where S − n S + n is given by the following spectral theorem ǫ being an infinitesimal positive constant. Equation (24) becomes where the Green's function g n,n is obtained by the solution of Eq. (23) g n,n = |M| n |M| (27) |M| n is the determinant obtained by replacing the n-th column of |M| by u.
To simplify the notations we put ω i = E i and ω = E in the following. By expressing we see that E i (i = 1, ..., N T ) are the poles of the Green's function. We can therefore rewrite g n,n as where f n (E i ) is given by Replacing Eq. (29) in Eq. (26) and making use of the following identity we obtain where n = 1, ..., N T . As < S z n > depends on the magnetizations of the neighboring layers via E i (i = 1, ..., N T ), we should solve by iteration the equations (32) written for all layers, namely for n = 1, ..., N T , to obtain the layer magnetizations at a given temperature T .
The critical temperature T c can be calculated in a self-consistent manner by iteration, letting all < S z n > tend to zero.
Let us show in Fig. 8 two examples of spin-wave spectrum, one without surface modes as in a simple cubic film and the other with surface localized modes as in body-centered cubic ferromagnetic case.
It is very important to note that acoustic surface localized spin waves lie below the bulk frequencies so that these low-lying energies will give larger integrands to the integral on the right-hand side of Eq. (32), making < S z n > to be smaller. The same effect explains the diminution of T c in thin films whenever low-lying surface spin waves exist in the spectrum. Figure 9 shows the results of the layer magnetizations for the first two layers in the films considered above with N T = 4.
Calculations for antiferromagnetic thin films with collinear spin configurations can be performed using Green's functions [36]. The physics is similar with strong effects of localized surface spin waves and a non-uniform spin contractions near the surface at zero temperature due to quantum fluctuations [37].

C. Frustrated Films
We showed above for a pedagogical purpose a detailed technique for using the Green's function method. In the case of frustrated thin films, the ground-state spin configurations are not only non collinear but also non uniform from the surface to the interior layers. In a class of helimagnets, the angle between neighboring spins is due to the competition between the NN and the NNN interactions. Bulk spin configurations of such helimagnets were discovered more than 50 years ago by Yoshimori [39] and Villain [40]. Some works have been done to investigate the low-temperature spin-wave behaviors [41][42][43] and the phase transition [44] in the bulk crystals. For surface effects in frustrated films, a number of our works have been recently done among which we can mention the case of a frustrated surface on a ferromagnetic substrate film [45], the fully frustrated antiferromagnetic facecentered cubic film [46], and very recently the helimagnetic thin films in zero field [47,48] and under an applied field [49].
The Green's function method for non collinear magnets has been developed for the bulk crystal [50]. We have extended this to the case of non collinear thin films in the works just mentioned. Since two spins S i and S j form an angle cos θ ij one can express the Hamiltonian in the local coordinates as follows [47]: The last term is an anisotropy added to facilitate a numerical convergence for ultra thin films at long-wave lengths since it is known that in 2D there is no ordering for isotropic Heisenberg spins at finite temperatures [51]. The determination of the angles in the ground state can be done either by minimizing the interaction energy with  respect to interaction parameters or by the so-called steepest descent method which has been proved to be very efficient [45,46]. Using their values, one can follow the different steps presented above for the collinear magnetic films, one then obtains a matrix which can be numerically diagonalized to get the spin-wave spectrum which is used in turn to calculate physical properties in the same manner as for the collinear case presented above. Let us show the case of a helimagnetic film. In the bulk, the turn angle in one direction is determined by the ratio between the antiferromagnetic NNN interaction J 2 (< 0) and the NN interaction J 1 . For the body-centered cubic lattice, one has cos θ = −J 1 /J 2 . The helimagnetic phase is stable for |J 2 |/J 1 > 1. Consider a film with the c axis perpendicular to the film surface. For simplicity, one supposes the turn angle along the c axis is due to J 2 . Because of the lack of neighbors, the spins on the surface and on the second layer have the turn angles strongly deviated from the bulk value [47]. The results calculated for various J 2 /J 1 are shown in Fig. 10 (right) for a film of N z = 8 layers. The values obtained are shown in Table I where one sees that the angles near the surface (2nd and 3rd columns) are very different from that of the bulk (last column).  Table 1: a strong rearrangement of spins near the surface is observed.
The spectrum at two temperatures is shown in Fig. 11 where the surface spin waves are indicated. The spin lengths at T = 0 of the different layers are shown in Fig. 12 as functions of J 2 /J 1 . When J 2 tends to -1, the spin configuration becomes ferromagnetic, the spin has the full length 1/2.
The layer magnetizations are shown in Fig. 13 where one notices the crossovers between them at low T . This is due to the competition between quantum fluctuations, which depends on the strength of antiferromagnetic interaction, and the thermal fluctuations which depends on the local coordinations.  As said earlier, Monte Carlo methods can be used in complicated systems where analytical methods cannot be efficiently used. Depending on the difficulty of the investigation, we should choose a suitable Monte Carlo technique. For a simple investigation to have a rough idea about physical properties of a given system, the standard Metropolis algorithm is sufficient [52,53]. It consists in calculating the energy E 1 of a spin, then changing its state and calculating its new energy E 2 . If E 2 < E 1 then the new state is accepted. If E 2 > E 1 the new state is accepted with a probability proportional to exp[−(E 2 − E 1 )/(k B T )]. One has to consider all spins of the system, and repeat the "update" over and over again with a large number of times to get thermal equilibrium before calculating statistical thermal averages of physical quantities such as energy, specific heat, magnetization, susceptibility, ... We need however more sophisticated methods if we wish to calculate critical exponents or to detect a first-order phase transition. For calculation of critical exponents, histogram techniques [54,55] are very precise: comparison with exact results shows an agreement often up to 3rd or 4th digit. To detect very weak first-order transitions, the Wang-Landau technique [56] combined with the finite-size scaling theory [57] is very efficient. We have used this technique to put an end to a 20-year-old controversy on the nature of the phase transition in Heisenberg and XY frustrated stacked triangular antiferromagnets [58,59].
To illustrate the efficiency of Monte Carlo simulations, let us show in Fig. 14 the layer magnetizations of the classical counterpart of the body-centered cubic helimagnetic film shown in section III C (figure taken from Ref. [47]). Though the surface magnetization is smaller than the magnetizations of interior layers, there is only a single phase transition. To see a surface transition, let us take the case of a frustrated surface of antiferromagnetic triangular lattice coated on a ferromagnetic film of the same lattice structure [45]. The in-plane surface interaction is J s < 0 and interior interaction is J > 0. This film has been shown to have a surface spin reconstruction as displayed in Fig. 15. We show an example where J s = −0.5J in Fig. 16. The left figure is from the Green's function method. As seen, the surface-layer magnetization is much smaller than the second-layer one. In addition there is a strong spin contraction at T = 0 for the surface layer. This is due to quantum fluctuations of the in-plane antiferromagnetic surface interaction J s . One sees that the surface becomes disordered at a temperature T 1 ≃ 0.2557 while the second layer remains ordered up to T 2 ≃ 1.522. Therefore, the system is partially disordered for temperatures between T 1 and T 2 . This result is very interesting because it confirms again the existence of the partial disorder in quantum spin systems observed earlier in the bulk [32,33]. Note that between T 1 and T 2 , the ordering of the second layer acts as an external field on the first layer, inducing therefore a small value of its magnetization. Results of Monte Carlo simulations of the classical model are shown on the right of Fig. 16 which have the same features as the quantum case. To close this paragraph we mention that the question of surface criticality has been a long-standing debate. On the one hand, pure theories would tell us that as long as the thickness is finite the correlation in this direction can be renormalized so that the nature of a phase transition in a thin film should be that of the corresponding 2D model. On the other hand, experimental observations and numerical simulations show deviations of critical exponents from 2D universality classes. The reader is referred to Refs. [60,61] for discussions on this subject.

E. Surface Reorientation
In this paragraph, we would like to show the case of thin films where a transition from an in-plane spin configuration to a perpendicular spin configuration is possible at a finite temperature. Such a reorientation occurs when there is a competition between a dipolar interaction which tends to align the spins in the film surface and a perpendicular anisotropy which is known to exist when the film thickness is very small [10,11]. Experimentally, it has been observed in a thin Fe film deposited on a Cu(100) substrate that the perpendicular spin configuration at low temperatures undergoes a transition to a planar spin configuration as the temperature (T ) is increased [62][63][64][65]. Theoretically, this problem has been studied by many people [66][67][68][69][70]. Let us consider a 2D surface for simplicity. The case of a film with a small thickness has very similar results [70]. The Hamiltonian includes three parts: a 6-state Potts model H p , a dipolar interaction H d and a perpendicular anisotropy H a : where σ i is a variable associated to the lattice site i. σ i is equal to 1, 2, 3, 4, 5 and 6 if the spin at that site lies along the ±x, ±y and ±z axes, respectively. J ij is the exchange interaction between NN at i and j. We will assume that (i) J ij = J s if i and j are on the same film surface (ii) J ij = J otherwise. For the dipolar interaction, we write where r i,j is the vector of modulus r i,j connecting the site i to the site j. One has r i,j ≡ r j − r i . S(σ i ) is a vector of modulus 1 pointing in the x direction if σ i = 1, in the −x direction if σ i = 2, etc. The perpendicular anisotropy is where A is a constant.
Using Monte Carlo simulations, we have established the phase diagram shown in Fig. 17 for two dipolar cutoff distances. Several remarks are in order: (i) in a small region above D = 0.1 (left figure) there is a transition from the in-plane to the perpendicular configuration when T increases from 0, (ii) this reorientation is a very strong first-order transition: the energy and magnetization are discontinuous (not shown) at the transition. Comparing to the Kagomé Ising case shown in Fig. 6  a first-order transition as also observed near phase frontiers in other systems such as in a frustrated body-centered cubic lattice [34].
To conclude this subsection, we mention that competing interactions determine frontiers between phases of different symmetries. Near these frontiers, we have seen many interesting phenomena such as reentrance, disorder lines, reorientation transition, ... when the temperature increases.

F. Spin Transport in Thin Films
The total resistivity stem from different kinds of diffusion processes in a crystal. Each contribution has in general a different temperature dependence. Let us summarize the most important contributions to the total resistivity ρ t (T ) at low temperatures in the following expression where A 1 , A 2 and A 3 are constants. The first term is T -independent, the second term proportional to T 2 represents the scattering of itinerant spins at low T by lattice spin-waves. Note that the resistivity caused by a Fermi liquid is also proportional to T 2 . The T 5 term corresponds to a low-T resistivity in metals. This is due to the scattering of itinerant electrons by phonons. At high T , metals however show a linear-T dependence. The logarithm term is the resistivity due to the quantum Kondo effect caused by a magnetic impurity at very low T . We are interested here in the spin resistivity ρ of magnetic materials. We have developed an algorithm which allows us to calculate the spin resistivity in various magnetically ordered systems [71-76, 76, 77, 77], in particular in thin films. Unlike the charge conductivity, studies of spin transport have been regular but not intensive until recently. The situation changes when the electron spin begins to play a central role in spin electronics, in particular with the discovery of the colossal magnetoresistance [2,3].
The main mechanism which governs the spin transport is the interaction between itinerant electron spins and localized spins of the lattice ions (s − d model). The spin-spin correlation has been shown to be responsible for the behavior of the spin resistivity [80][81][82]. Calculations were mostly done by mean-field approximation (see references in [83]). Our works mentioned above were the first to use intensive Monte Carlo simulations for investigating the spin transport. We also used a combination of the Boltzmann equation [15,84] and numerical data obtained by simulations [73,74]. The Hamiltonian includes three main terms: interaction between lattice spins H l , interaction between itinerant spins and lattice spins H r , and interaction between itinerant spins H m . We suppose where S i is the spin localized at lattice site i of Ising, XY or Heisenberg model, J i,j the exchange integral between the spin pair S i and S j which is not limited to the interaction between nearest-neighbors (NN). For H r we write where σ i is the spin of the i-th itinerant electron and I i,j denotes the interaction which depends on the distance between electron i and spin S j at lattice site j. For simplicity, we suppose the following interaction expression where r ij = | r i − r j |, I 0 and α are constants. We use a cut-off distance D 1 for the above interaction. Finally, for H m , we use where with K i,j being the interaction between electrons i and j. The system is under an electrical field ǫ which creates an electron current in one direction. In addition we include also a chemical potential which keeps electrons uniformly distributed in the system. Simulations have been carried out with the above Hamiltonian. The reader is referred to the original papers mentioned above for technical details. As expected, the spin resistivity reflects the spin-spin correlation of the system: ρ has the form of the magnetic susceptibility, namely it shows a peak at the phase transition. It is noted however that unlike the susceptibility which diverges at the transition, the spin resistivity is finite at T c due to the fact that only short-range correlations can affect the resistivity (see arguments given in [81]). Moreover, it is known that near the phase transition the system is in the critical-slowing-down regime. Therefore, care should be taken while simulating in the transition region. This point has been considered in our simulations by introducing the relaxation time [78].
To illustrate the efficiency of Monte Carlo simulations, we show in Fig. 18 the excellent agreement of our simulation [79] and experiments performed on MnTe [85]. The interactions and the crystalline parameters were taken from Ref. [86]. When a film has a surface phase transition at a low temperature in addition to the transition of the bulk at a higher temperature, one observes two peaks in the spin resistivity as shown in Ref. [73].

IV. CONCLUSION
To conclude this review let us discuss on the relation between theories and experiments, in particular on the difficulties encountered when one is confronted, on the one hand, with simplified theoretical pictures and hypotheses and, on the other hand, with insufficient experimental knowledge of what is really inside the material. We would like to emphasize on the importance of a sufficient theoretical background to understand experimental data measured on systems which are more complicated, less perfect than models used to describe them. Real systems have always impurities, defects, disorder, domains, ... However, as long as these imperfections are at extremely small amounts, they will not affect observed macroscopic quantities: theory tells us that each observable is a result from a statistical average over all microscopic states and over the space occupied by the material. Such an averaging will erase away rare events leaving only common characteristics of the system. Essential aspects can be thus understood from simple models if one includes correct ingredients based on physical arguments while constructing the model. One of the striking points shown above is the fact that without sophisticated calculations, we cannot discover hidden secrets of the nature such as the existence of disorder lines with and without dimension reduction, the extremely narrow reentrant region between two ordered phases, the coexistence of order and disorder of a system at equilibrium etc. These effects are from the competition between various interactions which are unavoidable in real materials. These interactions determine the boundaries between various phases of different symmetries in the space of physical parameters. Crossing a boundary the system will change its symmetry. Theory tells us that if the symmetry of one phase is not a subgroup of the other then the transition should be of first order or the two phases should be separated by a narrow reentrant phase. Without such a knowledge, we may overlook such fine effects while examining experimental data.
We have used frustrated thin films to illustrate various effects due to a combination of frustration and surface magnetism. We have seen that to understand when and why surface magnetization is small with respect to the bulk one we have to go through a microscopic mechanism to recognize that low-lying localized surface spin-wave modes when integrated in the calculation of the magnetization will indeed mathematically lower its value. Common effects observed in thin films such as surface reconstruction and surface disordering can be theoretically explained.
We should emphasize on the importance of a combination of Monte Carlo simulation and theory. We have seen for example that we can determine the exact transition temperature in the Kagomé model but the nature of the ordering is understood only with Monte Carlo simulation. This is not the only example: as theorists we need also hints and checks while constructing a theory, and as a simulator we need to understand what comes out from the computer since we often loose the track between inputs and outputs in a simulation. In this respect, a combination of numerical and theoretical methods is unavoidable.
I am grateful to my numerous former and present doctorate students with whom I shared uncountable wonderful moments of scientific discussions and from whom I draw my source of energy and joy.