Heliknoton in a film of cubic chiral magnet

1 Introduction

Cubic chiral magnets have attracted significant theoretical and experimental attention due to the vast diversity of experimentally observed topological magnetic solitons. These materials include various Si- and Ge-based alloys with B20-type crystal structures, such as Fe_1−xCo_x[1]; [2], [3]; [4]; [5]; [6], [7], and others [8]. The competition between the Heisenberg exchange interaction and the chiral Dzyaloshinskii-Moriya interaction [9]; [10] (DMI) stabilizes topological solitons in these materials. The most extensively studied types of magnetic solitons in these systems are magnetic skyrmions [11]; [12]. Skyrmions are vortex-like strings or tubes characterized by a topological index Q = −1. Experimental observations have revealed that clusters of skyrmion tubes can form complex three-dimensional superstructures, such as skyrmion braids [13]. Additionally, the skyrmion antiparticle, antiskyrmions, with topological charge Q = +1, have been observed in thin films of FeGe [14]. Skyrmion bags with arbitrary topological charge and skyrmions with chiral kinks are other solitons that have been reported [15]; [16]; [17]. Recently, the experimental observation of skyrmion bags with positive topological charge and their current-induced motion has been reported by Tang et al. [18].

Skyrmions, antiskyrmions, and skyrmion bags are two-dimensional (2D) topological solitons that are localized in the plane of the sample and confined by the free surfaces of the sample in the third dimension. The homotopy classification of these solitons is based on the continuous mapping between two spheres, ${\mathbb{S}}^2\to {\mathbb{S}}^2$. Where the pre-image is a sphere that is homeomorphic to the two-dimensional region where the skyrmion is localized, while the image is a sphere that corresponds to a parameter space of the magnetization unit field, given by n(r) = M(r)/M_s. For the magnetic texture in xy-plane the corresponding topological charge [19]:

$Q=\frac{1}{4\pi }\int \left(\mathbf{F}\cdot {\widehat{\mathbf{e}}}_z\right)\mathrm{d}x\mathrm{d}y,(1)$

where the vector field

$\mathbf{F}=\left(\begin{array}{c}\hfill \mathbf{n}\cdot \left[{\partial }_y\mathbf{n}\times {\partial }_z\mathbf{n}\right]\hfill \\ \hfill \mathbf{n}\cdot \left[{\partial }_z\mathbf{n}\times {\partial }_x\mathbf{n}\right]\hfill \\ \hfill \mathbf{n}\cdot \left[{\partial }_x\mathbf{n}\times {\partial }_y\mathbf{n}\right]\hfill \end{array}\right)(2)$

is the curvature vector of the vector field n and is frame-invariant, meaning it does not depend on the choice of coordinate frame [20,21]. For an unambiguous definition of Q, the direction of the z-axis should be chosen consistently with the magnetization n₀ around the region of texture localization such that ${\widehat{\mathbf{e}}}_z\cdot {\mathbf{n}}_0>0$, see details in Refs. [22]; [23].

Three-dimensional (3D) topological magnetic solitons belong to a distinct class of solutions that are localized in all three spatial dimensions. The classification of 3D magnetic solitons is based on the mapping ${\mathbb{S}}^3\to {\mathbb{S}}^2$, which was originally proposed by Hopf [24]. The Hopf index H, which is the 3D analog of the topological invariant Q, was derived by Whitehead [25] and is given by:

$H=-\frac{1}{16{\pi }^2}\int \left(\mathbf{F}\cdot \mathbf{V}\right)\mathrm{d}x\mathrm{d}y\mathrm{d}z,(3)$

where V is a vector potential satisfying ∇ ×V = F. For localized 3D magnetic textures without singularities, the Hopf index H is an integer that has the following meaning: two morphologically distinct magnetic textures with identical H can be continuously transformed into each other without the appearance of singularities, whereas continuous transformations between magnetic configurations with different H are impossible. Statically stable structures with non-zero Hopf index are commonly referred to as Hopf solitons or hopfions. It should be noted that the topological charge Q is zero for any hopfion cross-section. An example of a hopfion texture with H = 1 is depicted in the first column of Figure 1.

FIGURE 1

FIGURE 1. The images illustrate the spin texture inside the simulated box of the size 4L_D × 4L_D in the xy-plane and thickness of 2L_D along the z-axis. The first column of images correspond to the hopfion anzatz (5). The second column represents the hopfion ansatz after applying the spiralization (6) with the k-vector of the spiral parallel to the y-axis, |k| = 2π/L_D. The third column of images corresponds to the spin texture after the energy minimization assuming periodical boundary conditions in the xy-plane and free boundaries along the z-axis. The fourth column illustrates the relaxed spin texture after applying despiralization. The first row of images shows the isosurfaces n_z = 0. The images depicted in the second, third, and fourth rows illustrate the magnetization field in the middle planes. For illustrative purposes, these images are bounded by the blue box of the size 2.5L_D × 2.5L_D × 2L_D depicted in the top left image. The scale bar in the bottom right image corresponds to 1L_D.

Theoretical studies of hopfions in various models of magnetic crystals began in the 1970s and continue to the present day [26]; [27]; [28]; [29]; [30]; [31]; [32]; [33]; [34]. There is also some limited experimental evidence for direct observation of such magnetic textures [35]; [23]. An intriguing type of hopfion has been predicted in cubic chiral magnets [31], but has not been observed so far. This hopfion is embedded not in a collinear ferromagnetic state but in a helical spin spiral state. We refer to it as a heliknoton following the terminology used in Refs.[36]; [31].

The primary goal of this study is to investigate the stability of the heliknoton under realistic conditions, taking into account the demagnetizing field effects, and estimate its stability range in films of different thicknesses. Additionally, we explore how the contrast in theoretical Lorentz transmission electron microscopy (TEM) images changes with film thickness and suggest optimal conditions for experimental observations of the heliknoton. Furthermore, we estimate the stability of the heliknoton at finite temperatures using two complementary methods: harmonic transition state theory (HTST) [37]; [38] and stochastic spin dynamics [39].

2 Model

To estimate the stability of heliknoton in the films of finite thickness, we use a micromagnetic model, which contains the Heisenberg exchange term, the DMI, the Zeeman interaction, and the self-energy of the demagnetizing field [13]; [40]:

$\mathcal{E}=\underset{{\mathbb{R}}^3}{\int }\left[\mathcal{A}|\nabla \mathbf{n}{|}^2+\mathcal{D}\mathbf{n}\cdot \left(\nabla \times \mathbf{n}\right)-M_{\text{s}}\mathbf{n}\cdot \mathbf{B}+\frac{|\nabla \times \mathbf{A}{|}^2}{2{\mu }_0}\right]d\mathbf{r},(4)$

where n(r) is the magnetization unit vector field, M_s is the saturation magnetization, $\mathcal{A}$ is the exchange stiffness constant, $\mathcal{D}$ is the constant of isotropic bulk DMI, and μ₀ is the vacuum permeability. The magnetic field B(r) in Eq. 4 represents the sum of a homogeneous external magnetic field B_ext and the demagnetizing field produced by the sample itself, B = B_ext + ∇ ×A, where A(r) is the magnetic vector potential due to the presence of the magnetization field M(r) = M_sn. It should be noted that while the unit vector n is defined only within the volume of the sample, the vector potential A(r) is defined in the entire space of ${\mathbb{R}}^3$. As a result, the integration for the first three terms in Eq. 4 must be carried out only within the sample volume. For the present analysis, we will consider an infinite plate lying in the xy-plane and assume that the external magnetic field is parallel to the z-axis.

The static equilibrium heliknoton was obtained by numerically minimizing (4) with respect to the pair of fields n and A using a nonlinear conjugate gradient method (for details, see Ref. [13]). The simulated domain was discretized on a regular mesh with 256 × 256 × 181 nodes along the x, y, and z axes, respectively. To simulate an extended film, periodic boundary conditions were applied in the xy-plane. The results presented below for the micromagnetic calculations were obtained using GPU-accelerated software Excalibur [41]. Additionally, the results were verified with the publicly available software MuMax3 [42]. For definiteness, we used the material parameters for FeGe [43]: $\mathcal{A}=4.75$ pJm⁻¹, $\mathcal{D}=0.853$ mJm⁻², and M_s = 384 kAm⁻¹. To ensure that our results are generally applicable, we utilize reduced units for the external magnetic fields and distances throughout the paper. Specifically, we employ units concerning the cone saturation field $B_{\text{D}}={\mathcal{D}}^2/2M_{\text{s}}\mathcal{A}$ and the equilibrium period of the spin spiral at the ground state $L_{\text{D}}=4\pi \mathcal{A}/\mathcal{D}$ (70 nm for FeGe). As has been shown earlier [43], when accounting for demagnetizing fields, the critical field for cone saturation changes to B_c = B_D + μ₀M_s (682 mT for FeGe), and B_c = B_D (199 mT for FeGe) only if the demagnetizing field is neglected.

3 Results

3.1 Initial guess

The initial configuration for heliknoton can be obtained in two steps. First, we insert in the simulated domain the classical hedgehog-based ansatz [24]; [44] for a hopfion in the ferromagnetic background:

$\mathbf{n}=\frac{2{\sin }^{2}G}{r^2}\left(\begin{array}{c}\hfill xz\hfill \\ \hfill -yz\hfill \\ \hfill z^2\hfill \end{array}\right)+\frac{1}{r}\left(\begin{array}{c}\hfill y\sin 2G\hfill \\ \hfill x\sin 2G\hfill \\ \hfill r\cos 2G\hfill \end{array}\right),(5)$

where function G = G(r) with $r=\sqrt{x^2+y^2+z^2}$ describes the skyrmion-like profile: G (r = 0) = π and G (r → ∞) = 0. In our simulations we use G(r) = 2 arctan (exp (−2r)/r). The first column in Figure 1 illustrates the spin texture obtained with the ansatz (5). Note the opposite chirality of the spin texture in x and y directions.

In the second step, we implemented the spiralization procedure, which involves rotating all spins by an angle dependent on their coordinate along the wavevector k of the spiral. This transforms the ferromagnetic background into a helical one. Note that we applied this procedure to the entire volume, including the area containing the hopfion. For the case k‖e_y the spiralization procedure can be written as

${\mathbf{n}}^{\prime }=\left(\begin{array}{ccc}\hfill \cos \phi \hfill & \hfill 0\hfill & \hfill \sin \phi \hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill -\sin \phi \hfill & \hfill 0\hfill & \hfill \cos \phi \hfill \end{array}\right)\mathbf{n},\phi =\frac{2\pi }{L_{\text{D}}}y.(6)$

where n corresponds to the ansatz (5) and n′ is the ansatz for the hopfion in the helical background. The second column in Figure 1 illustrates the spin texture obtained after applying the spiralization (6) to the initial hopfion ansatz (5). The spin texture after energy minimization is shown in the third column of Figure 1. The quality of our ansatz is seen from the comparison of the spin textures before and after energy minimization. To better visualize the spin texture of the heliknoton, we utilize the despiralization procedure, which is illustrated in the last column of Figure 1. This procedure serves as the inverse of the spiralization operation (6). After the energy minimization or capturing a snapshot of the magnetic texture, we apply the transformation (6) with φ↦ − φ to obtain the despiralized spin texture. By using this procedure, one can gain a better understanding of the heliknoton’s spin texture.

3.2 Demagnetizing field effect

To estimate the demagnetizing field effect on the stability of heliknoton in a film of finite thickness, we performed systematic energy minimization at different external magnetic fields and film thickness starting with the ansatz (5)-(6). The results of these calculations performed with and without demagnetizing fields are presented in Figure 2. In these calculations, we keep the size of the simulated domain in xy-plane fixed, 4L_D × 4L_D, and vary only the film thickness between 0.5L_D and 3.5L_D. We have checked the interaction between heliknotons due to the periodic boundary conditions and found that it is relatively weak for this size of the simulation domain and does not significantly affect the heliknoton stability. Since the absolute values of the critical fields for the cases with and without demagnetizing fields are very different, we provide the critical fields in reduced units with respect to the saturation field of conical phase, B_c. This approach objectively estimates the range of heliknoton stability in both cases.

FIGURE 2

FIGURE 2. The diagram of stability for heliknoton in extended film calculated with and without demagnetizing field. For consistency, the external fields for both cases are given in reduced units with respect to the saturation field of the cone phase, B_c = B_D + μ₀M_s and B_c = B_D, for the case with and without demagnetizing filed, respectively. For B20-type FeGe, B_c = B_D + μ₀M_s = 0.682 T and L_D = 70 nm, see the top axis.

As follows from the diagram in Figure 2, taking into account the demagnetizing fields reduces the heliknoton stability range down to $\sim 0.3B_{\text{c}}$, which for the material parameters of FeGe corresponds to $\sim 200$ mT. This is a reasonably high range of fields accessible in most of the experimental setups.

At the film thickness below L_D, the heliknoton does not fit the size of the film and becomes unstable in both cases, with and without demagnetizing fields. With increasing thickness, the critical heliknoton fields tend to saturate in both cases. Near the thicknesses commensurate to the period of helical modulations, e.g., L/L_D = 2, the critical fields show minor picks. Similar but less pronounced pick is also seen at L/L_D = 3 for the case with demagnetizing field only. We attribute the non-monotonic behavior of the collapse field to the chiral surface twist effect [45]. This effect manifests as an additional twist of magnetization near the free edges of the sample, which can penetrate into the volume of the sample. In the case of thin films with L/L_D ≲ 3, these twists can distort the magnetic texture of the heliknoton. It is worth noting that the twist of magnetization near the free edges of the sample is also supported by the demagnetizing fields. The major difference is that the demagnetizing field does not favor a particular chirality of the surface twist, while the chiral surface twist gains energy only for the chirality favorable by the DMI. Thus, when we take into account the demagnetizing field, these two effects exhibit a cumulative effect. The most distinguishing is the behavior of the heliknoton in the thickness range between 1L_D and 1.5L_D. In this range, without demagnetizing fields, the solution remains stable only in the presence of the external field. Such anomaly in the stability of heliknoton in this case can be explained as the absence of cumulative effect from the chiral surface twist and demagnetizing fields. At realistic conditions, with demagnetizing fields, however, it is not the case, and the heliknoton is stable in the whole range of L > 1.1L_D even at zero magnetic fields. The thickness of $\sim 1.1L_{\text{D}}$ can be thought of as the lower bound for experimental observations of heliknotons.

3.3 Theoretical analysis of Lorentz TEM contrast

According to the diagram shown in Figure 2, there is no upper bound limit for the sample thickness, and even in a bulk crystal, a heliknoton can exist. However, the thickness of the sample plays a crucial role in experimental observation. In particular, the sample’s thickness is a significant factor for the applicability of Lorentz TEM. In our previous studies [13], we estimated the upper limit for the thickness of the FeGe sample to be approximately 300 nm. Above this thickness, the sample becomes no longer transparent for electrons. In the case of thicker samples beyond 300 nm, one must use a TEM setup with an acceleration voltage above 300 kV, the standard voltage for the most modern TEM instruments. Besides that, it is well known that for some magnetic textures, the Lorentz deflection forces can be completely canceled and thus give no contrast in non-tilted samples [46]; [47].

Since the heliknoton is the texture localized in all three dimensions, it is natural to expect that the contrast it produces in Lorentz TEM is thickness dependent. Figure 3 shows the Lorentz TEM contrast provided by the heliknoton in films of different thicknesses, which we calculated assuming material parameters for FeGe. For these calculations, we used a well-established method based on the phase object approximation [48]. For details of implementation, see Ref. [13].

FIGURE 3

FIGURE 3. (A) The theoretical Lorentz TEM images of heliknoton at zero external field in the plate of different thicknesses. The thicknesses in reduced units with respect to the period of helical modulations L_D are indicated on the top of each pair of over-focus and under-focus images. The defocus distance is 400 μm. (B) The over-focus and under-focus Lorentz TEM images of heliknoton in the film of thickness 2L_D calculated for different defocus distances. The scale bar in (A, B) corresponds to 70 nm.

As seen from Figure 3A, above the thickness of $\ge 2.5L_{\text{D}}$ ($\sim 170$ nm for FeGe), the characteristic features in both over-focus and under-focus images of heliknoton are hardly seen. The quality of the contrast can be improved by varying the defocus distance, as shown in Figure 3B. However, in real experiments this is not always possible, mainly due to diffraction from the crystal structure, which usually also contributes to the Lorentz TEM contrast. Based on a comparison of the theoretical Lorentz TEM images and the stability diagram, the optimal thickness range for observing heliknotons in a TEM experiment is between 1.5L_D and 2L_D (around 100 nm and 140 nm, respectively, for FeGe).

3.4 Minimum energy path calculations and HTST analysis

To estimate the stability of heliknoton at finite temperature within HTST, we use an effective atomistic spin Hamiltonian defined on a simple cubic lattice:

$E=-J\sum _{⟨i,j⟩}{\mathbf{n}}_i\cdot {\mathbf{n}}_j-\sum _{⟨i,j⟩}{\mathbf{D}}_{ij}\cdot \left[{\mathbf{n}}_i\times {\mathbf{n}}_j\right]-\mu \sum _i{\mathbf{B}}_{\text{ext}}\cdot {\mathbf{n}}_i,(7)$

were n_i is the unit vector along the magnetic moment at the lattice site i, J and $\mathbf{D}=D{\widehat{\mathbf{r}}}_{ij}$ are the Heisenberg exchange constant and Dzyaloshinskii-Moriya (DM) vector, respectively, ${\widehat{\mathbf{r}}}_{ij}$ is the unit vector along the segment connecting sites i and j, μ is the magnitude of the magnetic moment at each site. The symbol $⟨i,j⟩$ denotes the summation over the nearest neighbor pairs over all lattice sites.

Here we ignore the demagnetizing field effect and introduce characteristic parameters L_D = 2πJa/D and B_D = D²/(Jμ) for the atomistic model (7), where a is the lattice constant. The calculations for heliknoton energy barriers and lifetime provided below were performed at zero external magnetic field, B_ext = 0.

First, we apply HTST to estimate the lifetime of a particular magnetic state and the possible mechanisms of its collapse. Within the HTST, the rate of transition between states X and Y at temperature T is described by the Arrhenius law,

$k^{X\to Y}={\nu }^{X\to Y}\exp \left(-\frac{\mathrm{\Delta }{E}^{X\to Y}}{k_{\text{B}}T}\right),(8)$

where the energy barrier ΔE^X→Y can be computed knowing the minimal energy path (MEP) connecting X and Y states on the energy surface as the energy difference between the highest point along the MEP—The first-order saddle point (SP) on the energy surface of the system—And the minimum at X. The pre-exponential factor ν^X→Y incorporates dynamical ν_dyn and entropic ν_ent contributions to the transition rate

${\nu }^{X\to Y}=\frac{1}{2\pi }{\nu }_{\mathit{dyn}}{\nu }_{\mathit{ent}}.$

Within the harmonic approximation, the energy of the system in the vicinity of a stationary state is approximated by a quadratic form that allows the prefactors to be determined explicitly [49]. The quadratic forms for the minimum n^X and the saddle point n^SP are determined by the corresponding Hessian matrices ${\mathcal{H}}^X$ and ${\mathcal{H}}^{SP}$, respectively. The entropy of the states is expressed in terms of determinants of the matrices. If the energy is an invariant of some transformation, e.g., the energy of the hopfion is preserved during translations, then some eigenvalues of the Hessians are zero, and the corresponding modes are called zero modes. In the harmonic approximation, the contribution of zero-modes to the state n^X and saddle point n^SP which have Z^X and Z^SP, respectively, the zero-modes are estimated by their volumes V^X and V^SP. The entropy prefactor is given by

${\nu }_{\mathit{ent}}={\left(2\pi k_{B}T\right)}^{\frac{Z^X-Z^{SP}}{2}}\frac{V^{SP}}{V^X}\sqrt{\frac{\det {\mathcal{H}}^X}{|\det {\mathcal{H}}^{SP}|}},$

where the Hessians are restricted to subspaces consisting of non-zero modes.

The dynamical prefactor is expressed in terms of the negative eigenvalue ζ of the Hessian ${\mathcal{H}}^{\text{SP}}$ with the corresponding eigenvector e,

${\nu }_{\mathit{dyn}}=\sqrt{\frac{\mathbf{b}\cdot {\mathcal{H}}^{\text{SP}}\mathbf{b}}{|\zeta |}},b_i=\frac{\gamma \zeta }{\mu }{\mathbf{n}}_i^{\text{SP}}\times e_i.$

The index i numbers the lattice sites.

The heliknoton in bulk has at least two zero modes corresponding to translations in directions orthogonal to the helical axis. We eliminate one of the zero modes by pinning the spins at the boundaries x = const and y = const in such a way that the direction of the magnetic moments coincides with one of the helix. The pinned spins were not affected by the optimization procedures and were not included as degrees of freedom to the Hessians. On the surfaces z = const the boundary conditions are free, which allows us to simulate the escape through the boundary. Translation along the z-axis becomes a quasi-zero mode because of the effect of the boundary.

The computation of MEP was performed in Cartesian coordinates with constraints on the magnetic moments length taking into account by introduction of Lagrange multipliers, which allows us to avoid singularities and the necessity to transform coordinates while changing the map in the atlas [50]. The string method [51] with stable tangent estimate [52] was used for MEP calculation.

The most challenging part of the computation of the transition rate is the calculation of the determinants. When a standard LU or QR decomposition is used, the complexity of the problem is $\mathcal{O}(N^3)$, where N is the number of spins. Taking into account the block band structure of the Hessian matrix for close-range interactions in the absence of the demagnetization field, the complexity can be reduced to $\mathcal{O}(N^{7/3})$ [49]. Nevertheless, the problem remains time-consuming, especially for the domain of size 5L_D × 5L_D × 2L_D used in our calculations.

To find the optimal mesh density, we calculated the MEPs at different mesh densities of LD equal to.

To find the optimal mesh density, we calculated the MEPs at different mesh densities of L_D equal to 12a, 24a, 36a, and 48a. While the coarsest grid of 12a gives significantly different values of both the Heliknoton energy and the activation barrier, the results for the three densest grids (with 24a, 36a, and 48a) are only slightly different (see right inset in Figure 4A). Therefore, we used the discretization density L_D = 24a for the transition rate computation.

FIGURE 4

FIGURE 4. (A) The minimum energy paths for two distinct mechanisms for heliknoton collapse in the film: the escape through the free edges of the plate (blue curve) and the decay inside the film via the nucleation of a pair of Bloch points (red curve). The saddle points are indicated by hollow symbols. The right inset illustrates the convergence of the energies of the heliknoton, E_H, and the energy barrier for the decay, ΔE_D, to the micromagnetic limits with increasing mesh density, L_D/a. (B) The lifetime of heliknoton as a function of temperature. The escape mechanism (blue line) is less probable than the decay mechanism for all temperatures. Insets in (B) show the contributions of the exponent (activation barrier) in Arrhenius law to the lifetime. (C) Snapshots of the system along the minimum energy path for heliknoton decay. The spin texture after despiralization is represented by the isosurfaces m_z = 0. The number above each image corresponds to the reaction coordinate, r. Similar to Figure 1, the images show the volume confined by the blue box of the size 2.5L_D × 2.5L_D × 2L_D—The quarter of the whole simulated domain of the size 5L_D × 5L_D × 2L_D. (D) Snapshots of the system along the minimum energy path for heliknoton escape through the free surface. The field of view and notations are the same as in (C). In (C, D), the images in the middle correspond to the saddle points of corresponding MEPs.

The results of the HTST analysis are summarized in Figure 4. The MEPs presented in Figure 4A correspond to two distinct mechanisms for heliknoton collapse, which we refer to as decay and escape. The decay mechanism is similar to the collapse of hopfions in a model of a frustrated magnet [33]; [34]. In this case, the heliknoton collapses via the emergence of a pair of magnetic singularities, namely Bloch points of opposite signs. Representative images along the MEP corresponding to heliknoton decay are shown in Figure 4C.

As seen from the images in Figure 4D, when the heliknoton escapes through the boundary, magnetic singularities do not appear at any point along the MEP. Surprisingly, the escape mechanism for the heliknoton collapse is characterized by a higher energy barrier compared to heliknoton decay for sufficiently dense lattices. Thereby, unlike the case of hopfions in the model of frustrated magnets [32], the escape through the boundary for the heliknoton requires overcoming some activation barrier, which we attribute to the emergence of the chiral surface twist effect [53]. The latter suggests that the heliknoton in a thin plate of a chiral magnet may be more promising for experimental observations due to its finite energy barrier for escape.

The total transition rate for a metastable state is the sum of transition rates for all mechanisms of the collapse. The main contributions to the mean lifetime (inverse to the transition rate) of the heliknoton are the decay and escape mechanisms, with decay dominating:

$\tau =\frac{1}{k^{\mathit{decay}}+k^{\mathit{escape}}}\approx \frac{1}{k^{\mathit{decay}}}.$

The dependence of the heliknoton lifetime on temperature is shown in Figure 4B. The inverse to the transition rate for the decay mechanism gives an estimate of the lifetime of the heliknoton in thick films because the escape transition rate is inversely proportional to the film thickness.

Overall, our HTST analysis sheds light on the decay mechanisms of the heliknoton and highlights potential avenues for its experimental observation, especially at low temperatures. For instance, at T = 0.455J/k_B (or 0.35T_c), the lifetime for heliknoton is τ ∼ 10⁶ in dimensionless time units scaled by Jγμ⁻¹, see Figure 4B. For J = 112 meV and μ = 1030μ_B, we obtain the heliknoton lifetime of $\sim 1000$ seconds, which is quite reasonable for experimental observation. Assuming every model spin n_i corresponds to a cluster of 7³ atoms, where each atom has a magnetic moment of 3μ_B and a lattice constant of 5Å, the diameter of the observable heliknotons should be in the order of 100 nm. Note that the temperature of 0.35T_c for FeGe corresponds to $\sim 95$ K, which is accessible for the experiment performed at liquid nitrogen temperature.

3.5 Stochastic LLG dynamics

The stochastic LLG equation can be written as

$\frac{\partial {\mathbf{n}}_i}{\partial t}={\mathbf{n}}_i\times \left(\frac{1}{J}\frac{\partial E}{\partial {\mathbf{n}}_i}-{\mathbf{B}}_{\text{fluc}}^i\right)+\alpha {\mathbf{n}}_i\times \frac{\partial {\mathbf{n}}_i}{\partial t},(9)$

where t is a dimensionless time scaled by Jγμ⁻¹, with γ being the gyromagnetic ratio, α is the Gilbert damping parameter, ${\mathbf{B}}_{\text{fluc}}^i$ is the fluctuating field representing uncorrelated Gaussian white noise with a prefactor proportional to the temperature, T. For the numerical integration of Eq. 9, we use the semi-implicit method provided in Ref. [54] implemented in the publicly available software Magnoom [55]; [56]. For the chosen coupling parameters, we estimate the critical temperature, T_c ≃ 1.345J/k_B (see Ref. [57]). The simulations were performed at zero external magnetic field with α = 0.3 and a fixed time step of Δt = 0.01. We perform ten independent simulations on a domain with L_x = L_y = 4L_D, L_z = 2L_D and discretization density L_D/a = 32. We use periodic boundary conditions in the xy-plane. To increase the probability of heliknoton collapse in a reasonable time, the simulations were performed at an elevated temperature of T = 1.0J/k_B $(\sim 0.74T_{\text{c}})$. The total simulation time of each run is t = 10⁵ (10⁷ iterations).

To visualize the heliknoton isosurface n_z = 0 in the presence of strong thermal fluctuations, we first apply Fourier filtering as described in Ref. [57]. After that, we apply the despiralization procedure as introduced above, Figure 5. Among the few independent runs of the LLG simulations, we chose the two most representative examples illustrated in Figures 5A, B; Supplementary Movies S1, S2. In agreement with the results of HTST presented in the previous section, we observe two mechanisms of heliknoton collapse. In particular, Figure 5A shows the snapshots of the systems during the heliknoton collapse via the decay, accompanied by the Bloch points nucleation. Figure 5B illustrates the case when heliknoton collapses via escape through the top free surface of the plate. These are the only two mechanisms observed in stochastic LLG simulations. In the first case, the heliknoton shrinks until the tor hole completely disappears. At this moment, a pair of Bloch points of the opposite charge appears, and we see the transition from the heliknoton to the dipole string or toron [57]. The dipole string represents a non-stable state at these parameters and collapses through shrinking. In the second case, we observe the shifting of the heliknoton to the open surface as a result of its Brownian motion. The escaping happens if the distance between the heliknoton and the surface is smaller than the critical one. It is worth noting that despite the temporal fluctuations, the escape through the free surface represents a smooth transition. Remarkably, the heliknoton remains stable near the surface of the plate for a noticeably long time. That might indicate the existence of a metastable state resembling surface modulations. Details of the investigation of such localized states will be provided elsewhere.

FIGURE 5

FIGURE 5. Snapshots of the system in two different stochastic LLG simulations are shown. The set of images in row (A) corresponds to the heliknoton collapse through the nucleation of Bloch points, while images in row (B) correspond to its escaping through the free surface of the plate. The shown box has size 2L_D in each dimension.

Among ten independent simulations, we observed the collapse of the heliknoton in the times range between 2 × 10⁴ to 8 × 10⁴ in reduced time units. Using these results, we estimated the lifetime of the heliknoton at T = 1.0J/k_B to be approximately τ ∼ 10⁴. This result agrees with the findings of the HTST calculations presented in Figure 4B, which suggest a lifetime of τ ∼ 10³ for this temperature. The limited sampling (only ten independent runs) may account for the observed discrepancy. However, given the variation of the heliknoton’s lifetime by ten orders of magnitude for different temperatures, we consider the agreement between the simulations and HTST calculations to be satisfactory.

4 Conclusion

In this study, we investigated the stability and decay mechanisms of the heliknoton in realistic conditions using direct energy minimization, stochastic LLG simulations, and HTST calculations. Our results showed that the heliknoton is stable in a wide range of magnetic fields and plate thicknesses, as demonstrated by the stability diagram calculated at zero temperature. However, our calculations indicated that for plate thicknesses above 2.5L_D, the magnetic contrast in Lorentz TEM becomes too weak for reliable observations.

Based on our findings, we estimate the optimal plate thickness for heliknoton observations to be 2L_D ± 0.5L_D (140 ± 35 nm for FeGe parameters). Our analysis using the GNEB method revealed two main mechanisms of heliknoton decay: collapse via the formation of Bloch points and escape through the plate surface. The LLG simulations at finite temperatures supported these observations and showed reasonable agreement with the estimation of heliknoton lifetime, at least at elevated temperatures. Moreover, the HTST calculations allowed us to estimate heliknoton lifetimes at different temperatures. Overall, our results suggest that the heliknoton should be stable with a reasonably long lifetime at the standard measurement temperature of liquid nitrogen (T = 95 K) used in TEM experiments. We anticipate that these findings will stimulate further experimental observations of heliknotons in magnetic systems.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author contributions

VK, NK, and FR performed micromagnetic and LLG simulations, IL performed minimum energy path and lifetime calculations. NK wrote the first draft of the manuscript. SB and VU participated in the interpretation of the results, review and editing, supervision. All authors contributed to the manuscript writing. All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication. All authors contributed to the article and approved the submitted version.

Funding

VK acknowledges financial support from the Icelandic Research Fund (Grant No. 217750). This project has received funding from the European Research Council under the European Union’s Horizon 2020 Research and Innovation Programme (Grant No. 856538—Project “3D MAGiC”). NK and SB acknowledge financial support from the Deutsche Forschungsgemeinschaft through SPP 2137 “Skyrmionics”, Grant No. KI 2078/1-1 and Grant No. BL 444/16, respectively. FR acknowledges support from the Swedish Research Council.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2023.1201018/full#supplementary-material

References

1. Yu XZ, Onose Y, Kanazawa N, Park JH, Han JH, Matsui Y, et al. Real-space observation of a two-dimensional skyrmion crystal. Nature (2010) 465:901–4. doi:10.1038/nature09124

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Park HS, Yu X, Aizawa S, Tanigaki T, Akashi T, Takahashi Y, et al. Observation of the magnetic flux and three-dimensional structure of skyrmion lattices by electron holography. Nat Nanotechnology (2014) 9:337–42. doi:10.1038/nnano.2014.52

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Yu XZ, Kanazawa N, Onose Y, Kimoto K, Zhang WZ, Ishiwata S, et al. Near room-temperature formation of a skyrmion crystal in thin-films of the helimagnet FeGe. Nat Mater (2010) 10:106–9. doi:10.1038/nmat2916

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Kovács A, Caron J, Savchenko AS, Kiselev NS, Shibata K, Li Z-A, et al. Mapping the magnetization fine structure of a lattice of bloch-type skyrmions in an FeGe thin film. Appl Phys Lett (2017) 111:192410. doi:10.1063/1.5004394

CrossRef Full Text | Google Scholar

5. Du H, Zhao X, Rybakov FN, Borisov AB, Wang S, Tang J, et al. Interaction of individual skyrmions in a nanostructured cubic chiral magnet. Phys Rev Lett (2018) 120:197203. doi:10.1103/PhysRevLett.120.197203

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Yu X, Morikawa D, Yokouchi T, Shibata K, Kanazawa N, Kagawa F, et al. Aggregation and collapse dynamics of skyrmions in a non-equilibrium state. Nat Phys (2018) 14:832–6. doi:10.1038/s41567-018-0155-3

CrossRef Full Text | Google Scholar

7. Yu X, Kikkawa A, Morikawa D, Shibata K, Tokunaga Y, Taguchi Y, et al. Variation of skyrmion forms and their stability in mnsi thin plates. Phys Rev B (2015) 91:054411. doi:10.1103/PhysRevB.91.054411

CrossRef Full Text | Google Scholar

8. Shibata K, Yu XZ, Hara T, Morikawa D, Kanazawa N, Kimoto K, et al. Towards control of the size and helicity of skyrmions in helimagnetic alloys by spin–orbit coupling. Nat Nanotechnology (2013) 8:723–8. doi:10.1038/nnano.2013.174

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Dzyaloshinsky I A thermodynamic theory of “weak” ferromagnetism of antiferromagnetics. J Phys Chem Sol (1958) 4:241–55. doi:10.1016/0022-3697(58)90076-3

CrossRef Full Text | Google Scholar

10. Moriya T Anisotropic superexchange interaction and weak ferromagnetism. Phys Rev (1960) 120:91–8. doi:10.1103/PhysRev.120.91

CrossRef Full Text | Google Scholar

11. Bogdanov AN, Yablonskii DA Thermodynamically stable “vortices” in magnetically ordered crystals. the mixed state of magnets. Sov Phys JETP (1989) 68.

Google Scholar

12. Tokura Y, Kanazawa N Magnetic skyrmion materials. Chem Rev (2020) 121:2857–97. doi:10.1021/acs.chemrev.0c00297

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Zheng F, Rybakov FN, Kiselev NS, Song D, Kovács A, Du H, et al. Magnetic skyrmion braids. Nat Commun (2021) 12:5316. doi:10.1038/s41467-021-25389-7

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Zheng F, Kiselev NS, Yang L, Kuchkin VM, Rybakov FN, Blügel S, et al. Skyrmion–antiskyrmion pair creation and annihilation in a cubic chiral magnet. Nat Phys (2022) 18:863–8. doi:10.1038/s41567-022-01638-4

CrossRef Full Text | Google Scholar

15. Rybakov FN, Kiselev NS Chiral magnetic skyrmions with arbitrary topological charge. Phys Rev B (2019) 99:064437. doi:10.1103/PhysRevB.99.064437

CrossRef Full Text | Google Scholar

16. Foster D, Kind C, Ackerman PJ, Tai J-SB, Dennis MR, Smalyukh II Two-dimensional skyrmion bags in liquid crystals and ferromagnets. Nat Phys (2019) 15:655–9. doi:10.1038/s41567-019-0476-x

CrossRef Full Text | Google Scholar

17. Kuchkin VM, Barton-Singer B, Rybakov FN, Blügel S, Schroers BJ, Kiselev NS Magnetic skyrmions, chiral kinks, and holomorphic functions. Phys Rev B (2020) 102:144422. doi:10.1103/PhysRevB.102.144422

CrossRef Full Text | Google Scholar

18. Tang J, Wu Y, Wang W, Kong L, Lv B, Wei W, et al. Magnetic skyrmion bundles and their current-driven dynamics. Nat Nanotechnology (2021) 16:1086–91. doi:10.1038/s41565-021-00954-9

CrossRef Full Text | Google Scholar

19. Dubrovin BA, Fomenko AT, Novikov SP Modern geometry – methods and applications. Part II. The geometry and topology of manifolds. New York: Springer-Verlag (1985). doi:10.1007/978-1-4612-1100-6

CrossRef Full Text | Google Scholar

20. Aminov JA. Sources of curvature of a vector field. Mathematics USSR-Sbornik (1969) 9:199–211. doi:10.1070/sm1969v009n02abeh001128

CrossRef Full Text | Google Scholar

21. Aminov Y The geometry of vector fields. Amsterdam: Gordon & Breach (2000).

Google Scholar

22. Rybakov P Topological excitations in field theory models of superconductivity and magnetism. Ph.D. thesis. Stockholm: KTH Royal Institute of Technology (2021).

Google Scholar

23. Zheng F, Kiselev NS, Rybakov FN, Yang L, Blügel S, Dunin-Borkowski RE Hopfion rings in a cubic chiral magnet. Preprint Res Square (2023). doi:10.21203/rs.3.rs-2681064/v1

CrossRef Full Text | Google Scholar

24. Hopf H Über die abbildungen der dreidimensionalen sphäre auf die kugelfläche. Mathematische Annalen (1931) 104:637–65. doi:10.1007/bf01457962

CrossRef Full Text | Google Scholar

25. Whitehead JHC An expression of hopf's invariant as an integral. Proc Natl Acad Sci (1947) 33:117–23. doi:10.1073/pnas.33.5.117

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Dzyloshinskii IE, Ivanov BA Localized topological solitons in a ferromagnet. JETP Lett (1979) 29:540.

Google Scholar

27. Bogolubsky I Three-dimensional topological solitons in the lattice model of a magnet with competing interactions. Phys Lett A (1988) 126:511–4. doi:10.1016/0375-9601(88)90049-7

CrossRef Full Text | Google Scholar

28. Cooper NR Propagating magnetic vortex rings in ferromagnets. Phys Rev Lett (1999) 82:1554–7. doi:10.1103/PhysRevLett.82.1554

CrossRef Full Text | Google Scholar

29. Sutcliffe P Skyrmion knots in frustrated magnets. Phys Rev Lett (2017) 118:247203. doi:10.1103/physrevlett.118.247203

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Liu Y, Lake RK, Zang J Binding a hopfion in a chiral magnet nanodisk. Phys Rev B (2018) 98:174437. doi:10.1103/PhysRevB.98.174437

CrossRef Full Text | Google Scholar

31. Voinescu R, Tai J-SB, Smalyukh II Hopf solitons in helical and conical backgrounds of chiral magnetic solids. Phys Rev Lett (2020) 125:057201. doi:10.1103/PhysRevLett.125.057201

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Rybakov FN, Kiselev NS, Borisov AB, Döring L, Melcher C, Blügel S Magnetic hopfions in solids. APL Mater (2022) 10:111113. doi:10.1063/5.0099942

CrossRef Full Text | Google Scholar

33. Sallermann M, Jónsson H, Blügel S Stability of hopfions in bulk magnets with competing exchange interactions. Phys Rev B (2023) 107:104404. doi:10.1103/PhysRevB.107.104404

CrossRef Full Text | Google Scholar

34. Lobanov IS, Uzdin VM Lifetime, collapse, and escape paths for hopfions in bulk magnets with competing exchange interactions. Phys Rev B (2023) 107:104405. doi:10.1103/PhysRevB.107.104405

CrossRef Full Text | Google Scholar

35. Kent N, Reynolds N, Raftrey D, Campbell ITG, Virasawmy S, Dhuey S, et al. Creation and observation of hopfions in magnetic multilayer systems. Nat Commun (2021) 12:1562. doi:10.1038/s41467-021-21846-5

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Tai J-SB, Smalyukh II Three-dimensional crystals of adaptive knots. Science (2019) 365:1449–53. doi:10.1126/science.aay1638

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Bessarab PF, Uzdin VM, Jónsson H Harmonic transition-state theory of thermal spin transitions. Phys Rev B (2012) 85:184409. doi:10.1103/PhysRevB.85.184409

CrossRef Full Text | Google Scholar

38. Bessarab PF, Uzdin VM, Jónsson H Potential energy surfaces and rates of spin transitions. Z für Physikalische Chem (2013) 130708000310008. doi:10.1524/zpch.2013.0403

CrossRef Full Text | Google Scholar

39.D ter Haar, editor. Collected papers of L.D. Landau (pergamon). Chap. 18 - on the theory of the dispersion of magnetic permeability in ferromagnetic bodies (1965). p. 101–14. doi:10.1016/B978-0-08-010586-4.50023-7

CrossRef Full Text | Google Scholar

40. Di Fratta G, Muratov CB, Rybakov FN, Slastikov VV Variational principles of micromagnetics revisited. SIAM J Math Anal (2020) 52:3580–99. doi:10.1137/19M1261365

CrossRef Full Text | Google Scholar

41. Rybakov FN, Babaev E Excalibur software (2020). URL: http://quantumandclassical.com/excalibur/.

Google Scholar

42. Vansteenkiste A, Leliaert J, Dvornik M, Helsen M, Garcia-Sanchez F, Van Waeyenberge B The design and verification of mumax3. AIP Adv (2014) 4:107133. doi:10.1063/1.4899186

CrossRef Full Text | Google Scholar

43. Zheng F, Rybakov FN, Borisov AB, Song D, Wang S, Li Z-A, et al. Experimental observation of chiral magnetic bobbers in b20-type FeGe. Nat Nanotechnology (2018) 13:451–5. doi:10.1038/s41565-018-0093-3

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Skyrme THR A non-linear field theory. Proc R Soc Lond Ser A. Math Phys Sci (1961) 260:127–38. doi:10.1098/rspa.1961.0018

CrossRef Full Text | Google Scholar

45. Rybakov FN, Borisov AB, Blügel S, Kiselev NS New spiral state and skyrmion lattice in 3D model of chiral magnets. New J Phys. IOP Publishing (2016) 18 (4):045002. doi:10.1088/1367-2630/18/4/045002

CrossRef Full Text | Google Scholar

46. Denneulin T, Caron J, Hoffmann M, Lin M, Tan H, Kovács A, et al. Off-axis electron holography of néel-type skyrmions in multilayers of heavy metals and ferromagnets. Ultramicroscopy (2021) 220:113155. doi:10.1016/j.ultramic.2020.113155

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Denneulin T, Caron J, Müller-Caspary K, Boulle O, Kovács A, Dunin-Borkowski RE Visibility and apparent size of néel-type magnetic skyrmions in fresnel defocus images of multilayer films. Microsc Microanalysis (2021) 27:1356–65. doi:10.1017/s1431927621012927

CrossRef Full Text | Google Scholar

48. Graef MD 2. lorentz microscopy: Theoretical basis and image simulations. In: Experimental methods in the physical sciences. Elsevier (2001). p. 27–67. doi:10.1016/s1079-4042(01)80036-9

CrossRef Full Text | Google Scholar

49. Lobanov IS, Uzdin VM The lifetime of micron scale topological chiral magnetic states with atomic resolution. Comp Phys Commun (2021) 269:108136. doi:10.1016/j.cpc.2021.108136

CrossRef Full Text | Google Scholar

50. Lobanov IS, Potkina MN, Uzdin VM Stability and lifetimes of magnetic states of nano-and microstructures (brief review). JETP Lett (2021) 113:801–13. doi:10.1134/S0021364021120109

CrossRef Full Text | Google Scholar

51. E W, Ren W, Vanden-Eijnden E Simplified and improved string method for computing the minimum energy paths in barrier crossing events. J Chem Phys (2007) 126:164103. doi:10.1063/1.2720838

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Henkelman G, Jónsson H Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. J Chem Phys (2000) 113:9978–85. doi:10.1063/1.1323224

CrossRef Full Text | Google Scholar

53. Rybakov FN, Borisov AB, Blügel S, Kiselev NS New type of stable particlelike states in chiral magnets. Phys Rev Lett (2015) 115:117201. doi:10.1103/PhysRevLett.115.117201

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Mentink JH, Tretyakov MV, Fasolino A, Katsnelson MI, Rasing T Stable and fast semi-implicit integration of the stochastic landau–lifshitz equation. J Phys Condensed Matter (2010) 22:176001. doi:10.1088/0953-8984/22/17/176001

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Kiselev NS Magnoom software (2016). URL: https://github.com/n-s-kiselev/magnoom.

Google Scholar

56. Savchenko AS, Kuchkin VM, Rybakov FN, Blügel S, Kiselev NS Chiral standing spin waves in skyrmion lattice. APL Mater (2022) 10:071111. doi:10.1063/5.0097651

CrossRef Full Text | Google Scholar

57. Müller GP, Rybakov FN, Jónsson H, Blügel S, Kiselev NS Coupled quasimonopoles in chiral magnets. Phys Rev B (2020) 101:184405. doi:10.1103/physrevb.101.184405

CrossRef Full Text | Google Scholar