INTRODUCTION

Among the various structures in magnetic media investigated in the last decade, there has been a particular interest in vortex structures. The important role of vortices in the description of magnetic topological phase transitions has been noted by many researchers (see, e.g., [1–5]). Such structures are not only of academic interest. Vortices in chiral magnets (chiral skyrmions), predicted two decades ago [6, 7], can find important applications in running memory [8–10]. To date, two-dimensional vortices in ferro- and antiferromagnets have been studied in sufficient detail [6, 7, 11–15]. Note that, despite the years-long study of vortices in “large” (bulk) samples, three-dimensional vortex structures have been insufficiently studied theoretically and experimentally. However, a new experimental technique: off-axis electron holography, has recently appeared [16], which makes it possible to study three-dimensional localized structures. It was revealed that the two-dimensional chiral skyrmions predicted in [5, 6] have a three-dimensional structure [17]. Of the three-dimensional structures, only magnonic drops [10], magnetic “hedgehogs”, and spiral structures [18] have been theoretically described to date. Therefore, the study of three-dimensional structures by analytical methods is of great importance.

Such methods are generally applicable to models with a high degree of symmetry. Although the applicability of analytical methods is limited, their scientific value is not in doubt. They allow one to thoroughly investigate the structure of the nuclei of nonlinear structures and qualitatively take into account the influence of other interactions not included in the initial model. In addition, exact solutions are “seed” functions in computer simulations of energy density minimization, including magnetic fields, anisotropy energy, etc. Vortex structures in a ferromagnet are formed mainly by exchange interaction (Heisenberg model), and the objective of this work is to find three-dimensional vortices in this model.

We discuss new types of three-dimensional vortex structures in the Heisenberg model in this paper. The article is planned as follows. In the first section, in the equations of an isotropic magnet, we use a substitution that reduces the model to the pendulum equation and four equations with a simple geometric interpretation. The latter are reduced to two equations for a complex function \(S\left( {x,y,z} \right).\) We find the general solution of these equations, which depends on an arbitrary function. The simplest solution of these equations is discussed in Section 2. It is shown that it describes a new magnetic structure, which is comprised of two intersecting rectilinear vortex filaments, which change the topological charge after the intersection. At the end of the section, the experimental implementation of the structures found is discussed.

1 SIMPLIFYING THE EQUATIONS OF THE HEISENBERG MODEL

The Heisenberg model is based on using the Hamiltonian H of the exchange interaction,

$$H = \sum\limits_{i{\text{ > }}j} {{{J}_{{ij}}}{\mathbf{S}_{i}}{\mathbf{S}_{j}}} ,$$

where the summation is performed over all pairs {i, j} of different sites of the crystal accomodating ions with spins Si and Sj and the constants Jij characterize the exchange interaction between these ions. Taking into account the isotropic exchange interaction between nearest neighbors in the continuum limit, the density of the Hamiltonian of a ferromagnet has a simple form:

$$H = J~\nabla \mathbf{S}~\nabla \mathbf{S}.$$
(1)

For the ferromagnetic Heisenberg model, J > 0, and, for the antiferromagnetic model, J < 0 [19]. Note that decomposition (1) takes into account only quadratic terms in the spin gradients. In the general case, to describe even the exchange interaction in some systems, it is necessary to take into account in (1) terms of a higher order in \(\nabla S{\text{,}}\) including the biquadratic terms for the spin \(S \geqslant 1.\) In addition, we further investigate the case of zero magnetic field, when the ground state of the Heisenberg model is a uniformly magnetized ferromagnet without any structures.

Then, the three-dimensional Landau–Lifshitz equations can be written in the form

$$\left[ {\mathbf{n} \times \left( {\partial _{x}^{2} + \partial _{y}^{2} + \partial _{z}^{2}~} \right)\mathbf{n}} \right] = 0,\,\,\,\,{\mathbf{n}^{2}} = 1,$$
(2)

where n is the unit vector of magnetization S = M0n and M0 is the spontaneous magnetization. Hereinafter, \({{\Delta }} = \partial _{x}^{2} + \partial _{y}^{2} + \partial _{z}^{2}\) is the three-dimensional Laplace operator. The vector n is parameterized by the fields \({{\Theta }}\) and \({{\Phi :}}\)

$$\mathbf{n} = \left( {\cos {{\Phi }}\sin {{\Theta }},\sin {{\Phi }}\sin {{\Theta }},\cos {{\Theta }}} \right).$$

In these variables, Eq. (2) turns into a system of nonlinear differential equations:

$${{\Delta \Theta }} = \frac{1}{2}{{\left( {\nabla {{\Phi }}} \right)}^{2}}\sin 2{{\Theta ,}}\,\,\,\,\nabla \left[ {\left( {\nabla {{\Phi }}} \right){{{\sin }}^{2}}{{\Theta }}~} \right] = 0,$$
(3)

Equation (2) is invariant with respect to the group \(SO\left( 3 \right) \times SO\left( {3n} \right)\) of spin and spatial rotations. This symmetry makes it possible to find a wide class of exact solutions. The analytical solution of Eq. (2) is possible only in certain classes of solutions. To single out one of them, it is necessary to generalize the procedure proposed in [18] and assume that the field \({{\Theta }}\) locally depends on the auxiliary field: \({{\Theta }} = {{\Theta }}\left( {a\left[ {x,y,z} \right]} \right).\) Then, by direct calculations, it is easy to verify that the equations

$${{\Theta }}{\kern 1pt} {''}\left( a \right) = \frac{1}{2}{\text{sin}}2{{\Theta }}\left( a \right),$$
(4)
$${{\Delta }}a = {{\Delta \Phi }} = 0,\,\,\,{{\left( {\nabla a} \right)}^{2}} = {{\left( {\nabla {{\Phi }}} \right)}^{2}},\,\,\,\nabla a\nabla {{\Phi }} = 0$$
(5)

imply Eqs. (3). Indeed, after the substitution \({{\Theta }} = {{\Theta }}\left( {a\left[ {x,y,z} \right]} \right)\), the first equation in (3) becomes

$${{\Theta }}{\kern 1pt} {''}\left( a \right){{\left( {\nabla a} \right)}^{2}} - \frac{1}{2}{\text{sin}}2{{\Theta }}\left( a \right){{\left( {\nabla {{\Phi }}} \right)}^{2}} + {{\Theta }}{\kern 1pt} '\left( a \right){{\Delta }}a = 0,$$

which implies Eq. (4). The second equation in (3) after the substitution reduces to

$$2\cos {{\Theta }}~{{\Theta }}{\kern 1pt} '\left( a \right)\left( {\nabla a~\nabla {{\Phi }}} \right) + \sin {{\Theta }}~{{\Delta \Phi }} = 0.$$

As we will see below, such substitutions lead to a wide class of exact solutions of the nonintegrable model (2).

Note that Eqs. (5) have a simple geometric interpretation. They describe two surfaces: \(a\left( {x,y,z} \right) = 0\) and \({{\Phi }}\left( {x,y,z} \right) = 0,\) intersecting at each point, the normals to which, \(\nabla a\) and \(\nabla {{\Phi }}\), respectively, have unit length and are orthogonal and divergence-free:

$$\nabla \nabla a = 0,\,\,\,\,\nabla \nabla {{\Phi }} = 0.$$

Let us turn to solving Eqs. (4) and (5). We introduce a complex field \(S = a + i{{\Phi }}\) and write (5) as a system of two equations for the field S:

$$\left( {\nabla S} \right)\left( {\nabla S} \right) = 0;$$
(6)
$${{\Delta }}S = 0.$$
(7)

This system has a remarkable property of invariance with respect to arbitrary changes of the field S, which we will use below. It is easy to verify that, if the field \(S\) is a solution of this system, then an arbitrary function \(\tilde {S} = F\left( S \right)\) is also its solution. To solve system (6) and (7), we use a procedure proposed in [18]. We introduce a new complex field \(T\left( {x,y,z} \right)\) as

$${{S}_{{,y}}}\left( {x,y,z} \right) = T\left( {x,y,z} \right){{S}_{{,z}}}\left( {x,y,z} \right).$$
(8)

Here, the subscript denotes differentiation with respect to the corresponding variable:

$${{S}_{{,y}}}\left( {x,y,z} \right) = \frac{{\partial S\left( {x,y,z} \right)}}{{\partial y}}$$

etc. The substitution of (8) into (6) determines \({{S}_{{,x}}}\left( {x,y,z} \right){\text{:}}\)

$${{S}_{{,x}}}\left( {x,y,z} \right) = {\text{i}}\sqrt {1 + {{T}^{2}}\left( {x,y,z} \right)} {{S}_{{,z}}}\left( {x,y,z} \right).$$
(9)

The consistency condition for (8) and (9) gives a nonlinear equation for the field T:

$${{T}_{{,z}}} + T{{T}_{{,y}}} + {\text{i}}\sqrt {1 + {{T}^{2}}} {{T}_{{,x}}} = 0.$$

From the theory of first-order partial differential equations, it follows that the field \(T\) is determined by the implicit equation

$$G\left[ {{{H}_{1}},{{H}_{2}},{{H}_{3}}} \right] = 0,$$
(10)

where \({{H}_{1}},\) \({{H}_{2}},\) and \({{H}_{3}}\) are integrals of the characteristic system of equations for the coordinates \(x\left( t \right),\) \(y\left( t \right),\) \(z\left( t \right){\text{:}}\)

$$\begin{gathered} \frac{{\text{d}}}{{{\text{d}}t}}z\left[ t \right] = 1,\,\,\,\,\frac{{\text{d}}}{{{\text{d}}t}}y\left[ t \right] = T\left[ {x,y,z} \right], \\ \frac{{\text{d}}}{{{\text{d}}t}}x\left[ t \right] = {\text{i}}\sqrt {1 + {{T}^{2}}\left[ {x,y,z} \right]} ,\,\,\,\,\frac{{\text{d}}}{{{\text{d}}t}}T\left[ {x,y,z} \right] = 0. \\ \end{gathered} $$

The integrals have the form

$$\begin{gathered} {{H}_{1}} = T\left[ {x,y,z} \right],\,\,\,\,{{H}_{2}} = T\left[ {x,y,z} \right]z - y, \\ {{H}_{3}} = - {\text{i}}~\sqrt {1 + {{T}^{2}}\left[ {x,y,z} \right]} z + x. \\ \end{gathered} $$
(11)

From this wide class of solutions (10) and (11), we choose one with the function \(G\) in the form

$$G = {\text{i}}{{{\text{H}}}_{1}}{{H}_{2}} + {{H}_{3}}\sqrt {1 + H_{1}^{2}.} $$

Then, the field T satisfies the equation

$$ - {\text{i}}\left( {z + yT} \right) + x\sqrt {1 + {{T}^{2}}} = 0$$

and, therefore,

$$T = \frac{{ - yz + {\text{i}}x\sqrt {{{x}^{2}} + {{y}^{2}} + {{z}^{2}}} }}{{{{x}^{2}} + {{y}^{2}}}}.$$
(12)

It turns out that the field (12) satisfies not only Eq. (6), but also Eq. (7) (with the substitution \(S \to T\)). Then, according to the partition,

$$a + {\text{i}\Phi } = Q\ln T.$$

With a parameter Q taking integer values, we find the explicit form of the fields a and \({{\Phi :}}\)

$$\begin{gathered} {{\Phi }} = Q~{\text{arctan}}\left[ {\frac{{x\sqrt {{{x}^{2}} + {{y}^{2}} + {{z}^{2}}} }}{{ - yz}}} \right]; \\ a = Q\ln \frac{{\sqrt {{{x}^{2}} + {{z}^{2}}~} }}{{\sqrt {{{x}^{2}} + {{y}^{2}}} }}~, \\ \end{gathered} $$
(13)

which, as is easy to verify by direct calculations, satisfy Eqs. (5).

Let us choose as \({{\Theta }}\left( {a\left( {x,y,z} \right)} \right)\) a solution of Eq. (4) in the form of a soliton lattice:

$$\cos {{\Theta }} = {\text{sn}}~\left[ {\frac{a}{k},k} \right],\,\,\,\,0 < k < 1.$$
(14)

For \(k = 1\), expression (14) is simplified to

$${{\Theta }} = 2~{\text{arctan}}~\left[ {\exp ( - a)} \right].$$
(15)

Further on, we analyze formulas (13) and (15) of the new structure in an isotropic magnet.

2 ANALYSIS OF THE VORTEX STRUCTURE

Note first that substitution (4) and (5) and system (6) and (7) are a generalization of the instanton theory [20] in the two-dimensional case for an isotropic magnet. A wide class of solutions found in [20] is described by the simple formula

$$w = {\text{cot}}~\frac{{{\Theta }}}{2}\exp {\text{i}\Phi } = F\left[ z \right],\,\,\,z = x + {\text{i}}y$$

with an analytic function F. In the two-dimensional case, Eqs. (6) and (7) are identically satisfied with S = S(z), and the relationship with the field w, according to (15), is established by the simple formula

$$w = {\text{exp}}\left( S \right).$$

Despite the apparent simplicity, the expression for the field \({{\Phi }}\) (13) has a rich three-dimensional vortex structure. Previously investigated [6, 7, 11, 12] two-dimensional vortices comprised a straight filament with dependences \({{\Phi }} = {{\Phi }}\left( \varphi \right)\) and \({{\Theta }} = {{\Theta }}\left( r \right)\) in a polar coordinate system \(\left( {r,\varphi } \right).\) Formula (13) describes two intersecting vortex filaments. One of them is a rectilinear vortex filament located on the axis \(Oz\) and, in the polar coordinate system (\(x = r\cos \varphi ,\) \(y = ~r\sin \varphi \)) in the plane \(z = {\text{const}}\), has the form

$${{\Phi }}\left( {r,z,\varphi } \right) = Q~{\text{arctan}}\left( {\frac{{\sqrt {{{r}^{2}} + {{z}^{2}}} }}{{ - z}}{\text{cot}}~\varphi } \right)$$
(16)

with a topolgical charge

$$q = \frac{1}{{2\pi }}\oint\limits_{{\gamma }} {\nabla {{\Phi d}}{\mathbf{r}}} $$

equal to \(Q~{\text{sgn}}\left( z \right).\) Here, \(\gamma \) is an arbitrary closed contour in the plane \(z = {\text{const,}}\) enclosing the vortex center. Integration is performed counterclockwise. Formula (16) has the simplest form near the vortex center (\(r = 0\)):

$${{\Phi }}\left( {r,z,\varphi } \right) \to {\text{sgn}}\left( z \right)~Q\left( {\varphi - \frac{\pi }{2}} \right)$$

and describes a vortex for \(z > 0\) and an antivortex for \(z < 0\), with topological charges \( \pm Q\), respectively. From the boundary conditions

$$\begin{gathered} {{\Phi }} \to \pm \frac{{Q\pi }}{2}\,\,\,\,\left( {x \to \pm \infty } \right), \\ {{\Phi }} \to \pm \frac{Q}{2}{\text{arctan}}\left( {\frac{x}{{ - z}}} \right)\,\,\,\,\left( {y \to \pm \infty } \right) \\ \end{gathered} $$

and a numerical experiment, it follows that (13) describes two domain walls with nonexponential behavior at infinity. The centers of these domain walls (the locus of the points with the maximum value of the derivative) coincides, respectively, with the axes \(x = 0\) and \(y = 0.\) Such vortex formation by the intersection of domain walls was studied in detail for the sine-Gordon model in [21, 22].

The vortex structure is more evident in the distribution of the vector field

$$\mathbf{Z} = \left( {\cos {{\Phi }},\sin {{\Phi }}} \right).$$

This field in the planes \(z = - 100\) and \(z = 100\) is shown in Figs. 1 and 2, respectively. It can be seen from Fig. 1 that, at a constant \(y\), the vector Z rotates from \({{ - \pi } \mathord{\left/ {\vphantom {{ - \pi } 2}} \right. \kern-0em} 2}\) to \({\pi \mathord{\left/ {\vphantom {\pi 2}} \right. \kern-0em} 2}\) counterclockwise if y < 0 and clockwise if y > 0; therefore, after crossing the domain boundaries, they change sign. Expression (13) is not defined at \(x = 0,\) \(z = 0;\) therefore, in the plane \(y = {\text{const}}\), it describes the second straight vortex filament. In the coordinate system \(z = {{r}_{1}}\sin \varphi ,\) \(x = {{r}_{1}}\cos \varphi \), the vortex structure described by the expression

$${{\Phi }}\left( {{{r}_{1}},y,\varphi } \right) = - Q~{\text{arccot}}\frac{{{\text{tan}}~\varphi y}}{{\sqrt {r_{1}^{2} + {{y}^{2}}} }}$$

is a vortex (if \(y > 0\)) with a topological charge Q and an antivortex (if \(y < 0\)) with the same charge, which is also formed by the intersection of domain boundaries.

Fig. 1.
figure 1

Distribution of vector Z in the plane \(z = - 100.\)

Full size image
Fig. 2.
figure 2

Distribution of vector Z in the plane \(z = 100.\)

Full size image

The vortex structure is shown in Figs. 3 and 4, in which we clearly see jumps of the field \({{\Phi }}\) by \( - 2\pi \) (Fig. 3, \(y = - 100\)) and \(2\pi \) (Fig. 4, \(y = 100\)).

Fig. 3.
figure 3

Vortex structure of the field \({{\Phi }}\) in the plane \(xOz\) at \(y = - 100.\)

Full size image
Fig. 4.
figure 4

Vortex structure of the field \({{\Phi }}\) in the plane \(xOz\) at \(y = 100.\)

Full size image

The structure of the field \({{\Theta }}\) at \(k = 1\) in the standard polar system,

$${{\Theta }} = 2{\text{arccot}}{{\left( {\frac{{\sqrt {{{z}^{2}} + {{r}^{2}}{{{\cos }}^{2}}\varphi } }}{r}} \right)}^{Q}}$$

depends on three spatial variables.

Although, at \(z \to \pm \infty \), the field \({{\Theta }}\) tends to the ground state (\({{\Theta }} \to 0\)), at \(r \to \infty \), the azimuthal angle depends on the polar angle \(\varphi {\text{:}}\)

$${{\Theta }} \to 2{\text{arcot}}\left[ {{{{\left| {\sec \varphi } \right|}}^{Q}}} \right]\,\,\,\,{\text{at}}\,\,\,\,r \to \infty .$$

Therefore, the field n (2) has the following asymptotics: it is not singular at the origin and

$$\begin{gathered} \mathbf{n} \to \left( {0,0,Q} \right)\,\,\,\,{\text{at}}\,\,\,\,z \to \pm \infty , \\ n \to \left( {0,Q\frac{{2\cos \varphi }}{{1 + {{{\cos }}^{2}}\varphi }}, - Q\frac{{{{{\sin }}^{2}}\varphi }}{{1 + {{{\cos }}^{2}}\varphi }}} \right)\,\,\,\,{\text{at}}\,\,\,\,r \to \infty , \\ \mathbf{n} \to \left( {0,0,Q} \right)\,\,\,\,{\text{at}}\,\,\,\,r \to 0. \\ \end{gathered} $$
(17)

The energy density of the Heisenberg model,

$$e = \frac{1}{2}\left[ {{{{\left( {\nabla {{\Theta }}} \right)}}^{2}} + {{{\left( {\nabla {{\Phi }}} \right)}}^{2}}{{{\sin }}^{2}}{{\Theta }}} \right]JM_{0}^{2}$$

under constraints (4) and (5) has a simple form

$$e = {\text{cos}}{{{\text{h}}}^{{ - 2}}}a{{\left( {\nabla ~a} \right)}^{2}}~JM_{0}^{2}$$

and the total energy \(E\) is proportional to the size of the system, L:

$$E = \int {e~{{{\text{d}}}^{3}}\mathbf{r}} = 8\pi LJM_{0}^{2}.$$

This energy is significantly lower than the energy of a two-dimensional vortex [11, 12] in an easy-plane ferromagnet, which is proportional to \({\text{ln}}~\left[ {{L \mathord{\left/ {\vphantom {L {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}}} \right]\) (\({{r}_{0}}\) is the order of the lattice constant) per atomic layer.

The vortex structures predicted can nucleate in a cylinder with surface anisotropy, by thermal fluctuations, or an alternating magnetic field. In this case, the magnetization on the lateral surfaces takes the form (17).