Statistics of Gyrotropic Vortex Dynamics in Submicron Magnetic Disks

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.

local magnetization (), which is the magnetic dipole moment per unit volume, at position .A material is considered with saturation magnetization , 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 () as a function of position may tend to do two things: (1) () will have a strong tendency to point within the plane of the disk [3], if possible; (2) () 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 () 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 () 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 (), 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.= 30 nm, thickness = 10 nm, from the spin alignment relaxation scheme for the micromagnetics model, Section 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 () 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 changes sign is obvious as a circle in Figure 1.

Vortex charges
The magnetization profile () 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 = ± 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 () also has a vorticity charge = ± 1, which corresponds to the direction of rotation of () as one moves along a closed path encircling the core.The value = + 1 holds for the vortices described here, which are controlled by demagnetization effects ( being forced to follow the boundary).The value = − 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, (0) = ± = , where = ± 1 is the polarization charge.Because there is an energy barrier to flip the core polarization from = + 1 to = − 1, it offers yet another charge that could be useful for data storage and manipulation.

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 points perpendicular to the disk) defines the location of the vortex core, which we denote by position vector = (,), measured along the , Cartesian axes within the disk plane.Because the system is assumed to be thin, only two coordinates , are needed to locate the core.Further, the magnetization itself has little dependence on the coordinate () perpendicular to the disk.We take = (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 () takes place, while the demagnetization effects at the disk edge still try to maintain () parallel to the edge.The net result of the displacement is a slight increase in total energy.The vortex, as a quasiparticle, lives in some effective potential (), which has something approximating a parabolic form [6], with the minimum at the disk center, where is a force constant.This further implies an effective force back 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 elliptical shape [7], defined by principal axes and < : 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 , along 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.

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 , but a spatially varying direction, by writing () = (), where () is a unit vector.From the energetics expressed in terms of (), 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.

Energetics of a continuum nanomagnet
The system is governed by ferromagnetic exchange energy and interactions of generally with the demagnetization field (self-generated by ) and any possible externally applied field e .A continuum Hamiltonian for the system is where 0 is the permeability of free space, and 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 = 860 kA/m [12].The magnetization changes its direction over a length scale e called the exchange length.Exchange energy of the order / e 2 competes with demagnetization energy of the order For Py, e ≈ 5.3 nm.Exchange forces dominate over lengths less than e , but demagnetization dominates over larger lengths, allowing the () field to change its direction.At a boundary, the exchange effects are much less present, and demagnetization helps () to point parallel to the boundary, if possible.

The demagnetization field 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 = 0 ( + ).Gauss' law, ∇ .= 0, then becomes By assuming the demagnetization field comes from a scalar potential via = − ∇ Φ , a Poisson equation for the magnetostatics is obtained: Therefore, the magnetization () produces an effective magnetic charge density , which is the source in this Poisson equation.The solution for scalar potential Φ 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, ≪ , where is the radius (or the semi-major axis for an elliptical disk).In this case, both and do not depend on the vertical coordinate along the disk axis.Then, the problem can be solved by effective two-dimensional (2D) Green's functions [14].The components α = ,, of the demagnetization field can be expressed as 2D convolution integrals, where ( − ′) represents a tensor Green's function, and the integration is over 2D positions, for example now = (,).The evaluation of these integrals can be accelerated through the use of fast Fourier transforms [15].

Discretization and micromagnetics for simulations
For numerical solutions of the magnetization field () = (), it is necessary to partition the system into cells labeled by index for positions .We use a square grid of cells of individual size cell × cell × L, with cell = 2.0 nm, and disk thickness L = 5 nm and L = 10 nm.At the center of each cell is a unit direction vector , whose motion is to be followed.Each cell contains a magnetic dipole moment of magnitude μ = L cell 2 M s and direction .This micromagnetics approach [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 = 2, and magnetic fields have been brought to dimensionless forms, The presence of the factor cell 2 / e 2 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 cell should then be less than the exchange length.The micromagnetics approach, with the assumption that only the direction of is changing, makes the implicit assumption that demagnetization effects are a perturbation on exchange effects.Obviously, the Green's functions ( − ′) must also be brought to a discrete form to carry out the calculation of .
The Hamiltonian can be used to define the net magnetic inductions that act on each cell's magnetic dipole , according to where the dimensionless magnetic inductions are ( ) The first term involves a sum over the nearest neighbors () of cell ; 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 cell = 2.0 nm, and Py parameters, one has 0 2 ≈ 1.08 T and 0 ≈ 7.59 T. 0 gives the order of magnitude of the exchange fields; the demagnetization fields are considerably weaker.

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: (16) Note that this holds because / = is the spin angular momentum of the cell, whose time derivative is the torque, × .The equation can be written in terms of the dimensionless quantities, also defining a dimensionless time , The unit of time used for simulations is 0 ≡ ( 0 ) −1 .For Py parameters, it takes the value 0 ≈ 0.75 ps.Its reciprocal also defines a simulation frequency unit, 0 ≡ 0 ≈ 1.336 THz.
For static configurations, however, the time derivatives in Eq. ( 17) are zero.This implies that each dipole or its unit vector must align with the local field in that cell, .
Thus, an algorithm that iteratively points each along its current value of 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 ( , ), according to θ φ θ φ θ ˆ= (cos cos , cos sin , sin ).
In this notation, is referred to as the in-plane angle and is the out-of-plane angle, which can be positive or negative.The approximate in-plane structure of a vortex located at position = (,) in a disk can be expressed using the vorticity = + 1 as where (x i , y i ) is the 2D location of micromagnetics cell , and 0 = 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 = 0, that is, a planar vortex.However, the iteration will be such that all will remain zero, unless some small nonzero deviation is included.Therefore, small random values of 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 become insignificant (less than about one part in 10 8 ).
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 e (region where 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 reverses.Vortex structures for an elliptical nanodisk with = 60 nm, = 30 nm, = 10 nm, using cell size cell = 2.0 nm in the spin alignment relaxation scheme, including a Lagrange constraint on position. (a) Vortex is held at = 16 nm, = 0, resulting in total energy = 15.244 after 2000 iterations.(b) Vortex is held at = 0, = 16 nm, resulting in total energy = 16.001 after 3500 iterations.Compare Figure 1, where = = 0 and the energy is lower.

Effective potentials of a vortex in a nanodisk
Spin alignment relaxation can also be used to estimate the effective potential () for the vortex, by including a constraint on the vortex position = (, ).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 e along ± will displace the vortex along ±, 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 and .
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 .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 quasisingle-domain state (nearly uniform ()), or some other multi-domain state without a vortex.This is especially likely in the case of elliptic disks with a high aspect ratio ( ≪ ), where demagnetization will strongly favor 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 ( ≫ ), again, demagnetization will cause 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 = 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 for circular disks or even and for elliptical disks can be estimated quite accurately.In the example of Figure 3, one can observe that 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 and with < , we have found that for adequately large nanodisks and of sufficient size to stabilize the vortex, the ratio of force constants asymptotically approaches the relation, This has the correct limit for a circular disk, = .The relation is summarized by saying that the geometric ellipticity, b/, directly determines an energetic ellipticity, / .The energetic ellipticity can be seen to determine the shape of the elliptical vortex orbits at constant energy in the phase space.

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 , = 0].Here, we take that one step further and also include stochastic magnetic fields , () 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 are a superposition of deterministic effects (from ) and stochastic effects (from , ).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 is 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 , , 2 where is Boltzmann's constant and 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 = 0 or > 0.

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 = 0 LLG equations) that the steady-state dynamics of the core velocity = ˙ is described by a Thiele equation [22], The motion is governed by the gyrovector , 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 .In terms of the magnetization per unit area, 0 ≡ , and the gyromagnetic ratio = − 1.76 × 10 11 T −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, = + 1, while there are two core polarizations = ± 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 particlelike 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 = (,), 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 = with the Thiele equation, A vector identity is useful, The vortex velocity points in the plane of the disk, but is perpendicular to that plane, so ⋅ = 0.This gives the velocity as ( ) With = ( ˙, ˙), this is a pair of first-order differential equations, which can be directly integrated, starting from some initial vortex position (0) = ( 0 , 0 ).An elementary calculation gives elliptical motion, with instantaneous coordinates where the gyrotropic frequency is determined by the geometric mean of the force constants, ≡ : The negative square root is used, because a vortex with along + and a centrally directed force will move in the clockwise (or negative) direction in the 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, e ; especially, e 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 .
With e = 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 is found to be ( ) ( ) Dividing through by the constant, , this is the standard equation of an ellipse, with the semimajor and minor axes max , max , given by max max 2 2 , .
Their ratio is then The last approximate result in terms of the disk axes , 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), 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 = .This equivalent circular coordinate is useful for the analysis of vortex position statistics in elliptical disks.

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

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 10 12 over as many as 5.0 × 10 5 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 .

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 , , 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 are triplets of random deviates with zero mean and unit variance, for each site . 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 to get stochastic field components of the correct mean and variance (it does not need to be a Gaussian distribution) [13,7].

Vortex gyrotropic motion at zero temperature
In a circular nanodisk at zero temperature, with a radial force as 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 ) 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 is the geometric mean of the force constants along the principal axes.For a circular system, .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 2 /, that is, The frequency unit 0 = 0 −1 = 0 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, This is 0 = 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, With vorticity = + 1 assumed, the sign of is determined by the core polarization .This expression suggests using 0 ≡ / e as the unit of force constant and e as the unit of length.
Simulations can verify the frequency predictions from the Thiele equation.As shown below in some examples, the dimensionless periods of gyrotropic motion can be estimated precisely, in simulation time units ( = / 0 ).Dimensionless angular frequencies are then 2/ , which are given physical values by multiplying by 0 = 0 .Using Eq. ( 15), these can then be converted into units of 0 as follows: We use this below to convert the raw numerical data ( ) 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 inplane magnetization angle , to locate the four cells nearest to the vorticity center, , 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 ≈ only to a precision equal to the cell size.
It can be greatly improved by making a weighted average of the cell locations , using their squared out-of-plane components, which are largest in the vortex core, as the weighting function: |< 4 e ( ) = .( ) 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 changes smoothly.This algorithm even works very well for vortices moving in response to thermal fluctuations.

Circular nanodisks simulations
Some typical vortex motions in circular nanodisks of radius = 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 = (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, , 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 ≈ 5800, 3270, and 1872, respectively.From statics calculations of the effective potentials as described earlier, the corresponding raw force constants are k f cell / ≈ 0.033863, 0.120192 and 0.419143, respectively, using cell = 2.0 nm.Rescaling by a factor λ ex / cell = 5.3/2.0 converts them into e /, which appears in the Thiele theory expression (43).
For these and other similar vortex motion simulations with ranging from 2.0 to 20 nm, one can compare to the Thiele prediction by plotting the frequency versus / with units as suggested from Eq. (43) (see Figure 5).Note that for a given radius , the disk with the smallest has the largest frequency.The result is that , obtained from dynamics simulations, is very close to linearly related to /, 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.The damping = 0.02 was turned off after dimensionless time = 1000.The periods can be calculated accurately from the undamped motion.Graphs of () are of the same amplitudes and frequencies but shifted a quarter of a period.Gyrotropic frequency magnitudes versus geometric ellipticity /, for elliptic nanodisks of semi-major axis = 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 in 6, as expected from ∝ in Eq. (43).

Elliptical nanodisk simulations
Simulations for elliptic nanodisks [7] offer an even wider range of possibilities, because the variations with geometric ellipticity / 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 = 120 nm and = 10 nm.Both and were determined from the potentials derived by spin alignment with a position constraint.Their geometric mean , which determines gyrotropic frequencies, is also shown.The curves for these force constants do not go below a minimum value of /, where the vortex becomes unstable.
The corresponding gyrotropic frequencies for = 10 nm and also for = 5.0 nm are shown in Figure 7, versus /.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 in

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.

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, = 2 determines the energy scale and depends on the disk thickness.As an example, we consider a disk with = 60 nm, = 30 nm, and thickness = 5.0 nm, at temperature = 300 K.For Py parameters ( = 13 pJ/ m), the energy unit is = 130 zJ, while the thermal energy scale is = 4.14 zJ, which gives the dimensionless temperature, = = 0.032, for = 300 , = 5.0 nm. 2 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 ((), ()) 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 ≈ 3300).The period is somewhat longer than that found at zero temperature, ≈ 2970.
This softening of the mode with temperature is to be expected.In addition, the amplitudes of and motions are not equal, as expected from the elliptical disk shape.The gyrotropic orbital motion continues indefinitely; it was followed out to = 2.5 × 10 5 to get vortex statistics.
For comparison, an identical simulation but with the disk thickness increased by a factor of 2 to = 10 nm in shown in Figure 10, again starting the vortex from the disk center.The greater thickness approximately quadruples , but also doubles the gyrovector, thereby resulting in the frequency being double that for = 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 .

Thermal vortex motion as described by the Thiele equation
Next, we consider the statistical mechanics of the vortex core position = (,), based on an effective Lagrangian and Hamiltonian that give back the Thiele equation.The analysis [7] makes use of the general elliptic potential () 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 = ⋅ − .As P is determined by and , without any time derivatives, one can interpret this as a pair of constraint relations between components of and .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 = ˙= 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 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, .

Thermalized vortex probability distributions from the Thiele equation
One can assume that any of the coordinates, ,, , , as well as effective circular coordinate = ( , ), 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, 2 2 , in the limit of a circular disk.Using expression (50) for , with the energy shared equally between and motions (equipartition theorem for quadratic degrees of freedom) implies that each coordinate has a mean square value, For the systems we study, with b < a and < , this implies a wider range of motion for the 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 − , where = () −1 .Employing the circular symmetry for this coordinate, the probability () of finding the vortex core within some range centered at radius is proportional to the area 2 in a ring of radius , and the Boltzmann factor e − : By including a normalization constant, the unit normalized probability distribution function is easily found to be The root-mean-square radius implied from relation (54) can be verified with this probability function.One can also get the mean radius and the most probable radius: 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 = 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).= 30, and thicknesses = 5.0,10 nm.Data out to final time = 2.5 × 10 5 was used.The solid curves are the theory expression (57), using force constants = 0.1753 0 for = 5.0 nm and = 0.6832 0 for = 10 nm, as obtained from spin alignment calculations, with force constant unit 0 = / e .Using expressed in terms of both and , the probability to find the vortex core within some range and of the location (,) is (,) ∝ e − , where the normalized probability function is found to be This is a product of Gaussian distributions in each coordinate, (,) = ()(), with zero mean values, but variances given by The function () 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 1 2 2 for a particle, where is 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, is linearly proportional to thickness [see expression (25)], whereas tends to increase approximately with 2 [see expression (41) and also Section 4.2], making this mass nearly independent of L. Probably, is more strongly determined by the disk area, .In the case of circular disks, using the approximate expression (41) for = and the definition (25) of gives a quantitative result, Thus, the mass is determined primarily by the disk radius , and it does not depend on the material parameters such as the exchange stiffness or saturation magnetization.For a radius = 100 nm, the mass is 1.2 × 10 −22 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 = /, the mass is written equivalently as With depending only on disk radius or possibly area in the plane, and the gyrovector proportional to , this re-expresses that is also proportional to , as shown implicitly in

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, ≪ , the magnetization points primarily within the plane of the disk, except within the vortex core, and it has only weak dependence on the coordinate 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 , thereby allowing for the calculation of vortex potential () within the disk.For a vortex near the center of an elliptic disk, the force constants and for displacements along the principal axes a, b are found, with ≤ when ≤ .
However, the disk ellipticity / 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 multidomain 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 (,) by the dynamics of only two Cartesian components in the vortex core location, = ((), ()).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 = − /, which is confirmed in the dynamics simulations while using force constants from the Lagrangeconstrained static vortex structures.Generally, the zero-temperature gyrotropic frequencies are roughly proportional to / 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, = , 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 = 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 and depend on / , where is either or or , 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 .However, these noisy elliptical motions simply reflect the equipartition of energy into the two collective degrees of freedom available to the vortex (, ), with each receiving an average thermal energy of .The radial coordinate, in contrast, receives a full 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 () can be characterized by a mass proportional to the disk radius , 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 1 2 2 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 and spring constant , that is, = /

2. 4 .
The cell size is cell = 2 nm.Arrows show only the in-plane projection, ( , ).Blue line (red open) arrows indi- cate positive (negative) values of out-of-plane component .The core, where

Figure 3 .
Figure 3. Numerically determined vortex potentials, in units of the effective cell exchange constant = 2, for circu- lar Py nanodisks of thickness = 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 .
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 () are of the same amplitudes and frequencies but

Figure 5 .
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 .
Figure 6.Effective potential force constants versus geometric ellipticity /, for elliptic nanodisks of semi-major axis = 120 nm, and thickness = 10 nm.Results were found by Lagrange-constrained spin alignment relaxation.

Figure 7 .
Figure 7. Gyrotropic frequency magnitudes versus geometric ellipticity /, for elliptic nanodisks of semi-major axis = 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 in 6, as expected from ∝ in

Figure 6 ,
which is to be expected if the Thiele theory (43) is valid.The additional results for = 5.0 nm are included to demonstrate the dependence on disk thickness.With thicker disks having a greater restoring force and ∝ 2 , due to the extra area on the disk edges, the dependence of ∝ results is gyrotropic frequencies increasing roughly linearly in .The results can be presented in another view in Figure8, showing / versus ellipticity for different .One again gets a clear and quantitative verification of the Thiele theory result (43), seeing that / ∝ e / with the correct constant of proportionality.

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

Figure 9 .
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 = 5.0 nm.The vortex was initiated at the disk center, = = 0.This graph shows only 1/5 of the total data generated and used subsequently for analysis of vortex statis- tics, corresponding to hundreds of gyrotropic periods.

Figure 10 .
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, = 10 nm.Note the considerably smaller amplitude of gyrotropic oscilla- tions, and the much higher frequency.

For the simulations shown in Figures 9 and 10 ,
with = 60 nm, = 30 nm, = 300 K, position data out to = 2.5 × 10 5 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 Figure11.

Figure 11 .
Figure 11.The radial distribution of vortex core positions for the simulations in Figures 9 and 10, with = 60 nm,