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 solitons^{4,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 *E*_{0} an additional potential *E*_{ST} 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\u2212\Delta EkBT$, where Δ*E* is the height of the barrier that separates two (local or global) energy minima, and the prefactor *f*_{0} 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:

The conservative energy density is:

where *A* is the exchange constant, *K* is the magnetocrystalline energy constant, *μ*_{0} is the permeability of free space, **H**_{Z} 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 $\mu 0Ms22(Q\u22121)$, with $Q=2K\mu 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 **m**_{p} 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:

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=\u22121\mu 0Ms\delta E\delta m$.

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

where $L=\u2212\gamma \u2032m\xd7+\alpha m\xd7m\xd7$ is a linear operator that accounts for precession and damping, $\gamma \u2032=\gamma 01+\alpha 2$ and $\gamma 0=2.211\xd7105mA\u22c5s$. 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 ($\u2207E=0$) appear stationary ($m\u0307=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, $\rho lex=rQ\u22121$, with $lex=2A\mu 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=HZz\u0302Ms(Q\u22121)$ and current $\sigma =J\u210f\mu 0edMs2Q\u22121$, and the current has a profile *V*(*ρ*) in the form of a Heavyside step function $V(\rho )=H(\rho \u2212\rho *)$ where *ρ** is the nanocontact radius. A persistent soliton profile (for which $\Theta \u0307=0$) can be either an energy minimum or a transition state and satisfies the differential equation:^{8}

with boundary conditions

The solutions to Eqs. (5-6) are more easily found by imposing a third boundary condition Θ(*ρ* = 0) = Θ_{max} and solving for $\omega h=\sigma V(\rho )\alpha \u2212hZ$ as an eigenvalue problem. The amplitude *m*_{z}(*ρ* = 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 $\omega 0=\sigma \u2032\alpha $. 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 *m*_{z}(*ρ* = 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 currents^{10} and added an energy term without an associated effective field (we were motivated by a contribution from Fong^{11} 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 *M*_{s} = 8 × 10^{5}A/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 *m*_{z} = −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.

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.

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 (*I*_{sus} ≈ 1.15mA) below which the soliton cannot be sustained; at the other extreme, beyond a maximum value (*I*_{crit} ≈ 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.

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 **m**_{p} will not cause an increase of $EST$, and any possible changes of energy must occur elsewhere. Furthermore, in a purely symmetric scenario where **m**_{p} and **h**_{Z} 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 *m*_{z} 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 **m**_{p} 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).

## DATA AVAILABILITY

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