Homotopy transitions and 3D magnetic solitons

This work provides a concept for three-dimensional magnetic solitons based on mapping the homotopy path between various two-dimensional solutions onto the third spatial axis. The representative examples of statically stable configurations of that type in the model of an isotropic chiral magnet are provided. Various static and dynamic properties of such three-dimensional magnetic solitons are discussed in detail.


I. INTRODUCTION
Chiral magnets are a distinct class of materials, where the competition between the Heisenberg exchange and the chiral Dzyaloshinskii-Moria interaction 1,2 (DMI) gives rise to a wide diversity of magnetic solitons -localized magnetic textures with particle-like properties. The most representative example of magnetic solitons in chiral magnets are magnetic skyrmions which have been intensively studied both theoretically [3][4][5] and experimentally [6][7][8][9][10][11][12] in the last decades. In thick films and bulky samples, the magnetization vector field of chiral skyrmions resembles a filamentary structure composed of vortex-like tubes.
More recent studies have revealed a diversity of 3D solitons in chiral magnets beyond skyrmions. The most prominent examples are chiral bobbers 13,14 , skyrmion bags [15][16][17][18] , heliknotons 19 , and skyrmion braids 20 . Because of distinct topology, the above solitons possess different static, dynamic, and transport properties. In particular, the chiral bobbers are distinct by the presence of magnetic singularity -Bloch point, while other magnetic textures are characterized by smooth magnetic vector fields. Because of that, contrary to other solitons, the chiral bobber is not a topological soliton. The skyrmions and heliknotons, on the other hand, belong to different topological groups characterized by distinct topological invariants.
In this work, we present another type of magnetic solitons which belong to a topological group of skyrmions and are stabilized in bulk crystals of chiral magnets. Because of the unique properties of these solutions, it is reasonable to attribute them to a distinct class of magnetic solitons. Here we will refer to them as hybrid skyrmion tubes. The crosssections of well-studied skyrmion tubes usually represent nearly identical configurations with small additional modulations due to the braiding effect 20 or the presence of the noncollinear background 21 and/or free surfaces 22 . We show that besides such nearly homogeneous tubes, there are also solutions where the crosssections represent a continuous transformation between the skyrmions of different configurations but the same topological index -skyrmions of the same homotopy group. In a representative example, we illustrate the stability of such exotic spin textures and discuss their unique dynamic properties. Besides that, we discuss an important consequence of applying such a homotopical concept to truly 3D localized magnetic spin textures. In particular, we present the solution represent-ing a topologically trivial but statically stable and truly threedimensional soliton -3D chiral droplet.

II. MODEL
The model Hamiltonian for chiral magnet with bulk type DMI has the folowing form: where n = M/M s is normalized magnetization |n| = 1, A and D are micromagnetic constants of the Heisenberg exchange and the DMI, respectively. V m is the volume of the magnetic sample. The potential energy term U(n) includes an uniaxial (easy-axis or easy-plane) magnetic anisotropy, K u , and the interaction with the external magnetic field, B ext e z : Following the standard procedure 3-5 , the functional (1) can be written in a more suited for analysis, dimensionless form where E = E/2A is the reduced energy, u = K u /M s B D and h = B ext /B D are the dimensionless values of anisotropy constant and the strength of the external magnetic field, respectively. The integration in (3) is over the reduced volume The characteristic parameters L D = 4πA /D and B D = D 2 /2M s A are the period of spin spiral and the saturation field for isotropic case, K u = 0.
The energy functional (3) remains valid for the 2D (or quasi 2D) systems where the magnetization does not change along the z-axis and integration is carried out over dV = ldxdy, where l is the film thickness. The 2D model of the chiral magnet, besides ordinary axially symmetric π-skyrmions, provides a large variety of magnetic solitons e.g. skyrmion bags 15,16 and skyrmions with chiral kinks 17 . The homotopy classification of localized solutions in 2D is provided by the following invariant The solutions with identical integer index Q are called the solutions of one homotopy class. The latter implies that one can continuously transform the vector fields of such solutions into each other. For instance, let us assume there are two localized magnetic textures n 1 (x, y) and n 2 (x, y) of identical charge Q. Then, there is a vector field transformation, n(x, y; s), where s ∈ [0, 1], such that for n(x, y; 0) = n 1 (x, y) and n(x, y; 1) = n 2 (x, y). When n(x, y; s) is differentiable with respect to x, y, and s at any point such transformation can be called a homotopy path. Since for any n 1 (x, y) and n 2 (x, y), there is an infinite number of homotopy paths, there is no unique method to construct such paths. To make the search more definite, we consider only those homotopy paths that also satisfy the minimum energy path (MEP) criterion.
In particular, we use geodesic nudged elastic band (GNEB) method 23,24 implemented in Spirit code 25 . The details of the MEP calculation are provided in the Appendix A.

III. RESULTS
In Fig. 1a we provide representative example of the MEP between two skyrmions with Q = 1. We show only the spin texture corresponding to local minima and saddle points. Noticeably, the textures corresponding to the central minimum state (d) and saddle points (c), (e) contain a chiral kink 17 while another minimum state is a skyrmion bag free of kinks (b). The MEP presented in Fig. 1a satisfy the above criteria for the homotopy path. The parameter s can be associated with the reaction coordinate which has a meaning of the relative distance between the images (snapshots of the vector field) in the multidimensional parameter space.
To construct an initial configuration for the 3D skyrmion tube, one can take the stable 2D skyrmion configuration and place it at each xy-plane of the 3D simulated domain. Statically stable configuration of such homogeneous skyrmion tube corresponding to the 2D skyrmion in Fig. 1(b) is provided in Fig. 2(a). To construct the initial state for a nonhomogeneous 3D skyrmion tube, we use a mapping from the homotopy path to the third spatial axis, s → z. In other words, to create a 3D magnetic texture, we sequentially lay down the spin texture of the images from MEP on top of each other along with the z-axis. The statically stable spin configuration obtained by the energy minimization of that initial state is provided in Fig. 2(b) and its crosssections are shown in (c)-(h). The intermediate region resembling a knot on the isosurface of the skyrmion string in Fig. 2(b) is well localized and thus can be thought of as a soliton that is hosted by a skyrmion string. This nonhomogeneous 3D skyrmion tube is a representative example of the solutions to which we refer as hybrid skyrmion tubes. For conciseness, the localized intermediate region, we call a knot in the following. The choice of such terminology is justified by pure visual analogy and has nothing to do with the topological knots. Whether one can continuously unwind that knot without the appearance of Bloch points represents an intriguing question that, however, goes out of the scope of the present work.
At the strong magnetic field, h ≥ 1 − 2u, the hybrid skyrmion tubes represent a metastable state embedded in the ferromagnetic vacuum. For h < 1 − 2u, the vacuum for such states is the cone phase. The range of fields and anisotropies where hybrid skyrmion tubes remain stable depends on the particular configuration -the type of 2D skyrmions in the tube crosssections. For instance, the hybrid skyrmion tube depicted in Fig. 2 In general, the hybrid skyrmion tubes can be composed of a few intermediate states representing stable 2D skyrmions. Here we consider the solution composed of two intermediate states only. Such configurations are easier to handle because they satisfy periodical boundary conditions. It is worth noting, the stabilization of hybrid skyrmion tubes does not require the periodic boundary conditions along z-axis. The latter makes the experimental observation of hybrid skyrmion tubes in thick films of chiral magnets quite promising. In the case of free boundaries along z-axis, the tube has an additional surface twist modulations. For thick enough films (∼ 6L D = 420 nm for FeGe) These modulations do not decrease the tube stability. On the other hand the presence of the free boundaries the hybrid skyrmion tube depicted in Fig. 2 (b) can continuously unwind into the host skyrmion tube shown in Fig. 2 (a). This process is illustrated by Supplementary Movie 1, where each frame was obtained by a manual shift of the entire spin texture along the z-axis and following incomplete relaxation of the system. In the case of complete relaxation, the system converges to its initial state, which indicates the presence of the energy barrier for the knot to escape through the free boundary.
To demonstrate the unique dynamic properties of hybrid skyrmion tubes we perform micromagnetic simulations based on the Landau-Lifshitz-Gilbert equation 26 : where γ is the gyromagnetic ratio, α is the Gilbert damping, is the effective field. The last term in (5) is the Zhang-Li torque 27 because of the electric current: where the vector I = jµ B p(1+ξ 2 ) −1 (eM s ) −1 is proportional to the current density j, ξ is the degree of non-adiabaticity, p is the polarization of the spin current, µ B is the Bohr magneton and e is the electron charge. When the 3D soliton is embedded in the FM state, the current direction can be chosen arbitrary. On the contrary, when the soliton is embedded in non homogeneous vacuum, e.g., cone phase, it is essential to chose the current direction perpendicular to the q-vector of the cone to prevent the excitation of the background. To simplify the following discussions we consider the case of hybrid skyrmion tube in ferromagnetic vacuum.
Supplementary Movie 2 illustrates the dynamics of the hybrid skyrmion tube when the current is along the z-axis. The simulations were performed with Mumax 28 code. The example of Mumax script and the initial state files are provided in Supplementary Data. We used periodical boundary conditions in all spatial directions to simulate a bulk crystal. Under the electric current parallel to the skyrmion tube, the knot on the isosurfaces moves along the tube in the direction opposite to the current. In this case, the position of the host skyrmion tube in the xy-plane remains fixed.
In Supplementary Movie 3, we show the dynamics of the skyrmion tube when the current is perpendicular to the skyrmion tube, I e x . In this case, both the host skyrmion tube and the knot move. As a result, all three components of the knot velocity are non-zero.
To quantify the knot velocity, we follow the approach used in our previous work on 2D skyrmion dynamics 29 . Extending this approach to 3D textures the position of the knot R = R x , R y , R z can be defined as follows where N i j (r k ) = n z (r)dr i dr j is magnon density averaged in r k -plane, indices i, j, k ∈ {x, y, z}. Sigh ± accounts for the direction of solitons motion: along the basis vector e k (+) or in the opposite direction (−). The integer l k is the number of times the soliton crossed the corresponding domain boundary. By tracing the position R at each moment in time, we can estimate the instant velocity of the knot, v = dR/dt, and compare it with the results of a semi-analytical approach based on the method of collective coordinates suggested by Thiele 30 . Assuming the rigid motion of magnetic textures with velocity v, i.e. n(r,t) = n(r − vt), from (5) and (6), one can derive the Thiele equation: where the gyro-vector, G, and the dissipation tensor,Γ, have components G i and Γ i j defined as: where ε i jk is Levi-Civita symbol. The advantage of the Thiele approach is that the solution for the soliton velocity can be written explicitly. For instance, for I = Ie x , the solution of (8) for the uniform skyrmion tubes, Γ xz = Γ yz = Γ zz = 0 takes a form identical to the 2D case v z = 0.
For the habridized skyrmion tube all the entries of the dissipation tensorΓ have non-zero values and the solution of (8) takes the form It is worth noting two remarkable features of the solution (12). Firstly, there is no continuous transition between (12)  and (11) even when Γ xz , Γ yz Γ zz tend to zero. Secondly, in the most general case of 3D soliton, the components of gyrovector G can be thought of as 2D topological charges in corresponding x, y and z crosssections. For hybrid skyrmion tubes, however, similar to the 2D case, only the third component of the gyro-vector is non zero and proportional to 2D topological charge in xy-plane, G = 4πQL z e z .
To calculate the velocity components in (12), one has to find the components of the gyro-vector (9) and the dissipation tensor (10) for a particular magnetization distribution. For that, we use the magnetization vector field at the static equilibrium obtained by the numerical energy minimization of the functional (1). It is worth emphasizing that the soliton velocity obtained from the solutions of the Thiele equation (8) corresponds to the velocities in a steady-state only. Thus the velocities estimated from LLG simulations can be compared with the analytical solutions only when the transient state is over, and all components of the v reached the saturation.
In the most general case of 3D soliton, when all three components of the velocity are non-zero, one can introduce two deflection angles β 1 = arctan v y /v x and β 2 = π/2 − arctan v 2 x + v 2 y /v z . The angle β 1 is also known as skyrmion Hall angle. Because of that we refer to β 2 as the second skyrmion Hall angle. Table I demonstrates a good agreement between the values of β 1 and β 2 estimated from LLG simulations and calculated with the Thiele approach for parameters j = −5 · 10 8 e x A/m 2 , α = 0.01 and ξ = 0.05. Noticeably, for particular tube depicted in Fig. 2(b), the angle β 2 is larger than β 1 . As follows from the Thiele equation (8), when I = Ie z , in agreement with the micromagnetic simulations, only the z-component of the knot velocity is non zero, v z = −Iξ /α. The homotopy concept used for constructing hybrid skyrmion tubes can be extended further for solitons localized in all three dimensions. For that one can consider the homotopy transformation of any 2D skyrmion with Q = 0 into the FM phase. Here we examine such transformation on the example of chiral droplet shown in Fig. 3. The MEP in Fig. 3 is quite similar to that shown in Fig. 1. Two stable states representing the isolated soliton and the saturated FM state are separated by the saddle point. In this transformation, the collapse of the droplet goes by its shrinking, which has been previously discussed in Refs. 31,32 .
By mapping the obtained images of the corresponding homotopy path onto the spatial z-axis, we constructed the initial state of the 3D soliton. However, it turned out that such an initial guess does not lead to a stable configuration in the FM phase but instead becomes stable in the cone phase vacuum only. The isosurfaces and cross-sections of the stable configuration of the 3D chiral droplet obtained by numerical energy minimization are shown in Fig. 4. As seen in (a) the red color in the surface is missing, which means that by projecting the vector field of the 3D chiral droplet onto S 2 sphere, one can not cover it entirely. The latter means that both topological invariant (4) and Hopf invariant are zero in this case. The cross-sections in (b) and (c) show well-localization of the soliton. We estimated the characteristic size of the 3D chiral droplet equal to ∼ 1L D in all spatial directions.
The stability range for the 3D chiral droplet in terms of the magnetic field, h, and anisotropy, u, is shown in Fig. 5. Note, the phase transition lines between the phases are reproduced from Ref. 33 . As follows from the diagram, the 3D chiral droplet stability region is in the range where the skyrmion lattice is the lowest energy state, while the cone phase is a meta-stable state. The blurred edges at h 0.28 and u 0.15 denote that we did not succeed to identify the stability range below these values reliably. At these values of magnetic field and anisotropy, the helicoid competes in energy with the cone phase. A more precise estimation of the lower bound for the 3D chiral droplet stability requires much larger sizes of the simulated domain above the limit of our computational resources.
The 3D chiral droplet motion can be induced by applying different external stimuli. First, let us consider the motion of the 3D chiral droplet induced by the external magnetic field slightly tilted and rotating about the z-axis, h = 0.34(sin θ h cos 2πνt, sin θ h sin 2πνt, cos θ h ). We performed the micromagnetic simulations assuming the tilt angle θ h = 0.1 and the small frequency of rotation ν = 5 MHz which is close to a quasi-stationary process. Supplementary Movie 4 demonstrates such motion of the 3D chiral droplet consisting of two types of motion -the translation along the z-axis and the rotation about the z-axis. Remarkably, similar dynamics has been reported for heliknoton 19 and seems to be a common property of all 3D localized solitons hosted by the conical phase. The motion of the 3D chiral droplet free of its rotation can be induced by the cone phase rotation at the static external field. To demonstrate this, we performed the simulations where we manually rotated a small number of spins located far from the droplet. As illustrated in Supplementary Movie 5, the translation motion of the droplet along the z-axis, in this case, occurs without its rotation.
In the case of Zhang-Li torque, the gyro vector of the 3D droplet equals zero and the Thiele equation (8) provides a trivial solution, v = ξ I/α. On the other hand, due to the presence of the cone phase, one has to choose I ⊥ q to suppress the excitation of the cone phase itself. The analysis of the cone phase dynamics induced by I q, is provided in Appendix B.
One has to note that contrary to the Thiele equation, which predicts that the 3D chiral droplet should move without deflection, v −I, in numerical simulations, we still observe a slight deflection in the soliton motion. Since, in the case of the hybrid skyrmion tube, we got pretty well agreement between numerical simulations and the Thiele approach, we attribute this effect to the numerical artifact caused by the excitation of the cone phase, which in turn prevents reaching the regime of a steady motion. Supplementary Movie 6 shows the dynamics of the droplet caused by the in-plane electric current I e x . To prevent the excitation of the cone phase, in this simulation, we pinned the spins on the bottom plane of the simulated domain.

IV. CONCLUSIONS
In this work, we have discussed new types of 3D magnetic solitons stabilized in chiral magnets -hybrid skyrmion tube and 3D chiral droplet. We show that the magnetic textures of these solutions can be explained in terms of homotopy transitions between various 2D skyrmions.
We have studied the static and dynamic properties of new solitons. The parameters of the external magnetic field and magnetic anisotropy at which the presented solutions remain stabilized are different and characterized by the various states of the vacuum. hybrid skyrmion tubes can be stable in the FM state and cone phase, while the 3D chiral droplet is stable only in the conical phase surroundings. It is shown that the dynamics of the hybrid skyrmion tube can be induced by the electric current modeled by the Zhang-Li torque in the LLG equation. The electric current applied along the tube causes the motion of the knot on the skyrmion tube in the direction opposite to the current. The electric current applied perpendicular to the skyrmion tube leads to more complicated dynamics, which can be described in terms of two skyrmion Hall angles. The results of LLG simulations agree well with the analysis based on the Thiele equation.
For the case of the 3D chiral droplet, we show that its motion besides the electric current can be induced by the rotating The MEP in Fig. 1 (a) represents a homotopy path between two states (b) and (d). To get the MEP without singularities, one has to make a reasonable assumption for the initial path. The straightforward and often used approach based on the interpolation between two configurations by simple rotation of the spins is not appropriate for this purpose. It is easy to check that the topological index Q is often not conserved along the MEP with this approach. To avoid such a behavior of the GNEB solver, we have constructed the ansatz for the initial path using the energy minimization and drag-and-drop function implemented in Spirit. In particular, as seen in Fig. 1(b)-(e) the intermedeate states represents the elongation and twisting of the magnetic texture near the chiral kink. Starting with the state depicted in Fig. 1(b) and using the drag-and-drop op-tion we enforce the chiral kink to elongate and bend to form a state similar to that in (c) and (e). We intentionally stop the relaxation process before the system converges to one of two local minima as in Fig. 1(b) or (d). Then, we collect such unrelaxed snapshots of the spin texture with the different levels of elongation of the part containing the chiral kink and use them as an initial guess for MEP.
The calculation were perfomed at the following parameters: Appendix B: Cone phase rotational dynamics In terms of the spherical angles (Θ, Φ) the magnetization writes n = (sin Θ cos Φ, sin Θ sin Φ, cos Θ) T . Assuming that at t = 0, the magnetization profile is given by the cone phase Θ = Θ c , Φ = 2πz + φ 0 , we can find the analytic solution of the LLG equation with the current I = Ie z turned on at t > 0. The uniformity of the magnetic texture in the xy-plane leads to the following system of equations for (Θ, Φ): where τ = γB D t/4π 2 is the dimensionless time and i = 4π 2 I/γL D B D is the parameter of electric current. Although the system (B1) represents two coupled non-linear partial differential equations, its solution can be written for the case when Θ = Θ(τ). We found that strong currents can lead to the transition from the cone phase to the FM state. The critical current at which this happens is given by: The formula (B2) has sense only at ξ = α. In the case ξ = α, the cone phase can not be excited by the current. Below the critical regime, i.e. i ∈ (i − c , i + c ), the angle Θ monotonically changes from Θ c -the equilibrium cone angle without electric current to Θ I -the equilibrium cone angle in dynamical steady state in the presence of the current: For i / ∈ (i − c , i + c ) the angle Θ I equals to 0 or π depending on the sign of the current, i. At i ∈ (i − c , i + c ), the solution of (B1) for Θ can be written as: where a 1 = h − i(ξ − α)/2πα, a 2 = 1 − 2u. The solution for Φ function is given by: As follows from (B4), at τ → ∞, one has cos Θ = cos Θ I = a 1 /a 2 . In this case, Eq. (B5) describes the rotation of the cone phase with constant velocity Iξ /α which agrees with the Thiele prediction for topologically trivial configurations. No- ticeably, the dynamics described by the LLG equation reaches the dynamical steady state exponentially fast. Note, the solutions (B4), (B5) remain valid not only for the bulk systems but also for the films with free boundaries along z-axis. The latter remain valid for any h and u where the surface modulations 13 do not perturb the cone phase. This follows from the fact that found solutions satisfy the boundary conditions: ∂ z Θ = 0, ∂ z Φ = 2π on the free surface, z = const, for any τ.
To verify the found solutions, we compare them with the results of LLG simulations performed for different values of the current but fixed values of h = 0.34, u = 0.26, α = 0.01, and ξ = 0.05. The critical current values (B2) for these parameters are: i − c ≈ −0.22, i + c ≈ 1.29. We have performed simulations for two current values within the critical range: i 1 = −0.15, i 2 = 0.15, and for one out of this range i 3 = 2. The results are provided in Fig. 6. For all three cases we see a good agreement between numerical and analytical solutions at least for the chosen simulation time. At longer times, z-component of the magnetization tends to the limit value accordingly to (B3). For the cases shown in (a) and (b) these limit values are about 0.51 and 0.91, respectively. For the case shown in (c) corresponding to the current above the critical value, the magnetization tends to −1. As one can deduce, applying the current I q allows manipulating the cone phase angle and its dynamics in a controllable way.
The found analytical solutions besides the pure academic interests can be used also for testing the accuracy of numerical schemes for solving the LLG equation. In particular, the solutions of (B1) in a special case of zero damping, α = 0, can be written in a more compact form: The stability of LLG solvers typically requires the presence of non-zero damping. Therefore, the solutions (B6) can be useful for testing new LLG solvers which are free of this limitation. Both cases α = 0 and α = 0 can be generalized for the case of time-dependent currents I = I(t) easily. The detailed analysis of this case is out of the scope of the present study and will be provided elsewhere.