Key role of temperature in ferromagnetic Bloch point simulations

PHYSICAL REV[EW B 86, 094409 (2012)

Key role of temperature in ferromagnetic Bloch point simulations

K. M. Lebecki: D. Hinzke, and U. Nowak

Department of Physics, University of Konstanz, D-78457 KonS/(lnz, Germany

O. Chubykalo-Fesenko

Instituto de Ciencia de Materiales de Madrid, CSIC, Cantoblanco, E-28049 Madrid, Spain (Received 12 March 2012; revised manuscript received 20 August 2012; published 7 Septe\nber 2012)

Bloch points in permalloy cylinders are investigated using a micromagnetic framework, where thermal effects are included via the Landau-Lifshitz-Bloch equation of motion. We show that this enables micromagnetic modeling of a Bloch point avoiding the problem of singularities, which have been rep0l1ed in the literature so far. The details of the Bloch point which we reveal are compared with earlier analytic approximations describing its geometry and the magnetization drop in its center. The temperature dependence of characteristic parameters, like the Bloch point radius or the azimuthal inflow angle is given in the full temperature range.

00[: 10.11 03/PhysRevB.86.094409 PACS number(s): 75.40.Mg, 75.60.Jk, 75.60.Nt, 75.75.Fk

I. INTRODUCTION

Bloch points (BP) are small objects that have attracted attention recently as they play an essential role during the process of vOltex core (VC) switching.^I-3 Magnetic vortices themselves are of great interest for their possible application in new magnetic storage concepts4 or in sub- gigahertz wave generators.5 A magnetic vortex is basically composed of a rotational exterior and a narrow elongated interior-called a core-where the magnetization is point- ing out of the plane containing the vortex. A BP is the point, where the orientation (polarity) of the VC changes (see Fig. I), occurring, for example, during switching of the Vc.^I•3 In "classical" micromagnetism the magnetization length is kept constant.6 With that assumption classical micromagnetism is actually a zero-temperature approach, neglecting thermal fluctuations, and a BP is a singularity.7 In its close vicinity "any direction of the magnetization is present."g

There have already been several trials to investigate the BP by micromagnetic simulations'.3,8-11 All of these investigations were, however, performed in the framework of classical micromagnetism, where the BP is a singularity and serious problems appear when treating it in a numerical wayY·11 Criterion of reliable simulations is as following: The maximum angle between magnetization vectors in neighboring cells

~ipmax must remain small (typically ;S 30⁰; see Ref. 12).

When the BP is a singularity this condition cannot be fulfilled, since one always gets a large ~ipmax in the vicinity of the BP, up to 180⁰• One could try to circumvent this problem by exploring the limit of small cell sizes c -+ 0.⁹ This, however, does not remove the core of the problem: micromagnetism is a continuum theory, while the BP at T = 0 is discontinuous.

For finite temperatures, due to thermal fluctuations, there is no reason to assume a constant magnetization magnitude.

Following a mean-field approximation for atomistic spins interacting with a heat bath via the stochastic Landau-Lifshitz- Gilbert equation an improved dynamic equation of motion for the thermal macroscopic magnetization of the ferro- magnet was proposed,13 the Landau-Lifshitz-Bloch (LLB)

equation:

. Ms

M

= -9

M x Herr

+ 9

all M2 (M . Herr) M - 9 a.l M2 M Ms x (M x Herr),

Herr

=

Hd

+ ~V2M

^{_ (M2 -}

I)~.

/-toM; M; 2XII

Here, anisotropy effects are omitted (we focus our attention on permalloy, Py), M is the magnetization, 9 is the gyro- magnetic ratio, Herr is the effective field, all

=

a 2T 13Tc , a.l

=

^{a (1} - T 13 Te), 14 a is the Gilbert damping constant, Tc is the Curie temperature, Ms is the saturation magnetization at T

=

0, M

= IMI,

H" is the demagnetization field, A is the exchange constant, /-to is the vacuum permeability, Me(T) is the equilibrium magnetization [Me(T

=

0)

=

Ms], and XII is the longitudinal susceptibility.

Analytical theories predict that removal of the constraint M(r) = const leads to an important result: The BP is not a singularity anymore.^IO•15 As we will show in this paper this opens new possibilities to micromagnetic modeling of the BP. The use of the LLB equation removes the singularity of the problem, but it remains still a very difficult task to simulate reliably BP-related phenomena. This is because of the different length scales that have to be simultaneously considered: (comparatively) large sample, thin vortex core, and a BP having often a radius that is much smaller than the VC radius. As a result even though we used over 2 million cells, our low-temperature results were still dominated by discretization

.uncertainties. Although first large-scale simulations using the

LLB equation already exist,16.17 magnetostatic interactions were not treated exactly so far. We have implemented the LLB equation in the well-established OOMMF package. 18

II. DETAILS AND RESULTS

In our implementation we have followed the work of Kazantseva et al. 16 without taking thermal noise into account.

The full model of a BP should of course contain fluctuations.

1098-0121/2012/86(9)/094409(5) 094409-1 ©20 12 American Physical Society

First publ. in: Physical Review B ; 86 (2012), 9.- 094409

Konstanzer Online-Publikations-System (KOPS)

LEBECKl, HINZKE, NOWAK, AND CHUBYKALO-FESENKO PHYSICAL REVIEW B 86,094409 (2012)

(c)

FIG. I. (Color online) Permalloy cylinder containing two anti parallel vortices and a Bloch point (BP) in between them. (Left) Schema of magnetization flow (thick arrows) around the BP. The anti parallel vortex cores are marked in blue (lower section) and red (upper section), where the magnetization has primarily a z component. Close to the xy plane the magnetization has primarily rotational components with some inflow tendency described by the angle y > 90°, marked in gray. In other points in space it is a mixture of both these tendencies, approximated by us with model function MLH(r) [see Eq. (5)]. (Right) Cross sections through the sample (results of simulations for T = 600 K). Cross section (a) is in the xy plane, while (b) and (c) are in the xz plane. Color coding in part (b): red (upper) and blue (lower) according to M,(r) value. Color coding in parts (a) and (c): green (gray)-white according to magnetization length-here dark means M(r) ~ Me, while white means M(r) ~ O. Parts (a)-(c) show only the central region of the sample as used for the fitting procedure. Note that the small asymmetry visible in parts (b) and (c) is due to the even number of simulations cells.

Dynamical effects are namely in this case important as the BP is a highly unstable object.³ As already mentioned modeling of a BP is, however, complex anyhow. Thus, in the first approximation we consider a static BP. The BP can be artificially stabilized in numerical simulations (as described later on), but for that it is necessary to omit fluctuations from the modeling procedure.

As already mentioned monitoring of quality of simulations is especially important in the case o"f a BP. To control the validity of variable-M modeling one has to introduce a parameter slightly different from f:..rpmax. We propose to use a parameter _{f:..m max} = max(IMi - Mjl)IMe, where indices i,j enul11er\lte neighboring cells. Following the condition. f:..rpmax ;S 30° one concludes _{f:..m max} ;S 1/2.

We discuss two cases. First, we simulated a cylinder with 64-nm height, equal to its diameter. In this case the discretiza- tion was finer; see details below. This sample was for us a major source of results; if not specifically stated presented results refer to this sample. For a comparison we have also simulated a structure that much better resembles islands available in experiments-a disk with 200-nm diameter and 20-nm thick- ness. Here the discretization is coarser; we used it mainly to check the influence of the sample geometry. As an input for the LLB calculations one has to provide temperature-dependent functions describing material parameters-equilibrium mag- netization Me(T), exchange constant A(T), and parallel susceptibility XII(T). These functions can be evaluated from the mean-field approximation, taken from experiment, or in a more modern way evaluated from the multiscale approach.

In our case we have followed the latter option and followed Refs. 16,19, and 20. For simplicity we scaled the functions calculated earlier for FePt^l6 using the zero-temperature values for Py and its Curie temperature: A(O) = 13 x 1O-¹²J/m, Me(O)

=

0.86 X 10⁶ Aim, Tc

=

870 K. Crystalline anisotropy was neglected; Ci was set to 0.5 (value reasonable for static simulations). For the cylinder sample cubic discretization cells were used. For every temperature we started our simulations

with a sample containing two vortices with the same chirality and opposite polarity with large cell size c = 2 nm (stabilizing the BP is easier if the grid is coarser). After relaxing the system it was used as a starting configuration for simulations with c = I nm and again we have relaxed it. This was used as a starting configuration for our final simulations with c = 0.5 nm. For the disk sample a similar procedure was performed. Here the cells were, however, noncubic, with di- mension: 0.78125 x 0.78125 x 1.25 nm³. All our simulations resulted in magnetization structures qualitatively similar to those shown in Fig. I. The BP is an unstable object. As found out by Thiaville et al.⁹ due to introduction of discretization it is, however, possible to stabilize the BP position. To use this idea we had to use an even number of discretization cells, as the

"BP prefers to sit farthest from mesh points"g [see Figs. I (a) and I(c)].

Galkina et al.¹⁵ and later Elias and Verga 10 have predicted the following dependence of the Cartesian components of the magnetization vector in the close vicinity of a Bloch point centered in r

=

(0,0,0),

(

COS(¢ + y) sin

e)

MGE(r)=Me u(~)

sin(¢+y)sine , cose

(1)

where

e,

¢ refer to spherical coordinate angles. Here, the function u(x) describes the drop of the magnetization mag- nitude (see Fig. 2). It is a solution of complex differential equation; see Eq. (14) in Ref. 15. For fitting purposes we approximated it with a rational function having appropriate boundary properties:

x³

+

0.827x²

+

1.371x

v(x) = x³+O.827x²+2.37Ix+2x 1.371 (2) Both functions, together with their difference, are shown in Fig. 2. The length-scale factor, a, called Bloch point radius in the following, is a distance, where the magnetization length drops roughly by 50%. Analytical theories ^10.15 predict for

KEY ROLE OF TEMPERATURE IN FERROMAGNETIC ...

\.0 r===--'-~--c:::::::=====]

---~

~ _'-' f \ / ""--

;r' 0.5

,- ""l

^0.00;',

'

^\ _{\ .. /}

.

^{/ /}

^{-~-- I}

-0.01 0 2 4 6 0.0 I<..-_ ' - - _L - _ ' - - _ L - _ - ' - - _ L - - - - '

o

² 4 ^(j

Normalized distance from the HI' center, ria FIG. 2. (Color online) Theoretically predicted magnetization drop (normalized) close to the BP center as a function of the distance.

Horizontal axis is scaled similarly as in Ref. 15. Exact function, as got by numerical solving of differential equation u, is shown together with its approximation v [see Eq. (2)]. Difference U - v is shown enlarged in the inset.

this scaling factor the following dependence on the material parameters:

The second important parameter entering Eq. (I) is the

"azimuthal inflow angle" y. In the case of VCs pointing away from the BP this angle exceeds 90° (see Fig. I). This inward tilting reduces the demagnetization energy.7,IO Using an analytical approach for the case of M

=

constand an infinite sample Doring obtained y ~ 112°.7 Following his work Elias and Verga 10 considered a more general case of a variable M. When only the vicinity of the BP is taken into account they obtained y ~ 113°.

Obviously, Eq. (I) is not able to describe the magnetization behavior fUlther away from the BP. For small angles 8 and for larger distances it does not converge to the solution for a vortex core, where the

z

component of the magnetization, Mzvc, is a peaked function around the axis of the vortex (x,y)

=

(0,0).

The function Mzvc had to be incorporated into Eq. (I).

Among several VC profiles described in the literature21 the Gaussian profile, Mzvdx,y) ex exp(-rr2(x 2

+

y2)/8r~c),22 fits our simulation data best. Additionally, one should consider a monotonic dependence of the VC radius on the

z

distance from the BP center, rvdz) (see Fig. I). When searching for a simple analytical function rvdz) that could be used to model results of our simulations we found the best agreement for a dependence rvdz) ex IzI1f3. This is also consistent with vortices in films as shown in the Appendix.

In Eq. (I) the polar angle of the magnetization M(r), 8M

in spherical coordinates, and the polar angle of the vector r are equal: 8M = 8. To improve Eq. (I) by inclusion of the VC shape we thus introduce parametrization:

(4)

where K is a parameter describing the

z

variation of the VC radius. It is equal to the value of the VC radius power 3/2 at a distance along the

z

axis that is equal to the unit of length.

The improved equation that should be valid also for larger

PHYSICAL REVIEW B 86, 094409 (2012)

1.5

0

0

\.0

..

0

0

1.0

S

<

0

^!\m

..s

~

a

^{0.5 mrlX} _{0.5 ci'} ^':;J£

O.OF,..-II-

o

200 400 600 800 Temperature, T (K)

FIG. 3. (Color online) Results of our simulations as a function of the temperature in the range from 0 to Tc. (Left axis) Quality parameter, L:l.m_{m" .} The gray region corresponds to the suggested condition L:l.m_max

:s

1/2 (see text for more details). (Right axis) BP radius a, as returned by our fitting procedure (solid points; fitting errors are smaller than the points). The line shows the result of analytical theory; see Eq. (3).

distances from the BP is, hence,

(5)

where again

e,

1jJ, and

eM

refer to spherical coordinate angles.

Equations (4) and (5) contain three parameters: K, the BPradius a, and the inflow angle y. Using the function MLH(r) we have fitted the results of our simulations in a central region of the 4 x 4 x 4 nm cylinder sample, where K, a, and y were fitting parameters. Because simulations return the magnetization averaged over every discretization cell, in fitting we used three- dimensional convolutions of M_LH. The fitted function MLH agreed well with our simulations; the resulting parameters were not strongly dependent on the chosen fitting region. The results of our fitting procedure are shown in Figs. 3-5.

Figures I(a) and I(c) show clearly that indeed in our simulations the magnetization magnitude drops close to the BP center. On the right axis in Fig. 3 we show the

----

5 5

~ -;v;

~.

2

10

o

^{Cylinder: K}

.6

Disk: K

- - -Infinite film: ¢(T)

--- ---

OL-~--~~--~~--~--'--~~

o

200 400 600 800

Temperature (K)

FIG. 4. (Color online) Parameter determining the variation of the VC radius in the vicinity of the BP, K (points), as returned by our fitting procedure. Circles are for the cylinder sample, while triangles are for the disk. Fitting errors are smaller than the points. For comparison

~(T) from Eq. (AI) is shown.

LEBECKl, HINZKE, NOWAK, AND CHUBYKALO-FESENKO

0 0 0 0 0 0 00«)

100°

... 0

^{Cylinder Disk}

9()O L-.--,_--'-_~_..I.-_,----L_~_....I.-....J

o

200 400 600 800 Temperature (K)

FIG. 5. (Color online) Inflow angle y as returned by our fitting procedure. Circles are for the cylinder sample, while triangles are for the disk. Fitting errors are smaller than the points.

magnetization-drop-related parameter a together with the analytically expected au,(T). The temperature dependence of a,h(T) is dominated by the function XII(T), diverging at Tc and going to zero at low temperatures. a,h(T) and a agree quite well although there is some discrepancy for small temperatures.

We attribute this problem to the relatively large cell size [as compared to a,h(T); see below]. An appropriate description of the BP for these temperatures requires probably smaller cells. On the left axis of Fig. 3 we show the quality parameter L'l.mmax. One can see that for larger temperatures L'l.mmax is indeed smaller than the threshold 1/2 as suggested earlier in the text. The reason we put BP radius (physical quantity) and L'l.mmax (numerical quantity) on the same figure is because they are strongly related to each other. For larger temperatures a,h(T) increases and the magnetization drop in the BP center is more gradual, thus the magnetization changes between neighboring cells are smaller and L'l.mmax gets smaller. The threshold value 1/2 is crossed roughly for a temperature, where the simulation cell size (c

=

0.5 nm) is comparable to the theoretically predicted BP radius. One could thus conclude that BP simulations fulfilling the criterion L'l.mmax ;S 1/2 should meet the condition c ;S a,h(T). This is also consistent with our simulations for larger cell sizes. Fitting the results for the disk sample returned similar results for a.

The parameter K, together with HT) from Eq. (AI), is shown in Fig. 4-separately for the flat island (disk) and for the cylinder sample. For both considered geometries K has quite similar value. This can be a sign that the VC radius drop in the vicinity of the BP is not strongly affected by the sample surfaces. Actually, the values for K and HT) should not be directly compared with each other, as they are related to different geometries: K describes a VC that is close to a BP and measures how strongly the VC radius depends on the distance from the BP. In contrast, HT) is related to a single vortex in a film and measures how strongly its radius depends on the film thickness. Despite these differences there is some analogy between K and HT) and both have comparable value.

This analogy was for us a stimulus leading to Eq. (4). The fact that ~(T) < K might be a little bit surprising. This is because a VC on a surface of a sample containing a static BP can be wider as compared to a case of a VC alone in the same sample-in agreement with earlier reportsY

PHYSICAL REVIEW B 86, 094409 (2012)

The inflow angle y is shown in Fig. 5- again, separately for the disk island and for the cylinder sample. We see a remarkable difference between these two geometries: In the case of the disk sample y is much closer to 90°, which would mean no inflow at all. The angle y depends only on the dipolar interactions and through them on the sample shape.⁷•lo This explains the different values we got in our simulations. It explains also why our results differ from analytical theories: Elias et ai. got 113° close to the BP center and 112° (Doring's result) further away for a round sample.^1o Analytical calculations for flat structures, like our disk sample, are not present in the literature. What we know, however, is that y should not depend on the magnetization value itself.⁷•lo

This is consistent with our results: Despite large changes of the magnetization over the temperature range regarded here, Me(830K)/ Me (OK) ^~ 27%, we find small changes in y for both considered sample geometries.

III. CONCLUSIONS

In summary we have performed a series of simulations revealing details of a static Bloch point (BP) in permalloy in temperatures up to 95% of the Curie temperature. To avoid problems related to the singular character of the BP- affecting all simulations described in the literature so far-we propose to investigate the BP at elevated temperatures. This allowed us to answer an important question: how to perform reliable numerical modeling of BP-related phenomena, like, for example, vortex core reversal. The improvement of our simulations (in comparison to earlier reports) is because of two factors. Major one is the removal of the constraint

I

M

I

= const.

Another factor is the temperature dependence of the exchange length: For higher temperatures this length is slightly larger, thus the simulations are more reliable. Comparison with analytical theories confirmed the validity of our approach.

Performing reliable BP simulations at elevated temperatures is related to a proper choice of sufficiently small simulation cells. This might be difficult to fulfill when using the finite difference method. Among our results only simulations for larger temperatures fulfill the proposed quality condition. On the other hand, it might be easier to accomplish it in finite element simulations, where the cloSe vicinity of the BP can be densely discretized. The described approach can be applied to evaluate dynamical effects, like the vortex core reversal process. Helpful here might also be the proposed model function MLH. Another interesting issue revealed by our study is the relatively small temperature, where the application of the LLB approach turns out to be important. Up to now the use of the LLB equation was reported to be necessary only for temperatures close to the Curie temperature. 16.17

ACKNOWLEDGMENTS

The authors thank Michael Donahue (National Institute of Standards and Technology, Gaithersburg, MD) for valuable help during the implementation phase and Marek Gutowski (IF PAN, Warsaw, Poland) for suggestions regarding the application of convolutions. In the case of O.C. financial support from the Spanish Ministry of Science and Innova- tion under Grant No. FIS20 I 0-20979-C02-02 is gratefully

KEY ROLE OF TEMPERATURE IN FERROMAGNETIC ...

- - - R ',pr(D) - -R'(D)

Q 1 /

>-'

~ O~~--~--~--~--~--~--~~

() 5 10 15 20

Film thickness, D

FIG. 6. (Color online) VC radius as a function of the film thickness. All values are normalized similarly to Ref. 22. The function got by Feldtkeller, R'(D),22 is compared with our approximation, R;pr(D) (shown dashed).

acknowledged. The work in Konstanz was supported by the Kompetenznetz Funktionelle Nanostrukturen.

APPENDIX

We present here results from certain approximations to Feldtkel1er's theory of the Gaussian vortex profile.22 In films,

*kristof.lebecki@uni-konstanz.de

I K. Yamada, S. Kasai, Y. Nakatani, K. Kobayashi, H. Kohno, A. Thiaville, and T. Ono, Nat. Mater. 6, 270 (2007).

2B. Van Waeyenberge, A. Puzic, H. Stoll, K. W. Chou, T. Tyliszczak, R. Hertel, M. Fahnle, H. Bruckl, K. Rott, G. Reiss, I. Neudecker, D. Weiss, C. H. Back, and G. Schutz, Nature (London) 444, 461 (2006).

3R. Hertel, S. Gliga, M. Fiihnle, and C. M. Schneider, Phys. Rev.

Lett. 98, 11720 I (2007).

4K. Nakano, D. Chiba, N. Ohshima, S. Kasai, T. Sato, Y. Nakatani, K. Sekiguchi, K. Kobayashi, and T. Ono, Appl. Phys. Lett. 99, 262505 (20 II).

sR. Moriya, L. Thomas, M. Hayashi, Y. B. Bazaliy, C. Rettner, and S. S. P. Parkin, Nat. Phys. 4, 368 (2008).

6J. Miltat and M. J. Donahue, in Handbook of Magnetism and Advanced Magnetic Materials, Vol. 2, edited by H. Kronmiiller and S. S. P. Parkin (Wiley-Interscience, Chichester, 2007), pp.

742-764; R. Skomski, J. Phys.: Condens. Matter 15, R841 (2003).

7w.

Doring, J. Appl. Phys. 39, 1006 (1968).

8R. Hertel and J. Kirschner, J. Magn. Magn. Mater. 278, L291 (2004).

9 A. Thiaville, 1. M. Garcia, R. Dittrich, J. Miltat, and T. Schrefi, Phys. Rev. B 67, 094410 (2003).

lOR. G. Elias and A. Verga, Eur. Phys. J. B 82, 159 (2011).

PHYSICAL REVIEW B 86, 094409 (2012)

the VC radius has a monotonic dependence on the film thickness R'(D) (see solid line in Fig. 6). The function R'(D) is obtained by minimizing the total energy

* ⁺

'~~

=

^{0, where}

fJ =

Jr/4R', ^EA is the exchange energy, and

ES is the magnetostatic energy.22 For thicker films, where (DfJ)-1 ~ 0 _, one can expand the Taylor series - ~ _lJ.oM,D

,if,

_dfj

=

3lT~'/2 ^(f3b),[1

+

O(fj~)l. Leaving only the first expansion term leads to a final result (see Fig. 6),

R' (D)

=

_Jr____ D¹/ 3

=

~2/3 D1/ 3

(

5/2 A ) 1/3

apr 6fJ-o M; , (AI)

where we have introduced a proportionality factor ~ = (Jr5/2 A/6fJ-o M

;y/2

with the unit of length.

Despite the differences between the main subject of our paper (the Bloch point) and the topic covered by this Appendix, Eq. (A I) is an important result. It helped us to formulate Eq. (4), where the VC radius close to the BP has similar power dependence. Thus, there is some analogy between factors ~

and K; we show them both in Fig. 4.

liS. Gliga, Y. Liu, and R. Hertel, J. Phys.: Conf. Ser. 303, 012005 (2011).

12H. Fangohr, T. Fischbacher, M. Franchin, G. Bordignon, J.

Generowicz, A. Knittel, M. Walter, and M. Albert, NMAG User Manual, Chap. II [http://bit.ly/yWyQXb].

I3D. A. Garanin, Phys. Rev. B 55, 3050 (1997).

14We consider here only the case T ~ Te .

ISE. G. Galkina, B. A. Ivanov, and V. A. Stephanovich, J. Magn.

Magn. Mater. 118, 373 (1993).

16N. Kazantseva, D. Hinzke, U. Nowak, R. W. Chantrell, U. Atxitia, and O. Chubykalo-Fesenko, Phys. Rev. B 77, 184428 (2008).

17 O. Chubykalo-Fesenko, U. Nowak, R. W. Chantrell, and D. Garanin, Phys. Rev. B 74, 094436 (2006).

18M. J. Donahue and D. G. Porter, OOMMFUser's Guide, Release 1.2a4pre, Tech. Rep. 6376 [http://math.nist.gov/oommfj (NIST, Gaithersburg, 1999).

19p' Asselin, R. F. L. Evans, J. Barker, R. W. Chantrell, R. Yanes, O. Chubykalo-Fesenko, D. Hinzke, and U. Nowak, Phys. Rev. B 82,054415 (2010).

20U. Atxitia, D. Hinzke, O. Chubykalo-Fesenko, U. NOWak, H. Kachkachi, O. N. Mryasov, R. F. Evans, and R. W. Chantrell, Phys. Rev. B 82, 134440 (20 I 0).

21 D. Altbir, J. Escrig, P. Landeros, F. S. Amaral, and M. Bahiana, Nanotechnology 18,485707 (2007).

22E. Feldtkeller and H. Thomas, Phys. kond. Materie 4, 8 (1965).