InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Physics » "Vortex Structures in Fluid Dynamic Problems", book edited by Hector Perez-de-Tejada, ISBN 978-953-51-2944-8, Print ISBN 978-953-51-2943-1, Published: March 1, 2017 under CC BY 3.0 license. © The Author(s).

# Statistics of Gyrotropic Vortex Dynamics in Submicron Magnetic Disks

By Gary Matthew Wysin
DOI: 10.5772/64849

Article top

## Overview

Figure 1. The magnetization field M→(r)=Msm^(r) of a vortex centered in a nanodisk with principal axes a = 60 nm, b=30 nm, thickness L=10 nm, from the spin alignment relaxation scheme for the micromagnetics model, Section 2.4. The cell size is acell = 2 nm. Arrows show only the in-plane projection, (mix,miy). Blue line (red open) arrows indicate positive (negative) values of out-of-plane component miz. The core, where miz is larger, appears as a hole in this projection. Even though the system is elliptical, note that the core region remains close to circular.

Figure 2. Vortex structures for an elliptical nanodisk with a=60 nm, b=30 nm, L=10 nm, using cell size acell=2.0 nm in the spin alignment relaxation scheme, including a Lagrange constraint on position. (a) Vortex is held at X=16 nm, Y=0, resulting in total energy E=15.244J after 2000 iterations. (b) Vortex is held at X=0, Y=16 nm, resulting in total energy E=16.001J after 3500 iterations. Compare Figure 1, where X=Y=0 and the energy is lower.

Figure 3. Numerically determined vortex potentials, in units of the effective cell exchange constant J=2AL, for circular Py nanodisks of thickness L=4.0 nm and indicated radii, as found from the Lagrange-constrained spin alignment relaxation. The vortex becomes unstable towards escaping from the disk in the regions of downward curvature.

Figure 4. Motions for one component of vortex position in circular nanodisks from RK4 integration of the LLG equations (shifted for clarity). The damping α=0.02 was turned off after dimensionless time τ=1000. The periods can be calculated accurately from the undamped motion. Graphs of Y(τ) are of the same amplitudes and frequencies but shifted a quarter of a period.

Figure 5. Vortex gyrotropic frequency magnitudes from RK4 (dynamics) simulations for circular nanodisks, with thicknesses ranging from L = 2.0–20 nm, and indicated disk radii, versus force constants scaled by disk thickness, obtained from Lagrange-constrained spin alignment (static) simulations. The dashed line is the theoretical result (43) from the Thiele equation.

Figure 6. Effective potential force constants versus geometric ellipticity b/a, for elliptic nanodisks of semi-major axis a=120 nm, and thickness L=10 nm. Results were found by Lagrange-constrained spin alignment relaxation.

Figure 7. Gyrotropic frequency magnitudes versus geometric ellipticity b/a, for elliptic nanodisks of semi-major axis a=120 nm, and thicknesses L = 5.0 and 10 nm. Results were found by simulations of the LLG equations using RK4 integration. For L = 10 nm, compare the similar shape of the curve of k¯ in Figure 6, as expected from ωG∝k¯ in Eq. (43).

Figure 8. Gyrotropic frequency magnitudes (from dynamics) scaled by mean force constants (from statics), versus geometric ellipticity b/a, for elliptic nanodisks with a=120 nm and two different thicknesses. The results confirm the predictions from the Thiele theory, dashed lines from Eq. (43), using exchange length λex=5.3 nm, with no adjustable parameters.

Figure 9. Spontaneous vortex core motion caused by thermal fluctuations, as found by H2 integration of the LLG-Langevin equations (21) for a nanodisk with thickness L=5.0 nm. The vortex was initiated at the disk center, X=Y=0. This graph shows only 1/5 of the total data generated and used subsequently for analysis of vortex statistics, corresponding to hundreds of gyrotropic periods.

Figure 10. Spontaneous vortex core motion caused by thermal fluctuations, for a nanodisk simulation identical to that in Figure 9, but with double the thickness, L=10 nm. Note the considerably smaller amplitude of gyrotropic oscillations, and the much higher frequency.

Figure 11. The radial distribution of vortex core positions for the simulations in Figures 9 and 10, with a=60 nm, b=30, and thicknesses L=5.0,10 nm. Data out to final time τ=2.5×105 was used. The solid curves are the theory expression (57), using force constants k¯=0.1753 k0 for L=5.0 nm and k¯=0.6832 k0 for L=10 nm, as obtained from spin alignment calculations, with force constant unit k0=A/λex.

Figure 12. Distributions of vortex core coordinates for the LLG-Langevin simulation in Figure 9 with a=60 nm, b=30, and thickness L=5.0 nm. The solid curves are the theory expressions from Eq. (59) based on the Thiele theory for vortex motion, using force constants kx=0.1156 k0 and ky=0.2657 k0 from spin alignment relaxation, where k0=A/λex.

# Statistics of Gyrotropic Vortex Dynamics in Submicron Magnetic Disks

Gary Matthew Wysin
Show details

## Abstract

Topological vortex excitations in thin magnetic nanodisks have attracted a lot of attention because of their dynamic stability and various charge-like properties, which make them suitable objects for data storage. They also have a natural gyrotropic orbital motion that can be described rather well by an approximate Thiele gyrotropic equation for the magnetization dynamics. The gyrotropic oscillation makes them available as a basis for natural oscillators at close to gigahertz frequencies. This gyrotropic motion is excited naturally even by thermal fluctuations. In addition, the gyrotropic oscillation frequency can be affected by external perturbations, which allows possibilities for the design of nanoscale detectors. The vortex moves in an effective potential, strongly determined by the shape anisotropy of the magnetic disk, which then determines the force appearing in the Thiele equation of motion. The motion of an individual vortex within a disk of circular or elliptical shape is considered theoretically, including stochastic thermal effects together with the deterministic gyrotropic effects. From simulations of the motion at different parameter values, a picture of the typical vortex position and velocity distribution within the disk is developed and compared with what is expected from the Thiele equation.

Keywords: magnetic vortex, topological charge, vortex potential, magnetic resonance, magnetic dots

## 1. Introduction: vortices in thin submicron magnetic disks

A cylindrically shaped thin disk of soft ferromagnetic material with a radius a on the order of 100 nm to a few microns and a thickness La  on the order of 10–50 nm provides an interesting system for the study of vortices [1, 2]. A magnetic configuration is described by its local magnetization M(r), which is the magnetic dipole moment per unit volume, at position r. A material is considered with saturation magnetization Ms, which is the magnitude when the medium is completely magnetized along some axis. Due to the demagnetization effect, which is responsible for shape anisotropy, any such magnetic system tends to avoid the formation of magnetic poles on the surfaces, if possible, which would raise the total energy. For a thin circular disk, the local magnetization M(r) as a function of position r may tend to do two things: (1) M(r) will have a strong tendency to point within the plane of the disk [3], if possible; (2) M(r) may then follow the curved circular boundary at the disk edge, thereby completely avoiding the generation of any poles on the edge. This prevents any strong magnetic field lines from passing through the space surrounding the disk edge.

Within the disk, the forces of ferromagnetic exchange cause M(r) also to have a circular structure and remain close to the disk plane. At the disk center, which is a singular point, remaining in the disk plane is impossible, and M(r) then points perpendicular to the plane of the disk, forming tiny north/south poles on opposite faces at the disk center. The resulting circular form of M(r), together with its central poles in a core region, is a magnetic vortex. It is a type of magnetic excitation that is topologically stable and acts in many ways like a particle, when exposed to forces.

#### Figure 1.

The magnetization field M(r)=Msm^(r) of a vortex centered in a nanodisk with principal axes a = 60 nm, b=30 nm, thickness L=10 nm, from the spin alignment relaxation scheme for the micromagnetics model, Section 2.4. The cell size is acell = 2 nm. Arrows show only the in-plane projection, (mix,miy). Blue line (red open) arrows indicate positive (negative) values of out-of-plane component miz. The core, where miz is larger, appears as a hole in this projection. Even though the system is elliptical, note that the core region remains close to circular.

One can also consider deviations from circular symmetry, such as in elliptic nanodisks, where magnetic vortex dynamics has been studied by measuring their radio frequency oscillations [4] and even by direct electrical contact [5] to a nanodisk. An example of a magnetic vortex centered in a thin elliptical nanodisk is shown in Figure 1. It has been obtained from a numerical relaxation algorithm [6], see Section 2.4 below. Although there is a tendency for M(r) to follow the boundary, one sees that the exchange forces are more dominating, and especially the core structure of the vortex retains a circular shape. The locus where the perpendicular component Mz changes sign is obvious as a circle in Figure 1.

### 1.1. Vortex charges

The magnetization profile M(r) of a vortex may point either in a counterclockwise (CCW) or clockwise (CW) direction around the disk. This twofold degeneracy is associated with its circulation charge c=±1, which is also referred to as chirality. It provides one topologically stable geometric property that could be used for data storage in a vortex, if it can be reliably controlled and detected.

In principle, a vortex profile M(r) also has a vorticity charge q=±1, which corresponds to the direction of rotation of M(r) as one moves along a closed path encircling the core. The value q=+1 holds for the vortices described here, which are controlled by demagnetization effects (M being forced to follow the boundary). The value q=1, known as an antivortex, would only be energetically stable if demagnetization effects were not present. The limit of zero thickness would eliminate the relevant demagnetization and make antivortices energetically possible.

The magnetization at the vortex core can take one of two values, M(0)=±Msz^=pMsz^, where p=±1 is the polarization charge. Because there is an energy barrier to flip the core polarization from p=+1 to p=1, it offers yet another charge that could be useful for data storage and manipulation.

### 1.2. Vortex potential and forces

The above-described magnetic vortex will have its minimum energy when it is centered in the disk. The location of the poles (where M points perpendicular to the disk) defines the location of the vortex core, which we denote by position vector R=(X,Y), measured along the x,y Cartesian axes within the disk plane. Because the system is assumed to be thin, only two coordinates X,Y are needed to locate the core. Further, the magnetization itself has little dependence on the coordinate (z) perpendicular to the disk. We take R=(0,0) for the vortex at the center of the disk. It is possible to imagine that the vortex core becomes slightly displaced from the disk center. In that case, a slight deformation of the vortex structure M(r) takes place, while the demagnetization effects at the disk edge still try to maintain M(r) parallel to the edge. The net result of the displacement is a slight increase in total energy. The vortex, as a quasi-particle, lives in some effective potential U(R), which has something approximating a parabolic form [6], with the minimum at the disk center,

 U(R)≈12kFR2, (1)

where kF is a force constant. This further implies an effective force F back towards the disk center, according to the gradient of the potential,

 F=−∇→U(R)≈−kFR. (2)

For a circular disk, the potential is circularly symmetric, and then small displacements lead to a circularly symmetric Hooke’s law type of force. It is also possible to consider magnetic vortices in a cylindrical disk of elliptical shape [7], defined by principal axes a and b< a:

 x2a2+y2b2=1. (3)

This situation leads correspondingly to a modification of the potential also to an approximately elliptic form [8],

 U(R)≈12(kxX2+kyY2). (4)

The parabolic functional form now has separate force constants kx, ky along the two principal axes. It gives a force,

 F=(−kxX,−kyY). (5)

While a vortex in a nanodisk experiences a force directed roughly towards the disk center, its motion tends to be in an orbital sense, which is the gyrotropic oscillation mode [9, 10]. This is discussed further in Section 3 on dynamics. Before coming to that, we begin by a quantitative description of the calculation of vortex structures.

## 2. Analysis of quasi-stationary vortices in a nanodisk

The theoretical analysis is based on the statics and dynamics of the magnetization field, which is now assumed to keep a uniform magnitude Ms, but a spatially varying direction, by writing M(r)=Msm^(r), where m^(r) is a unit vector. From the energetics expressed in terms of m^(r), equations for the vortex structure and motion can be developed. See Ref. [11] for a general discussion of the calculation of magnetic vortex structures and properties.

### 2.1. Energetics of a continuum nanomagnet

The system is governed by ferromagnetic exchange energy and interactions of M generally with the demagnetization field HM (self-generated by M) and any possible externally applied field Hext. A continuum Hamiltonian for the system is

 H=∫dV{A∇m→⋅∇m→−μ0(H→M+12H→ext)⋅M→}, (6)

where μ0 is the permeability of free space, and A is the exchange stiffness. One commonly used material is Permalloy-79 (Py, 79% nickel, 21% iron), with exchange stiffness about 13 pJ/m and saturation magnetization Ms=860 kA/m [12]. The magnetization changes its direction over a length scale λex called the exchange length. Exchange energy of the order A/λex2 competes with demagnetization energy of the order 12μ0Ms2. Equating these terms gives the definition of the length scale,

 λex=2Aμ0Ms2. (7)

For Py, λex5.3 nm. Exchange forces dominate over lengths less than λex, but demagnetization dominates over larger lengths, allowing the M(r) field to change its direction. At a boundary, the exchange effects are much less present, and demagnetization helps M(r) to point parallel to the boundary, if possible.

### 2.2. The demagnetization field H→M in a thin magnetic film

The demagnetization field is determined by the global configuration of the magnetization of the system; it is derived from considerations of magnetostatics (see Ref. [13] for the details of the approach used here). In the absence of an external applied magnetic field, one has magnetic induction B=μ0(HM+M). Gauss’ law, .B=0, then becomes

 ∇→⋅H→M=−∇→⋅M→. (8)

By assuming the demagnetization field comes from a scalar potential via HM=ΦM, a Poisson equation for the magnetostatics is obtained:

 −∇2ΦM=ρM,  ρM≡−∇→⋅M→. (9)

Therefore, the magnetization M(r) produces an effective magnetic charge density ρM, which is the source in this Poisson equation. The solution for scalar potential ΦM can be obtained by various numerical methods. Generally, we have used a scheme based on discretization of the system (introduction of a spatial grid), together with appropriate Green’s functions for the Poisson equation. In addition, it is extremely helpful to use the approximation that the disk is very thin, L a, where a is the radius (or the semi-major axis for an elliptical disk). In this case, both M and HM do not depend on the vertical coordinate z along the disk axis. Then, the problem can be solved by effective two-dimensional (2D) Green’s functions [14]. The components α=x,y,z of the demagnetization field can be expressed as 2D convolution integrals,

 HαM(r)=∫d2r′∑β=x,y,zGαβ(r−r′)Mβ(r′) (10)

where Gαβ(rr) represents a tensor Green’s function, and the integration is over 2D positions, for example now r=(x,y). The evaluation of these integrals can be accelerated through the use of fast Fourier transforms [15].

### 2.3. Discretization and micromagnetics for simulations

For numerical solutions of the magnetization field M(r)=Msm^(r), it is necessary to partition the system into cells labeled by index i for positions ri. We use a square grid of cells of individual size acell × acell × L, with acell=2.0 nm, and disk thickness L = 5 nm and L = 10 nm. At the center of each cell is a unit direction vector m^i, whose motion is to be followed. Each cell contains a magnetic dipole moment μi of magnitude μ = Lacell2Ms and direction m^i. This micromagnetics approach [16, 17] then represents the original continuum system, but with a discretized 2D micromagnetics Hamiltonian,

 H=−J[∑(i,j)m^i⋅m^j+acell2λex2∑i(H˜iext+12H˜iM)⋅m^i], (11)

where the effective exchange constant and energy scale between nearest-neighbor cells is J=2AL, and magnetic fields have been brought to dimensionless forms,

 H˜iext=H→iext/Ms,  H˜iM=H→iM/Ms. (12)

The presence of the factor acell2/λex2 gives the relative strength of demagnetization effects compared with exchange effects. For the micromagnetics approach to be valid, this factor should be much less than 1. The transverse cell size acell should then be less than the exchange length. The micromagnetics approach, with the assumption that only the direction of M is changing, makes the implicit assumption that demagnetization effects are a perturbation on exchange effects. Obviously, the Green’s functions Gαβ(rr) must also be brought to a discrete form to carry out the calculation of H˜iM.

The Hamiltonian can be used to define the net magnetic inductions that act on each cell’s magnetic dipole μi, according to

 B→i=−δHδμ→i=Jμb→i, (13)

where the dimensionless magnetic inductions are

 b→i=∑j∈z(i)m^j+acell2λex2(H˜iext+H˜iM). (14)

The first term involves a sum over the nearest neighbors z(i) of cell i; it is the exchange field. The second term represents the combination of external and demagnetization fields. The effective strength of magnetic inductions is indicated by the unit we use for simulations,

 B0≡J/μ=μ0Msλex2/acell2. (15)

In the results presented here with acell=2.0 nm, and Py parameters, one has μ0Ms21.08 T and B07.59T. B0 gives the order of magnitude of the exchange fields; the demagnetization fields are considerably weaker.

### 2.4. Static vortex configurations from a relaxation scheme

Static vortex configurations are derived as the stationary solutions of the dynamic equations of motion. At zero temperature, the undamped dynamic equation of motion is a simple torque equation for each magnetic dipole, which interacts with its local net field:

 dμ→idt=γμ→i×B→i. (16)

Note that this holds because μi/γ=Si is the spin angular momentum of the cell, whose time derivative is the torque, μi×Bi. The equation can be written in terms of the dimensionless quantities, also defining a dimensionless time τ,

 dm^idτ=m^i×b→i, τ≡γB0t. (17)

The unit of time used for simulations is t0(γB0)1. For Py parameters, it takes the value t00.75 ps. Its reciprocal also defines a simulation frequency unit, f0γB01.336 THz.

For static configurations, however, the time derivatives in Eq. (17) are zero. This implies that each dipole μi or its unit vector m^i must align with the local field in that cell, bi.

Thus, an algorithm that iteratively points each m^i along its current value of bi will tend to move the system towards a static configuration. We call this approach a spin alignment relaxation scheme [18]. To carry it out, some initial state must be chosen from which to begin the iteration. Assume that the direction vectors are defined in terms of spherical planar angles (ϕi,θi), according to

 m^i=(cosθicosϕi,cosθisinϕi,sinθi). (18)

In this notation, ϕi is referred to as the in-plane angle and θi is the out-of-plane angle, which can be positive or negative. The approximate in-plane structure of a vortex located at position R=(X,Y) in a disk can be expressed using the vorticity q=+1 as

 ϕi=q tan−1(yi−Yxi−X)+ϕ0 (19)

where (xi, yi) is the 2D location of micromagnetics cell i, and ϕ0=cπ2 depends on the vortex circulation charge. There is not a corresponding analytic form for the out-of-plane component. Instead, one can start with all θi=0, that is, a planar vortex. However, the iteration will be such that all θi will remain zero, unless some small nonzero deviation is included. Therefore, small random values of θi can be used for the initial state. A nonzero out-of-plane component will then grow naturally as the system relaxes into a vortex state. The process is repeated until the changes in the m^i become insignificant (less than about one part in 108).

For a circular or elliptic disk, if the vortex is initiated away from the center, as the spin alignment relaxation proceeds, the vortex will be found to both develop an out-of-plane component and also move to the disk center. Spin relaxation is an energy minimization algorithm; the system moves to its nearest minimum energy state, which is that configuration centered in the disk. A profile of a vortex obtained this way in an elliptic disk is shown in Figure 2. The projection of the dipoles onto the disk plane is shown. Note that there is a core region with a radius of the order of λex (region where M has significant out-of-plane components, appearing as a hole in the diagram). Interestingly, the core tends to keep a reasonably circular form, as seen by the locus of points where the sign of Mz reverses.

### Figure 2.

Vortex structures for an elliptical nanodisk with a=60 nm, b=30 nm, L=10 nm, using cell size acell=2.0 nm in the spin alignment relaxation scheme, including a Lagrange constraint on position. (a) Vortex is held at X=16 nm, Y=0, resulting in total energy E=15.244J after 2000 iterations. (b) Vortex is held at X=0, Y=16 nm, resulting in total energy E=16.001J after 3500 iterations. Compare Figure 1, where X=Y=0 and the energy is lower.

### 2.5. Effective potentials of a vortex in a nanodisk

Spin alignment relaxation can also be used to estimate the effective potential U(R) for the vortex, by including a constraint on the vortex position R=(X,Y). The effective potential is the system energy less any constraint energy, for a chosen vortex position. A constraint on vortex position can be enforced with the use of Lagrange’s undetermined multipliers [6]. Physically, a vortex can be shifted away from the disk center, by the application of a magnetic field within the disk plane. A uniform field Hext along ±x will displace the vortex along ±y, and vice versa, with the sign determined by the vortex chirality. Buchanan et al. [8] were able to map out the vortex potential energy numerically using the field to move the vortex to different equilibrium positions. This gives one way to obtain the effective force constants kx and ky.

Rather than using a uniform applied field, it is possible to imagine the application of a spatially varying field, which primarily acts on the core region of the vortex. These fictitious extra fields are the undetermined Lagrange multipliers; they are determined through course of the calculation. Simultaneously, another constraint is applied that ensures unit length for the direction vectors m^i. The fictitious fields exert torques on the cells in the core region, which hold the vortex in the desired location, without significantly changing the overall vortex structure. Thus, a quasi-static vortex structure can be obtained numerically, for whatever position is desired, within reason. The approach works best for a vortex near the disk center. For the same elliptical disk of Figure 1, the vortex has been relaxed by this scheme to positions 16 nm from the center, in Figure 2. Note that the energy is higher for a displacement along the shorter axis of the ellipse [8].

The work here considers stable vortex states. It should be kept in mind that for some parameters or disk sizes, the vortex could become unstable towards the formation of a lower energy quasi-single-domain state (nearly uniform M(r)), or some other multi-domain state without a vortex. This is especially likely in the case of elliptic disks with a high aspect ratio (b a), where demagnetization will strongly favor M aligning with the longer axis [7]. The vortex state will also become unstable in a circular disk if it is too thin, which minimizes the demagnetization forces from the circular edge, which usually stabilize a vortex. Also, if the disk is too thick (L ≫  a), again, demagnetization will cause M to approximately align with the long axis and a vortex will not be stable.

Typical vortex potentials obtained from Lagrange-constrained spin alignment for circular nanodisks are shown in Figure 3, for various radii with fixed thickness L=4.0 nm. The minimum energy region is close to parabolic form; however, as the vortex is placed closer to the edge, a lack of stability (downward curvature) appears. Using the interior region of the potential, the effective force constants kF for circular disks or even kx and ky for elliptical disks can be estimated quite accurately. In the example of Figure 3, one can observe that kF becomes smaller for the larger radii disks. See Refs. [6, 7] for further details.

In elliptical disks [4, 7, 8], the force constant for displacement along the longer disk axis is found to be weaker than that along the shorter disk axis; see Figures 1 and 2. Thus, the potential acquires an elliptical shape that is determined by the original geometrical shape of the disk. For a disk with semi-major axes a and b with b< a, we have found that for adequately large nanodisks and b of sufficient size to stabilize the vortex, the ratio of force constants asymptotically approaches the relation,

 kx/ky≈b/a. (20)

This has the correct limit for a circular disk, kx=ky. The relation is summarized by saying that the geometric ellipticity, b/a, directly determines an energetic ellipticity, kx/ky. The energetic ellipticity can be seen to determine the shape of the elliptical vortex orbits at constant energy in the phase space.

### Figure 3.

Numerically determined vortex potentials, in units of the effective cell exchange constant J=2AL, for circular Py nanodisks of thickness L=4.0 nm and indicated radii, as found from the Lagrange-constrained spin alignment relaxation. The vortex becomes unstable towards escaping from the disk in the regions of downward curvature.

## 3. Magnetic dynamics and the Landau-Lifshitz-Gilbert-Langevin equations

The dynamics described by Eq. (16) or its dimensionless form in Eq. (17) is not completely realistic, because it does not include the effects of damping nor of temperature and its statistical fluctuations. Both the damping and thermal effects could be quite large on a vortex. When only damping with a dimensionless parameter α is included, the well-known Landau-Lifshitz-Gilbert (LLG) dynamic equation [19, 20] is obtained [Eq. (21) but with all bs,i=0]. Here, we take that one step further and also include stochastic magnetic fields bs,i(t) that represent the effects of temperature. This leads to a Langevin equation derived from the LLG equation [21], for an individual micromagnetics cell,

 dm^idτ=m^i×(b→i+b→s,i)−αm^i×[m^i×(b→i+b→s,i)]. (21)

The changes in m^i are a superposition of deterministic effects (from bi) and stochastic effects (from bs,i). The stochastic fields act to bring the system to thermal equilibrium. That takes place provided their correlations follow the fluctuation–dissipation (FD) theorem, which can be written for this problem in the dimensionless quantities as (site index i is suppressed)

 〈bsλ(τ) bsλ′(τ′)〉=2α T δλλ′ δ(τ−τ′). (22)

Here, δλλ′ is a Kronecker delta and the indices λ,λ′ refer to any of the Cartesian coordinates; δ(ττ) is a Dirac delta function. The dimensionless temperature is thermal energy scaled by J,

 T≡kTJ=kT2AL, (23)

where k is Boltzmann’s constant and T is the absolute temperature. The FD relation indicates how the stochastic magnetic fields move energy into and out of the system, in random processes that nevertheless can be quantitatively measured. The stochastic fields are included only when a study of temperature effects in real time is desired. They can be set to zero if the zero-temperature dynamics is of interest, producing the LLG equation. Below we use solutions of Eq. (21) obtained appropriately for the type of system under study, be it T=0 or T>0.

### 3.1. The Thiele equation for vortex core motion

Magnetic excitations such as domain walls and vortices do not obey Newtonian dynamics. Instead, it can be shown from magnetic torque considerations (i.e., analysis of the T=0 LLG equations) that the steady-state dynamics of the core velocity V=R˙ is described by a Thiele equation [22],

 F+G×V=0. (24)

The motion is governed by the gyrovector G, which for vortices is a vector that points perpendicular to the disk plane, in a direction determined by the magnetization at the vortex core. Consider a material with saturation magnetization Ms. In terms of the magnetization per unit area, m0MsL, and the gyromagnetic ratio γ=1.76×1011 T−1 s−1 of an electron (its magnetic moment divided by its angular momentum), the gyrovector of a vortex is

 G=Gz^=2πpqm0z^/γ. (25)

For the vortices in a disk, q=+1, while there are two core polarizations p=±1 possible. The gyrovector points perpendicular to the disk in two possible directions. A solution of the Thiele equation then gives a description of the motion of a vortex, provided it remains as a particle-like stable object under the dynamic environment it is found in. A general review of vortex motion obeying a Thiele equation, even including an intrinsic mass, is given in [11].

Here, we suppose that a vortex is moving within a nanodisk of elliptical shape, at position R=(X,Y), with the force in Eq. (5) acting on it. One finds that it makes an elliptical orbital motion, whose gyrotropic frequency can be estimated from the Thiele equation [7]. A solution for the vortex velocity is obtained quickly by taking the cross product of G=Gz^ with the Thiele equation,

 G×F+G×(G×V)=0. (26)

A vector identity is useful,

 G×(G×V)=(G⋅V)G−(G⋅G)V. (27)

The vortex velocity points in the plane of the disk, but G is perpendicular to that plane, so GV=0. This gives the velocity as

 V=G×FG2=1G(kyY,−kxX). (28)

With V=(X˙,Y˙), this is a pair of first-order differential equations, which can be directly integrated, starting from some initial vortex position R(0)=(X0,Y0). An elementary calculation gives elliptical motion, with instantaneous coordinates

X(t)=X0cosωGt+Y0kyGωGsinωGt
 Y(t)=Y0cosωGt−X0kxGωGsinωGt (29)

where the gyrotropic frequency is determined by the geometric mean of the force constants, k¯kxky:

 ωG=−kxkyG=−k¯G. (30)

The negative square root is used, because a vortex with G along +z^ and a centrally directed force will move in the clockwise (or negative) direction in the xy plane. This result applies even when the vortex equilibrium position is displaced from the disk center by an applied magnetic field, using the effective force constants at that displaced location [8]. Buchanan et al. [8] found that the experimentally measured vortex oscillation frequency can be controlled by the application of an in-plane field, Hext; especially, Hext along the short (or hard) axis of the ellipse displaces the vortex on the long axis, where its frequency increases substantially due to position-dependent increases in both force constants with X.

With Hext=0, one can find the shape of the elliptical orbit and compare with the shape of the nanodisk. The vortex in undamped motion must move along an equipotential centered in the disk. The orbital energy U is found to be

 U=12(kxX2+kyY2)=12(kxX0+kyY0)2. (31)

Dividing through by the constant, U, this is the standard equation of an ellipse, with the semi-major and minor axes Xmax,Ymax, given by

 A=2Ukx,  B=2Uky. (32)

Their ratio is then

 BA=kxky≈ba. (33)

The last approximate result in terms of the disk axes a,b was obtained by using relation (20), valid only in the limit of larger ellipses. Thus, the shape of the vortex orbit is nearly the same as the shape of the nanodisk. The energetic ellipticity (not to be confused with the eccentricity),

 e≡kxky, (34)

determines the ratio of the orbital axes. Indeed, the potential can be brought to a circular form, with a new coordinate ρ:

 U(ρ→)=12k¯ρ→2,  ρ→≡(eX,1eY). (35)

Then, it is possible to show that the velocity follows a typical expression for circular motion,

 ρ→˙=(ρ˙x,ρ˙y)=ω→G×ρ→. (36)

where ωG=ωGz^. This equivalent circular coordinate ρ is useful for the analysis of vortex position statistics in elliptical disks.

### 3.2. Numerical methods for magnetization dynamics

The analysis of vortex motion via the Thiele equation is expected to be approximate. Numerical simulations can be used to give a more complete and reliable description of the dynamics. We require the time evolution from Eq. (21) solved either for zero temperature or finite temperature. These results are generated for Py parameters, based on the exchange length of λex=5.3 nm, together with a micromagnetics cell size of acell=2.0 nm.

#### 3.2.1. Zero temperature: fourth-order Runge-Kutta

At zero temperature, a stable integration scheme is the well-known fourth-order Runge-Kutta (RK4) scheme, which we have used. A time step in dimensionless simulation units of Δτ = 0.04 is sufficient to insure good energy conserving dynamics (at zero damping), resulting in energy conservation to one part in 1012 over as many as 5.0×105 time steps, in systems with up to 4000 cells. To insure this high precision control of the energy, it is essential to evaluate the demagnetization field continuously during every substep of the RK4 algorithm. In the zero temperature simulations used to estimate gyrotropic frequencies, the initial state of the dynamics is a vortex obtained by spin alignment relaxation to a desired position. It is also helpful to run the time evolution initially with some weak damping (α0.02) for a limited time, followed by energy conserving dynamics (α=0). The inclusion of damping for a short interval helps to remove any spin wave oscillations that may be generated by a less than perfect initial vortex state. The subsequent energy conserving dynamics then gives precise estimates of the frequencies ωG.

#### 3.2.2. Finite temperature: Langevin dynamics via second-order Heun method

For finite temperatures, the Eq. (21) has been solved effectively by a second-order Heun method (H2) [21, 23]. This scheme is equivalent to a two-stage predictor-corrector algorithm, where the predictor stage is an Euler step and the corrector stage is the trapezoid rule. Both stages use the same random fields bs,i, which are produced by a random number generator. Any of the Cartesian components of these fields are to be random deviates with a zero mean value and a variance that must depend on both the dimensionless time step and temperature according to

 σs=2αT Δτ. (37)

This is a result of the FD theorem Eq. (22), and it is used to replace the stochastic fields integrated over a time step, by the relation

 ∫τnτn+Δτdτ b→s,i(τ)→σsw→i. (38)

The vectors wi are triplets of random deviates with zero mean and unit variance, for each site i. In usual programming, there are standard random number generators, which return a uniform deviate from 0 to 1, with a variance found to be 1/12. These can be shifted into the range from 0.5 to +0.5 and then rescaled by 12σs to get stochastic field components of the correct mean and variance (it does not need to be a Gaussian distribution) [13, 7].

## 4. Vortex gyrotropic motion at zero temperature

In a circular nanodisk at zero temperature, with a radial force F as in Eq. (2), the analysis from the Thiele equation (24) shows that the vortex velocity always points along the azimuthal direction:

 V=Gz^×FG2=−γkFR2πpqLMsϕ^. (39)

The minus sign indicates that a vortex with positive gyrovector (along z^) will move in the clockwise or negative sense, in uniform circular motion, and oppositely for those with negative gyrovector. More generally, for elliptic nanodisks, the predicted gyrotropic frequency obtained from the Thiele equation is

 ωG=−k¯G=−γk¯2πpqLMs, (40)

where k¯ is the geometric mean of the force constants along the principal axes. For a circular system, k¯kF. These results depend strongly on the force constants, which can be estimated from the static vortex configurations. It has been found [9, 24, 13] that for sufficiently large circular disks far from any stability limits of the vortex, the force constants are very roughly proportional to L2/R, that is,

 kF≈14L2RAacell2≈0.878 μ0Ms2L2R. (41)

The frequency unit f0=t01=γB0 used in the simulations depends on the cell size, which is inconvenient for comparison with experiment. Thus, it is important to convert the results to a commonly used frequency unit,

 ω0≡μ04πγMs. (42)

This is ω0=γMs in the centimeter-gram-second (CGS) system of units. With the help of definition (7) for the exchange length, expression (40) for gyrotropic frequency can be written,

 pωG=−(μ04πγMs)(k¯λexA)(λexL). (43)

With vorticity q=+1 assumed, the sign of ωG is determined by the core polarization p. This expression suggests using k0A/λex as the unit of force constant and λex as the unit of length.

Simulations can verify the frequency predictions from the Thiele equation. As shown below in some examples, the dimensionless periods τG of gyrotropic motion can be estimated precisely, in simulation time units (τ=t/t0). Dimensionless angular frequencies are then 2π/τG, which are given physical values by multiplying by f0=γB0. Using Eq. (15), these can then be converted into units of ω0 as follows:

 ωG=2πτGf0=2πτGγμ0Msλex2acell2=2πτGλex2acell2 4π ω0. (44)

We use this below to convert the raw numerical data (τG) into frequencies in ω0 units.

Of course, to get precise estimates of the frequency, the vortex must be instantaneously located to high precision. That is a two step process. The first step is to use the singularity in the in-plane magnetization angle ϕ, to locate the four cells nearest to the vorticity center, Rv, defined implicitly according to the relation

 ∇→×∇→ϕ(r)=2πz^δ(r−Rv). (45)

For the micromagnetics square grid, the vorticity center falls between the four cells that have a net 2π circulation in ϕ. This gives the location RRv only to a precision equal to the cell size. It can be greatly improved by making a weighted average of the cell locations ri, using their squared out-of-plane components, which are largest in the vortex core, as the weighting function:

 R=∑|ri−Rv|<4λex(miz)2 ri∑|ri−Rv|<4λex(miz)2. (46)

For better efficiency, the sum is restricted to cells within four exchange lengths of the vorticity center. This avoids using useless data from the core of interest (e.g., spin wave oscillations near the edge of the disk should be ignored). As the vortex moves, the resulting estimate for R changes smoothly. This algorithm even works very well for vortices moving in response to thermal fluctuations.

### 4.1. Circular nanodisks simulations

Some typical vortex motions in circular nanodisks of radius a=120 nm are presented in Figure 4, as obtained from integration of the LLG equations by the RK4 scheme. The initial states came from Lagrange-constrained spin alignment to the initial position R=(4.0,0) nm. A weak damping with parameter α=0.02 was included but turned off at dimensionless time τ=1000. The remaining evolution was used to estimate the periods, τG, which are then converted using Eq. (44).

For the motions displayed in Figure 4, the dimensionless periods for L = 5 nm, 10 nm, and 20 nm are τG5800, 3270, and 1872, respectively. From statics calculations of the effective potentials as described earlier, the corresponding raw force constants are kf acell /A0.033863, 0.120192 and 0.419143, respectively, using acell = 2.0 nm. Rescaling by a factor λex/acell = 5.3/2.0 converts them into kFλex/A, which appears in the Thiele theory expression (43). For these and other similar vortex motion simulations with L ranging from 2.0 to 20 nm, one can compare to the Thiele prediction by plotting the frequency ωG versus kF/L with units as suggested from Eq. (43) (see Figure 5). Note that for a given radius a, the disk with the smallest L has the largest frequency. The result is that ωG, obtained from dynamics simulations, is very close to linearly related to kF/L, obtained from static simulations, with a unit slope for these units. This gives a strong verification of the Thiele equation being applicable to vortex motion in nanodisks where the vortex is stable. Note that all simulations here used a reasonably small vortex orbital radius of about 4.0 nm, avoiding having the vortex core approach the disk edge, which would tend to destabilize the vortex.

### Figure 4.

Motions for one component of vortex position in circular nanodisks from RK4 integration of the LLG equations (shifted for clarity). The damping α=0.02 was turned off after dimensionless time τ=1000. The periods can be calculated accurately from the undamped motion. Graphs of Y(τ) are of the same amplitudes and frequencies but shifted a quarter of a period.

### Figure 5.

Vortex gyrotropic frequency magnitudes from RK4 (dynamics) simulations for circular nanodisks, with thicknesses ranging from L = 2.0–20 nm, and indicated disk radii, versus force constants scaled by disk thickness, obtained from Lagrange-constrained spin alignment (static) simulations. The dashed line is the theoretical result (43) from the Thiele equation.

### Figure 6.

Effective potential force constants versus geometric ellipticity b/a, for elliptic nanodisks of semi-major axis a=120 nm, and thickness L=10 nm. Results were found by Lagrange-constrained spin alignment relaxation.

### Figure 7.

Gyrotropic frequency magnitudes versus geometric ellipticity b/a, for elliptic nanodisks of semi-major axis a=120 nm, and thicknesses L = 5.0 and 10 nm. Results were found by simulations of the LLG equations using RK4 integration. For L = 10 nm, compare the similar shape of the curve of k¯ in Figure 6, as expected from ωGk¯ in Eq. (43).

### 4.2. Elliptical nanodisk simulations

Simulations for elliptic nanodisks [7] offer an even wider range of possibilities, because the variations with geometric ellipticity b/a can be studied. For instance, the variation of the effective potential force constants has a behavior like that in Figure 6, for the particular case a=120 nm and L=10 nm. Both kx and ky were determined from the potentials derived by spin alignment with a position constraint. Their geometric mean k¯, which determines gyrotropic frequencies, is also shown. The curves for these force constants do not go below a minimum value of b/a, where the vortex becomes unstable.

The corresponding gyrotropic frequencies ωG for L=10 nm and also for L=5.0 nm are shown in Figure 7, versus b/a. These were obtained from simulations the same as those described for circular nanodisks. Note that the shapes of these curves are very similar to the curves of k¯ in Figure 6, which is to be expected if the Thiele theory (43) is valid. The additional results for L=5.0 nm are included to demonstrate the dependence on disk thickness. With thicker disks having a greater restoring force and k¯L2, due to the extra area on the disk edges, the dependence of GL results is gyrotropic frequencies increasing roughly linearly in L. The results can be presented in another view in Figure 8, showing ωG/k¯ versus ellipticity for different L. One again gets a clear and quantitative verification of the Thiele theory result (43), seeing that ωG/k¯λex/L with the correct constant of proportionality.

### Figure 8.

Gyrotropic frequency magnitudes (from dynamics) scaled by mean force constants (from statics), versus geometric ellipticity b/a, for elliptic nanodisks with a=120 nm and two different thicknesses. The results confirm the predictions from the Thiele theory, dashed lines from Eq. (43), using exchange length λex=5.3 nm, with no adjustable parameters.

## 5. Spontaneous gyrotropic motion from thermal fluctuations

Now we consider the effects of temperature on a vortex. Specifically, the temperature effectively acts as a bath of random magnetic fields that exchange torques and energy with the vortex. Even though that exchange is somewhat random, one sees that it is able to spontaneously initiate the organized gyrotropic motion of the vortex. That motion proceeds over a noisy background of spin waves. Even so, it is readily apparent and persistent. Here, we show typical time evolutions, and then later discuss statistical properties.

### 5.1. Simulation of a vortex initially at disk center

A vortex that has been relaxed to its minimum energy configuration (e.g., by the spin alignment scheme) is situated in the disk center, whether it be circular or elliptical. This assumes that a quasi-single-domain state is not lower in energy. Then, in the absence of any external forces or forces due to a thermal environment, it would statically remain centered in the disk and exhibit no dynamics. However, Machado et al. [25] noticed that finite temperature micromagnetics simulations demonstrate the spontaneous motion of the vortex, even if it starts in it minimum energy location. This is rather surprising, although it is really not much different than any spin wave mode from being excited thermally in an equilibrium system with temperature. From the point of view of statistical mechanics, any excitable modes (i.e., independent degrees of freedom) should share equally in available thermal energy, and because the energy present in the vortex gyrotropic motion is quite small, rather large orbital motions can develop solely due to the effects of temperature.

In the numerical solution [13] of the magnetic Langevin equation (21), the dimensionless temperature is required. For the simulation units being used, J=2AL determines the energy scale and depends on the disk thickness. As an example, we consider a disk with a=60 nm, b=30 nm, and thickness L=5.0 nm, at temperature T=300 K. For Py parameters (A=13 pJ/m), the energy unit is J=130 zJ, while the thermal energy scale is kT=4.14 zJ, which gives the dimensionless temperature,

 T=kT2AL=0.032, forT=300K,L=5.0nm. (47)

This was used to determine the variance of the random magnetic fields, Eq. (37), together with a damping parameter α=0.02. A dimensionless time step Δτ=0.01 for the second-order Heun method was used. The resulting vortex core coordinates (X(τ),Y(τ)) are displayed in Figure 9, out to a time of τ=50,000. From Figure 9, a clockwise orbital motion takes place, together with a noisy background, and there are about 15 complete orbits for τ<50,000 (period τG3300). The period is somewhat longer than that found at zero temperature, τG2970. This softening of the mode with temperature is to be expected. In addition, the amplitudes of X and Y motions are not equal, as expected from the elliptical disk shape. The gyrotropic orbital motion continues indefinitely; it was followed out to τ=2.5×105 to get vortex statistics.

For comparison, an identical simulation but with the disk thickness increased by a factor of 2 to L=10 nm in shown in Figure 10, again starting the vortex from the disk center. The greater thickness approximately quadruples k¯, but also doubles the gyrovector, thereby resulting in the frequency being double that for L=5.0 nm. It is also clearly apparent that the amplitude of the thermally induced motion is reduced in the thicker nanodisk (the graphs have different vertical scales). These differences then are primarily driven by the modifications to the force constants and to G.

### Figure 9.

Spontaneous vortex core motion caused by thermal fluctuations, as found by H2 integration of the LLG-Langevin equations (21) for a nanodisk with thickness L=5.0 nm. The vortex was initiated at the disk center, X=Y=0. This graph shows only 1/5 of the total data generated and used subsequently for analysis of vortex statistics, corresponding to hundreds of gyrotropic periods.

### Figure 10.

Spontaneous vortex core motion caused by thermal fluctuations, for a nanodisk simulation identical to that in Figure 9, but with double the thickness, L=10 nm. Note the considerably smaller amplitude of gyrotropic oscillations, and the much higher frequency.

### 5.2. Thermal vortex motion as described by the Thiele equation

Next, we consider the statistical mechanics of the vortex core position R=(X,Y), based on an effective Lagrangian and Hamiltonian that give back the Thiele equation. The analysis [7] makes use of the general elliptic potential U(R) in Eq. (4). It is straightforward to check that a Lagrangian whose Euler-Lagrange variations gives back the Thiele equation is [13]

 L=−12G(XY˙−YX˙)−12(kxX2+kyY2) (48)

This is a particular choice of gauge and this Lagrangian is not unique (see Ref. [26] for a different choice). To transform to the associated Hamiltonian, one finds the canonical momentum for this symmetric gauge,

 P=∂L∂V=12(GY,−GX). (49)

This shows that the Lagrangian can be expressed as L=PVU. As P is determined by X and Y, without any time derivatives, one can interpret this as a pair of constraint relations between components of P and R. It means that of four original coordinates plus momenta, only two are independent.

The Hamiltonian is obtained in the usual way,

 H=P⋅V−L=U=12(kxX2+kyY2)=12k¯ρ2. (50)

Curiously, this has no momenta present. This strange situation seems to imply that there is no dynamics, because the Hamilton equations of motion are

 P=−∂H∂R,  X=∂H∂P. (51)

That would give V=R˙=0, which is clearly wrong. This singular situation comes about because of the constraint (49) between momentum and position components. In order to get a true dynamics, one needs to rewrite the Hamiltonian half as a potential part and half as a kinetic part, that is,

 H=14(kxX2+kyY2)+14(2G)2(kxPy2+kyPx2). (52)

This is exactly equal to H in Eq. (50), but now it does give back the Thiele equation when its time dynamics is found from Eq. (51). Because of the constraint, the vortex motion depends on only two independent coordinates, or degrees of freedom, rather than four. For the purposes of statistical mechanics, then, the thermalized vortex motion contains an average energy, H=2×12kT.

### 5.3. Thermalized vortex probability distributions from the Thiele equation

One can assume that any of the coordinates, X,Y,PX,PY, as well as effective circular coordinate ρ=(ρx,ρy), obey a Boltzmann distribution, whose parameters are determined by the average energy,

 〈H〉=kT. (53)

This directly gives the mean squared effective circular radius for an elliptic disk,

 〈ρ2〉=〈2H/k¯〉=2kT/k¯. (54)

This becomes the usual mean squared radius, ρ2R2, in the limit of a circular disk. Using expression (50) for H, with the energy shared equally between X and Y motions (equipartition theorem for quadratic degrees of freedom) implies that each coordinate has a mean square value,

 〈X2〉=kT/kx, 〈Y2〉=kT/ky. (55)

For the systems we study, with b < a and kx<ky, this implies a wider range of motion for the X coordinate, as could obviously be expected. These relations for the mean square values indicate the importance of the force constants for describing the statistical distribution of vortex position.

Now consider determining the probability distributions for the vortex core location. The Hamiltonian is circularly symmetric when expressed in terms of the square of the effective circular coordinate ρ. We can suppose that each possible location has a probability determined from a Boltzmann factor, eβH, where β=(kT)1. Employing the circular symmetry for this coordinate, the probability p(ρ)dρ of finding the vortex core within some range dρ centered at radius ρ is proportional to the area 2πρdρ in a ring of radius ρ, and the Boltzmann factor eβH:

 p(ρ)dρ∝2πρdρe−βH=2πρdρe−12βk¯ρ2. (56)

By including a normalization constant, the unit normalized probability distribution function is easily found to be

 p(ρ)=βk¯ρ e−12k¯ρ2. (57)

The root-mean-square radius ρrms=2kT/k¯ implied from relation (54) can be verified with this probability function. One can also get the mean radius and the most probable radius:

 〈ρ〉=12πkT/k¯,  ρmax=kT/k¯. (58)

For the simulations shown in Figures 9 and 10, with a=60 nm, b=30 nm, T=300 K, position data out to τ=2.5×105 was used to find histograms of vortex core position, and thereby get the radial probability distribution to compare with Eq. (57). The results are shown in Figure 11. To compare with theory, the force constants from spin alignment relaxations were used (see the Figure 11 caption). Note also that as the gyrotropic frequency is considerably larger for L=10. nm, those data correspond to many more orbits of the vortex, equivalent to a more complete averaging. Even so, the distributions for both thicknesses follow very closely to the expected form that depends on the validity of the Thiele equation, with no adjustable parameters (see Figure 12).

### Figure 11.

The radial distribution of vortex core positions for the simulations in Figures 9 and 10, with a=60 nm, b=30, and thicknesses L=5.0,10 nm. Data out to final time τ=2.5×105 was used. The solid curves are the theory expression (57), using force constants k¯=0.1753k0 for L=5.0 nm and k¯=0.6832k0 for L=10 nm, as obtained from spin alignment calculations, with force constant unit k0=A/λex.

Using H expressed in terms of both X and Y, the probability to find the vortex core within some range dX and dY of the location (X,Y) is p(X,Y)dXdYeβHdXdY, where the normalized probability function is found to be

 p(X,Y)=βkx2πe−12βkxX2βky2πe−12βkyY2. (59)

This is a product of Gaussian distributions in each coordinate, p(X,Y)=p(X)p(Y), with zero mean values, but variances given by

 σx=kT/kx,  σy=kT/ky. (60)

The distributions p(X) and p(Y) found from the simulation data of Fig. 9 are shown in Fig. 12, and compare closely to the theoretical expression (59).

### Figure 12.

Distributions of vortex core coordinates for the LLG-Langevin simulation in Figure 9 with a=60 nm, b=30, and thickness L=5.0 nm. The solid curves are the theory expressions from Eq. (59) based on the Thiele theory for vortex motion, using force constants kx=0.1156k0 and ky=0.2657k0 from spin alignment relaxation, where k0=A/λex.

Clearly one could also find the corresponding distributions of the momentum components by similar reasoning.

Instead of looking at the momentum components, we can equivalently calculate a theoretical speed distribution for the vortex core [13]. This is simplest if we use the effective circular coordinate ρ, and consider that fact that its velocity in Eq. (36) implies a speed u|ρ˙| given by

 u=|ωG|ρ. (61)

As u is proportional to ρ, so are their probability distributions. If g(u) is the desired speed probability distribution, then conservation of probability states that

 g(u)du=p(ρ)dρ=p(u/|ωG|)du/|ωG|. (62)

Thus, the speed distribution is derived from the effective circular coordinate distribution by

 g(u)=|ωG|−1 p(u/|ωG|). (63)

With |ωG|=k¯/G, one obtains

 g(u)=βG2k¯u e−12βG2u2/k¯=2uurms2 e−u2/urms2. (64)

This depends on the root-mean-square speed, determined from ρrms,

 urms=|ωG|ρrms=2kTk¯/G. (65)

The function g(u) is a Maxwellian speed distribution similar to that for an ideal gas. One could consider the factor in the exponent as depending on a kinetic energy term 12mGu2 for a particle, where mG is some mass associated with that particle in gyrotropic motion. From Eq. (64), one can read off the value needed for this mass,

 mG=G2/k¯. (66)

This curious result gives a kind of effective mass that depends on the potential experienced by the vortex. Thus, it should not be consider an intrinsic vortex mass. Generally, G is linearly proportional to thickness L [see expression (25)], whereas k¯ tends to increase approximately with L2 [see expression (41) and also Section 4.2], making this mass nearly independent of L. Probably, mG is more strongly determined by the disk area, πab. In the case of circular disks, using the approximate expression (41) for k¯=kF and the definition (25) of G gives a quantitative result,

 mG≈(2π)2a0.878μ0γ2. (67)

Thus, the mass is determined primarily by the disk radius a, and it does not depend on the material parameters such as the exchange stiffness or saturation magnetization. For a radius a=100 nm, the mass is 1.2×1022 kg, an extremely small value. Even so, the mass can be taken to represent how a vortex responds dynamically to the potential. With the gyrotropic frequency given by |ωG|=k¯/G, the mass is written equivalently as

 mG=G/|ωG|∝L/|ωG|. (68)

With mG depending only on disk radius or possibly area in the xy plane, and the gyrovector proportional to L, this re-expresses that |ωG| is also proportional to L, as shown implicitly in Figures 5 and 7.

## 6. Summary and interpretation of results

This chapter has provided an overview of some methods for finding the static, dynamic, and statistical properties of vortex excitations in thin nanodisks of soft magnetic material. By assuming the thickness is much less than the principal radius, L ≪ a, the magnetization points primarily within the plane of the disk, except within the vortex core, and it has only weak dependence on the coordinate z perpendicular to the plane. This allowed for the transformation to an equivalent 2D problem, which has been studied here using a form of micromagnetics, converting the continuum problem to one on a square grid.

The Lagrange-constrained spin alignment scheme was used to find static vortex energies while securing the vortex in a desired location R, thereby allowing for the calculation of vortex potential U(R) within the disk. For a vortex near the center of an elliptic disk, the force constants kx and ky for displacements along the principal axes a, b are found, with kxky when ba. However, the disk ellipticity b/a must be above a lower limiting value for a vortex to stable; a very narrow disk will prefer the formation of a quasi-single-domain state, or even a multi-domain state, but not a vortex. A vortex energetically prefers a displacement along the longer axis of the disk; that is consistent with the shape of its elliptic orbits, which have the same ellipticity as the disk itself [see Eq. (33)].

The vortex gyrotropic orbits can be described very well through the use of the Thiele equation (24), which replaces the dynamics of the many degrees of freedom in the magnetization field M(r,t) by the dynamics of only two Cartesian components in the vortex core location, R=(X(t),Y(t)). This works best for a vortex near the disk center, where it is unlikely to be destabilized by deformations caused by the boundaries. For zero temperature, the dynamics from RK4 integration of the LLG equations is completely consistent with that from the Thiele equation. The Thiele equation predicts the vortex gyrotropic frequencies to be ωG=k¯/G, which is confirmed in the dynamics simulations while using force constants from the Lagrange-constrained static vortex structures. Generally, the zero-temperature gyrotropic frequencies are roughly proportional to L/a with only a weak dependence on disk ellipticity, as can be concluded from the results in Figure 7. The frequencies are determined by the geometric mean force constant, k¯=kxky, which shows why knowledge of the vortex potential is important for this problem.

Thermal effects for nonzero temperature have been included by introducing a Langevin equation (21) that results from including stochastic magnetic fields into the LLG equation. This Langevin equation gives the time evolution in the presence of thermal fluctuations. Solved numerically using a second-order Heun algorithm, a vortex initially at the disk center (the minimum energy point) will spontaneously undergo gyrotropic orbital motion, on top of a noisy spin wave background. The orbital motion takes place at a slightly lower frequency compared with its motion for T=0, because the presence of spin waves weakens the exchange stiffness of the system. The resulting distribution of vortex position can be predicted using an effective Lagrangian and Hamiltonian that result from the Thiele equation. That Hamiltonian can be expressed in a form in Eq. (50) containing only a potential energy. This then shows that the distributions (and variances) of effective radial coordinate ρ and Cartesian coordinates X and Y depend on kT/kF, where kF is either k¯ or kx or ky, respectively [see Eqs. (57) and (60)]. Surprisingly, large vortex rms displacements on the order of 1–10 nm can result, with the larger values taking place in the weaker potentials of thinner disks (Figure 11) and in disks with larger radii a. However, these noisy elliptical motions simply reflect the equipartition of energy into the two collective degrees of freedom available to the vortex (X,Y), with each receiving an average thermal energy of 12kT. The radial coordinate, in contrast, receives a full kT of energy on average. The theoretical probability distributions are confirmed in simulations provided the time evolution averages over a large number of gyrotropic orbits.

A vortex speed distribution can also be derived from the position distribution, essentially because the momentum and position coordinates of a vortex are not independent. The speed distribution g(u) can be characterized by a mass mG proportional to the disk radius a, but independent of material properties. The mass has the sense that as the vortex position fluctuates, it has some Maxwellian speed distribution, with a kinetic energy 12mGu2 that enters in the Boltzmann factor. This is in contrast to the Thiele equation, which has been used here with no intrinsic mass term. Indeed, the vortex gyrotropic frequency is the same as that for a corresponding 2D harmonic oscillator of mass mG and spring constant k¯, that is, ωG=k¯/mG

## Acknowledgements

Portions of this work benefited substantially from discussions with Afranio Pereira and Winder Moura-Melo at the Universidade Federal de Viçosa, Minas Gerais, Brazil, and Wagner Figueiredo at the Universidade Federal de Santa Catarina, Florianópolis, Brazil, and from use of computation facilities at both universities.

## References

1 - N.A. Usov and S.E. Peschany, Magnetization curling in a fine cylindrical particle. J. Magn. Magn. Mater. 118, 290 (1993).
2 - J. Raabe, R. Pulwey, S. Sattler, T. Schweinbock, J. Zweck and D. Weiss, Magnetization pattern of ferromagnetic nanodisks. J. Appl. Phys. 88, 4437 (2000).
3 - G. Gioia and R.D. James, Micromagnetics of very thin films. Proc. R. Soc. London Ser. A 453, 213 (1997).
4 - K.S. Buchanan, P.E. Roy, F.Y. Fradkin, K. Yu Guslienko, M. Grimsditch, S.D. Bader and V. Novosad, Vortex dynamics in patterned ferromagnetic ellipses. J. Appl. Phys. 99, 08C707 (2006).
5 - D. Ruzmetov and V. Chandrasekhar, The Dynamics of magnetic vortex states in a single permalloy nanoparticle. J. Magn. Magn. Mater. 320, 47 (2008).
6 - G.M. Wysin, Vortex-in-nanodot potentials in thin circular magnetic dots. J. Phys.: Condens. Matter 22, 376002 (2010).
7 - G.M. Wysin, Vortex dynamics in thin elliptic ferromagnetic nanodisks. Low Temp. Phys. (Fizika Nizkih Temperatur), 41(10), 1009–1023 (2015).
8 - K.S. Buchanan, P.E. Roy, M. Grimsditch, F.Y. Fradkin, K.Yu Guslienko, S.D. Bader and V. Novosad, Magnetic-field tunability of the vortex translational mode in micron-sized permalloy ellipses: Experiment and micro-magnetic modeling. Phys. Rev. B 74, 064404 (2006).
9 - K.Yu Guslienko, B.A. Ivanov, V. Novosad, Y. Otani, H. Shima and K. Fukamichi, Eigenfrequencies of vortex state excitations in magnetic submicron-size disks. J. Appl. Phys. 91, 8037 (2002).
10 - K.Yu Guslienko, X.F. Han, D.J. Keavney, R. Divan and S.D. Bader, Magnetic vortex core dynamics in cylindrical ferromagnetic dots. Phys. Rev. Lett. 96, 067205 (2006).
11 - G.M. Wysin, Magnetic Excitations & Geometric Confinement: Theory and Simulations, ISBN: 978-0-7503-1074-1 (Institute of Physics Expanding Physics ebook, 2015).
12 - R.P. Cowburn, D.K. Koltsov, A.O. Adeyeye, M.E. Welland and D.M. Tricker, Single-domain circular nanomagnets. Phys. Rev. Lett. 83, 1042 (1999).
13 - G.M. Wysin and W. Figueiredo, Thermal vortex dynamics in thin circular ferromagnetic nanodisks. Phys. Rev. B 86, 104421 (2012).
14 - Z. Huang, High accuracy numerical method of thin-film problems in micromagnetics. J. Comp. Math. 21(1), 33 (2003).
15 - J. Sasaki and F. Matsubara, Circular phase of a two-dimensional ferromagnet with dipolar interactions. J. Phys. Soc. Jpn. 66, 2138 (1997).
16 - C.J. García-Cervera, “Magnetic Domains and Magnetic Domain Walls,” Ph.D. thesis, New York University, New York (1999).
17 - C.J. García-Cervera, Z. Gimbutas and E. Weinan, Accurate numerical methods for micromagnetics simulations with general geometries. J. Comp. Phys. 184, 37 (2003).
18 - G.M. Wysin, Magnetic vortex mass in two-dimensional easy-plane magnets. Phys. Rev. B 54, 15156 (1996).
19 - L.D. Landau and E.M. Lifshitz, On the theory of the dispersion of magnetic permeability in ferromagnetic bodies. Phys. Z. Sowjet. 8, 153 (1935).
20 - F.H. de Leeuw, R. van den Doel and U. Enz, Dynamic properties of magnetic domain walls and magnetic bubbles. Rep. Prog. Phys. 43, 44 (1980).
21 - J.L. García-Palacios and F.J. Lázaro, Langevin-dynamics study of the dynamical properties of small magnetic particles. Phys. Rev. B 58, 14937 (1998).
22 - A.A. Thiele, Steady-state motion of magnetic domains. Phys. Rev. Lett. 30, 230 (1973).
23 - U. Nowak, in Annual Reviews of Computational Physics IX, p. 105, edited by D. Stauffer (World Scientific, Singapore, 2000).
24 - B.A. Ivanov and C.E. Zaspel, Gyrotropic mode frequency of vortex-state permalloy disks. J. Appl. Phys. 95, 7444 (2004).
25 - T.S. Machado, T.G. Rappoport and L.C. Sampaio, Vortex core magnetization dynamics induced by thermal excitation. Appl. Phys. Lett. 100, 112404 (2012).
26 - B.A. Ivanov, E.G. Galkina and A.Yu Galkin, Quantum dynamics of vortices in small magnetic particles. Low Temp. Phys. 36, 747 (2010).