Magnetic vortex dynamics in the non-circular potential of a thin elliptic ferromagnetic nanodisk with applied fields

Abstract: Spontaneous vortex motion in thin ferromagnetic nanodisks of elliptical shape is dominated by a natural gyrotropic orbital part, whose resonance frequency ωG = k/G depends on a force constant and gyrovector charge, both of which change with the disk size and shape and applied in-plane or out-of-plane fields. The system is analyzed via a dynamic Thiele equation and also using numerical simulations of the Landau-Lifshitz-Gilbert (LLG) equations for thin systems, including temperature via stochastic fields in a Langevin equation for the spin dynamics. A vortex is found to move in an elliptical potential with two principal axis force constants kx and ky, whose ratio determines the eccentricity of the vortex motion, and whose geometric mean k = √ kxky determines the frequency. The force constants can be estimated from the energy of quasi-static vortex configurations or from an analysis of the gyrotropic orbits. kx and ky get modified either by an applied field perpendicular to the plane or by an in-plane applied field that changes the vortex equilibrium location. Notably, an out-ofplane field also changes the vortex gyrovector G, which directly influences ωG. The vortex position and velocity distributions in thermal equilibrium are found to be Boltzmann distributions in appropriate coordinates, characterized by the force constants.


Introduction
Magnetic vortices in thin ferromagnetic disks of sub-micron size offer an interesting system for the study of collective dynamics of fundamental excitations [1].A single vortex centered in a circular disk can be the absolute minimum energy state or it can be metastable, separated from a nearby quasi-singledomain state by a weak energy barrier.A vortex experiences a restoring force F = −k F R dependent on its displacement R = (X, Y) from the center, mostly caused by demagnetization effects from weak pole formation on the disk edges, where k F is a force constant [2].The response to this force is the gyrotropic orbital vortex motion at a frequency ω G , which can be detected in resonance experiments [3].
While much work has been developed for disks with a circular boundary, less studied are disks with an elliptical boundary [4].Quite generally in physical problems, deviation of a circular system into one with elliptical symmetry leads to interesting modifications, due to the breaking of the circular symmetry.We consider an elliptical edge characterized by semi-major radius a along the x-axis and a semi-minor radius b along the y-axis, for a disk of thickness L with L a.In a circular disk, the in-plane angle of local magnetization in a vortex state is determined by vorticity charge q = ±1, and a chirality or circulation charge c = ±1, via a relation The gyrotropic charge is G = 2πqp, or gyrovector G = Gẑ, where p = ±1 is the polarization direction (magnetization along ±z) of the vortex core.Simulations show that the vortex structure itself is not squeezed along the narrow direction of the ellipse.Rather, the vortex retains close to a circular shape, but experiences a modified potential.When taken as an assumption, this is the rigid vortex approximation.In numerical simulations, it need not hold precisely.Regardless of that, the deviation of the disk edge from circular symmetry is found to introduce two non-equivalent force constants k x and k y , corresponding to the principal axes of the ellipse.The force constants change with the shape of the disk [5], until reaching a high in-plane aspect ratio b/a 1, where the vortex becomes unstable and a nearly uniform quasi-single-domain state or possibly a two-domain state emerges, similar to those found [6] in circular nanodisks.Here only the vortex state is considered.
With a magnetic field applied in the plane of the disk, the vortex equilibrium position will be displaced away from the disk center, perpendicular to the field in a direction depending on chirality c.In an elliptic disk, displacements along the two principal directions are non-equivalent.Buchanan et al. [7] have noted that a field H ext y applied along the shorter (y) axis, that shifts the vortex minimum position along the long (x) axis, results in an increase in its gyrotropic frequency.To the contrary, a field H ext x along the long axis, shifting the vortex minimum position along the short axis, does not significantly change the frequency.We confirm these results, showing how the vortex effective potential and force constants are modified by the displaced vortex equilibrium location.The gyrotropic resonant frequency is then seen to shift, without any modification of the gyrotropic charge or gyrovector G.If a field is applied instead perpendicular to the disk plane (z-axis), there will be two non-degenerate vortex gyrotropic modes, as has been seen in resonance experiments [8] and micromagnetics [9] for circular disks.Here further analysis of this effect is given, and we find how a perpendicular applied field modifies both the force constants and gyrotropic frequencies.The gyrovector is shifted with an out-of-plane applied field H ext z , such that the resonant frequency changes nearly linearly with H ext z .A field pointing out-of-plane in the same direction as the vortex core magnetization increases ω G .
In this article vortex motion in elliptic disks is considered, as obtained from two-dimensional micromagnetics simulations, and from analysis of the Thiele equation [10] for magnetization dynamics of a collective excitation such as a domain wall or vortex.The Thiele equation analysis depends directly on the force or the effective potential that the vortex moves in.This analysis is considered first for the zero temperature motion as obtained from Landau-Lifshitz-Gilbert (LLG) equations.The studies verify that the Thiele equation gives a good description of the motion and can predict the gyrotropic frequencies, based on the force constants, even in the presence of applied fields.
Quasi-static vortex structures and their energies can be used to estimate the force constants.The dynamic motion itself can also be used to estimate the force constants, especially for vortices displaced by an in-plane applied field.In the first set of studies presented here, some of the behaviors of the force constants and the gyrovector with disk geometry and applied field are discussed.Note that the gyrovector only changes significantly for an out-of-plane applied field.In a second set of studies, the stochastic effects due to finite temperature are included, by using the Langevin-LLG equations for the micromagnetics.The simulations can be compared with the vortex statistics expected from applying the principle of equipartition, not to the numerous spin degrees of freedom, but rather, to only the two degrees of freedom for the position of the vortex core.That is, good agreement is found for the vortex position statistics, based on a theory with only two degrees of freedom, when analyzing simulation data for the 2N degrees of freedom represented in the dynamics of N micromagnetics cells for an elliptical nanodisk of magnetic material.

The System, Energetics and Dynamic Equations
The magnetic medium is assumed be of thickness L along the z-axis and have an elliptical boundary in the xy-plane, The vortex magnetization structure is not strongly modified by the boundary, however, the vortex experiences a non-circular effective potential U(R) caused by the boundary.R = (X, Y) represents the vortex core location, measured from the disk center.For slight deviations from its equilibrium position, the potential experienced by a vortex is found to be of elliptical form, The in-plane aspect ratio or ellipticity b/a ≤ 1 controls the properties of the potential in which the vortex moves, which is represented in terms of the force constants k x and k y .When the ratio b/a becomes too small, a vortex will be destabilized, and some other vortex-free state such as a quasisingle-domain state will be prefered.
The underlying dynamics is that of the local magnetization M(r) = M s m(r).Analyzed numerically in the micromagnetics approximation, the magnitude is fixed at the saturation value M s and only the direction m(r) is changing.A continuum energy function for the system includes isotropic exchange, and demagnetization (H M ) and applied (H ext ) fields: The exchange is characterized by the exchange stiffness A ex in units of J/m.Its competition with magnetostatic energy due to demagnetization effects leads to the exchange length, that sets a length scale for the problem: The magnetization tends to preserve its direction over this length scale.Two-dimensional (2D) micromagnetics is based on the idea of using a 2D grid with cells not larger than this scale, so that the spatial variations in m can be correctly described.For Permalloy with A ex ≈ 13 pJ/m and M s ≈ 860 kA/m, and µ 0 M s ≈ 1.08 T, that gives λ ex ≈ 5.3 nm.In these micromagnetics simulations we have used cells of size a cell × a cell × L, with a cell = 2.0 nm, so that weak changes in magnetization direction will be included.This 2D analysis is based on the assumption that there is little dependence of m on the z-coordinate, through the thickness of the disk.That should be true for thin disks.The numerical simulations keep track of cell magnetic dipoles µ i = µm i where i labels the cells, and µ ≡ La 2 cell M s is their fixed magnitude.Neighboring cells have an effective exchange constant, J = 2A ex L. The demagnetization fields H M i are produced as a result of the current state of the m i .Their calculation is based on magnetostatics theory for an isolated thin 2D system, using effective Green's functions appropriate for the thin disk problem [11].The calculations of H M i can be accelerated somewhat through the use of fast Fourier transforms applied to the defining convolution integrals.
For the discretized 2D system, the dynamic equations of motion resulting from (4), and including an additional damping term with dimensionless parameter α, are the Landau-Lifshitz-Gilbert (LLG) torque equations, This includes the gyromagnetic ratio γ and the effective local magnetic induction B i acting on a cell, Of course, there are contributions to B i due to exchange fields, demagnetization fields, and the externally applied field.For numerics, we measure B i in a basic unit B 0 ≡ J/µ = 2A ex /(a 2 cell M s ), defining dimensionless fields as b i = B i /B 0 .With cell size a cell = 2.0 nm, this is B 0 ≈ 7.59 T for the simulations described here using Permalloy parameters.Then the time is measured in units t 0 ≡ (γB 0 ) −1 , leading to dimensionless time variable τ = t/t 0 .Thus the dynamics follows the dimensionless equations, For the numerical simulations, the magnetization unit vectors m i (t) are evolved forward by some updating procedure.The method used depends on whether static or dynamic results are desired.Static or quasi-static vortex structures were used to get force constants, as an example.Dynamic simulation is necessary to obtain the gyrotropic frequencies.

Quasi-Static Vortex Properties
For finding quasi-static or relaxed structures, a local spin alignment relaxation scheme has been used, iteratively pointing each magnetic moment to align with its effective field, until some convergence is reached.The vortex energy (same as total system energy) was evaluated for different positions, which were enforced by a Lagrange constraint [12] on the vortex core, explained briefly here.The constraint is included by adding additional terms to the energy functional, creating the modified energy functional on the 2D grid with N cells, Here the set of α i at N grid sites are Lagrange multipliers used to enforce ultimately a fixed length m = 1 for the reduced magnetization vectors.The vector λ = (λ x , λ y ) with only xy components is a pair of Lagrange multipliers for enforcing the position constraint.λ is a fictitious field applied only to a limited set of grid sites N c = 24 symmetrically surrounding the vortex core.One can note that the in-plane magnetization structure of a vortex is highly symmetrical around its core, with the in-plane components of m pointing in opposite directions on opposing sides of the vortex.This is the reason for including the last term in the expression for Λ.The requirements that Λ be an extremum with respect to variations in the set of α i and λ x , λ y leads to a set of equations that can be solved iteratively for the cell magnetizations m i and the constraining field.This works for a desired core position exactly between four nearest grid cells.A slight modification to this scheme allows for securing the position offset arbitrarily from these symmetrical positions, see Ref. [12] for further details.These increase quickly with disk thickness, due to stronger pole density at the disk edge.
For elliptical nanodisks without applied fields, Figure 1 shows results for the vortex force constants obtained by this scheme, for vortices near the center of elliptical disks.The semi-major radius is a = 120 nm, while the semi-minor radius b takes on a range of values, corresponding to ellipses of different shapes.The force constants were estimated by using the energy change for a displacement of ∆X = 4.0 nm or ∆Y = 4.0 nm away from the disk center, where the vortex energy is the minimum value, U 0 .Assuming the potential in Eq. ( 3), the force constants are obtained quasi-statically by expressions, The results in Figure 1 show some interesting features.First, the potential is stiffer for vortex motion along the shorter (y) direction.Thus, k y ≥ k x , where the equality holds only in the circular limit.The vortex moves much more freely along the longer (x) axis.Secondly, for ellipses with a higher in-plane aspect ratio (i.e., smaller b/a), k x reduces slightly while k y increases more rapidly.At the same time, the geometric mean force constant k remains nearly constant.Eventually all of the force constants tend towards zero for small enough b/a, where the vortex is destabilized.Finally, also note that the force constants increase with the thickness of the disk, more than linearly with L. The disk with greater thickness have a much stronger demagnetization effect, which leads to a much stronger restoring force on the vortex.

About Finding the Vortex Location
The vortex core position R = (X, Y) can be determined with a precision much smaller than the numerical grid.This is done by first locating the set of four cells that surround the vorticity center or vortex core, among which the change in in-plane angle φ changes by 2π as expected from Eq. ( 1).Then, using a set of the cells within about 4 exchange lengths from that preliminary position estimate r v , an improved estimate is found from an average weighted by the squares of out-of-plane scaled magnetization components m z i .This uses the fact that the magnetization tilts out of the disk plane at the vortex center, with m z i decaying away towards its boundary value over a distance on the order of the exchange length.We use an expression to estimate the position, Each r i is the center position of a cell, with the sum restricted to the core region.Especially for zero-temperature simulations this weighted location gives a very smoothly changing vortex position, even when following the dynamics.It is verified by comparison with the time-dependent plots of the magnetization as it evolves in the simulations.It is used below for the comparison of simulations with the Thiele theory for vortex core motion, and also in the study of vortex position statistics in Section 7.

Thiele Equation Analysis
The results found for vortex dynamics based on numerics can be analyzed in light of the Thiele collective coordinate equation for a localized magnetic excitation.If the effective force F = − ∇U(R) is acting on the vortex core, then the Thiele equation for the core velocity V = Ṙ predicts the motion by This depends on the topological charge or gyrovector G of the vortex, which is determined by the vorticity charge q = +1 (antivorticity with q = −1 is not considered here), the out-of-plane core polarization p = ±1, and the magnetic dipole moment per unit area, m 0 = LM s , The potential in (3) is assumed, which depends on force constants k x and k y .With the gyrovector having only an out-of-plane component, the equations of motion are equivalent to those of an elliptical oscillator, e.g., Starting at location (X 0 , Y 0 ) at time t = 0, the solution is that of elliptical motion, The geometric mean of the force constants, k, determines the gyrotropic frequency ω G , which can also be taken as a vector perpendicular to the plane, ω G = ω G ẑ.The minus sign is included in (18) to indicate clockwise motion in the xy-plane when the gyrovector has a positive z-component.While the mean force constant determines ω G , the ratio of those force constants controls the shape of the orbit.Considering using Y 0 = 0, one gets the ratio of maximum displacements on the two axes (orbital ellipticity e, or ratio of semi-minor to semi-major axes) to be Thus, the magnetic dynamics leads to vortex elliptical motion, whose ellipticity is directly related to the square root of the force constant ratio.Simulations of static vortex structure that lead to k x and k y , such as in Figure 1, show that to a good approximation, e ≈ b/a for adequately large nanodisks.
For some analysis, a stretching of the coordinate system into a new variable is useful, because it returns the potential to a circular symmetry: For the same reason, the vortex core motion then takes a simple form, The variable, ρ and especially its magnitude is convenient for analysis of vortex position statistics.

Vortex Gyrotropic Frequencies
For zero temperature dynamics, fourth order Runge-Kutta (RK4) scheme has been used to get the time evolution.For finite temperature dynamics, additional stochastic fields are included into the equations (8), and the resulting Langevin-LLG equations (also known as the stochastic LLG equations) can be evolved forward in time using a second order Heun method.The chosen temperature T determines the relative strength of the stochastic magnetic fields.See Refs.[2,5] for further details.
At zero temperature, the validity of the Thiele analysis is confirmed by simulating vortices starting with a small displacement (4.0 nm) from the disk center, and evolving the undamped LLG equations forward in time with an RK4 scheme.The motion can be followed over 10 to 20 periods, from which the frequency is measured.The frequencies obtained dynamically are found to be directly proportional to the mean force constants k obtained from statics.Results versus disk ellipticity are summarized in a compact form in Figure 2. The frequencies follow closely the prediction (18) of the Thiele equation, which for q = 1 can be transformed to a form: The RHS contains λ ex as length unit, a frequency unit ω 0 ≡ µ 0 4π γM s , and the force constant unit k 0 ≡ A ex /λ ex .Thus, the Thiele equation predicts in these units: This is confirmed in the simulations for various disk geometries to within a few percent, except for small ellipticity, for which there is limited vortex stability, see Figure 2.
a=120 nm, L=10 nm a=120 nm, L=5 nm a=60 nm, L=10 nm a=60 nm, L=5 nm a=30 nm, L=10 nm a=30 nm, L=5 nm From simulations of LLG eqns.by RK4 time evolution Figure 2. Summary of vortex gyrotropic frequencies in elliptical nanodisks, using p = −1 and frequency unit ω 0 ≡ µ 0 4π γM s , force constant unit k 0 ≡ A ex /λ ex and λ ex as the unit of length.These data fall very close to the prediction of the Thiele equation, which is the unit value, Lω G / k = λ ex ω 0 /k 0 .

Vortex in an Out-of-plane Applied Field
Next, consider an applied field H ext z , or in dimensionless simulation units, include a nonzero field b ext z = µ 0 H ext z /B 0 .Circular magnetic disks with this field orientation have been considered by Loubens et al. [8] and more recently by Fried et al. [9], where it was found that the two opposite vortex polarizations become energy-and frequency-split by the field.Here we consider a vortex with q = p = +1, and positive (negative) values of b ext z correspond to the applied field pointing in the same (opposite) direction as the magnetization in the vortex core region.The effect on the vortex effective potential for a system with b/a = 0.5 is shown in Figure 3, for applied fields b ext z = 0, ±0.05.The total system energy is plotted as a function of the vortex displacement from the disk center.Generally, the total energy is reduced with an applied field, and due to the symmetry, the disk center remains the location of the minimum.A positive field causes the larger reduction in total energy, as more of the magnetization is strongly aligned to the applied field.Another example for the same system, but with b ext z = ±0.15, is shown in Figure 4. Obviously, an even greater field causes a larger downward energy shift.More importantly, the force constants are also modified by b ext z , although this effect is difficult to see in the plots of U(R).Using Eq. ( 10), the results for k x , k y and k are shown in Figure 5.This clearly shows how all of these are maximized at zero field, and tend towards k x = k y = k → 0 at an upper positive field limit, where the vortex is destabilized.Similarly, the vortex will be destabilized by a strong enough negative field, however, this takes place partly because the vortex core region in that case will acquire a very short radius (the core will be oppositely polarized to its surroundings).Note that the gyrotropic frequency ω G = k/G would be diminished by positive or negative b ext z , if the gyrovector were constant.However, that is not the case, and the gyrotropic frequencies versus b ext z does not have the shape of the k curves.
. Vortex force constants versus the dimensionless out-of-plane applied field, in an elliptical system of indicated parameters, obtained by applying Eq. 10 to potentials like those in Figs. 3 and 4. The force constants go to zero near b ext z =≈ 0.7 (a very strong field, equivalent to B ext z ≈ 5.3 T in Permalloy) as the magnetization becomes uniformly out-of-plane and the vortex is destabilized.
Without applied fields, the dimensionless gyrovector γG/m 0 = 2πpq represents the total steradians of a unit sphere covered by the magnetization direction of the vortex, which is half of the unit sphere.G is also given by a formula involving the out-of-plane reduced magnetization at the core, m z (0) = p and at infinite radius, where the scale is determined by the zero-field continuum gyrovector value, Once a field is applied along the z-axis, the boundary value m z (∞) will be modified, which directly leads to a modification of G. Considering the case p = +1, G is reduced (increased) for positive (negative) b ext z , compared to its value at zero field.For a large enough system, the vortex core region is small compared to the rest of the area.Then an approximate expression for the value at large radius is to use the average over the whole system, m z (∞) ≈ m z .This gives a rough estimate of the gyrovector, A plot showing the behaviors of m z and the resulting G for an elliptical system with a = 120 nm, b = 60 nm, is displayed in Figure 6.There results close to a linear dependence of G on the out-of-plane field.
. From a set of zero-temperature LLG simulations, the gyrotropic periods τ G (in simulation units t 0 on right axis) and frequencies ω G (in units of ω 0 = µ 0 4π γM s , left axis) for an elliptical disk as a function of the out-of-plane field.Solid symbols connected by dotted lines are derived from the numerical simulations.Open symbols with solid lines result from the theory expression (18) using the estimates of force constants and G as in Figure 6.Note that the field dependence of ω G is nearly linear at small field values, but exhibits a strong nonlinearity at the upper destabilizing field.
Using LLG simulations of vortex motions, starting from small displacements X = 4.0 nm, many orbits can be followed, from which the periods t G can be determined, then leading to the gyrotropic frequencies ω G = 2π/t G .Note that in dimensionless units the period is τ G = t G /t 0 , where t 0 = (γB 0 ) −1 is the simulation time unit.For the example system with a = 120 nm, b = 60 nm, L = 20 nm, the resulting periods and frequencies are shown in Figure 7 over a range of out-of-plane applied field.There is very good agreement between the frequencies from the simulations, and those derived from the Thiele theory expression (18), using the gyrovector G that varies with field in Figure 6.One finds a rather strong effect of the field on ω G .The frequency exhibits a nearly linear dependence on b ext z at weaker fields, and even for negative fields, where there is little destabilization of the vortex.At the higher positive field values, the frequency deviates from the linear behavior, as the vortex becomes unstable.The weak-field linear behavior is consistent with that found by Loubens et al. [8], while the nonlinear dependence at higher field strength may be explained by core deformation, according to Fried et al. [9].

Effect of an In-plane Applied Field
A magnetic field applied within the plane of the disk will displace the equilibrium vortex position away from the disk center, along a line perpendicular to that field; the direction depends on the vortex chirality c or twist direction, Eq. ( 1).Doing simulations of the quasi-static vortex potentials, results such as those in Fig. 8 are found.A field applied along the y direction, for positive chirality c, displaces the vortex in the −x direction, all the more so for stronger field.Similar potential curves were obtained by Buchanan et al. [7] using what they called a "static approach" in micromagnetics.A field applied on the whole system was used to move the vortex around to different equilibrium positions.Then the system energy for the vortex relaxed through micromagnetics at different equilibrium positions gives the potential curves.
The Lagrange-constrained approach used here appears to give equivalent potential curves, and has the advantage that this iterative scheme is faster than running a micromagnetics relaxation for every desired point on the potential curve.However, one could argue that the scheme applied by Buchanan et al. is conceptually simpler and numerically more straightforward.
The asterisks in Figure 8 indicate the minima of the different potential curves.The resulting equilibrium vortex locations, as functions of the applied field, are those shown in Figure 9.The vortex appears to become unstable when reaching some edge region of the disk, which is about the same distance for the two disk thicknesses tested here. .From the vortex potentials such as those in Fig. 8, the equilibrium vortex locations as functions of the in-plane applied field.The vortex becomes unstable if pushed too close to the edge, hence the curves end at specific field strengths.A larger maximum field is needed to eject the vortex in the thicker disk.
Buchanan et al. [7] have noticed further that in addition to this displacement, there is an upward shift in the gyrotropic frequency, resulting in close to a 100% increase, as long as the vortex has been shifted on the longer axis of the ellipse.Here we give further analysis to this effect; we find that the shifted location changes the force constants, but k still primarily determines ω G , because for the most part G is unaffected by an in-plane field.
Initially a vortex is relaxed with damping, in the presence of an in-plane field, b ext y along the shorter ellipse axis.This produces some equilibrium location (X eq , 0) away from the center, on the longer axis, and a corresponding minimum energy E eq .Then, another simulation is done without damping, starting from a nearby position, which results in gyrotropic motion around location (X eq , 0), at some higher energy E, whose period t G is measured.The resulting orbital shape (X(τ), Y(τ)) has a semi-major axis A path and semi-minor axis B path .Fitting the energy difference E − E eq to expression (3) for U(R), gives a dynamic evaluation of the force constants by Some typical results for k x , k y , and the resulting k are indicated in Figure 10, for different disk thicknesses.At least at weaker field magnitude, the force constants increase with b ext y .This is the primary cause of an increasing gyrotropic frequency.It appears that there is a greater tendency for k x to increase rather than k y .As the vortex is pushed into the narrower end of the disk, it experiences stronger demagnetization fields and a resulting stiffer potential. .From the vortex potentials such as those in Figure 8, the vortex force constants obtained dynamically using expressions (28), as functions of the in-plane applied field.The vortex becomes unstable if pushed too close to the disk edge; for example, see the lowest curve in Figure 8.Hence, these curves drop off and then terminate at the limiting field that expels the vortex from the disk. .Vortex gyrotropic frequencies (for p = −1) under the presence of an inplane magnetic field that displaces the equilibrium position.ω G found directly from the periods in simulations (filled symbols with dotted lines) is compared with the Thiele prediction, Eq. (18), using k evaluated from the orbital shape via Eq.( 28) and assuming fixed G = G 0 (open symbols with solid lines).Beyond the ends of the curves, the vortex is ejected from the disk, see also Figure 9.
The gyrotropic frequency ω G = 2π/t G then is compared with the Thiele prediction, ω G = − k/G, using their geometric mean value k, and assuming the zero-field gyrovector value, G = G 0 .Results are shown in Figure 11, as functions of the dimensionless applied field b ext y , up to a limit where the vortex is forced out the edge of the disk.The increase of ω G with applied field is seen to be significant, consistent with the results of Buchanan et al. [7].The largest frequency change takes place as the vortex is forced to move close to the disk edge.A weakening in the effect takes place near the limiting value of b ext y , which is most simply explained by the reduction of k when the vortex is near the disk edge.

Thermally Induced Spontaneous Motion
It was pointed out by Machado et al. [13] and studied further by Wysin and Figueiredo [2] that gyrotropic motion can be spontaneously generated at finite temperature.The motion can self-generate even for a vortex initially at the disk center (no applied field is considered here).The amplitude of the motion is determined by equipartition, which can be analyzed under the assumption that the vortex core obeys the Thiele equation.This means that the vortex is considered to possess only two primary degrees of freedom, being the two Cartesian coordinates (X, Y) of its core position.A exercise shows that a Lagrangian that leads back to the Thiele equation with elliptic potential is [5]   Using the resulting canonical momentum P, this implies a Hamiltonian that can be expressed as purely potential energy, E = P • Ṙ − L = U(X, Y) = 1 2 k ρ 2 .Then assuming equipartition for only two degrees of freedom, each with average energy of 1  2 k B T , with k B is Boltzmann's constant, the thermally averaged Hamiltonian takes the value, E = k B T ≡ β −1 .This gives a prediction for the distribution of vortex core radial position using the effective circular coordinate ρ defined in (20) as The most probable effective radius is given by ρ max = (β k) −1/2 , which demonstrates how weaker restoring force leads to a wider distribution.Simulations can be used to test these expectations, solving the Langevin-LLG equations by a second order Heun method [11].The integration was done out to time τ = 2.5×10 5 , with a weak damping constant α = 0.02, starting with a vortex at the disk center.The vortex motion initiates spontaneously due to thermal fluctuations, then proceeds in a noisy gyrotropic orbit for many periods, from which a histogram of ρ can be calculated.Some typical results without an applied field are shown in Figure 12 for disks with a = 60 nm, b = 48 nm at T = 300 K using Permalloy parameters.There is reasonable agreement between the simulation and the Thiele theory, however, averaging over at least 100 orbits is required.The distribution of vortex velocity can also be found to follow a Boltzmann form.The application of an out-of-plane field can be seen to modify the thermally generated vortex radial distribution.For instance, for an elliptical Permalloy disk with a = 60 nm, b = 48 nm, L = 10 nm, the zero-field mean force constant is k = 1.632 × 10 −3 N/m, as used in Figure 12.With the field b ext z = 0.30 applied, quasi-static vortex relaxation shows that the force constant is reduced to k = 1.304 × 10 −3 N/m.This is similar to the force constant reductions as displayed in Figure 5.A slight change in p(ρ) results, as shown in Figure 13, where the zero-field and non-zero field cases are plotted.The slightly weaker potential in the presence of the field allows the vortex to explore a larger area of the disk, for fixed temperature.If a much stronger field were applied, it could be possible to weaken the potential sufficiently for the vortex to destabilize or exit the disk.

Figure 1 .
Figure1.Vortex force constants versus disk ellipticity, from static vortices obtained by spin alignment relaxation.The quantity A ex /λ ex serves as the unit of k x , k y and k ≡ k x k y .These increase quickly with disk thickness, due to stronger pole density at the disk edge.

0Figure 3 .
Figure3.Vortex effective potential curves for an elliptical nanodisk with b/a = 0.5, at the indicated applied fields.The vortex configurations were found quasi-statically, using a Lagrange constraint on the vortex position.The energy unit is J = 2A ex L. This type of result is used to make estimates of the force constants, which are clearly different for x and y displacements.

Figure 4 .
Figure 4.More vortex effective potential curves for an elliptical nanodisk with b/a = 0.5, at the indicated applied fields, like those in Figure 3.

2 GFigure 6 .
Figure 6.For a moderately sized elliptical disk, the magnitude of the gyrovector versus out-of-plane applied field, as estimated from the average out-of-plane reduced magnetization via Eq.(27).Note that G → 0 at the same field value (b ext z ≈ 0.7) as where k → 0, compare Figure5.This occurs where m z → 1 and the vortex out-of-plane component becomes indistinguishable from the uniform background.

Figure 8 .
Figure 8. Example of the effects of an in-plane applied field along the shorter side of the ellipse, on a vortex with chirality c = +1, see Eq. (1).The different curves are labeled by values of b ext y .The vortex potentials are shifted such that the minima (shown as asterisks) move perpendicular to the field, and the effective force constants are modified.

Figure 9
Figure 9. From the vortex potentials such as those in Fig.8, the equilibrium vortex locations as functions of the in-plane applied field.The vortex becomes unstable if pushed too close to the edge, hence the curves end at specific field strengths.A larger maximum field is needed to eject the vortex in the thicker disk.

Figure 10
Figure 10.From the vortex potentials such as those in Figure8, the vortex force constants obtained dynamically using expressions (28), as functions of the in-plane applied field.The vortex becomes unstable if pushed too close to the disk edge; for example, see the lowest curve in Figure8.Hence, these curves drop off and then terminate at the limiting field that expels the vortex from the disk.
Figure 11.Vortex gyrotropic frequencies (for p = −1) under the presence of an inplane magnetic field that displaces the equilibrium position.ω G found directly from the periods in simulations (filled symbols with dotted lines) is compared with the Thiele prediction, Eq. (18), using k evaluated from the orbital shape via Eq.(28) and assuming fixed G = G 0 (open symbols with solid lines).Beyond the ends of the curves, the vortex is ejected from the disk, see also Figure9.

Figure 13 .
Figure 13.Vortex core radial position distributions for Permalloy parameters in a disk at T = 300 K. Symbols are from Langevin LLG simulations to time τ = 2.5×10 5 .Curves are theory expression Eq. (30) with k = 1.632 × 10 −3 N/m for b ext z = 0 and k = 1.304 × 10 −3 N/m for b ext z = 0.30, from quasi-static relaxed vortex calculations.The distribution is shifted outward slightly due to the reduction in force constant caused by the field, see also Figure 5. )