Nonlinear magnetic vortex dynamics in a circular nanodot excited by spin-polarized current

We investigate analytically and numerically nonlinear vortex spin torque oscillator dynamics in a circular magnetic nanodot induced by a spin-polarized current perpendicular to the dot plane. We use a generalized nonlinear Thiele equation including spin-torque term by Slonczewski for describing the nanosize vortex core transient and steady orbit motions and analyze nonlinear contributions to all forces in this equation. Blue shift of the nano-oscillator frequency increasing the current is explained by a combination of the exchange, magnetostatic, and Zeeman energy contributions to the frequency nonlinear coefficient. Applicability and limitations of the standard nonlinear nano-oscillator model are discussed.


Background
Spin torque microwave nano-oscillators (STNO) are intensively studied nowadays. STNO is a nanosize device consisting of several layers of ferromagnetic materials separated by nonmagnetic layers. A dc current passes through one ferromagnetic layer (reference layer) and thus being polarized. Then, it enters to an active magnetic layer (so-called free layer) and interacts with the magnetization causing its high-frequency oscillations due to the spin angular momentum transfer. These oscillation frequencies can be tuned by changing the applied dc current and external magnetic field [1][2][3] that makes STNO being promising candidates for spin transfer magnetic random access memory and frequency-tunable nanoscale microwave generators with extremely narrow linewidth [4]. The magnetization in the free layer can form a vortex configuration that possesses a periodical circular motion driven by spin transfer torque [1,[5][6][7][8][9][10][11]. For practical applications of such nanoscale devices, some challenges have to be overcome, e.g., enhancing the STNO output power. So, from a fundamental point of view as well as for practical applications, the physics of STNO magnetization dynamics has to be well understood.
In the present paper, we focus on the magnetic vortex dynamics in a thin circular nanodot representing a free layer of nanopillar (see inset of Figure 1). Circular nanodots made of soft magnetic material have a vortex state of magnetization as the ground state for certain dot radii R and thickness L. The vortex state is characterized by in-plane curling magnetization and a nanosize vortex core with out-of-plane magnetization. Since the vortex state of magnetization was discovered as the ground state of patterned magnetic dots, the dynamics of vortices have attracted considerable attention. Being displaced from its equilibrium position in the dot center, the vortex core reveals sub-GHz frequency oscillations with a narrow linewidth [2,7,12]. The oscillations of the vortex core are governed by a competition of the gyroforce, Gilbert damping force, spin transfer torque, and restoring force. The restoring force is determined by the vortex confinement in a nanodot. Vortex core oscillations with small amplitude can be well described in the linear regime, but for increasing of the STNO output power, a large-amplitude motion has to be excited. In the regime of large-amplitude spin transfer-induced vortex gyration, it is important to take into account nonlinear contributions to all the forces acting on the moving vortex. The analytical description and micromagnetic simulations of the magnetic field and spin transfer-induced vortex dynamics in the nonlinear regime have been proposed by several groups [12][13][14][15][16][17][18][19][20][21][22], but the results are still contradictory. It is unclear to what extent a standard nonlinear oscillator model [13] is applicable to the vortex STNO, how to calculate the nonlinear parameters, and how the parameters depend on the nanodot sizes.
In this paper, we show that a generalized Thiele approach [23] is adequate to describe the magnetic vortex motion in the nonlinear regime and calculate the nanosize vortex core transient and steady orbit dynamics in circular nanodots excited by spin-polarized current via spin angular momentum transfer effect.

Analytical method
We apply the Landau-Lifshitz-Gilbert (LLG) equation of motion of the free layer magnetization M s is the saturation magnetization, γ > 0 is the gyromagnetic ratio, H eff is the effective field, and α G is the Gilbert damping. We use a spin angular momentum transfer torque in the form suggested by Slonczewski [24], τ s = σJm × (m × P), where σ = ℏη/(2|e|LM s ), η is the current spin polarization (η ≅ 0.2 for FeNi), e is the electron charge, P is direction of the reference layer magnetization, and J is the dc current density. The current is flowing perpendicularly to the layers of nanopillar and we assume P ¼ Pẑ. The free layer (dot) radius is R and thickness is L.
We apply Thiele's approach [23] for the magnetic vortex motion in circular nanodot (inset of Figure 1). We assume that magnetization distribution can be characterized by a position of its center X(t) that can vary with time and, therefore, the magnetization as a function of the coordinates r and X(t) can be written as m(r,t) = m(r,X(t)). Then, we can rewrite the LLG equation as a generalized Thiele equation for X(t): where W is the total magnetic energy, α,β = x,y, and ∂ α = ∂/∂X α . The components of the gyrotensorĜ, damping tensorD, and the spin-torque force can be expressed as follows [16]: We assume that the dot is thin enough and m does not depend on z-coordinate. The magnetization m(x,y) has the components m x þ im y ¼ 2w= 1 þ w w ð Þ and Þexpressed via a complex function w ζ; ζ À Á [25]. Inside the vortex core, the vortex configuration is described as a topological soliton, w ζ; ζ is an analytic function. Outside the vortex core region, the magnetization distribution is w ζ; ζ For describing the vortex dynamics, we use two-vortex ansatz (TVA, no side surface charges induced in the course of motion) with function f(ζ) being written as f ζ where C is the vortex chirality, ζ = (x + iy)/R, s = s x + is y , s = X/R, c = R c /R, and R c is the vortex core radius. The total micromagnetic energy Equation 1 including volume W v m and surface W s m magnetostatic energy, exchange W ex energy, and Zeeman W Z energy of the nanodot with a displaced magnetic vortex is a functional of magnetization distribution W[m(r, t)]. Using m = m(r, X(t)) and integrating over-the-dot volume and surface, the energy W can be expressed as a function of X within TVA [16]. The Zeeman We introduce a time-dependent vortex orbit radius and phase by s = u exp(iΦ). The gyroforce in Equation 1 is determined by the gyrovector G ¼ Gẑ, where G = G z = G xy . The functions G(s) and W(s) depend only on u = |s| due to a circular symmetry of the dot. G(0) = 2πpM s L/γ, where p is the vortex core polarity. The damping forceD _ X and spin-torque force F ST are functions not only on u = |s| but also on direction of s. Nonlinear Equation 1 can be written for the circular dot in oscillator-like form where It is assumed here that F ST (s) = a(u)(z × s) [14], where a is proportional to the CPP current density J and a (0) = πRLM s σJ.
To solve Equation 3, we need to answer the following questions: (1) can we decompose the functions W(s), G (s), D αβ (s), and F ST (s) in the power series of u = |s| and keep only several low-power terms? and (2) what is the accuracy of such truncated series accounting that u = |s| can reach values of 0.5 to 0.6 for a typical vortex STNO? Some of these functions may be nonanalytical functions of u = |s|. If the answer to the first question is yes, then we should decompose W(s) up to u 4 , F ST (s) up to u 3 , and G(s), D αβ (s) up to u 2 -terms to get a cubical equation of the vortex motion. The series decomposition of G(s) does not contain u 2 -term; it contains only small c 2 u 2 -term, , although G(u) essentially decreases at large u, when the vortex core is close to be expelled from the dot [16]. The result of power decomposition of the total energy density w u and the coefficients are JCRς M s and There is an additional contribution to κ/2, 2(L e /R) 2 , due to the face magnetic charges essential for the nanodots with small R [27]. The contribution is positive and can be accounted by calculating dependence of the equilibrium vortex core radius (c) on the vortex displacement. This dependence with high accuracy at cu < < 1 can be described by the function Here, c(0) is the equilibrium vortex core radius at s = 0, for instance c(0) = 0.12 (R c = 12 nm) for the nanodot thickness L = 7 nm.
The last term in Equation 3 prevents its reducing to a nonlinear oscillator equation similar to the one used for the description of saturated STNO in [13]. Calculation within TVA yields the decomposition d n s ð Þ ¼ d 0 n þ d 1 n s x s y , where d 0 n ¼ 0 , i.e., the term containing d n (s) ≈ α G u 2 <<1 can be neglected. Then, substituting s = u exp(iΦ) to Equation 3, we get the system of coupled equations Equation 3 and the system (6) are different from the system of equations of the nonlinear oscillator approach [13]. Equations 6 are reduced to the autonomous oscilla- Þ only if the conditions d 2 < < 1 and dχ < < ω G are satisfied and we define the positive/negative damping parameters [13] as Γ + (u) = d(u)ω G (u) and Γ − (u) = χ(u). We note that reducing the Thiele equation (1) to a nonlinear oscillator equation [13] is possible only for axially symmetric nanodot, when the functions W(s), G(s), d(s) and χ(s) depend only on u = |s| and the additional conditions d n < < 1, d 2 < < 1, and dχ < < ω G are satisfied. The nonlinear oscillator model [13] cannot be applied for other nanodot (free layer) shapes, i.e., elliptical, square, etc., whereas the generalized Thiele equation (1) has no such restrictions.

Numerical method
We have simulated the vortex motion in a single permalloy (Fe 20 Ni 80 alloy, Py) circular nanodot under the influence of a spin-polarized dc current flowing through it. Micromagnetic simulations of the spin-torqueinduced magnetization dynamics in this system were carried out with the micromagnetic simulation package MicroMagus (General Numerics Research Lab, Jena, Germany) [28]. This package solves numerically the LLG equation of the magnetization motion using the optimized version of the adaptive (i.e., with the time step control) Runge-Kutta method. Thermal fluctuations have been neglected in our modeling, so that the simulated dynamics corresponds to T = 0. Material parameters for Py are as follows: exchange stiffness constant A = 10 −6 erg/cm, saturation magnetization M s = 800 G, and the damping constant used in the LLG equation α G = 0.01. Permalloy dot with the radius R = 100 nm and thickness L = 5, 7, and 10 nm was discretized in-plane into 100 × 100 cells. No additional discretization was performed in the direction perpendicular to the dot plane, so that the discretization cell size was 2 × 2 × L nm 3 . In order to obtain the vortex core with a desired polarity (spin polarization direction of dc current and vortex core polarity should have opposite directions in order to ensure the steady-state vortex precession) and to displace the vortex core from its equilibrium position in the nanodot center, we have initially applied a short magnetic field pulse with the out-of-plane projection of 200 Oe, the in-plane projection H x = 10 Oe, and the duration Δt = 3 ns. Simulations were carried out for the physical time t = 200 to 3,000 ns depending on the applied dc current because for currents close to the threshold current J c1 , the time for establishing the vortex steady-state precession regime was much larger than for higher currents (see Equation 8 below).

Results and discussion
Calculated analytically, the vortex core steady orbit radius in circular dot u 0 (J) as a function of current J is compared with the simulations (see Figure 1). There is no fitting except only taking the critical current J c1 value from simulations. Agreement is quite good, confirming that all the nonlinear parameters of the vortex motion were accounted correctly. The steady orbit radius u 0 (J) allows finding the STNO generation frequency ω G J ð Þ ¼ ω 0 J ð Þ þ ω 1 u 2 0 J ð Þ, which increases approximately linearly with J increasing up to the second critical current value J c2 when the steady oscillation state becomes unstable (see Figure 2). The instability is related with the vortex core polarity reversal reaching a core critical velocity or the vortex core expelling from the dot increasing the current density J [12,16]. We simulated the values of J c2 = 2.7, 5.0, and 10.2 MA/cm 2 for the dot thickness L = 5, 7, and 10 nm, respectively. The calculated STNO frequency is 15 to 20% higher than the simulated one due to overestimation of ω ′ 0 within TVA for β =0.1. The calculated nonlinear frequency part is very close to the simulated one, except the vicinity of J c2 , where the analytical model fails.
Our comparison of the calculated dependences u 0 (J) and ω G (J) with simulations is principally different from the comparison conducted in a paper [19], where the authors compared Equations 5 and 7 with their simulations fitting the model-dependent nonlinear coefficients N and λ from the same simulations. One can compare Figures 1  and 2 with the results by Grimaldi et al. [20], where the authors had no success in explaining their experimental dependences u 0 (J) and ω G (J) by a reasonable model. The realistic theoretical nonlinear frequency parameter N for Py dots with L = 5 nm and R = 250 nm should be larger than 0.11 that the authors of [21] used. N = 0.25 can be calculated from pure magnetostatic energy in the limit β → 0 (inset of Figure 2). Accounting all the energy contributions in Equation 4 yields N = 0.36, which is closer to the fitted experimental value N = 0.50.
The system (6) can be solved analytically in nonlinear case. Its solution describing transient vortex dynamics is where u(0) is the initial vortex core displacement and 1=τ þ J ð Þ ¼ 2d 0 ω ′ 0 J=J c1 −1 ð Þ is the inverse relaxation time for J > J c1 (order of 100 ns). u t; J c1 ð Þ∝1= ffiffi t p at t → ∞ and J = J c1 . If J < J c1 , the orbit radius u(t, J) decreases exponentially to 0 with the relaxation time The divergence of the relaxation times τ ± at J = J c1 allows considering a breaking symmetry secondorder phase transition from the equilibrium value u 0 = 0 to finite u 0 J ð Þ∝ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi J=J c1 −1 p defined by Equation 7. Equations 7 and 8 represent mean-field approximation to the problem and are valid not too very close to the value of J = J c1 , where thermal fluctuations are important [13,21]. Equation 8 describes approaching of the vortex orbit radius to a steady value u 0 (J) under influence of dc spinpolarized current. We note that the corresponding relaxation time is determined by only linear parameters, whereas the orbit radius (7) depends on ratio of the nonlinear and linear model parameters. The solution of Equation 8 is plotted in Figure 3 as a function of time along with micromagnetic simulations for circular Py dot with thickness L = 7 nm and radius R = 100 nm. The vortex was excited by in-plane field pulse during approximately the first 5 ns, and then the vortex core approached the stationary orbit of radius u 0 (J). We estimated u(0) after the pulse as u(0) = 0.1 and plotted the solid lines without any  fitting except using the simulated value of the critical current J c1 . Overall agreement of the calculations by Equation 8 and simulations is quite good, especially for large times t ≥ 3τ + , although the calculated relaxation time τ + is smaller than the simulated one due to overestimation of ω ′ 0 within TVA. The typical simulated ratio J c2 /J c1 ≈ 1.5; therefore, minimal τ + ≈ 20 to 30 ns. But the transient time of saturation of u(t, J) is about of 100 ns and can reach several microseconds at J/J c1 < 1.1. The simulated value of λ = 0.83, whereas the analytic theory based on TVA yields the close value of λ(J c1 ) = 0.81.
Typical experiments on the vortex excitations in nanopillars are conducted at room temperature T = 300 K without initial field pulse, i.e., a thermal level u(0) should be sufficient to start vortex core motion to a steady orbit. To find the thermal amplitude of u(0), we use the well-known relation between static susceptibility of the systemχ and magnetization fluctuations M 2 The in-plane components are , and M = ξM s s, where ξ = 2/3 within TVA [26]. This leads to the simple relation u 2 for interpretation of the experiments. u T (0) ≈ 0.05 (5 nm in absolute units) for the dot made of permalloy with L = 7 nm and R = 100 nm.
The nonlinear frequency coefficient N(β, R, J) = κ′(β, R, J)/ κ(β, R, J) is positive (because of κ, κ′ >0 for typical dot parameters), and it is a strong function of the dot geometrical sizes L and R and a weak function of J. For the dot radii R > > L e , N(β, R, 0) ≈ 0.21 − 0.25 (the magnetostatic limit, see inset of Figure 2). If R > > L e and β → 0, then N (β, R, 0) ≈ 0.25 [14]. For the realistic sizes of free layer in a nanopillar (R is about 100 nm and L = 3 to 10 nm), this coefficient is essentially larger due to finite β and exchange contribution, and it can be of order of 1. The exchange nonlinear contribution κ′ ex is important for R < 300 nm. However, the authors of [19][20][21] did not consider it at all. Note that N(0.089, 300 nm, 0) ≈ 0.5 recently measured [29] is two times larger than 0.25. The authors of [19] suggested to use an additional term~u 6 in the magnetic energy fitting the nonlinear frequency due to accounting a u 4 -contribution (N = 0.26) that is too small based on [14], while the nonlinear coefficient N(β, R) calculated by Equation 5 for the parameters of Py dots (L = 4.8 nm, R = 275 nm) [19] is equal to 0.38. Moreover, the authors of [19] did not account that, for a high value of the vortex amplitude u = 0.6 to 0.7, the contribution of nonlinear gyrovector G(u) ∝ c 2 u 2 to the vortex frequency is more important than the u 6 -magnetic energy term. The gyrovector G(u) decreases essentially for such a large u resulting in the nonlinear frequency increase. The TVA calculations based on Equation 5 lead to the small nonlinear Oe energy contribution κ′ Oe , whereas Dussaux et al. [19] stated that κ′ Oe is more important than the magnetostatic nonlinear contribution.

Conclusions
We demonstrated that the generalized Thiele equation of motion (1) with the nonlinear coefficients (2) considered beyond the rigid vortex approximation can be successfully used for quantitative description of the nonlinear vortex STNO dynamics excited by spin-polarized current in a circular nanodot. We calculated the nonlinear parameters governing the vortex core large-amplitude oscillations and showed that the analytical two-vortex model can predict the parameters, which are in good agreement with the ones simulated numerically. The Thiele approach and the energy dissipation approach [12,19] are equivalent because they are grounded on the same LLG equation of magnetization motion. The limits of applicability of the nonlinear oscillator approach developed for saturated nanodots [13] to vortex STNO dynamics are established. The calculated and simulated dependences of the vortex core orbit radius u(t) and phase Φ(t) can be used as a starting point to consider the transient dynamics of synchronization of two coupled vortex ST nano-oscillators in laterally located circular nanopillars [30] or square nanodots with circular nanocontacts [31] calculated recently.