Open access peer-reviewed chapter

Statistics of Gyrotropic Vortex Dynamics in Submicron Magnetic Disks

By Gary Matthew Wysin

Submitted: March 21st 2016Reviewed: July 8th 2016Published: March 1st 2017

DOI: 10.5772/64849

Downloaded: 1081


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.


  • 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 rmay 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 fieldM(r)=Msm^(r)of a vortex centered in a nanodisk with principal axesa= 60 nm,b=30nm, thicknessL=10nm, from the spin alignment relaxation scheme for the micromagnetics model, Section 2.4. The cell size isacell= 2 nm. Arrows show only the in-plane projection,(mix,miy). Blue line (red open) arrows indicate positive (negative) values of out-of-plane componentmiz. The core, wheremizis 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 Mzchanges 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 chargec=±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 chargeq=±1, which corresponds to the direction of rotationof M(r)as one moves along a closed path encircling the core. The value q=+1holds for the vortices described here, which are controlled by demagnetization effects (Mbeing 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=±1is the polarization charge. Because there is an energy barrier to flip the core polarization from p=+1to 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 Mpoints perpendicular to the disk) defines the location of the vortex core, which we denote by position vector R=(X,Y), measured along the x,yCartesian axes within the disk plane. Because the system is assumed to be thin, only two coordinates X,Yare 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,


where kFis a force constant. This further implies an effective forceFback towards the disk center, according to the gradient of the potential,


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 ellipticalshape [7], defined by principal axes aand b<a:


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


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


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 Mgenerally with the demagnetization field HM(self-generated by M) and any possible externally applied field Hext. A continuum Hamiltonian for the system is


where μ0is the permeability of free space, and Ais 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=860kA/m [12]. The magnetization changes its direction over a length scale λexcalled the exchange length. Exchange energy of the order A/λex2competes with demagnetization energy of the order 12μ0Ms2. Equating these terms gives the definition of the length scale,


For Py, λex5.3nm. 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 HMin 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


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


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 ΦMcan 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, La, where ais the radius (or the semi-major axis for an elliptical disk). In this case, both Mand HMdo not depend on the vertical coordinate zalong the disk axis. Then, the problem can be solved by effective two-dimensional (2D) Green’s functions [14]. The components α=x,y,zof the demagnetization field can be expressed as 2D convolution integrals,


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 ifor positions ri. We use a square grid of cells of individual size acell× acell× L, with acell=2.0nm, 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 μiof magnitude μ = Lacell2Ms and direction m^i. This micromagneticsapproach [16, 17] then represents the original continuum system, but with a discretized 2D micromagnetics Hamiltonian,


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


The presence of the factor acell2/λex2gives 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 acellshould then be less than the exchange length. The micromagnetics approach, with the assumption that only the direction of Mis 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


where the dimensionless magnetic inductions are


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,


In the results presented here with acell=2.0nm, and Py parameters, one has μ0Ms21.08T and B07.59T. B0gives 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:


Note that this holds because μi/γ=Siis 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 τ,


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

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

Thus, an algorithm that iteratively points each m^ialong its current value of biwill tend to move the system towards a static configuration. We call this approach a spin alignmentrelaxation 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


In this notation, ϕiis referred to as the in-plane angleand θiis the out-of-planeangle, 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=+1as

ϕi=q tan1(yiYxiX)+ϕ0E19

where (xi, yi) is the 2D location of micromagnetics cell i, and ϕ0=cπ2depends 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 θiwill remain zero, unless some small nonzero deviation is included. Therefore, small random values of θican 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^ibecome 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 Mhas 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 Mzreverses.

Figure 2.

Vortex structures for an elliptical nanodisk witha=60nm,b=30nm,L=10nm, using cell sizeacell=2.0nm in the spin alignment relaxation scheme, including a Lagrange constraint on position. (a) Vortex is held atX=16nm,Y=0, resulting in total energyE=15.244Jafter 2000 iterations. (b) Vortex is held atX=0,Y=16nm, resulting in total energyE=16.001Jafter 3500 iterations. CompareFigure 1, whereX=Y=0and 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 Hextalong ±xwill 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 kxand 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 (ba), where demagnetization will strongly favor Maligning 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 Mto 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.0nm. 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 kFfor circular disks or even kxand kyfor elliptical disks can be estimated quite accurately. In the example of Figure 3, one can observe that kFbecomes 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 aand bwith b< a, we have found that for adequately large nanodisks and bof sufficient size to stabilize the vortex, the ratio of force constants asymptotically approaches the relation,


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 constantJ=2AL, for circular Py nanodisks of thicknessL=4.0nm 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,


The changes in m^iare 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 iis suppressed)


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,


where kis Boltzmann’s constant and Tis 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=0or 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=0LLG equations) that the steady-state dynamics of the core velocity V=R˙is described by a Thiele equation[22],


The motion is governed by the gyrovectorG, 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×1011T−1 s−1 of an electron (its magnetic moment divided by its angular momentum), the gyrovector of a vortex is


For the vortices in a disk, q=+1, while there are two core polarizations p=±1possible. 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,


A vector identity is useful,


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


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


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


The negative square root is used, because a vortex with Galong +z^and a centrally directed force will move in the clockwise (or negative) direction in the xyplane. 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, Hextalong 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 Uis found to be


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


Their ratio is then


The last approximate result in terms of the disk axes a,bwas 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),


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


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


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.3nm, together with a micromagnetics cell size of acell=2.0nm.

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 1012over as many as 5.0×105time 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


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


The vectors wiare 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.5to +0.5and then rescaled by 12σsto 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 Fas in Eq. (2), the analysis from the Thiele equation (24) shows that the vortex velocity always points along the azimuthal direction:


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


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,


The frequency unit f0=t01=γB0used 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,


This is ω0=γMsin 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,


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

Simulations can verify the frequency predictions from the Thiele equation. As shown below in some examples, the dimensionless periods τGof 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 ω0as follows:


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

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


For the micromagnetics square grid, the vorticity center falls between the four cells that have a net 2πcirculation in ϕ. This gives the location RRvonly 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:


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 Rchanges 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=120nm 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.02was 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 kfacell/A0.033863, 0.120192and 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 Lranging from 2.0 to 20 nm, one can compare to the Thiele prediction by plotting the frequency ωGversus kF/Lwith units as suggested from Eq. (43) (see Figure 5). Note that for a given radius a, the disk with the smallest Lhas the largest frequency. The result is that ωG, obtained from dynamicssimulations, is very close to linearly related to kF/L, obtained from staticsimulations, 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.02was turned off after dimensionless timeτ=1000. The periods can be calculated accurately from the undamped motion. Graphs ofY(τ)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 fromL= 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 ellipticityb/a, for elliptic nanodisks of semi-major axisa=120nm, and thicknessL=10nm. Results were found by Lagrange-constrained spin alignment relaxation.

Figure 7.

Gyrotropic frequency magnitudes versus geometric ellipticityb/a, for elliptic nanodisks of semi-major axisa=120 nm, and thicknessesL = 5.0 and 10 nm. Results were found by simulations of the LLG equations using RK4 integration. ForL = 10 nm, compare the similar shape of the curve ofk¯inFigure 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/acan 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=120nm and L=10nm. Both kxand kywere 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 ωGfor L=10nm and also for L=5.0nm 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.0nm 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 GLresults 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/Lwith the correct constant of proportionality.

Figure 8.

Gyrotropic frequency magnitudes (from dynamics) scaled by mean force constants (from statics), versus geometric ellipticityb/a, for elliptic nanodisks witha=120nm and two different thicknesses. The results confirm the predictions from the Thiele theory, dashed lines from Eq. (43), using exchange lengthλex=5.3nm, 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=2ALdetermines the energy scale and depends on the disk thickness. As an example, we consider a disk with a=60nm, b=30nm, and thickness L=5.0nm, at temperature T=300K. For Py parameters (A=13pJ/m), the energy unit is J=130zJ, while the thermal energy scale is kT=4.14zJ, which gives the dimensionless temperature,


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.01for 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 Xand Ymotions are not equal, as expected from the elliptical disk shape. The gyrotropic orbital motion continues indefinitely; it was followed out to τ=2.5×105to get vortex statistics.

For comparison, an identical simulation but with the disk thickness increased by a factor of 2 to L=10nm 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.0nm. 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 thicknessL=5.0nm. 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 inFigure 9, but with double the thickness,L=10nm. 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]


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,


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

The Hamiltonian is obtained in the usual way,


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


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,


This is exactly equal to Hin 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,


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


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 Xand Ymotions (equipartition theorem for quadratic degrees of freedom) implies that each coordinate has a mean square value,


For the systems we study, with b< aand kx<ky, this implies a wider range of motion for the Xcoordinate, 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:


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


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:


For the simulations shown in Figures 9 and 10, with a=60nm, b=30nm, T=300K, position data out to τ=2.5×105was 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 inFigures 9and10, witha=60nm,b=30, and thicknessesL=5.0,10nm. Data out to final timeτ=2.5×105was used. The solid curves are the theory expression (57), using force constantsk¯=0.1753k0forL=5.0nm andk¯=0.6832k0forL=10nm, as obtained from spin alignment calculations, with force constant unitk0=A/λex.

Using Hexpressed in terms of both Xand Y, the probability to find the vortex core within some range dXand dYof the location (X,Y)is p(X,Y)dXdYeβHdXdY, where the normalized probability function is found to be


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


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 inFigure 9witha=60nm,b=30, and thicknessL=5.0nm. The solid curves are the theory expressions from Eq. (59) based on the Thiele theory for vortex motion, using force constantskx=0.1156k0andky=0.2657k0from spin alignment relaxation, wherek0=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


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


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


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


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


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 12mGu2for a particle, where mGis some mass associated with that particle in gyrotropic motion. From Eq. (64), one can read off the value needed for this mass,


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, Gis 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, mGis more strongly determined by the disk area, πab. In the case of circular disks, using the approximate expression (41) for k¯=kFand the definition (25) of Ggives a quantitative result,


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=100nm, the mass is 1.2×1022kg, 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


With mGdepending only on disk radius or possibly area in the xyplane, 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 zperpendicular 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 kxand kyfor displacements along the principal axes a, bare found, with kxkywhen ba. However, the disk ellipticity b/amust 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 dynamicssimulations while using force constants from the Lagrange-constrained staticvortex structures. Generally, the zero-temperature gyrotropic frequencies are roughly proportional to L/awith 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 spontaneouslyundergo 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 Xand Ydepend on kT/kF, where kFis either k¯or kxor 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 kTof 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 mGproportional 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 12mGu2that 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 mGand spring constant k¯, that is, ωG=k¯/mG



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.

© 2017 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Gary Matthew Wysin (March 1st 2017). Statistics of Gyrotropic Vortex Dynamics in Submicron Magnetic Disks, Vortex Structures in Fluid Dynamic Problems, Hector Perez-de-Tejada, IntechOpen, DOI: 10.5772/64849. Available from:

chapter statistics

1081total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Two-Dimensional Solitons and Vortices in Linear and Nonlinear Lattice Potentials

By Jianhua Zeng and Boris A. Malomed

Related Book

First chapter

Vortex Structures in Ultra-Cold Atomic Gases

By Nick Verhelst and Jacques Tempere

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us