We present overdamped micromagnetic simulations of the time evolution of magnetic droplet solitons which are formed in a thin film with perpendicular magnetic anisotropy by injecting a spin polarized current through a circular nanocontact. The overdamped dynamics help explore the effective energy landscape of these structures and permit identification of stationary states that are either energy extrema or saddle configurations. Our micromagnetic simulations start with configurations that are numerical solutions of a one-dimensional model where the magnetization depends only on the radial distance to the center of the nanocontact. We verify that these solutions correspond either to energy minima or saddle states, and use them to estimate thermal activation barriers for various current strengths. From the two-dimensional micromagnetic simulations, we identify a new persistent state which does not appear in the one-dimensional model.

Droplet solitons are dynamic states in thin-film materials having a strong magnetic anisotropy perpendicular to the plane of the film.1 Since they were first predicted by Kosevich et al.,2 numerous analytical and theoretical papers have been published on the properties of these objects. Most notably, Hoefer et al.3 introduced a theoretical model to predict that droplet solitons can be sustained in spin-valves with a nano-contact that injects a spin polarized current. Since then, experiments have been designed to detect and observe magnetic droplet solitons4,5 and study many of their properties.

An understanding of the thermal stability of such textures is a critical issue for many practical applications, and a variety of mechanisms have been proposed to estimate soliton lifetimes.6,7 The rotational symmetry of these systems can be used to define an effective energy E by adding to the conservative (i.e., precession only with no damping or spin transfer torque) energy E0 an additional potential EST that accounts for the work needed to rotate the magnetization against the spin transfer torque.8 The rate f of creation of annihilation of droplet solitons will obey an Arrhenius type law, f=f0eΔEkBT, where ΔE is the height of the barrier that separates two (local or global) energy minima, and the prefactor f0 is (roughly speaking) an attempt frequency of the order of GHz.

These energy terms are obtained after integrating corresponding energy densities over the volume of the sample:

E=E0+EST=VE0+ESTd3r.
(1)

The conservative energy density is:

E0(m)=Am2Kμ0Ms22mz2μ0MsHZm
(2)

where A is the exchange constant, K is the magnetocrystalline energy constant, μ0 is the permeability of free space, HZ is the external applied field. The terms on the right hand are, respectively, exchange, effective anisotropy and Zeeman energies. The term in parenthesis can be rewritten as μ0Ms22(Q1), with Q=2Kμ0Ms2. Inside the parenthesis, the first term corresponds to the crystalline anisotropy and the second to the magnetostatic energy of a thin extended film.

A current injected perpendicular to the film after passing by a fixed polarizing layer with magnetization mp will act on m and favor anti-parallel alignments between the layers. To account for this, a spin-torque induced pseudo-potential can be defined as follows:

EST=JedV(r)αmm̂pν=0JedV(r)αln1+νmm̂pνν0.
(3)

where α is the Gilbert damping constant, ν is an anisotropy parameter, J is the value of the current density, e is the electron’s charge, d is the thickness of the film, and V(r) describes the spatial profile of the current. The normalized effective field is obtained from the functional derivative of the energy density Heff=1μ0MsδEδm.

To recover the appropiate magnetization dynamics, the Landau-Lifshitz-Gilbert equation can be rewritten as follows:

ṁ=LHeffγ0μ0Msm×mEST
(4)

where L=γm×+αm×m× is a linear operator that accounts for precession and damping, γ=γ01+α2 and γ0=2.211×105mAs. The second term can be dropped by switching to a reference frame that rotates at a constant frequency. In this way, the critical points of the energy landscape (E=0) appear stationary (ṁ=0) and correspond either to energy minima or to transition states.

These stationary states can be found numerically in a one-dimensional model where profiles are assumed to depend only on the radial distance to the center of the nanocontact, r. In the 1D model, droplet solitons profiles are described by the angle Θ that the magnetization makes with the normal to the film as a function of the rescaled radial distance to the center of the nanocontact, ρlex=rQ1, with lex=2Aμ0Ms2. The magnetization is assumed to precess with a spatially uniform frequency ω and the azimuthal angle Φ is constant across the sample (i.e. ∇Φ = 0). The droplet’s frequency depends on the value of the scaled field hZ=HZẑMs(Q1) and current σ=Jμ0edMs2Q1, and the current has a profile V(ρ) in the form of a Heavyside step function V(ρ)=H(ρρ*) where ρ* is the nanocontact radius. A persistent soliton profile (for which Θ̇=0) can be either an energy minimum or a transition state and satisfies the differential equation:8 

d2Θdρ2+1ρΘρ12sin2Θ+σV(ρ)αhZsinΘ=0,
(5)

with boundary conditions

Θρρ=0=0limΘρ=0.
(6)

The solutions to Eqs. (5-6) are more easily found by imposing a third boundary condition Θ(ρ = 0) = Θmax and solving for ωh=σV(ρ)αhZ as an eigenvalue problem. The amplitude mz(ρ = 0) = cosΘmax can be used to associate solutions found at different nanocontact radii ρ*. In the absence of fields and for infinite nanocontact radii, the frequency of precession of the soliton reduces to ω0=σα. A value of ω0 for solutions of (5-6) at finite radii is assigned by selecting the value ω0 which of the solution that for infinite nano-contact radius has the same amplitude at the origin mz(ρ = 0).

To examine if the 1D model is robust when transformed to two dimensions we do micromagnetic simulations in OOMMF.9 To do this, we have used the OOMMF extension that incorporates spin polarized currents10 and added an energy term without an associated effective field (we were motivated by a contribution from Fong11 who uses the opposite approach: he incorporates additional STT dynamical terms using additional fields but does not incorporate an energy term). As a result, the dynamics are unchanged but we are now able to keep track of the the pseudo-potential in the total energy. These extensions were introduced in Ref. 8 and are provided there as supplementary material with the necessary problem specification files. Here, we use those extensions and repeat the procedure introduced there for different values of the current, since only one value of the current was calculated in that work.

The procedure, as explained in Ref. 8 is as follows. We simulate a disk with exchange constant A = 13pJ/m, saturation magnetization Ms = 8 × 105A/m, and to simplify the algebra we select an artificial value of K so that Q = 2. As initial configurations we create magnetization files from the numerical solution of Eqs. (5)-(6) (shown in Fig.1). We let the simulation run to reach one of two states: for small currents the soliton vanishes and the system relaxes to the uniform state (with mz = −1), for large currents the amplitude of the soliton increases towards a stable soliton with Θmax → 1. By tuning the current and following a bracketing-and-bisection procedure we numerically search for the value of the current where the relaxation behavior changes. After this we have two values of the current, I+ and I, which are very close to each other but the system relaxes to a high amplitude soliton for I+ and decays to the uniform state for I. This provides an estimate of the current for which the initial state is a saddle point of the energy landscape.

FIG. 1.

Soliton profiles for nanocontact radius ρ* = 6, and different amplitudes Θmax. These represent numerical solutions of Eqs. (5)-(6). The yellow region represents that in which the spin current flows.

FIG. 1.

Soliton profiles for nanocontact radius ρ* = 6, and different amplitudes Θmax. These represent numerical solutions of Eqs. (5)-(6). The yellow region represents that in which the spin current flows.

Close modal

A typical result is provided in Fig. 2(a) which shows the time evolution of the total energy. As can be seen, the energy changes very slowly at the beginning of the simulation for both current values, indicating that in the vicinity of the initial configuration the energy surface is almost flat. However, the final state is different for a very small change in the current, indicating that the initial state corresponds to the transition point between the basins of attraction of each of the final states. Subtracting the energies provides the values of the thermal excitation needed to overcome the barrier.

FIG. 2.

The structure of this graphic array is as follows: (left) Short-term (right) and long-term overdamped micromagnetic simulations for a nanocontact of diameter 68nm (ρ* = 6). (Top) Time evolution of the energy. (Middle) Average magnetization components for mx, my and (bottom) mz. The images between the top and middle panels are snapshots of magnetization configurations at the given times which are marked with points in panels (c-f). The major features of individual panels are as follows: (Panel a) Early relaxation of the total energy as a function of time. For comparison, the sum of conservative energy terms E0=Eex+EKeff and of the spin-torque energy term EST is also shown. This simulation starts with the configuration with amplitude Θmax = 0.7π (central inset). To the right, the current flowing through the nanocontact is I+ = 2.55993 mA and the magnetization relaxes to a high amplitude soliton (right inset); to the left the current is I = 2.55991 mA and the magnetization relaxes to the uniform state (left inset). The points in the graph of the total energy E indicate the times at which the configurations shown in the insets were taken. The activation barrier for creation of soliton Δ+ is obtained by subtracting the total energy of the uniform state from that of the initial state; and for annihilation, subtracting the total energy of the high amplitude soliton Δ+. (Panel b) Long time evolution of the energy terms. The inset shows a soliton configuration where the in-plane magnetization component winds along the soliton perimeter, resulting in a small energy increase δ. (Panels c and d) Average magnetization components as a function of time. To collect the data summarized in Fig. 3, the precessional periods of the saddle, stable soliton, and wounded soliton have been obtained from these images as shown with double-headed horizontal arrows.

FIG. 2.

The structure of this graphic array is as follows: (left) Short-term (right) and long-term overdamped micromagnetic simulations for a nanocontact of diameter 68nm (ρ* = 6). (Top) Time evolution of the energy. (Middle) Average magnetization components for mx, my and (bottom) mz. The images between the top and middle panels are snapshots of magnetization configurations at the given times which are marked with points in panels (c-f). The major features of individual panels are as follows: (Panel a) Early relaxation of the total energy as a function of time. For comparison, the sum of conservative energy terms E0=Eex+EKeff and of the spin-torque energy term EST is also shown. This simulation starts with the configuration with amplitude Θmax = 0.7π (central inset). To the right, the current flowing through the nanocontact is I+ = 2.55993 mA and the magnetization relaxes to a high amplitude soliton (right inset); to the left the current is I = 2.55991 mA and the magnetization relaxes to the uniform state (left inset). The points in the graph of the total energy E indicate the times at which the configurations shown in the insets were taken. The activation barrier for creation of soliton Δ+ is obtained by subtracting the total energy of the uniform state from that of the initial state; and for annihilation, subtracting the total energy of the high amplitude soliton Δ+. (Panel b) Long time evolution of the energy terms. The inset shows a soliton configuration where the in-plane magnetization component winds along the soliton perimeter, resulting in a small energy increase δ. (Panels c and d) Average magnetization components as a function of time. To collect the data summarized in Fig. 3, the precessional periods of the saddle, stable soliton, and wounded soliton have been obtained from these images as shown with double-headed horizontal arrows.

Close modal

Repeating this procedure we obtain the result summarized in Fig. 3, which confirms that the key features of the 1D model remain valid in two dimensions. There is a minimum current (Isus ≈ 1.15mA) below which the soliton cannot be sustained; at the other extreme, beyond a maximum value (Icrit ≈ 4.98mA) the uniform state is unstable. These two currents define the bistability region for which both states coexist. If we plot the observed frequency vs. applied current (inset of Fig. 3), our 2D simulations agree well with predictions from the 1D model. However, the activation barriers predicted by the 1D model are smaller than those obtained from the two-dimensional simulations, which likely arises from exchange energy contributions that result from variations of the azimuthal angle Φ, which the 1D model assumes to be spatially uniform. The soliton configuration shown in Fig. 2(a) presents a modulation of the in-plane component that supports this interpretation.

FIG. 3.

Main figure: activation barriers from overdamped OOMMF simulations. Inset: frequency as a function of applied current obtained from OOMMF simulations and compared to predictions from the 1D model.

FIG. 3.

Main figure: activation barriers from overdamped OOMMF simulations. Inset: frequency as a function of applied current obtained from OOMMF simulations and compared to predictions from the 1D model.

Close modal

We now discuss a feature that our 2D simulations in OOMMF revealed and which does not appear in a 1D model where Φ is taken to be uniform over the sample. As we have argued in Ref. 8, the second term in Eq. (4) is perpendicular to the field associated with EST; hence, rotations of the magnetization about mp will not cause an increase of EST, and any possible changes of energy must occur elsewhere. Furthermore, in a purely symmetric scenario where mp and hZ are perpendicular to the plane, the anisotropy and Zeeman energy should remain constant during precession if the soliton profile Θ(ρ) does not change.

Now consider a snapshot of a configuration of the magnetization in which Θ depends only on ρ, but Φ depends on the in-plane angular coordinate even as the precession across the sample occurs at the same rate, i.e. Θ = Θ(ρ) and Φ = ωt + Φ(t = 0, ϕ, ρ). Because Φ is no longer spatially uniform, there is a component of the exchange-induced field which is parallel to the direction of precession. When this occurs, the second term in (4) is responsible for a net change in the energy of the system. As we discuss next, based on Fig. 2 this unexpected situation is what our 2D simulations have revealed.

In the procedure outlined above, most simulations were stopped once it was possible to determine whether the system relaxed into the uniform state or the stable soliton. The left side of Fig. 2 shows the typical behavior of a system that relaxes into a stable soliton. If the simulations are allowed to continue, some unexpected features of the dynamics become apparent. As can be seen in the right side of Fig. 2, the total and conservative energy increase (panel b), while the average in-plane magnetization (panel d) continues to oscillate with the same frequency but with decreasing amplitude, forming a bottle-shaped envelope. The out-of-plane magnetization component (panel f) stays constant, as does the spin-torque energy. Because the anisotropy depends solely on mz and does not change, the increase in energy can only be attributed to the exchange term.

A look into the magnetization configurations confirms this. The size of the soliton is similar at short and long timescales, but at longer times Φ is no longer spatially uniform. Nevertheless, all spins seems to require the same time to complete their precessional cycles (see animation in supplementary material). Consequently these objects share the same frequency but have different in-plane magnetization averages. A comparison of vectors occupying the white bands in both images (where the magnetization is in-plane) shows that for long times the in-plane magnetization has considerable winding along that band, which will necessarily result in an increase of exchange energy.

It is worth mentioning that deformations of the perimeter of the soliton have been observed in simulations for strong currents when mp is tilted away from the normal to the plane.12 Here the soliton boundary deforms with multiple bulges. Similarly, we have previously predicted that a radial dependence Φ(ρ) should appear if the spin-torque has nonzero anisotropy parameter ν.8 The phenomenon described here is different from both of those situations: the perimeter remains close to circular but the azimuthal angle now depends on the angular coordinate Φ(ϕ).

Since these two separate configurations have similar frequencies and out-of-plane magnetization components, distinguishing them experimentally might prove challenging. Notice, however, that this new state has a slightly larger energy than the “unwound” soliton which indicates that a slightly lower thermal excitation might be needed for decay into the uniform state.

To summarize, we have confirmed that the 1D model predictions are robust when compared to more realistic 2D micromagnetic simulations. We have also found new dynamic states that require a 2D description. Further studies could explore the dynamic behavior of these objects in the underdamped regime.

Three short movies are provided as supplementary material with the following processes: relaxation into the uniform state, relaxation into the stable soliton, and precession at long times.

This research was supported in part by the US National Science Foundation Grant No. DMR 1610416 (D.L.S.) and the National Science Centre Poland under OPUS funding grant No. 2019/33/B/ST5/02013 (G.D.C.-O).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
F.
Macià
and
A. D.
Kent
, “
Magnetic droplet solitons
,”
J. Appl. Phys.
128
,
100901
(
2020
).
2.
A. M.
Kosevich
,
B. A.
Ivanov
, and
A. S.
Kovalev
, “
Magnetic solitons
,”
Phys. Rep.
194
,
117
238
(
1990
).
3.
M. A.
Hoefer
,
T. J.
Silva
, and
M. W.
Keller
, “
Theory for a dissipative droplet soliton excited by a spin torque nanocontact
,”
Phys. Rev. B
82
,
054432
(
2010
).
4.
S. M.
Mohseni
,
S. R.
Sani
,
J.
Persson
,
T. N. A.
Nguyen
,
S.
Chung
,
Y.
Pogoryelov
,
P. K.
Muduli
,
E.
Iacocca
,
A.
Eklund
,
R. K.
Dumas
,
S.
Bonetti
,
A.
Deac
,
M. A.
Hoefer
, and
J.
Åkerman
, “
Spin torque–generated magnetic droplet solitons
,”
Science
339
,
1295
1298
(
2013
).
5.
D.
Backes
,
F.
Macià
,
S.
Bonetti
,
R.
Kukreja
,
H.
Ohldag
, and
A. D.
Kent
, “
Direct observation of a localized magnetic soliton in a spin-transfer nanocontact
,”
Phys. Rev. Lett.
115
,
127205
(
2015
).
6.
R. O.
Moore
and
M. A.
Hoefer
, “
Stochastic ejection of nanocontact droplet solitons via drift instability
,”
Phys. Rev. B
100
,
014402
(
2019
).
7.
P.
Wills
,
E.
Iacocca
, and
M. A.
Hoefer
, “
Deterministic drift instability and stochastic thermal perturbations of magnetic dissipative droplet solitons
,”
Phys. Rev. B
93
,
144408
(
2016
).
8.
G. D.
Chaves-O’Flynn
and
D. L.
Stein
, “
Thermal activation barriers for creation and annihilation of magnetic droplet solitons in the presence of spin transfer torque
,”
Phys. Rev. B
101
,
184421
(
2020
).
9.
M. J.
Donahue
and
D.
Porter
, OOMMF User’s Guide, Version 1.0 Interagency Report No. NISTIR 6376,
1999
.
10.
J.
Xiao
,
A.
Zangwill
, and
M. D.
Stiles
, “
Boltzmann test of Slonczewski’s theory of spin-transfer torque
,”
Phys. Rev. B
70
,
172405
(
2004
).
11.
X.
Fong
, OOMMF_xf_stt [Source Code] https://github.com/seeder-research/OOMMF_xf_stt,
2015
.
12.
D.
Xiao
,
V.
Tiberkevich
,
Y. H.
Liu
,
Y. W.
Liu
,
S. M.
Mohseni
,
S.
Chung
,
M.
Ahlberg
,
A. N.
Slavin
,
J.
Åkerman
, and
Y.
Zhou
, “
Parametric autoexcitation of magnetic droplet soliton perimeter modes
,”
Phys. Rev. B
95
,
024106
(
2017
).

Supplementary Material