Metadynamics study of the temperature dependence of magnetic anisotropy and spin-reorientation transitions in ultrathin films

TL;DR

This study uses metadynamics simulations to explore how magnetic anisotropy and spin orientations in ultrathin ferromagnetic films change with temperature, providing insights into phase transitions and magnetic behavior.

Contribution

It introduces a systematic simulation approach combining first-principles parameters to analyze temperature-dependent magnetic phenomena in ultrathin films.

Findings

Reproduces the power-law behavior of magnetic anisotropy energy.

Identifies spin-reorientation transitions consistent with experimental data.

Demonstrates the effectiveness of metadynamics in magnetic studies.

Abstract

We employ metadynamics simulations to calculate the free energy landscape of thin ferromagnetic films and perform a systematic study of the temperature dependence of magnetic anisotropy and of the spin-reorientation transitions. By using a simple spin model we recover the well-known power-law behavior of the magnetic anisotropy energy against magnetization and present a rather detailed analysis of the spin-reorientation transitions in ultrathin films. Based on tensorial exchange interactions and anisotropy parameters derived from first-principles calculations we perform simulations for Fe double layers deposited on Au(001) and W(110). In case of Fe$_2$W(110) our simulations display an out-of-plane to in-plane spin-reorientation transition in agreement with experiments.

Figures

Tables

Table 1. Table 1: Calculated layer dependent magnetic anisotropy parameters (in units of meV) for the Fe 2 Au(001) layers. The Fe layer at the interface with Au is denoted by I, the one at the surface by S. Negative/positive values of the anisotropies prefer in-plane/normal-to-plane orientation of the magnetization.

layer	$\lambda$	$J_p^{zz}-J_p^{xx}$
I	$-0.097$	$0.181$
S	$0.360$	$-0.314$

Table 2. Table 2: Calculated layer dependent magnetic anisotropy parameters (in units of meV) for the Fe 2 W(110) layers. The Fe layer at the interface with W is denoted by I, the one at the surface by S. The notations x , y and z stand for the ( 1 ​ 1 ¯ ​ 0 ) 1 ¯ 1 0 (1\overline{1}0) , ( 001 ) 001 (001) and ( 110 ) 110 (110) directions, respectively.

layer	${\lambda }_x$	${\lambda }_y$	$J_p^{zz}-J_p^{xx}$	$J_p^{zz}-J_p^{yy}$
I	0.611	0.261	– 0.603	0.138
S	– 0.055	– 0.137	0.377	0.106

Equations

$$H=-\frac{1}{2}\sum _{⟨i,j⟩}\left(J{\mathbf{s}}_i{\mathbf{s}}_j-d{s}_{zi}{s}_{zj}\right)-\sum _i{\lambda }_i{s}_{zi}^2,$$

$$V_{\mathrm{b}}(\eta ,t+\tau )$$

$$V_{\mathrm{G}}(\eta -{\eta }_{\mathrm{act}})$$

$$P({\mathbf{s}}_i\to {\mathbf{s}}_i^{\prime })=\text{mbox}min\left\{1,e^{-\beta \left[E({\mathbf{s}}^{\prime })+V_{\mathrm{b}}(\eta ({\mathbf{s}}^{\prime }))-E(\mathbf{s})-V_{\mathrm{b}}(\eta (\mathbf{s}))\right]}\right\},$$

$$V_{\mathrm{b}}(\eta ,t+\tau )$$

$$E=$$

$${\left|\frac{{\partial }^{2}E}{\partial {\vartheta }_i{\vartheta }_j}\right|}_{{\vartheta }_i=0,\pi /2}=0,$$

$$\left(J+d-\frac{{\lambda }_1^{\mathrm{l}}}{2}\right)\left(J+d-\frac{{\lambda }_2^{\mathrm{l}}}{2}\right)=(d-J{)}^2$$

$$\left(J-2d+\frac{{\lambda }_1^{\mathrm{u}}}{2}\right)\left(J-2d+\frac{{\lambda }_2^{\mathrm{u}}}{2}\right)=J^2.$$

$$F=\frac{4J}{2}{\mathbf{M}}^2-\frac{4d}{2}{M^z}^2-k_{\mathrm{B}}T\ln (Z_1)-k_{\mathrm{B}}T\ln (Z_2)$$

$$Z_i$$

$$k_{\mathrm{B}}{T}_z$$

$$k_{\mathrm{B}}{T}_x$$

$$H=-\frac{1}{2}\sum _{p,q=1}^2\sum _{i,j}{\mathbf{s}}_{pi}^T{\mathcal{J}}_{pi,qj}{\mathbf{s}}_{qj}-\sum _{p=1}^2\sum _i{\lambda }_p({\mathbf{s}}_{pi}\hat{\mathbf{z}}{)}^2,$$

$${\mathcal{J}}_p=\frac{1}{2}\sum _{q=1,2}\sum _j{\mathcal{J}}_{p0,qj},$$

$$H$$

Full text

I introduction

Since Néel’s seminal paper in 1954 Néel (1954) considerable interest has been focused on the magnetism of thin films and multilayers. Magnetic anisotropy plays a key role in several phenomena important for technological applications. In a magnetic data storage device the information is stored by controlling the magnetic orientation of a small magnetic domain that is retained by magnetic anisotropy. In the early implementation of magnetic recording the magnetization of the bits were parallel with the plane of the film. Application of materials with perpendicular magnetic anisotropy (PMA) triggered an order of magnitude increase of the storage density. The first realization of a perpendicular magnetic recording occurred more than a decade ago ichi Iwasaki (2009), recent reviews on PMA can be found in Refs. Tudu and Tiwari (2017); Dieny and Chshiev (2017). In order to further increase the storage density, the grain size in the recording medium should be decreased which requires a high uniaxial magnetic anisotropy of the thin film. Due to the large magnetic anisotropy of the recording media the field produced by the write head might no longer be sufficient to overcome the barrier to switch the magnetization. To circumvent this issue a heat assisted magnetic recording (HAMR) is proposed Hamann et al. (2004); McDaniel (2005); Stipe et al. (2010). In HAMR the magnetic anisotropy is decreased by temporarily heating the domain storing the information.

The temperature dependence of magnetic anisotropy of thin films has been investigated both experimentally Okamoto et al. (2002); Thiele et al. (2002); Vaz et al. (2008); Fu et al. (2016) and theoretically Staunton et al. (2004); Mryasov et al. (2005); Buruzs et al. (2007); Asselin et al. (2010). The magnetic anisotropy energy (MAE) at finite temperature is usually defined as the difference between the free energy of the in-plane and that of the normal-to-plane magnetized system. Magnetic simulations provide different tools for sampling the complex free-energy surfaces. One branch of such schemes is formed by the adaptive biasing potential methods such as the Wang–Landau algorithm Wang and Landau (2001), umbrella sampling Procacci and Schettino (2006) and metadynamics Laio and Parrinello (2002). In metadynamics a biasing potential is constructed as a sum of Gaussians centered along the trajectory in the space of the collective variables Laio and Parrinello (2002). In well-tempered metadynamics the smooth convergence of the biasing potential is guaranteed by changing adaptively the height of the Gaussians Barducci et al. (2008). This algorithm is proved to converge to the exact free energy Dama et al. (2014).

In this work we perform a systematic study of the temperature dependence of magnetic anisotropy and spin-reorientation transitions (SRT) by using metadynamics. In Section II we outline the main features of metadynamics simulations with the aim at studying the free energy landscape of a thin ferromagnetic film. In Section III we first present a model study of the temperature dependence of magnetic anisotropy and a rather detailed analysis of the SRT in ultrathin films. Based on tensorial exchange interactions and anisotropy parameters derived from first-principle calculations we then present simulations on Fe bilayers deposited on Au(001) and W(110) and, finally, we summarize our results.

II Details of the metadynamics simulations

The magnetic properties of thin films of transition metals are often described by classical spin models Nowak (2007). In most part of this work we choose a simple Heisenberg model to describe the magnetic properties of an ultrathin films with uniaxial anisotropy and anisotropic exchange interactions:

[TABLE]

where $\mathbf{s}_{i}$ is a unit vector representing the direction of the atomic magnetic moment at site $i$, only nearest neighbors are considered in the first sum on the right-hand side with isotropic exchange coupling $J$ and an anisotropic part $d$, while $\lambda_{i}$ are the uniaxial anisotropy constants. More complex spin models will be presented and used only in Sections III.C and D in context of Fe bilayers on Au(001) and W(110).

The free energy is sampled along an appropriately chosen collective variable (CV) labeled by $\eta$. For our present study we chose the $z$ (normal-to-plane) component of the normalized magnetization, $\eta=M_{z}/M$, as the collective variable, where $M_{z}=\sum_{i}s_{zi}$ and $M=|\sum_{i}\mathbf{s}_{i}|$. The key quantity in metadynamics is the bias potential $V_{\mathrm{b}}(\eta)$ added to the energy of the system. Although in most of its applications metadynamics is implemented in molecular dynamics, there are examples where it is successfully used in Monte Carlo simulations Marini et al. (2008); Crespo et al. (2010) as well. If a Monte Carlo step (MCS) is interpreted as a time step the bias potential will be time dependent as well. After every $\tau$ MCS a Gaussan potential centered at the actual value of the CV, $\eta_{\mathrm{act}}$, is added to the bias potential:

[TABLE]

where $\sigma$ and $w$ are the width and the height of the Gaussian, respectively. In well-tempered metadynamics Barducci et al. (2008) the height of the Gaussian is chosen to change with the time $\tau$.

In our metadynamics simulations we applied a simple Metropolis algorithm Metropolis et al. (1953) with the probability of a random change of the spin at site $i$ $\mathbf{s}_{i}\rightarrow\mathbf{s}_{i}^{\prime}$,

[TABLE]

where $\beta$ is the inverse temperature and $E(\mathbf{s})$ is the energy of the spin configuration given by Eq. (1). After a predefined number of Monte Carlo steps the biasing potential is updated by adding a Gaussian centered at the actual value of the CV with the height of $w=w_{0}e^{-\frac{V_{\mathrm{b}}(\eta)}{k_{\mathrm{B}}T_{\mathrm{m}}}}$, where $T_{\mathrm{m}}$ is an appropriately chosen temperature as it is explained in the procedure of well-tempered metadynamics Barducci et al. (2008). In equilibrium, i.e. when the bias potential becomes stationary, the free energy $F(T)$ of the system is identified with the negative of the bias potential, $F(T)=-V_{\mathrm{b}}(\eta(T))$, where $\eta(T)$ stands for the equilibrium value of the CV Dama et al. (2014).

The values of the CV chosen for our model must be within the interval $[-1,1]$ and the free energy has a discontinuity at the boundaries which can not be accurately reproduced by a sum of finite-width Gaussians as it is detailed in Refs. Laio and Gervasio (2008); Crespo et al. (2010). In order to eliminate this problem, the procedure proposed by Crespo et al. Crespo et al. (2010) has been modified in the following manner. Whenever the bias potential is updated, an extra Gaussian with the same width and height is added out of the physically relevant interval of the CV:

[TABLE]

where $V_{\mathrm{G}}(\eta)$ is the Gaussian potential given Eq. (3). This scheme clearly makes the bias potential continuous at $\eta=\pm 1$. It should be noted that $V_{\mathrm{b}}(\eta)$ does not go smoothly to zero in the nonphysical region, but this part of the CV is never sampled during the simulation. In order to explore the free energy surface along the CV multiple walkers metadynamics Raiteri et al. (2006) was applied. The simulations were done simultaneously on typically four replicas each contributing equally to the growth of a joint bias potential.

III Results and discussions

III.1 Temperature dependence of the magnetic anisotropy energy

In order to validate metadynamics for the study of finite temperature magnetism we first investigated the temperature dependence of the MAE, $K(T)$ of a monolayer. In case of on-site uniaxial anisotropy, the magnetic anisotropy energy should exhibit a $K(T)\propto M^{3}(T)$ scaling at low temperature as predicted by Callen and Callen Callen and Callen (1966). If the magnetic anisotropy comes also from the exchange coupling, the low temperature behaviour of the MAE will be similar to the case of on-site uniaxial anisotropy, but at higher temperature the exponent will differ from three Staunton et al. (2004); Mryasov et al. (2005).

The isotropic exchange couplings in ferromagnetic systems are closely related to the Curie temperatures. For ultrathin films of transition metals the Curie temperature is few hundreds of K and the corresponding effective exchange coupling $J$ is few tens of meV. The uniaxial anisotropy constant for hcp Co is 70 $\mu$eV Baberschke (2001) and for a broad scale of thin films on different substrates it is in the range of 10-200 $\mu$eV Johnson et al. (1996). According to the above experimental values, as compared to the effective isotropic coupling, the uniaxial anisotropy constant ($\lambda/J$) and the anisotropy of the exchange coupling ($d/J$) have been chosen between $0.001$ and $0.01$ for the subsequent simulations.

The first simulation was performed for a monolayer containing ferromagnetic nearest neighbour exchange coupling and uniaxial on-site anisotropy with easy axis perpendicular to the plane, $(d=0,\;\lambda>0)$, see Eq. (1). The ground state of the system is ferromagnetic with a normal-to-plane orientation. The bias potential has a quadratic dependence on the CV as it is shown in the inset of Fig. 1. This parabolic behaviour is retained in the whole temperature range below the paramagnetic phase transition. The free energy has a maximum at $\eta=0$ referring to the in-plane configuration and it has minima at $\eta=\pm 1$ representing out-of-plane magnetic orientations. The difference between these two extrema is defined as the MAE. Numerically more efficiently, $K(T)$ can be obtained as the second order coefficient of a symmetric parabola fitted to the bias potential. This is plotted in Fig. 1. As the temperature is increasing the curvature of the free energy (bias potential) as a function of CV is gradually decreasing and it tends to zero above the Curie temperature. The Curie temperature is identified as the temperature corresponding to the maximum of the specific heat. Although the Curie temperature scales with the system size, it should be chosen compatible with the size of the system for which the MAE is calculated.

The magnetic anisotropy in Fig. 1 is almost linearly decreasing with the temperature similarly to the results obtained by using constrained Monte Carlo simulations Asselin et al. (2010) for uniaxial anisotropy. The non-zero value of the magnetic anisotropy above the Curie temperature is the consequence of the finite size of the system. In Fig. 2 the MAE is plotted against the magnetization on a log-log mesh. As can be seen, at low temperatures the results show excellent agreement with the scaling behavior predicted by Callen and Callen.

If the uniaxial anisotropy $\lambda$ is removed from the model, Eq. (1), and anisotropic exchange $d<0$ is introduced, the scaling behaviour of the anisotropy energy will be different as shown in Fig. 3. At low temperatures the system behaves as in the case of uniaxial on-site anisotropy, but at higher temperatures the exponent $\gamma$ in the relationship $K(T)\propto M(T)^{\gamma}$ changes from three to two. Such a behavior of the temperature dependence of the MAE was explored in earlier experimental Okamoto et al. (2002) and theoretical studies Staunton et al. (2004); Mryasov et al. (2005); Deák et al. (2014) for FePt alloys.

III.2 Spin-reorientation transitions

The interplay of different type of anisotropies often leads to a reorientation of the magnetization direction. The temperature driven spin-reorientation transition in thin films is usually explained by the competition of the uniaxial on-site anisotropy and the shape anisotropy Pappas et al. (1990); Fruchart et al. (1997); Sellmann et al. (2000); Ślęzak et al. (2010). For planar systems the shape anisotropy due to the magnetic dipolar interaction always prefers in-plane magnetization, while the on-site anisotropy of a magnetic overlayer frequently prefers a normal-to-plane orientation. The shape anisotropy due to the anisotropic exchange interaction which is the consequence of the spin-orbit coupling may also prefer both directions.

In the model given in Eq. (1) the two competing anisotropies are the on-site uniaxial anisotropies $\lambda_{i}$ and the anisotropy of the exchange coupling $d$. Considering a single square lattice, in case of $\lambda>2d$ the ground state is a normal-to-plane ferromagnetic. If $\lambda-2d$ is not too large, a temperature induced normal-to-plane to in-plane SRT can occur. In the inset of Fig. 4 the bias potentials for a monolayer with $\lambda=0.05375\,J$ and $d=0.025\,J$ are shown for different temperatures as obtained from metadynamics simulations. At low temperatures the maxima of the bias potentials — the minima of the free energy — correspond to $\eta=M_{\mathrm{z}}=\pm 1$, i.e. to a normal-to-plane configuration. As the temperature is increasing the curvature of the bias potential changes sign and the minimum of the free energy moves to $\eta=M_{\mathrm{z}}=0$, i.e. to in-plane magnetic orientation. The magnetic anisotropy energy $K(T)$ in Fig. 4 is zero at the transition temperature $T_{\mathrm{r}}$. It is worthwhile to mention that if the magnetization turns into the plane the system will have a gap-less magnetic excitation spectrum and long range magnetic order will no longer exist according to the Mermin-Wagner theorem. However, the magnetic anisotropy energy can still be defined as the free-energy difference between the normal-to-plane and in-plane magnetic orientations. The bias potential $V_{\text{b}}$ shown in the inset of Fig. 4 demonstrates a first order phase transition. Moschel and Usadel Moschel and Usadel (1995) using MC simulations and Fridman et al. Yu. A. Fridman (2002) applying a Hubbard-operator technique also confirmed that a monolayer exhibits first order SRT.

In a case of a bilayer our simple model results in a more feature-rich phase diagram where both first order and second order SRT can occur. A mean-field analysis of a very similar model has been performed almost two decades ago Udvardi et al. (2001) and here we recall some of the results of this study. As a model system we consider a bilayer on an fcc(001) surface with nearest neighbor interactions $J$ and $d$, and on-site anisotropy parameters $\lambda_{1}$ and $\lambda_{2}$. At zero temperature supposing uniform magnetization within each monolayer, but different orientations in the two monolayers the energy of the system can be written as:

[TABLE]

where $\vartheta_{i}$ is the polar angle with respect to the axis perpendicular to the layers ($z$). In the case of uniform in-plane $(\vartheta_{1}=\vartheta_{2}=\pi/2)$ and a normal-to-plane $(\vartheta_{1}=\vartheta_{2}=0)$ orientations the energy has an extremum. The energies of these two particular configurations coincide if $4d=\lambda_{1}+\lambda_{2}$ defining a line in the $\{\lambda_{1},\lambda_{2}\}$ parameter space. In the vicinity of this line a canted magnetic configuration exists. The boundaries of the region of the canted states can be obtained from the stability condition:

[TABLE]

yielding the lower boundary line,

[TABLE]

and the upper boundary line,

[TABLE]

Below the line given by Eq. (10), $\lambda_{1}+\lambda_{2}<\lambda^{\mathrm{l}}_{1}+\lambda^{\mathrm{l}}_{2}$, the ground state is in-plane ferromagnetic and above the line given by Eq. (11), $\lambda_{1}+\lambda_{2}>\lambda^{\mathrm{u}}_{1}+\lambda^{\mathrm{u}}_{2}$, it is normal-to-plane ferromagnetic.

At finite temperature the mean-field free energy of the double layer can be expressed as:

[TABLE]

where

[TABLE]

$I_{0}(x)$ is the modified Bessel function of the first kind, while $M_{x}$ and $M_{z}$ are the $x$ and $z$ component of the magnetization of the bilayer, respectively. As was shown in Ref. Udvardi et al. (2001) the magnetization can go to zero either via an in-plane or via a normal-to-plane direction at temperatures, $T_{x}$ and $T_{z}$, respectively, the higher of which can obviously be associated with the mean-field estimation of the Curie temperature $T_{\mathrm{C}}$. Minimizing the free-energy with respect to the magnetization of the system with the constraint $M_{z}=0$ or $M_{x}=0$ and using a high temperature expansion yields the following expressions for $T_{x}$ and $T_{z}$ to first order in $\lambda_{1}$ and $\lambda_{2}$:

[TABLE]

An out-of-plane to in-plane SRT can occur only when the ground-state magnetization is out of plane and $T_{z}<T_{x}=T_{\mathrm{C}}$. In the case of a reversed SRT the ground state magnetization has to be in-plane (or canted) and $T_{x}<T_{z}=T_{\mathrm{C}}$. In the parameter space $\{\lambda_{1},\lambda_{2}\}$ the region where SRT can occur are bounded by the line defined by Eq. (10) and by the line where $T_{x}=T_{z}$: $\lambda_{1}+\lambda_{2}=\tfrac{40}{3}d$.

We performed metadynamics Monte Carlo simulations to explore the phase diagram of a model bilayer. Although the anisotropy parameters $\lambda_{i}$ and $d$ can take both positive and negative values, in order to keep the MC simulations tractable, our investigations were restricted to the positive quarter of the parameter space $\{\lambda_{1}/d,\lambda_{2}/d\}$. The phase diagram for $d=0.005\,J$ is shown in Fig. 5. In this case, the region where canted ground states exist determined by Eqs. (10) and (11) is extremely narrow. The area where a normal-to-plane to in-plane SRT occurs provided by the metadynamics simulations (colored region) is considerably narrower then the corresponding area predicted by the mean field theory (bounded by the two solid red lines). The coloring clearly demonstrates that the reorientation temperature $T_{\mathrm{r}}$ gradually approaches the Curie temperature as the uniaxial anisotropy constants are increasing, while parallel to the lines $\lambda_{1}+\lambda_{2}=\mathrm{const.}$ it is almost constant. If the uniaxial anisotropy is further increased the system keeps its normal-to-plane ferromagnetic order till the ferromagnetic-paramagnetic phase transition.

Increasing the two-site anisotropy $d$ the area of canted ground states on the phase diagram becomes wider. In the case of $d=0.05\,J$, the lower and upper boundary of the canted region are indicated by the dotted and dashed lines in Fig. 5, respectively. For further investigations we choose two points in the phase diagram: A ($\lambda_{1}/d=8.66$, $\lambda_{2}/d=0$) representing a canted ground state, however, lying in the vicinity of the upper boundary line of this region (dashed line in Fig. 5) and B ($\lambda_{1}/d=\lambda_{2}/d=4.33$) corresponding to a normal-to-plane ferromagnetic ground state. For the first choice of ($\lambda_{1},\lambda_{2}$) the magnetization of the system continuously turns into the plane as the temperature is increasing and, considering the normal-to-plane component of the magnetization as order parameter, the system undergoes a second order SRT. This is demonstrated in Fig. 6, where the bias potentials of a $2\times 64\times 64$ lattice are shown as the function of the CV at different temperatures close to the SRT. Below the reorientation temperature, $T_{\mathrm{r}}/T_{\mathrm{C}}\sim 0.45$, the magnitude of the maximum position of the bias potential (minimum position of the free energy), $\eta_{\mathrm{max}}$, decreases continuously with increasing temperature, while at the in-plane magnetization $\eta=0$ there is a minimum in the bias potential. Above the reorientation transition temperature the in-plane configuration belongs to the maximum of the bias potential (minimum of the free energy), that means the order parameter is identical to zero.

If the uniaxial anisotropy parameters are the same for both layers, $\lambda_{1}=\lambda_{2}$, no canted ground state exists for the bilayer, therefore, the mean-field description of temperature dependent magnetism is analogous with that of the monolayer. The results of metadynamics simulations show, however, some different features for the bilayer and the monolayer. According to Fig. 4 the SRT for the monolayer is discontinuous and the normal-to-plane and in-plane phases can not coexist. The bias potentials for the bilayer with anisotropy parameters $\lambda_{1}=\lambda_{2}=4.33\,d$ are shown in Fig. 7. Below the reorientation temperature the bias potential has maxima at $\eta=\pm 1$ which correspond to a normal-to-plane average magnetization. As the temperature is increasing a local maximum of the bias potential evolves at $\eta=0$ referring to in-plane magnetization. Further increasing the temperature the local maximum at $\eta=0$ becomes the global maximum. The spin-reorientation transition is, therefore, of first order as in the case of the monolayer but the phases with in-plane and normal-to-plane magnetization can coexist.

III.3 Fe2Au(001)

Over the past three decades, thin iron films deposited on the surface of gold have been the subject of extensive investigations, especially in context of low-dimensional magnetism, see e.g. Ref. Wilgocka-Ślęzak et al. (2010) and references therein. An Fe monolayer grown on Au (001) has often been referred as a prototypical two-dimensional ferromagnet. The film FenAu(001) exhibits a normal-to-plane magnetic ground state for $n\leq 2$ and they undergo a thickness driven spin reorientation when the thickness of the Fe film reaches three monolayers Wilgocka-Ślęzak et al. (2010). While the driving force of this spin reorientation is the magnetostatic shape anisotropy, it is worth to study the temperature dependence of the spin-orbit induced MAE by using the metadynamics simulations introduced in this work. In this Section we present such a study for Fe2Au(001).

For the simulations we used the following spin Hamiltonian,

[TABLE]

where $p$ and $q$ denote layers, $i$ and $j$ stand for Fe atoms within each layer, $\hat{\mathbf{z}}$ is a unit vector parallel to the $z$ axis, the $\mathcal{J}_{ij}$ is a $3\times 3$ matrix of exchange interactions and the sum in the first term is not restricted to the nearest neighbours only. The trace of the tensor $\mathcal{J}_{ij}$ can be identified as three times the isotropic exchange coupling $J_{ij}$, while the symmetric and anti-symmetric part of the tensor correspond to the pseudo-dipolar and Dzyaloshinsky-Moriya interactions, respectively Udvardi et al. (2003). In order to determine the exchange tensors we applied the relativistic extension of the torque method Udvardi et al. (2003) implemented in the framework of the Screened Korringa-Kohn-Rostoker (SKKR) method Szunyogh et al. (1995). Since the (001) surface of fcc Au fits almost perfectly to the (001) surface of the bcc Fe (the lattice mismatch is less than 0.6 %) we used two-dimensional translational symmetry for the whole system using the lattice constant of Au (2.87 Å). The Fe-Fe inter-layer distance has been chosen to be the same as the bulk value (1.44 Å) and the Fe-Au inter-layer distance was 1.6 Å.

The calculated spin model parameters were then used in Monte Carlo and metadynamics simulations. In order to reduce finite size effects, the Curie temperature of the system has been determined from the intersection of the Binder cumulants yielding $T_{\mathrm{C}}\simeq 380\,\mathrm{K}$ in good agreement with the experiments Dürr et al. (1989). In order to characterize the anisotropy of the exchange tensors the lattice sum of the exchange couplings has been introduced:

[TABLE]

where $\mathcal{J}_{p0,qj}$ is the coupling tensor between an arbitrary site [math] in layer $p$ and site $j$ in the layer $q$. Due to the C4v symmetry of the lattice $\mathcal{J}_{p}$ is a diagonal matrix with identical $J_{p}^{xx}$ and $J_{p}^{yy}$ elements.

The layer dependent uniaxial anisotropy constants $\lambda_{p}$ and the anisotropy of exchange couplings $J_{p}^{zz}-J_{p}^{xx}$ are summarized in Table 1. Interestingly, the on-site anisotropies and the exchange anisotropies have different signs in both the interface (I) and the surface (S) layer, and they also change sign between the two layers. Nevertheless, in both layers the positive contributions dominate, resulting in an overall normal-to-plane magnetic ground state for the bilayer. The temperature dependent MAE obtained from metadynamics simulations is plotted in Fig. 8 at low temperatures as a function of the magnetization. Similar to the model simulations, see Figs. 2 and 3, the MAE follows the regular $K\propto M^{3}$ rule Callen and Callen (1966).

III.4 Fe2W(110)

Ultrathin Fe films epitaxially grown on W(110) have been studied intensively Sander (1999); Meyerheim et al. (2001) due to their peculiar magnetic properties, such as in- and out-of-plane anisotropy Hauschild et al. (1998), spin reorientation Weber et al. (1997); Ślęzak et al. (2010), and domain wall formation Heide et al. (2008). In this Section we consider the double layer system Fe2W(110). The magnetic ground state of this system strongly depends on the size and shape of the double-layer areas in the experiments von Bergmann et al. (2006); Ślęzak et al. (2010). Fe DL stripes exhibit a periodic magnetic structure with alternating out-of-plane domains separated by 180∘ walls Elmers et al. (1999). For larger DL islands there is a normal-to-plane ferromagnetic order at low temperature Weber et al. (1997), which turns into the $(1\overline{1}0)$ in-plane direction at higher temperature Dunlavy and Venus (2004).

As for Fe2Au(001), the electronic structure of Fe2W(110) was determined self-consistently via the SKKR method and the relativistic torque method was employed to find the exchange tensors and anisotropy parameters. Since a DL of Fe grows pseudomorphically on W(110)Gradmann and Waller (1982), two-dimensional translational symmetry is applied throughout the whole system with the lattice constant of bcc bulk W ($a_{\mathrm{W}}=3.16\,$Å). According experimental Santos et al. (2016) and theoretical Qian and Hübner (1999) studies there is a considerable inward relaxation of the Fe layers due to the large lattice mismatch between Fe and W. Following Ref. Heide, 2006, the Fe-W and Fe-Fe layer distances were chosen as $2.01\,$Å and $1.71\,$Å, respectively. In good agreement with previous calculations Qian and Hübner (1999); Heide et al. (2008), we obtained $2.18\,\mu_{\mathrm{B}}$ and $2.73\,\mu_{\mathrm{B}}$ for the spin-magnetic moments of Fe in the surface and the interface layer, respectively.

We employed a spin model similar to that we used for Fe2Au(001) layer, but because of the $C_{\mathrm{2v}}$ symmetry of the system bi-axial anisotropy applies,

[TABLE]

where $\hat{\mathbf{x}}$ and $\hat{\mathbf{y}}$ are unit vectors parallel to the $(1\overline{1}0)$ and $(100)$ in-plane directions, respectively. The layer-wise on-site and exchange anisotropy parameters, as explained in the case of Fe2Au(001), are summarized in Table 2. The anisotropy of the exchange couplings in the interface layer prefers the in-plane $(1\overline{1}0)$ direction which is partially compensated by the contribution from the surface layer. On the contrary, the on-site anisotropy of the interface layer clearly prefers the $(110)$ direction for the magnetization. The MAE calculated as the difference between the energy of the system magnetized in the $(1\overline{1}0)$ in-plane direction and parallel to the normal-to-plane $(110)$ direction, $E_{1\overline{1}0}-E_{110}=0.330\,\mathrm{meV}$, as well as the MAE related to the $(001)$ and $(110)$ directions, $E_{001}-E_{110}=0.368\,\mathrm{meV}$, imply indeed a normal-to-plane magnetic orientation in the ground state as also found in Refs. Qian and Hübner (1999); Heide et al. (2008).

According to susceptibility measurements Dunlavy and Venus (2004) the Curie temperature strongly depends on the Fe coverage and in case of 1.8 monolayer of Fe $T_{\mathrm{C}}=455\,K$ was measured. Our simulations on a perfect DL of Fe resulted in a Curie temperature of 520 K, in relatively good agreement with the experiment. In our metadynamics MC simulations the normal-to-plane component of the normalized magnetization was chosen again as the collective variable. In Fig. 9 the magnetic anisotropy energy defined as the difference of the free energy between the $(1\overline{1}0)$ in-plane orientation and the (110) normal-to-plane orientation is depicted for a wide range below $T_{\mathrm{C}}$. As can be inferred from this figure, the MAE changes sign at $T_{\rm r}=0.64\,T_{\mathrm{C}}$ indicating a SRT from the normal-to-plane to in-plane direction. The driving force of the spin reorientation is most probably a competition between the exchange anisotropy and the on-site anisotropy, since these contributions to the MAE exhibit different temperature dependence.

IV Summary

We introduced metadynamics as combined with Monte-Carlo simulations to study the thermal equilibrium of magnetic systems and demonstrated that the method can be applied to the temperature dependence of magnetic anisotropy of thin films. In particular, we reproduced the power-law scaling of the magnetic anisotropy vs. magnetization proposed by Callen and Callen Callen and Callen (1966) as for systems with on-site uniaxial anisotropy the simulations provided an exponent of three, whereas in case of dominating exchange anisotropy the exponent of two has been obtained in the high-temperature regime.

We applied the method to explore spin-reorientation transitions in thin films. By using a simple spin model, first we performed a detailed analysis of the SRT for a monolayer and a double-layer. For double-layers we have shown that, by setting appropriate model parameters, both first and second order SRT can occur as it is predicted within the mean field theory. Then we considered two kinds of iron double-layer systems with perpendicular magnetic anisotropy where we set up a more complex spin-model containing tensorial exchange interactions calculated from first-principles methods. In case of Fe2Au(001) the MAE followed the usual $M^{3}$ power law and no SRT was observed. In case of Fe2W(110) the MAE showed a more complex temperature dependence and our simulations reproduced the normal-to-plane to in-plane SRT seen in experiments Weber et al. (1997); Ślęzak et al. (2010). One of the future challenges for the simulations based on ab initio spin models is posed by exploring the effect of the Dzyaloshinsky-Moriya interactions on the temperature dependence of the magnetic anisotropy in thin films proposed recently Rózsa et al. (2017).

Acknowledgements.

The authors are grateful for the financial support by the Hungarian National Scientific Research Fund (NKFIH) under project No. K115575, as well as by the BME Nanotechnology FIKP grant of EMMI (BME FIKP-NAT).