Study on the critical properties of thin magnetic films using the clock model

During the last 40 years, physics of materials of nanometric size have attracted a wide interest. This is due to important applications in industry [1, 2]. An example is the so-called giant magneto-resistance (GMR) used in data storage devices, magnetic sensors... [3–6]. The experimental researches are often immediately used for industrial applications without waiting for a full theoretical understanding. In parallel to these experimental developments, much theoretical effort has also been devoted to the search of physical mechanisms lying behind new properties found in nanoscale objects such as ultrathin films, ultrafine particles, quantum dots, spintronic devices, etc. This effort aimed not only at providing explanations for experimental observations but also at predicting new effects for future experiments [7, 8]. The physics of two-dimensional (2D) systems is very exciting. Some of those 2D systems can be exactly solved: one famous example is the Ising model on the square lattice [9]. In thin films, there is a crossover from 2D to three-dimensional (3D) behavior as the film thickness increases [10]. Another interesting result is the absence of long-range ordering at finite temperatures for the continuous spin models (XY and Heisenberg models) in 2D [11]. In general, 3D systems for any spin models cannot be unfortunately solved. However, several methods in the theory of phase transitions and critical phenomena can be used to calculate the critical behaviors of these systems [12].

One of the most fascinating tasks of statistical physics is the study of the phase transition in systems of interacting particles. In particular, different kinds of transition from one phase to another have been efficiently studied by theory of renormalization group, high–low-temperature expansions [12], exact methods [13], numerical simulations [14], etc. Most phase transitions belong to either the first-order type or the second-order type. Specially, the Kosterlitz–Thouless (KT) transition is observed in superfluid systems which can be described by the 2D XY model [15].

On the other hand, q-state ferromagnetic clock model is a discrete version of the XY model, and is also known as the vector Potts model that consists of 2D planar spins with q different orientations. At the limiting cases, for $q=2$, the clock model reduces to the Ising model, while for $q\geqslant 6$, it is XY-like model. However, for $q\geqslant 4$, ones believed that a Berezinskii–Kosterlitz–Thouless (BKT) phase should occur between the ferromagnetic phase and the paramagnetic phase [15, 16]. This model is of great interest, because there is a controversy on the number of phase transitions and the character of the phase transitions [17–23]. In this paper we investigate the critical properties of thin magnetic film systems using the q-state clock model. For carrying out this purpose, we shall use Monte Carlo (MC) simulations with highly accurate Wang-Landau algorithm (flat histogram technique) [24]. We introduce a modified definition of the order parameter. Thus, the critical exponents can be determined for the transitions at both low and high temperature.

The paper is organized as follows. Section 2 is devoted to the description of the model and the simulation technique. Section 3 presented the simulation results. Concluding remarks are given in section 4.

We consider the clock model on a thin film made from the ferromagnetic square lattices. The size of the film is $\Omega =N\times N\times {{N}_{z}}$ with ${{N}_{z}}$ being the film thickness. The spin is considered as a planar vector of magnitude $\sigma =1$ with two components, $\boldsymbol{\sigma}=\left({{\sigma}^{x}},{{\sigma}^{y}}\right)$. The Hamiltonian is given by

$$\begin{eqnarray}&&H=-J\underset{\langle i,j\rangle}{\sum}\cos \left({{\theta}_{i}}-{{\theta}_{j}}\right)-{{h}_{q}}\underset{i}{\sum}\cos \left(q{{\theta}_{i}}\right),\end{eqnarray} \tag{ 1 }$$

where the sum ${{{\sum}^{}}_{\langle i,j\rangle}}$ runs over all the nearest-neighbor (NN) pairs. $J=1$ (in units of energy) is the ferromagnetic exchange coupling and ${{h}_{q}}$ being qth order symmetry-breaking field. ${{\theta}_{i}}\in 2\pi$ is the angle of spin $\boldsymbol{\sigma}_{{i}}$ makes with the x axis.

In order to investigate the nature of the phase transition of this system, we use the standard MC method and the Wang–Landau technique of simulation. Wang and Landau have proposed an MC algorithm for classical statistical models which allowed to study systems with difficultly accessed microscopic states [24]. The algorithm uses a random walk in the energy space to get an accurate estimate for the density of states (DOS) $g(E)$, which is defined as the number of spin configurations for any given $E$. This method is based on the fact that a flat energy histogram $P(E)$ is produced if the probability for the transition to a state of energy $E$ is proportional to $g{{(E)}^{-1}}$. We summarize how this algorithm is implied here. At the beginning of the simulation, the density of states $g(E)$ is unknown so all densities are set to unity, $g(E)=1$. We start a random walk in energy space $(E)$ by choosing a site randomly and flipping its spin with a transition probability

$$\begin{eqnarray}&&p\left(E\to {{E}^{\prime}}\right)=\min \left(\frac{g\left({{E}^{{}}}\right)}{g\left({{E}^{\prime}}^{{}}\right)},\,1\right),\end{eqnarray} \tag{ 2 }$$

where $E$ is the energy of the current state and ${{E}^{'}}$ is the energy of the proposed new state. Each time an energy level $E$ is visited, the DOS is modified by a modification factor $f>0$ whether the spin is flipped or not, i.e. $g(E)\to g(E)f$. Initial value of the modification factor $f$ can be as large as ${{\text{e}}^{1}}\approx 2.718\,2818$. A histogram $P(E)$ records the number of times a state of energy $E$ is visited. Each time the energy histogram satisfies a certain 'flatness' criterion, the histogram $P(E)$ is then reset to zero, and the modification factor is reduced, typically to the square root of the previous factor, to produce a finer estimate of $g(E)$. The reduction process of the modification factor $f$ is repeated several times until a final value ${{f}_{\text{final}}}$ which is close enough to one. The histogram is considered as flat if

$$\begin{eqnarray}&&P(E)\geqslant x\%\langle P(E)\rangle \end{eqnarray} \tag{ 3 }$$

for all energies, where x% is chosen between 90% and 95% and $\langle P(E)\rangle$ is the average histogram.

The thermal average of a thermodynamic quantity $A$ can be evaluated by

$$\begin{eqnarray}&&{{\langle A\rangle}_{\text{T}}}=\frac{1}{Z}\underset{E}{\sum}g(E)A\exp \left(-\frac{E}{{{k}_{\text{B}}}T}\right),\end{eqnarray} \tag{ 4 }$$

in which

$$\begin{eqnarray}&&Z=\underset{E}{\sum}g(E)\exp \left(-\frac{E}{{{k}_{\text{B}}}T}\right).\end{eqnarray} \tag{ 5 }$$

Thermal averages of physical quantities are thus calculated as continuous functions of T, now the results should be valid over a much wider range of temperature.

In MC simulations, one calculates the magnetic susceptibility defined by

$$\begin{eqnarray}&&\chi =\frac{1}{{{k}_{\text{B}}}T}(\langle {{m}^{2}}\rangle -{{\langle m\rangle}^{2}})\,,\end{eqnarray} \tag{ 6 }$$

$m$ being the magnetization (also called the order parameter). It has following form:

$$\begin{eqnarray}&&m=\frac{1}{\Omega}\sqrt{{{\left(\underset{i}{\sum}\cos \left({{\theta}_{i}}\right)\right)}^{2}}+{{\left(\underset{i}{\sum}\sin \left({{\theta}_{i}}\right)\right)}^{2}}}.\end{eqnarray} \tag{ 7 }$$

For the q-state Potts model, the value of spins is taken from a finite set of integers $\left\{1,\,2\,,\ldots,\,\,q\right\}$. Therefore, the order parameter is defined by other equation

$$\begin{eqnarray}&&{{m}_{\text{P}}}=\frac{q{{M}^{\max}}-1}{q-1},\end{eqnarray} \tag{ 8 }$$

with

$$\begin{eqnarray}&&{{M}^{\max}}=\frac{\max \left({{M}_{1}},{{M}_{2}},\,\ldots,\,{{M}_{\text{q}}}\right)}{\Omega},\end{eqnarray} \tag{ 9 }$$

where ${{M}_{i}}$ is the number of spins which have a same value $i\left(i=1,\,2,\,\ldots,\,q\right)$.

In order to study the critical behaviors, we use the scaling relations given by

$$\begin{eqnarray}&&V_{1}^{\max}\propto {{N}^{1/\nu}},\,\,V_{2}^{\max}\propto {{N}^{1/\nu}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{{\chi}^{\max}}\propto {{N}^{\gamma /\nu}},\end{eqnarray} \tag{ 11 }$$

$\nu$ and $\gamma$ being the critical exponents. ${{V}_{n}}$ is nth order cumulant of the order parameter which is written in the form

$$\begin{eqnarray}&&{{V}_{n}}=\frac{\partial \,\ln \left({{O}^{n}}\right)}{\partial \left(1/{{k}_{\text{B}}}T\right)}=\langle E\rangle -\frac{\langle {{O}^{n}}E\rangle}{\langle {{O}^{n}}\rangle}.\end{eqnarray} \tag{ 12 }$$

Here, the order parameter $O$ is replaced by either $m$ in equation (7) or ${{m}_{\text{P}}}$ in equation (8). In this paper, we only estimate $\nu$ from $V_{2}^{\max}$, with this value we calculate $\gamma$ from ${{\chi}^{\max}}$.

Let us discuss on a specific case of $q=4$, there exists two transitions which separate the system into three phases: ferromagnetic phase, BKT phase and the paramagnetic phase [22]. These phase transitions are characterized by non-universal critical exponents. The system undergoes from ferromagnetic phase to BKT phase passing through the first transition at low temperature (ferromagnetic-BKT transition), and then it passes through the second transition at high temperature (BKT-paramagnetic transition). From the MC simulation results, although the plot of the specific heat clearly showed a peak for the ferromagnetic-BKT transition at low temperature, ones could not estimate the critical exponents $\nu$ and $\gamma$ because there was no peak in the plots of susceptibility and nth order cumulant. In other words, the definitions of the order parameter (7) and (8) were not usable for calculating the critical exponents of the ferromagnetic-BKT transition at low temperature (the first transition).

Now, we suggest a new method for the calculation of the order parameter by modifying the definition of ${{M}_{i}}$ in equation (9), i.e. ${{M}_{i}}$ is the total number of the spin vectors within ith circular sector ${{S}_{i}}$ with central angle ${{\varphi}_{i}}$ $\,\left(i=1,\,2,\,\ldots,\,q\right)$. The central angle satisfies the condition $\sum\nolimits_{i=1}^{q}{{\phi}_{i}}=2\pi$ and it is defined by

$$\begin{eqnarray}&&{{\phi}_{i}}=\left({{\Phi}_{i}}-\frac{\pi}{q},{{\Phi}_{i}}+\frac{\pi}{q}\right]\,,\end{eqnarray} \tag{ 13 }$$

where ${{\Phi}_{i}}$ are the angles corresponding to the maximum of anisotropic term in the Hamiltonian (1), i.e. $\cos \left(q\theta \right)=\cos \left(\Phi \right)=1$ or ${{\Phi}_{i}}=2\pi \,n/q$ with $\,n=0,\,1,\ldots,\,q-1$. Then formula (13) becomes

$$\begin{eqnarray}&&{{\phi}_{i}}=\left(\frac{\pi (2i-3)}{q},\frac{\pi (2i-1)}{q}\right],\,i=1,\,2,\,\ldots,\,q.\end{eqnarray} \tag{ 14 }$$

For example, since the system has four-states ($q=4$), we obtain

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}{{\phi}_{1}}=\left(-\pi /4,\pi /4\right],\quad {{\phi}_{2}}=\left(\pi /4,3\pi /4\right],\\ \quad {{\phi}_{3}}=\left(3\pi /4,5\pi /4\right],\quad {{\phi}_{4}}=\left(5\pi /4,-\pi /4\right].\end{array}\end{eqnarray}$$

In this section, we consider a special case of ${{N}_{z}}=1$ for comparing with previous results. We shall present the results obtained by MC simulations using the Wang-Landau algorithm. Note that the order parameter is obtained from equations (8), (9) and (14). The xy linear sizes N  =  20, 30, ..., 80 have been used in our simulations. Errors shown in the following have been estimated using statistical errors, which are very small and fitting errors given by fitting software.

Let us consider a specific case of $q=2$ with a strong field ${{h}_{q}}$ for testing the new method of the calculations. It is well know that the system becomes 2D Ising-like model in limit ${{h}_{q}}\to \infty$ with a single phase transition. We calculate the universal critical exponents $\nu$ and $\gamma$ by following the expressions (10)–(12). We show in figures 1 and 2 the maximum of 2nd order cumulant and the maximum of susceptibility versus lattice size in the ln–ln scale, respectively. If we use the equations (10) and (11) to fit these lines, i.e. without correction to scaling, we obtain $1/\nu$ and $\gamma /\nu$ from the slopes of the remarkably straight lines. For ${{h}_{q}}=0.1$, we have $1/\nu =0.991\,63\pm 0.006$ which yields $\nu =1.0084$ and $\gamma /\nu =1.734\,76\pm 0.008$ yielding $\gamma =1.7494$. These results are in excellent agreement with the exact results ${{\nu}_{\text{2D}}}=1$ and ${{\gamma}_{\text{2D}}}=1.75$. Thus, we can conclude that our calculation method is applicable for the clock model. We summarize in table 1 the estimative values of $\nu$ and $\gamma$ for several values of symmetry-breaking field ${{h}_{q}}=0.01,\,0.05,0.1$ and 0.5. The critical exponents are slightly decrease with increasing ${{h}_{q}}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. The critical exponent for $q=2$ with several ${{h}_{q}}=0.01,\,0.05,0.1$ and 0.5.

${{h}_{q}}$	$1/\nu$	$\gamma /\nu$	$\nu$	$\gamma$
0.01	0.864 99  ±  0.006	1.640 55  ±  0.013	1.156 08	1.896 61
0.05	0.971 01  ±  0.005	1.737 89  ±  0.002	1.029 86	1.789 78
0.10	0.991 63  ±  0.006	1.734 76  ±  0.008	1.008 44	1.749 40
0.50	1.005 66  ±  0.003	1.749 01  ±  0.006	0.994 37	1.739 17

For $q=3$, there is still a single transition at high temperature. If ${{h}_{q}}>0.01$, the critical exponents are most likely to that of Ising-like case. As see in table 2, the critical exponents have an unusual value for ${{h}_{q}}=0.01$. This transition is mixed of the two transitions: ferromagnetic-BKT transition and BKT-paramagnetic transition. It indicates that the BKT phase starts to occur from $q=3$ with ${{h}_{q}}$ is small enough.

Table 2. The critical exponent for $q=3$ with several ${{h}_{q}}=0.01,\,\,0.03,\,0.05$ and 0.1.

${{h}_{q}}$	$1/\nu$	$\gamma /\nu$	$\nu$	$\gamma$
0.01	0.604 46  ±  0.003	1.149 02  ±  0.018	1.654 37	1.900 90
0.03	0.656 67  ±  0.002	1.793 59  ±  0.020	1.522 83	2.731 34
0.05	0.699 56  ±  0.004	1.774 09  ±  0.020	1.429 47	2.536 01
0.10	0.852 91  ±  0.018	1.860 17  ±  0.013	1.172 46	2.180 97

For $q=4$ and ${{h}_{q}}=0.05$, we show in figure 3 the peak of the susceptibilities located around the critical temperature $T\approx 1.05$. This is the second transition at high temperature (BKT-paramagnetic transition). The maximum of susceptibility increases with increasing the system size that is a signature of the phase transition. Note that, the susceptibilities ${{\chi}_{\text{m}}}$ are obtained by using equations (6) and (7). Therefore, our result is in good agreement with the previous result (see figure 1 in [22]). Besides, the curves of ${{\chi}_{\text{m}}}$ have a prominence located around $T\approx 0.75$ which is expected to be the first transition at low temperature (ferromagnetic-BKT transition). Of course, ones could not obtain the critical exponents of this transition because there is no peak in the plots of the susceptibility and the second order cumulant. In contrast, if the set of equations (6), (8), (9) and (14) are applied for calculating the magnetization, then we can obtain the expected susceptibilities which are presented in figure 4. It clearly shows two peaks corresponding to the two transitions: one is BKT-paramagnetic transition at high temperature $T\to 1.0$ and the other one is the ferromagnetic-BKT transition at low temperature $T\to 0.78$. The critical exponents $\nu$ and $\gamma$ are given in tables 3 and 4 for the first and the second phase transition, respectively. These results are again in good agreement with the experimental data [23].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. The critical exponent of the ferromagnetic-BKT transition with $q=4$ and several ${{h}_{q}}=0.01,\,\,0.05,\,0.1$ and 0.5.

${{h}_{q}}$	$1/\nu$	$\gamma /\nu$	$\nu$	$\gamma$
0.01	0.481 76  ±  0.008	1.115 49  ±  0.017	2.075 72	2.315 44
0.05	0.524 28  ±  0.006	1.316 56  ±  0.021	1.907 37	2.511 17
0.10	0.581 61  ±  0.007	1.389 87  ±  0.015	1.719 36	2.389 69
0.50	0.647 42  ±  0.005	1.531 46  ±  0.006	1.518 65	2.365 48

Table 4. The critical exponent of the BKT-paramagnetic transition with $q=4$ and several ${{h}_{q}}=0.01,\,\,0.05,\,0.1$ and 0.5.

${{h}_{q}}$	$1/\nu$	$\gamma /\nu$	$\nu$	$\gamma$
0.01	0.599 74  ±  0.007	1.607 60  ±  0.011	1.667 39	2.680 49
0.05	0.588 63  ±  0.009	1.607 78  ±  0.011	1.698 86	2.731 39
0.10	0.593 90  ±  0.009	1.604 07  ±  0.010	1.683 79	2.700 91
0.50	0.658 48  ±  0.006	1.529 24  ±  0.008	1.518 65	2.322 38

Before closing this section, let us discuss on the behavior of the phase transition for $q>4$. The critical exponents do not depend neither on the number of state nor symmetry-breaking field. From the simulation results, we obtain the critical exponents $\nu \approx 2$ and $\gamma \approx 3$ for the first transition, $\nu \approx 1.65$ and $\gamma \approx 2.6$ for the second transition.

We have studied a system, namely the q-state clock model on the thin film. The order parameter was calculated by using the new definitions. Therefore, we could easily determine the maxima of the susceptibility and nth order cumulant of order parameter, though the system size was very small. For the simulations, we considered the 2D case in order to clarify the point whether or not there was an existence of the BKT phase with varying q state and h_q field.

From the results obtained by very highly accurate flat histogram technique shown above, we conclude that the existence of the BKT phase strongly depends on both q and h_q. In the case q  =  3, the BKT phase starts to occur when ${{h}_{q}}\leqslant 0.01$. While it starts to occur when ${{h}_{q}}\leqslant 0.5$ for q  =  4. Note that, this phase is always observed for the five-state and six-state [17].

For the BKT-paramagnetic transition at high temperature, the critical exponents slightly change with varying q state and h_q due to the statistical errors. For the first transition at low temperature (ferromagnetic-BKT transition), the non-universal critical exponents strongly depend on h_q when q  <  5 while these have a few changes if otherwise. Our results confirmed the validation of the previous expectation [22] and the experimental data [23].

This work was supported by the NAFOSTED (National Foundation for Science and Technology Development, Vietnam), Grant No. 103.02–2011.55.