Spin kinetic Monte Carlo method for nanoferromagnetism and magnetization dynamics of nanomagnets with large magnetic anisotropy

Abstract

The Kinetic Monte Carlo (KMC) method based on the transition-state theory, powerful and famous for simulating atomic epitaxial growth of thin films and nanostructures, was used recently to simulate the nanoferromagnetism and magnetization dynamics of nanomagnets with giant magnetic anisotropy. We present a brief introduction to the KMC method and show how to reformulate it for nanoscale spin systems. Large enough magnetic anisotropy, observed experimentally and shown theoretically in terms of first-principle calculation, is not only essential to stabilize spin orientation but also necessary in making the transition-state barriers during spin reversals for spin KMC simulation. We show two applications of the spin KMC method to monatomic spin chains and spin-polarized-current controlled composite nanomagnets with giant magnetic anisotropy. This spin KMC method can be applied to other anisotropic nanomagnets and composite nanomagnets as long as their magnetic anisotropy energies are large enough.

Graphical abstract

Keywords

nanomagnet / spin / Monte Carlo method / information science

1 Introduction

With the development of nanoscience and information technology, it is highly desirable to explore the nanomagnetism and magnetization stability of nanomagnets. There is a nonzero equilibrium spin magnetization in a bulk ferromagnet, but one should observe a zero average magnetization in a nanoscale ferromagnet in equilibrium because of Mermin-Wagner theorem [1, 2]. Here by nanomagnet we mean that its size is tens of nanometers or smaller in all the three space directions. Rigorously speaking, the average magnetization of a nanoferromagnet is always changing with the time. The rate of the change varies enormously, depending on temperature, sample size, materials parameters, and so on. Generally speaking, the magnetization can be reversed through thermal fluctuations (and quantum tunnelling effects), or by applying magnetic fields or spin-polarized currents.

Gambardello et al. [3] synthesized one-dimensional (1-D) monatomic Co chains by self-assembly epitaxial techniques and showed that the magnetic anisotropy energy can reach to 9 meV per atom [4]. The accurate first- principle calculation indicated that much larger magnetic anisotropy may be obtained in near future. In these systems the magnetic anisotropy energy can be compared with or even exceed the interacting energy between moments. Giant magnetic anisotropy energy can stabilize the spin orientation and therefore is useful to applying such nanomagnets to information science and technology. It is highly desirable to investigate the effect of the giant magnetic anisotropy [5−9].

We are interested in the dynamical properties of these nanomagnets. Recently, we used the Kinetic Monte Carlo method (KMC) [10−12], powerful and famous for simulating atomic growth of epitaxial films and nanostructures, to simulate the dynamical spin properties of nanomagnets with giant magnetic anisotropy [13, 14]. Experimental nanoferromagnetism at low temperatures is explained naturally and convincingly, which in turn confirms that our spin KMC method is reliable. Here we review the spin KMC method and its applications to such nanomagnets.

The remaining part of this article is organized as follows. In the next section we briefly introduce the KMC method based on the transition-state theory and show how to reformulate it to simulate spin dynamics of nanomagnets with large enough magnetic anisotropy. In the third section we show how to apply the spin KMC method to simulate the intrinsic nanoferromagnetism of monatomic spin chains. In the fourth section we show how to simulate the magnetization reversal induced by spin-polarized current by the spin KMC method. In the fifth section we present a concluding remark on the further development of the spin KMC method and its possible applications to other nanoscale spin systems.

2 Kinetic Monte Carlo method for spin systems

The KMC method, different from other Monte Carlo methods [15, 16], is based on the transition-state theory [17, 18]. We use the term “process” to denote the change from a starting state to a destination state. All the transition states between the starting state and the destination state must be considered. Usually there exists a transition-state energy barrier, ∆E, for a process. It must be overcome to go from the starting state to the destination state within a KMC step. At finite temperature T, a process may happen with a rate defined in terms of an Arrhenius law [19]:

\begin{array}{}\text{(1)}& R=R_0{\text{e}}^{-\text{Δ}E/k_{\text{B}}T}\end{array}

where k_B is the Boltzmann constant and R₀ the attempt frequency.

We use Fig. 1 to give more details. The ball moves along the line. There are three equilibrium sites: A, B, and B′. The ball at site A must climb over peak C (or C′) to arrive at site B (or B′). Meanwhile it must overcome an energy barrier E_{\text{A}}^{\text{C}} (or E_{\text{A}}^{{\text{C}}^{\prime }}). Within a KMC step, the ball can (1) move to site B, (2) move to site B′, or (3) remain at site A. For the process from site A to B, we have the rate

\begin{array}{}\text{(2)}& R_{\text{A}\to \text{B}}=R_0^{\text{A}}{\text{e}}^{-E_{\text{A}}^{\text{C}}/(k_{\text{B}}T)}\end{array}

where R_0^{\text{A}} is determined by the local potential shape around site A. For a KMC simulation, we define the time sequence as t₀, t₁, t₂, t₃, ···, where the relationship t_n = t_n−1 + ∆t is satisfied for n > 0. Supposing the ball is at site A at time t_n, we have three possibilities for the ball:

(1) the ball moves to site B with probability

\begin{array}{}\text{(2)}& P_{n1}=\text{Δ}t\cdot R_{\text{A}\to \text{B}}\end{array}

(2) the ball moves to site B′ with probability

\begin{array}{}\text{(2)}& \begin{array}{l}{P}_{n2}=\text{Δ}t\cdot R_{\text{A}\to {\text{B}}^{\prime }}=R_0^A{\text{e}}^{-E_{\text{A}}^{{\text{C}}^{\prime }}/(k_{\text{B}}T)}\hfill \end{array}\end{array}

(3) the ball remains at site A with probability

\begin{array}{}\text{(2)}& P_{n0}=1-P_{n1}-P_{n2}\end{array}

Here the condition {\displaystyle \sum _{i=1}{P}_{ni}⩽1} must be satisfied by adjusting down the time interval ∆t. At the next step, the ball must be at one of the three sites A, B, and B′. We produce a random number, ρ, in the region [0, 1] and determine which site it will be at in the following way:

(1) if 0 ≤ ρ < P_n1, then the ball will be at site B;

(2) if P_n1 ≤ ρ < P_n1 + P_n2, then the ball will be at site B′;

(3) if P_n1 + P_n2 ≤ ρ ≤ 1, then the ball will be still at site A.

If the ball is already at site B or B′, we have only two possibilities instead of three for the next KMC step. If we have N balls, the above step must be repeated for N times within every KMC step so that every ball is selected once on average.

Now we turn to spin systems. Figure 2 shows a typical spin chain with inter-spin ferromagnetic exchange interaction. A spin alone can takes any direction, but all the spins tend to orient in one single direction because of the inter-spin interaction. A strong uniaxial magnetic anisotropy with the z easy axis is supposed. Hence, there is an energy favored axis for all the spins. In order to reformulate the KMC method for spin systems, we need two key ingredients: (1) equilibrium directions, and (2) transition-state barriers and transition rates. When the anisotropy is strong enough, a spin, labelled as S_i, naturally has two degenerate equilibrium directions: +z and −z. Then the two equilibrium directions for each spin correspond to equilibrium sites A, B, and B′ for the ball discussed above. Although it is quantum in essence, the spin can be described by a classical vector as long as the spin value is high enough, the temperature is high enough, and/or the anisotropy is strong enough. We use θ_i to describe the angle of spin S_i from its original equilibrium direction. θ_i takes values in the region [0,π]. “0” is the original state and “π” is the reversed state. Other states with θ_i satisfying 0 < θ_i < π are transition states. If the anisotropy is defined as -k_u{(S_i^z)}^2{{\displaystyle \cos }}^2{\theta }_i (k_u > 0), there is a transition-state barrier E_r = k_uS² for the reversal of an isolate spin. Because of the inter-spin interaction, the transition-state barrier E_r may be changed. It can also be changed by applying a magnetic field. It even disappears for some special processes. With the transition-state barrier E_r, we can calculate the reversal rate

\begin{array}{}\text{(3)}& R_r=R_{r0}{\text{e}}^{-E_r/(k_{\text{B}}T)}\end{array}

where R_r0 is the attempt frequency for spin freedom, usually taking values from 10⁹ to 10¹²/s.

With the equilibrium directions and transition rates, we can do our spin KMC simulation for such a spin system consisting of N spins with a uniaxial anisotropy k_u. We can do it in the following way:

(1) We define a KMC time sequence as t₀, t₁, t₂, t₃, ··· in the same way as the above and then construct the original spin configuration, {S_0i}.

(2) For step n (n=1, 2, 3, ···) we can calculate a reversal rate R_ni and a reversal probability P_ni =∆t · R_ni for each spin S_(n−1)i. Every spin is taken to make one try on average for reversal. It can reverse with probability P_ni or remains unchanged with 1−P_ni. This should be repeated N times if we have N spins in total. The updated spin configuration {S_ni} will be used for the next step.

(3) Adding n by 1 and replacing n + 1 with n formally, repeat step (2) unless the condition for stoping is satisfied.

(4) Calculate average magnetization m by averaging all the spins

\begin{array}{}\text{(4)}& m=\frac{1}{N}{\displaystyle \sum _{i=1}^N{S}_{n_{f}i}^z}\end{array}

where n_f denotes the sequence number for the last step. In order to reduce possible errors, we should repeat the above run independently many times to make the average magnetization stable enough. It is helpful to keep P_ni substantially smaller than 1.

3 Intrinsic nanoferromagnetism of nanomagnets with giant magnetic anisotropy

Now we discuss the application of the spin KMC method to the monatomic Co spin chains [8]. In our simulation we keep the temperature higher than 5 K in order to use the classical description of spins [13, 14]. The Co magnetic moment on average is approximately 4 µ_B, which is quite large, and has its easy axis in the z direction. Hence the spin Hamiltonian can be written as [13]

\begin{array}{}\text{(5)}& H=-{\displaystyle \sum _{i\mathrm{,}j}{J}_{ij}{\widehat{S}}_i\cdot {\widehat{S}}_j-k_u}{\displaystyle \sum _i{(S_i^z)}^2}-\mu H{\displaystyle \sum _i{S}_i^z}\end{array}

where the first term describes the Heisenberg inter-spin ferromagnetic exchange interaction, the second one the giant magnetic anisotropy energy, and the third one the Zeeman energy in the presence of a magnetic field along the z direction. Using the spin angle θ_i defined above, we obtain the Hamiltonian of the ith spin [13]:

\begin{array}{}\text{(6)}& \begin{array}{ll}{H}_i=\hfill & -{\displaystyle \sum _j}{J}_{ij}^s{s}_i{s}_j\cos {\theta }_i-k_u^s{{\displaystyle \cos }}^2{\theta }_i-{\mu }^{s}H{s}_i\cos {\theta }_i\hfill \end{array}\end{array}

where the s_i is the reduced spin variable, taking one of +1 or −1, and the superscript “s” for the constants means the spin value S (4 in the present case) has been absorbed into them. The variable s_i reminds us of the Ising model [20], but our original spin angle can change continuously. The two values of s_i is caused by the strong anisotropy and the KMC method together. We can derive the energy of transition states:

\begin{array}{}\text{(7)}& E_i=k_u^s{{\displaystyle \sin }}^2{\theta }_i-h_i(\cos {\theta }_i-1)\end{array}

where h_i is defined by

\begin{array}{}\text{(8)}& h_i=\left({\displaystyle \sum _j{J}_{ij}^s{s}_j+{\mu }^{s}H}\right)s_i\end{array}

If |h_i/2k_u^s| is less than 1, the energy E_i as a function of cos θ_i has a maximum and therefore the rate for the spin reversal is

\begin{array}{}\text{(9)}& R_i=R_{r0}{\text{e}}^{-\text{Δ}{E}_i/(k_{\text{B}}T)}\end{array}

where the barrier is equivalent to

\begin{array}{}\text{(10)}& \text{Δ}{E}_i=(2k_u^s+h_i{)}^2/(4k_u^s)\end{array}

When the condition |h_i/2k_u^s|⩾1 is satisfied, there is no barrier and we use the Glauber method [21] to determine the rate, with the prefactor kept the same.

Our KMC simulation was done with the above equilibrium directions and transition rates. We take N = 80, set R₀ to 7 × 10⁹, and use field sweeping rate 132 T/s with sweeping interval 0.1 T. The inter-spin interaction is limited to the nearest neighboring sites. We use the time interval ∆t = 7.576 × 10⁻⁴s and therefore we have ∆t · R_r0 = 5 303 030. We repeat the KMC simulation for 200 times to reduce random errors. Our simulated results for both 10 and 45 K are shown in Fig. 3. A hysteresis loop of average magnetization against the magnetic field is observed for 10 K. It is clear that the spin chain of 80 atoms has ferromagnetism at low temperatures such as 10 K, but the ferromagnetism crossovers to a paramagnetic phase at high enough temperatures such as 45 K. A coercive field H_c can be defined. Its nonzero value reflects the ferromagnetism. Our systematic KMC simulations indicate that H_c decreases with the temperature. Upon reducing the field sweeping rate down to zero, we obtain H_c = 0. This implies that the ferromagnetism is a special phenomenon for nanomagnets with large enough magnetic anisotropy.

4 Current controlled magnetization reversal of nanomagnets with giant magnetic anisotropy

Figure 4 shows our composite nanomagnet consisting of an array of fixed spins with giant magnetic anisotropy and a thin, narrow metallic stripe into which a spin-polarized current can be injected. The spin polarization p of the injected current is defined as p = (n_↑−n_↓)/(n_↑+ n_↓), where n_↓ and n_↑ are the number of down spins and that of up spins of the injected electrons over an appropriate time interval. Generally speaking, such a system is very complicated and we should treat the charge and spin freedoms at the same time. However, we are interested in small nanomagnets and suppose the current density is not so large that there is at most one electron in the region of the fixed spins at a time. A moving electron cannot meet another one in the region within a KMC step as long as the time interval ∆t is small enough. Spin wave excitations can be neglected due to the large enough magnetic anisotropy. Therefore, we can use a pure spin model to describe the magnetization dynamics in such composite nanomagnets.

The Hamiltonian of the fixed spins and the injected spins can be expressed as

\begin{array}{}\text{(11)}& \begin{array}{ll}H=\hfill & -{\displaystyle \sum _{i\mathrm{,}j}}{J}_{ij}{\widehat{S}}_i^{\text{B}}\cdot {\widehat{S}}_j^{\text{B}}-{\displaystyle \sum _i}{k}_u^{\text{B}}{(S_i^{\text{B}z})}^2\hfill \\ \hfill & -{\displaystyle \sum _{i\mathrm{,}j}}{J}_{ij}^{\text{'}}{\widehat{S}}_i^{\text{B}}\cdot {\widehat{S}}_j^{\text{A}}-{\displaystyle \sum _i}{k}_u^{\text{A}}{(S_i^{\text{A}z})}^2\hfill \end{array}\end{array}

where {\widehat{S}}_i^{\text{B}} describes the fixed spin at site i with anisotropy k_u^{\text{B}}, and {\widehat{S}}_i^{\text{A}} describes the electron spin of injected current with anisotropy k_u^{\text{A}}. J_ij (J_ii = 0) is the inter-spin Heisenberg exchange interaction and J_{ij}^{\text{'}} the Heisenberg exchange interaction between the fixed spin and the moving spin. For the present case we need two spin angles, θ_i and {\theta }_i^{\text{'}}, to describe the orientations of the two spins with respect to the corresponding equilibrium directions. The Hamiltonians of the fixed spin and the moving spin at site i can be written as [14]

\begin{array}{}\text{(12)}& \begin{array}{l}{H}_i=k_u^{Bs}{{\displaystyle \sin }}^2{\theta }_i-{\displaystyle \sum _j}(J_{ij}^s{s}_j^B+J_{ii}^{\text{'}s}{s}_i^A)s_i^B\cos {\theta }_i\hfill \end{array}\end{array}

and

\begin{array}{}\text{(13)}& H_i^{\text{'}}=k_u^{As}{{\displaystyle \sin }}^2{\theta }_i^{\text{'}}-J_{ii}^{\text{'}s}{s}_i^B{s}_i^A\cos {\theta }_i^{\text{'}}\end{array}

where the parameters with superscript “s” correspond to the above-defined ones with some S and S ² factors absorbed. The energies are similar to those in the last section, and can be expressed in terms of a simple function, f(x) = u(1−x²)−vx, with x being cos θ_i and cos {\theta }_i^{\text{'}}. If the condition is satisfied, f(x) gets to its maximum f_max = u + v²/4u and a transition-state barrier can be calculated by ∆E = f_max−f(0) = (v + 2u)²/4u. If not, we use the Glauber method [21] as we do in the last section. As a result, we obtain barriers \text{Δ}{E}_i^s and \text{Δ}{E}_i^{\text{'}s} of the two spins at site i for the transition states described by θ_i and {\theta }_i^{\text{'}}. Corresponding rates R_i and R_i^{\text{'}} can be calculated by Eq. (9), and the probabilities are equivalent to ∆t · R_i and \text{Δ}t\cdot R_i^{\text{'}}, respectively. We have N fixed spins and N′ moving spins. Because N′ can be 0 or 1 for our simulation, the total number of spins, N ^t = N + N′, is changing.

The spin KMC simulation was done for N = 20, with the movement of the injected spin supposed to be not affected by the fixed spins [14]. The simulated results of average magnetizations against applied currents are shown in Fig. 5. The magnetic anisotropy energies for the fixed spin and the moving spin are k_u^{\text{B}s}=16.8 and k_u^{\text{A}s}=12 meV. The exchange constant between fixed spins is J^s=11.2 meV and that between the fixed spin and the moving spin is J′^s=8 meV. The simulation temperature is 30 K. We observe hysteresis loops for the magnetization against applied currents for p = 1 and p = 0.75. There are substantially fluctuations for the magnetization in the case of p = 0.75. A coercive current I_c can be defined in the same way as we do for the field-driven hysteresis loops in the last section. Further simulations show that I_c decreases with the temperature and increases with the anisotropies. The time evolution of the magnetizations can be studied systematically. The magnetization m can be well described by a simple time function:

\begin{array}{}\text{(14)}& m=m_{\text{s}}-(1+m_{\text{s}})\exp \left[-\frac{{(t/t_0)}^2}{b+t/t_0}\right]\end{array}

where m_s is the saturated magnetization and the parameters t₀ and b are constants. Such simulation can be done with other starting magnetizations. Whatever the starting magnetization, the same saturated magnetization is obtained. It depends on the temperature, applied spin-polarized current, and materials parameters.

5 Conclusion

We present a brief introduction to the usual KMC method and show how to reformulate it for simulating the intrinsic and current-controlled magnetization dynamics of nanoscale spin systems with large enough magnetic anisotropy. The large enough magnetic anisotropy is not only essential to stabilize spin orientation but also necessary in making the transition-state barriers of spin reversals for spin KMC simulation. We show two applications of the spin KMC method to nanomagnets with giant magnetic anisotropy. This spin KMC method can be applied to other anisotropic nanomagnets and composite nanomagnets as long as their magnetic anisotropy energies are large enough [22].

References

[1] Landau L. D. and Lifshitz E. M., Statistical Physics, London: Pergmon, 1959, 5: 482

[2] Mermin N. D. and Wagner H., Phys. Rev. Lett., 1966, 17: 1133

[3] Gambardella P., Dallmeyer A., Maiti K., Malagoli M. C., Eberhardt W., Kern K., and Carbone C., Nature, 2002, 416: 301

[4] Gambardella P., Rusponi S., Veronese M., Dhesi S. S., Grazioli C., Dallmeyer A., Cabria I., Zeller R., Dederichs P. H., Kern K., Carbone C., and Brune H., Science, 2003, 300: 1130

[5] Wernsdorfer W., Orozco E. B., Hasselbach K., Benoit A., Barbara B., Demoncy N., Loiseau A., Pascard H., and Mailly D., Phys. Rev. Lett., 1997, 78: 1791

[6] Back C. H., Allenspach R., Weber W., Parkin S. S. P., Weller D., Garwin E. L., and Siegmann H. C., Science, 1999, 285: 864

[7] Koch R. H., Grinstein G., Keefe G. A., Lu Yu, Trouilloud P. L., Gallagher W. J., and Parkin S. S. P., Phys. Rev. Lett., 2000, 84: 5419

[8] Wernsdorfer W., Clerac R., Coulon C., Lecren L., and Miyasaka H., Phys. Rev. Lett., 2005, 95: 237203

[9] Curilef S., del Pino L. A., and Orellana P., Phys. Rev. B, 2005, 72: 224410

[10] Witten T. A. and Sander L. M., Phys. Rev. Lett., 1981, 47: 1400

[11] Zhang Z. and Lagally M. G., Science, 1997, 276: 377

[12] Liu B. G., Wu J., Wang E. G., and Zhang Z., Phys. Rev. Lett., 1999, 83: 1195

[13] Li Y. and Liu B. G., Phys. Rev. B, 2006, 73: 174418

[14] Li Y. and Liu B. G., Phys. Rev. Lett., 2006, 96: 217201

[15] Metropolis N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H., and Teller E., J. Chem. Phys., 1953, 21: 1087

[16] Binder K. and Heermann D. W., Monte Carlo Simulation in Statistical Physics: An Introduction, 3rd ed., Berlin: Springer, 1997

[17] Gold V., Pure. Appl. Chem., 1979, 51: 1725

[18] Laidler K. J. and King M. C., J. Phys. Chem., 1983, 87: 2657

[19] Kirby R. D., Shen J. X., Hardy R. J., and Sellmyer D. J., Phys. Rev. B, 1994, 49: 10810

[20] Ising E., Z. Phys., 1925, 31: 253

[21] Glauber R. J., J. Math. Phys., 1963, 4: 294

[22] Zhang K. C. and Liu B. G., to be published

Rights & permissions

Higher Education Press and Springer-Verlag 2007