Stability of Biskyrmions in Centrosymmetric Magnetic Films

TL;DR

This paper analyzes the stability of biskyrmions in centrosymmetric magnetic films, considering various interactions and lattice effects, and identifies conditions under which they can be stabilized by magnetic fields.

Contribution

It provides an analytical and numerical study of biskyrmion stability, including the effects of film thickness, interactions, and lattice discreteness, revealing conditions for their stabilization.

Findings

Biskyrmions with separation d of order λ can be stabilized by magnetic fields.

Shape of biskyrmions approximates Belavin-Polyakov solutions when d is small.

Conditions for experimental creation of biskyrmions are discussed.

Abstract

Motivated by the observation of biskyrmions in centrosymmetric magnetic films (Yu et al. Nature Communications 2014, Wang et al. Advanced Materials 2016), we investigate analytically and numerically the stability of biskyrmions in films of finite thickness, taking into account the nearest-neighbor exchange interaction, perpendicular magnetic anisotropy (PMA), dipole-dipole interaction (DDI), and the discreteness of the atomic lattice. The biskyrmion is characterized by the topological charge $Q=2$, the spatial scale $\lambda$, and another independent length $d$ that can be interpreted as a separation of two $Q=1$ skyrmions inside a $Q=2$ topological defect in the background of uniform magnetization. We find that biskyrmions with $d$ of order $\lambda$ can be stabilized by the magnetic field within a certain range of the ratio of PMA to DDI in a film having a sufficient number of atomic layers $N_{z}$. The shape of biskyrmions has been obtained by the numerical minimization of the energy of interacting spins in a $1000\times1000\times N_{z}$ atomic lattice. It is close to the exact solution of Belavin-Polyakov model when $d$ is below the width of the ferromagnetic domain wall. We compute the magnetic moment of a biskyrmion and discuss ways of creating biskyrmions in experiment.

Equations

$$E_{ex}=\frac{J}{2}\int dxdy\left(\nabla s_{\alpha }\cdot \nabla s_{\alpha }\right),$$

$$Q=\int \frac{dxdy}{4\pi }\mathbf{s}\cdot \frac{\partial \mathbf{s}}{\partial x}\times \frac{\partial \mathbf{s}}{\partial y}$$

$$\mathbf{s}\times {\nabla }^2\mathbf{s}=0.$$

$$E_{ex}=4\pi J|Q|.$$

$${\omega }_1=\frac{2s_x}{1+s_z},{\omega }_2=\frac{2s_y}{1+s_z}.$$

$$\frac{\partial {\omega }_1}{\partial x}=\pm \frac{\partial {\omega }_2}{\partial y},\frac{\partial {\omega }_1}{\partial y}=\mp \frac{\partial {\omega }_2}{\partial x}$$

$$\omega (z)=e^{i\gamma }\frac{(z-d/2)(z+d/2)}{l^2}=e^{i\gamma }\frac{[z^2-(d/2{)}^2]}{l^2}.$$

$$s_x=\frac{\mathrm{Re}(\omega )}{1+|\omega {|}^2/4},s_y=\frac{\mathrm{Im}(\omega )}{1+|\omega {|}^2/4},s_z=\frac{1-|\omega {|}^2/4}{1+|\omega {|}^2/4}.$$

$$s_x$$

$$s_y$$

$$s_z$$

$$s_z=s_{z,\mathrm{sad}}=\frac{{\lambda }^4-(d/2{)}^4}{{\lambda }^4+(d/2{)}^4}.$$

$$\mathbf{s}(\mathbf{r})=\left\{\frac{2{\lambda }^2{r}^2\cos (\gamma \pm 2\varphi )}{{\lambda }^4+r^4},\frac{2{\lambda }^2{r}^2\sin (\gamma \pm 2\varphi )}{{\lambda }^4+r^4},\frac{{\lambda }^4-r^4}{{\lambda }^4+r^4}\right\},$$

$$\omega =e^{i\gamma }{\left(z/l\right)}^{|Q|},\omega =e^{i\gamma }{\left(z^{\ast }/l\right)}^{|Q|},$$

$$\mathbf{s}(\mathbf{r})$$

$$I_1=S_z=2\pi \int rdr\left[s_z(r)+1\right]=\frac{2{\pi }^2{\lambda }^2}{|Q|\sin \left(\pi /|Q|\right)}$$

$$I_2=2\pi \int rdr{\left[s_z(r)+1\right]}^2=\frac{4{\pi }^2{\lambda }^2\left(|Q|-1\right)}{Q^2\sin \left(\pi /|Q|\right)},$$

$$E_{ex}=\frac{J}{2}\int dxdy\frac{|d\omega /dz{|}^2}{(1+|\omega {|}^2/4{)}^2}.$$

$$S_z=\int \frac{dxdy}{a^2}(s_z+1),$$

$$s_z+1=\frac{2}{1+|\omega {|}^2/4}.$$

$$|\omega {|}^2/4=u^2-2pu\cos (2\varphi )+p^2,u\equiv \frac{r^2}{{\lambda }^2},p\equiv \frac{d^2}{4{\lambda }^2}$$

$$S_z={\left(\frac{\lambda }{a}\right)}^{2}f(p),$$

$$f(p)={\int }_0^{\infty }{\int }_0^{2\pi }\frac{dud\varphi }{u^2-2pu\cos (2\varphi )+p^2+1}.$$

$$f(p)=2\pi {\int }_0^{\infty }\frac{du}{{\left[(u^2-p^2+1{)}^2+4p^2\right]}^{1/2}}$$

$$S_z={\left(\frac{\lambda }{a}\right)}^2\{\begin{array}{ll}{\pi }^2,& d=0\\ 16\pi {\left(\lambda /d\right)}^2\ln (d/\lambda ),& d\gg \lambda .\end{array}$$

$$S_z=8\pi {\left(\frac{\tilde{\lambda }}{a}\right)}^2\ln \left(\frac{1.05d}{\tilde{\lambda }}\right)$$

$$S_z^{\{1\}}=4\pi {\left(\frac{\tilde{\lambda }}{a}\right)}^2\ln \left(\frac{\Lambda }{\tilde{\lambda }}\right),$$

$${\mathcal{H}}_{ex}=-\frac{1}{2}\sum _{ij}{J}_{ij}{\mathbf{s}}_i\cdot {\mathbf{s}}_j,$$

$$E_{ex}=\frac{J}{2}\int dxdy\left[{\left({\partial }_x\mathbf{s}\right)}^2+{\left({\partial }_y\mathbf{s}\right)}^2\right]$$

$$\delta E_{ex}=-\frac{J{a}^2}{24}\int dxdy\left[{\left({\partial }_x^2\mathbf{s}\right)}^2+{\left({\partial }_y^2\mathbf{s}\right)}^2\right].$$

$$\delta E_{ex}=-\frac{{\pi }^{2}J}{3}\frac{Q^2-1}{\sin (\pi /|Q|)}{\left(\frac{a}{\lambda }\right)}^2.$$

Full text

I Introduction

Magnetic skyrmions in 2D films represent a very active field of research due to their potential for topologically protected data storage and information processing at the nanoscale Nagaosa2013 ; Zhang2015 ; Klaui2016 ; Leonov-NJP2016 ; Hoffmann-PhysRep2017 ; Fert-Nature2017 . Skyrmions were initially introduced in high energy physics as nonlinear field models of elementary particles SkyrmePRC58 ; Polyakov-book ; Manton-book . They entered condensed matter physics after it was realized that Skyrme’s theory described topological defects in ferro- and antiferromagnetic films BelPolJETP75 ; Lectures ; Brown-book . Similar topology leads to skyrmions in Bose-Einstein condensates AlkStoNat01 , quantum Hall effect SonKarKivPRB93 ; StonePRB93 , anomalous Hall effect YeKimPRL99 , and liquid crystals WriMerRMR89 .

Skyrmions are characterized by the topological charge $Q=\pm 1,\pm 2,...$. In a $2D$ exchange model of a continuous spin field the conservation of $Q$ is provided by topolgy: Different $Q$ arise from different homotopy classes of the mapping of the three-component fixed-length spin field onto the $2D$ plane. Similar topological properties are possessed by the magnetic bubbles that have been intensively studied in 1970s MS-bubbles ; ODell . They were in effect cylindrical domains surrounded by domain walls of thickness $\delta$ that is small compared to the radius of the domain $R$. In typical ferromagnets $\delta\sim 10$-$100$nm, so the bubbles of 1970’s were at least of a micron size or greater. With the emergence of nanoscience and nanoscale measuring techniques, the experimentalists have been able to observe topological defects in $2D$ magnetic films of size smaller than $\delta$. This is when the field of “skyrmionics” took off in condensed matter physics. Unlike the bubbles, nanoscale skyrmions are much closer to the topological defects described by the Skyrme model.

Research on magnetic skyrmions has focused on their stability, dynamics, and various symmetry properties. Perpendicular magnetic anisotropy (PMA), dipole-dipole interaction (DDI), magnetic field, and confined geometry can stabilize significantly large magnetic bubbles IvanovPRB06 ; IvanovPRB09 ; Moutafis-PRB2009 ; Ezawa-PRL2010 ; Makhfudz-PRL2012 . For small skyrmions, violation of the scale invariance by the crystal lattice with a finite atomic spacing $a$ leads to a stronger violation of the conservation of $Q$. In a pure exchange model, ferromagnetic skyrmions of size $\lambda$ collapse CCG-PRB2012 on a time scale proportional to $(\lambda/a)^{5}$. Their stability requires other than Heisenberg exchange coupling, strong random field or random anisotropy, or, most commonly, a non-centrosymmetric system with large Dzyaloshinskii-Moriya interaction (DMI) AbanovPRB98 ; Bogdanov-Nature2006 ; Heinze-Nature2011 ; Leonov-NatCom2015 ; Chen-APL2015 ; Boulle-NatNano2016 ; Lin-PRB2016 ; Leonov-NJP2016 ; EC-DG-PRL2018 ; EC-DG-NJP2018 .

To date, the presence of stable biskyrmions (see the computer-generated image in Fig. 1) has been independently reported at least in two centrosymmetric films of finite thickness: the La2-2xSr1+2xMn2O7 manganite Yu-2014 and the (Mn1-xNix)65Ga35 half Heusler alloy Zhang2016 . Given that stability of $Q=1$ skyrmions in films of large lateral dimensions normally require DMI, these findings are quite amazing and call for a theoretical analysis. Previously we have shown GCZ-EPL2017 that clusters with $Q>1$ are naturally generated due to the presence of Bloch lines in domain walls when labyrinth domains are destroyed by the magnetic field in a centrosymmetric magnetic film of finite thickness. We also observed that $Q=2$ biskyrmions were generated by slow relaxation of the system at $T=0$ starting from the disordered spin state or, equivalently, by slow cooling to $T=0$ from high temperature. A more detailed numerical investigation of biskyrmions, including the current-induced dynamics, was performed in Ref. XZhang2017, within a 2D frustrated micromagnetic model. Biskyrmions arising from frustrated Heisenberg exchange have also been reported in studies of triangular spin lattices Leonov-NatCom2015 and in Ginzburg-Landau theory of skyrmions Lin-PRB2016 . Metastable biskyrmion configurations have been observed in Landau-Lifshitz dynamics of a frustrated bilayer film Nowak2017 . However, the study of the separation of skyrmions in a biskyrmion, the shape of the biskyrmion, and its stability on the applied magnetic field, on the strength of the DDI, and on the film thickness in centrosymmetric systems has been absent so far. This, together with the above-mentioned experimental findings, provided motivation for our work.

In the past, biskyrmions have been intensively studied in nuclear physics in a hope that they would provide a model of a deuteron D1 ; D2 . The Belavin-Polyakov (BP) pure exchange model BelPolJETP75 in 2D contains an exact solution with $Q=2$ characterized by an arbitrary spatial dimension $\lambda$, and another arbitrary parameter $d$ that can be visually interpreted as a separation of two $Q=1$ skyrmions in a $Q=2$ topological defect (see Fig. 1), although the non-linearity of the model makes such interpretation meaningful only at a large separation. At $d=0$ the defect possesses symmetry with respect to the rotation in the $xy$-plane, which gets broken by any $d\neq 0$. We derive analytical formulas for the spin components and the magnetic moment of a BP biskyrmion for arbitrary $\lambda$ and $d$.

In real systems, the magnetic biskyrmions are more complicated as they are formed by a number of competing interactions. In this paper we show that a biskyrmion spin configuration naturally arises from the energy minimization in a centrosymmetric film of finite thickness in a lattice model that contains nearest-neighbor exchange interaction, perpendicular magnetic anisotropy (PMA), and dipole-dipole interaction (DDI). In accordance with the experimental findings, we find that stable biskyrmions exist within a certain range of parameters in films of sufficient thickness.

One interesting observation that follows from our studies is that the DDI always favors a biskyrmion with a finite $d$ over a $Q=2$ topological defect that does not split into $Q=1$ skyrmions. Another interesting observation is that, although biskyrmions are stabilized by interactions other than ferromagnetic exchange, sufficiently small biskyrmions are always close to the BP shape. On changing the parameters and the external magnetic field one can change the shape of the $Q=2$ topological defect from a BP biskyrmion to a thin-wall biskyrmion bubble.

This paper is organized as follows. Analytical formulas for the spin components, the magnetic moment of a $Q=2$ skyrmion in a continuous spin-field BP exchange model, as well as the effect of the discreteness of the lattice are derived in Section II. Numerical results on biskyrmions in a lattice model of a centrosymmetric film of finite thickness with the nearest-neighbor exchange interaction, PMA, DDI, and the external magnetic field are presented in Section III. Our results and suggestions for experiments are discussed in Section IV.

II Skyrmions and biskyrmions in the 2D exchange model

II.1 Spin field in a Belavin-Polyakov (BP) biskyrmion

We begin with a 2D exchange model with the energy

[TABLE]

where summation over spin components $\alpha=x,y,z$ is assumed. Here $J$ is the exchange constant and ${\bf s}$ is a three-component fixed-length spin field, ${\bf s}^{2}=1$. All spin-field configurations ${\bf s}(x,y)$ are divided into homotopy classes characterized by the topological charge BelPolJETP75

[TABLE]

that takes quantized values $Q=0,\pm 1,\pm 2,...$. The extremal spin-field configurations satisfy

[TABLE]

Throughout this paper we choose uniform magnetization, ${\bf s}=(0,0,-1)$, in the negative z-direction at infinity.

The absolute energy minimum inside each homotopy class is given by

[TABLE]

The corresponding solutions with $Q>0$ are called skyrmions, while solutions with $Q<0$ are called antiskyrmions. They have the simplest form BelPolJETP75 ; Lectures in terms of a complex variable $\omega=\omega_{1}+i\omega_{2}$ with

[TABLE]

Eq. (3) then reduces to $\partial\omega/\partial x=\mp i\partial\omega/\partial y$ or, equivalently, to

[TABLE]

that corresponds to the Cauchy-Riemann conditions for the complex function $\omega$. They are satisfied by any analytical function $\omega(z)$, where $z=x+iy$.

The minimum-energy solutions for a skyrmion with $Q=1$ and an antiskyrmion with $Q=-1$ are given by $\omega=z/l$ and $\omega=z^{*}/l$, respectively, with $z^{*}$ being the complex conjugate of $z$. The minimum-energy biskyrmion with $Q=2$ corresponds to

[TABLE]

It is parametrized by the chirality angle $\gamma$ and two lengths: the parameter $l$ that roughly describes the size of the biskyrmion, and another parameter $d$ that can be visually interpreted as the separation of two $Q=1$ skyrmions in a biskyrmion, see Fig. 1. Notice that the nonlinearity of Eq. (3) does not support this interpretation per se, rather Eq. (7) suggests that a biskyrmion is a product of two $Q=1$ skyrmions.

In terms of $\omega$ of Eq. (5) the components of ${\bf s}$ are given by

[TABLE]

These formulas allow one to obtain the spatial dependence of the components of the spin field in the BP biskyrmion:

[TABLE]

where we have replaced $l$ with $\lambda$ satisfying $\lambda^{4}=4l^{4}$ to have fewer numerical factors in the formulas.

Fig. 1 provides visualization of the spin field given by the above equations with $\gamma=\pi/2$ (Bloch-type biskyrmion). At the saddle point $x=y=0$ one has

[TABLE]

For $d\gg\lambda$, the BP biskyrmion becomes a superposition of two $Q=1$ BP skyrmions with the size $\tilde{\lambda}=\lambda^{2}/d\ll\lambda$.

For a biantiskyrmion Eq. (7) with $z\rightarrow z^{*}$ yields the expression similar to Eqs. (9)-(11) with the sign of the $xy$ terms changed.

At $d=0$ Eqs. (9)-(11) acquire a simple form in the polar coordinates, ${\bf r}=(r,\phi)$,

[TABLE]

with a plus sign for the $Q=2$ skyrmion and a minus sign for the $Q=-2$ antiskyrmion.

In a particular case of a compact single-centered $|Q|=2$ topological defect, that does not split into $Q=\pm 1$ skyrmions or antiskyrmions and is given by

[TABLE]

respectively, one arrives at a more general expression that is valid for an arbitrary positive or negative integer $Q$:

[TABLE]

To determine the size $\lambda$ of such a skyrmion from the numerically computed spin field, one can use the integrals that can be obtained with the help of Eq. (15):

[TABLE]

for $|Q|\geq 2$ and

[TABLE]

etc. Note that $I_{1}$ is the total spin, $S_{z}$, of the topological defect (see below). If the values of $\lambda_{n}$ obtained through different $I_{n}$ are close to each other, the shape of the skyrmion is close to the BP shape. For $Q=\pm 1$ the first integral logarithmically diverges and one has to use $I_{n}$ with higher powers of $n$ to characterize the shape of the bubble (see, e.g., Eq. (8) of Ref. AP-2018, ).

The energy of any spin configuration can be computed with the help of Eq. (1) or by noticing its equivalent form

[TABLE]

Substitution of Eqs. (7) or (14) into Eq. (18) and integration with $dxdy=rdrd\phi$ and $z=x+iy=re^{i\phi}$reproduces Eq. (4).

II.2 Magnetic moment of the BP biskyrmion

The magnetic moment of the topological defect equals $g\mu_{B}S_{z}$, where $g$ is the gyromagnetic factor, $\mu_{B}$ is the Bohr magneton and $S_{z}$ is the total spin of the defect defined as the difference between the spin of the $2D$ plane with and without the defect. For the boundary condition $s_{z}=-1$ at infinity one has

[TABLE]

where $a$ is the lattice constant and $s_{z}$ is given by the last of Eqs. (8), which yields

[TABLE]

For a biskyrmion given by Eq. (7) with $l^{2}=\lambda^{2}/2$, switching to polar coordinates, one has

[TABLE]

Substituting this into Eqs. (20) and (19) one obtains

[TABLE]

where

[TABLE]

Integration over the angle yields

[TABLE]

For an arbitrary $p$ this integral can be expressed via special functions in a rather cumbersome way that we do not provide here. The function $f(p)$ plotted numerically is shown in Fig. 2. Its asymptotic behavior can be computed analytically:

[TABLE]

Numerical evaluation gives a more accurate but close result, $\ln(1.1d/\lambda)$, for the logarithmic cutoff in the case of $d\gg\lambda$.

According to Fig. 2 and Eqs. (22), (25), at a fixed $\lambda$ the magnetic moment of the biskyrmion rapidly decreases on increasing $d$. One should notice, however, that in this limit the total spin of the biskyrmion can be written as

[TABLE]

in terms of the effective size $\tilde{\lambda}=\lambda^{2}/d\ll\lambda$ of the $Q=1$ skyrmion in a biskyrmion. Note that the magnetic moment of an isolated $Q=1$ skyrmion is given by quantum

[TABLE]

where $\Lambda$ is a long-distance cutoff determined by the lateral dimension of the film, $L$, or by $\delta_{H}=\sqrt{J/H}$, whichever is shorter. Thus, by order of magnitude, $S_{z}$ remains the same regardless of the separation of $Q=1$ skyrmions in the biskyrmion and similar to the magnetic moment of an isolated $Q=1$ skyrmion. As we shall see the energy minimum of a biskyrmion in a real magnetic film is realized at $d\sim\lambda$. In this case $S_{z}\propto d^{2}\sim\lambda^{2}$, that is, the magnetic moment is proportional to the area occupied by the biskyrmion.

II.3 Lattice-discreteness correction to the energy

The scale invariance of the skyrmion energy is broken by the discreteness of the lattice. The exchange energy of the atomic spins considered as classical spin vectors with $\left|\mathbf{s}_{i}\right|=1$ has the form

[TABLE]

where $J_{ij}=J$ for the nearest neighbors and zero otherwise. We use the simple cubic lattice. The corresponding correction to the energy can be obtained by calculating this lattice sum using the BP solution. For $\lambda\gtrsim a$, expanding the dot products to the lowest second order in spatial derivatives of the spin field, one obtains the well-known result

[TABLE]

that is equivalent to Eq. (1) and leads to Eq. (4). Expansion to the fourth order yields the energy correction

[TABLE]

For a singled-centered skyrmion with an arbitrary $Q$, substituting Eq. (15) into Eq. (29), we obtain with the help of computer algebra

[TABLE]

The previously obtained result in Ref. CCG-PRB2012, for $Q=\pm 1$,

[TABLE]

follows from this formula in the limit of $|Q|\rightarrow 1$.

III Biskyrmions in a magnetic film with perpendicular magnetic anisotropy

III.1 Lattice model and dipolar field

In the numerical work we study the lattice model of a ferromagnetic film of finite thickness with the energy given by the sum over lattice sites $i,j$

[TABLE]

Here ${\cal H}_{ex}$ is given by Eq. (28), $D$ is the easy-axis PMA constant, and $H\equiv g\mu_{B}SB$, with $S$ being the value of the atomic spin and $B$ being the induction of the applied magnetic field. In the DDI part of the energy

[TABLE]

where $\mathbf{r}_{ij}\equiv\mathbf{r}_{i}-\mathbf{r}_{j}$ is the displacement vector between the lattice sites and $\alpha,\beta=x,y,z$ denote cartesian components. The parameter $E_{D}=\mu_{0}M_{0}^{2}a^{3}/(4\pi)$ defines the strength of the DDI, with $M_{0}=g\mu_{B}S/a^{3}$ being the magnetization for our lattice model and $\mu_{0}$ being the magnetic permeability of vacuum.

The ratio of the PMI and DDI is given by the dimensionless parameter $\beta\equiv D/(4\pi E_{D})$. For $\beta>1$, the energy of the uniform state with spins directed along the $z$-axis is lower than that of the state with spins lying in the film’s plane. For $\beta<1$, the state with spins in the plane has a lower energy. The most interesting practical case is $\beta\sim 1$ realized in many materials in which there is a considerable compensation of the effects of the PMA and DDI.

In most materials the exchange interaction is much stronger than all other interactions of the spins. Consequently in a centrosymmetric system with magnetic anisotropy and nonsingular defects, such as skyrmions, the magnetization can only smoothly rotate on a large spatial scale gauged by two lengths: the domain-wall width, $\delta=a\sqrt{J/D}\gg a$, and another length, $\delta_{H}=a\sqrt{J/H}\gg a$, generated by the field. To describe such non-uniform spin states one needs to consider a system with macroscopically large number of spins. In the lattice model this can lead to impractically large computation times. Besides, the systems with interactions that are weak compared to the ferromagnetic exchange are magnetically soft because the energy of the non-uniform structures per spin is small. This makes the convergence of the energy-minimization routine very slow. The good news is that in such a case the atomic-scale spatial resolution is excessive.

To speed up the computation for structures that are much larger than the atomic spacing $a$ one can rescale the problem to another lattice constant $b>a$ by first rewriting the energy in the continuous approximation and then discretizing it again. The rescaled model with the parameters

[TABLE]

has smaller number of mesh points $N_{\alpha}^{\prime}=(a/b)N_{\alpha}$ and smaller mismatch between $J$ and other parameters. This provides faster convergence. After the computation with the rescaled model is completed, one obtains the results for the original system by rescaling the parameters back. Note that the domain-wall width and the field related length are the same in the original and rescaled models,

[TABLE]

The same is valid for all spatial structures such as magnetic bubbles, etc. The above-mentioned rescaling is only important for the study of films of large lateral dimensions. In a confined geometry, such as, e.g., nanotracks, one can perform computation at the atomic level.

The use of strong anisotropy and strong DDI allows one to work with a relatively small system and to have a reasonably fast convergence GCZ-EPL2017 . In this case, however, the shape of the topological defects is closer to that of the thin-wall skyrmion bubbles than to the BP skyrmions. To obtain the latter one has to work with the skyrmion size $\lambda$ satisfying $a\lesssim\lambda\lesssim\delta$ that requires a smaller anisotropy constant. To investigate biskyrmions, which is the purpose of this paper, we use $D/J=0.001$. This requires a large system size, $N_{x}\times N_{y}=1000\times 1000$, and results in longer computation times.

An important parameter controlling the DDI is the film thickness represented by $N_{z}$ in the units of the atomic spacing $a$. For thin films that are studied here, the magnetization inside the film is nearly constant along the direction perpendicular to the film. Thus one can make the problem effectively two-dimensional by introducing the effective DDI between the columns of parallel spins, considered as effective spins of the 2D model. This greatly speeds up the computation. To this end, for the simple cubic lattice, one can write the dipolar coupling, Eq. (33), as $\varPhi_{ij,\alpha\beta}=\phi_{\alpha\beta}\left(n_{x},n_{y},n_{z}\right)$, where $n_{x}\equiv i_{x}-j_{x}$ etc., are the distances on the lattice and

[TABLE]

The effective DDI is defined by

[TABLE]

Using the symmetry, one can express this result in the form with only one summation,

[TABLE]

that is used in the computations.

The effective DDI (that can be pre-computed) has different forms in different ranges of the distance $r$. At $r\apprge aN_{z}$ it scales as the interaction of magnetic dipoles $1/r^{3}$, while at $r\lesssim aN_{z}$ it goes as $1/r$ that corresponds to the interaction of magnetic charges at the surface of the film. Numerical results in Fig. 3 show that at large distances the effective DDI is stronger in films of finite thickness than in pure 2D systems. This has an important effect on the stability of biskyrmions.

In the computations, the dipolar field from the downward spins outside the system was added to obtain the results that are valid for an infinite system. This field was computed as the field of the plate of a very large size magnetized downward plus the dipolar field created by the finite system under the consideration magnetized upward.

III.2 Energy landscape of rigid-shape biskyrmions

Interactions other than exchange deform BP skyrmions and biskyrmions. As the result, the exchange energy increases. In the numerical work this increase provides the measure of the shape distortion. Also, one can consider the energy $\Delta E$ of the state with skyrmions with respect to that of the uniform state with all spins down. This energy is smaller than $4\pi J|Q|$ due to the contributions of the PMA and DDI.

We start with exploring the energy landscape, $E(\lambda,d$), in the presence of biskyrmions, assuming the rigid BP shape given by Eqs. (9)-(11) and numerically calculating the energy with the help of Eq. (32). This is valid when the exchange is much greater than all other interactions and when the biskyrmion size is small compared to the domain wall width: $\lambda,d\lesssim\delta$. In fact, it is more convenient to parametrize $\lambda$ as $\lambda^{2}=\tilde{\lambda}(\tilde{\lambda}+d)$, where $\tilde{\lambda}$ is the actual size of an individual $Q=1$ skyrmion in the limit of a large separation $d$.

Whereas in the pure exchange model the energy of biskyrmions is independent of $\lambda$ and $d$, other interactions (here PMA, DDI, and Zeeman), with the applied field $H$ as a control parameter, break this invariance of the skyrmion energy and select the values of $\lambda$ and $d$ that provide the energy minimum. The energy minimum of the topological defect with $Q=2$ always corresponds to a biskyrmion with $d\neq 0$.

To the contrary, we find that for a bi-antiskymion the energy always has a minimum at $d=0$, thus, there should be only $d=0$ single-centered antiskyrmions with $Q=-2$. This confirms the findings of Ref. GCZ-EPL2017, , see its Fig. 9 that shows different kinds of magnetic bubbles obtained by the relaxation from a random spin state: All objects with $Q=-2$ are spatially symmetric antiskyrmions, whereas all objects with $Q=2$ have a finite separation, $d\neq 0$.

Fig. 4 (upper panel) shows an example of the biskyrmion energy landscape for the system of size $N_{x}\times N_{y}\times N_{z}=1000\times 1000\times 100$ for which most of the computations have been done. We used $D/J=0.001$ and $\beta=1$. For $H/J=-0.00015$ there is a local energy minimum at $\lambda/a=27.1$ and $d/a=41.5$. The difference $\Delta E$ from the energy of the uniformly magnetized state with spins down is given in the units of $4\pi J$. For a weakly distorted BP biskyrmion it is close to 2. As the parameter $\lambda$ decreases, the energy goes down due the lattice-discreteness correction, see Sec. II.3. On the other hand, for the biantiskyrmion shown in Fig. 5 (upper panel), for the same parameters as above, there is no energy minimum at all, see lower panel in Fig. 4. Thus the biantiskyrmion evolves to a $Q=-2$ antiskyrmion shown in Fig. 5 (lower panel) It should be noted that the computed energy landscape is the same for $\phi=0$ and $\phi=\pi/2$. For smaller $|H|$, there is an energy minimum of the biantiskyrmion at $d=0$ and $\lambda>0$.

Whereas the rigid biskyrmion approximation provides a good qualitative description, it is not completely accurate since the balance of different interactions is very subtle and small deformations of the shape have a significant effect on the energy. More accurate results can be obtained by the numerical energy minimization described in the next section.

III.3 Numerical energy minimization and results

In this section we compute minimum-energy configurations of spins. The numerical method DCP-PRB2013 combines sequential rotations of spins ${\bf s}_{i}$ towards the direction of the local effective field, ${\bf H}_{{\rm eff},i}=-\partial{\cal H}/\partial{\bf s}_{i}$, with the probability $\alpha$, and the energy-conserving spin flips (so-called overrelaxation), ${\bf s}_{i}\to 2({\bf s}_{i}\cdot{\bf H}_{{\rm eff},i}){\bf H}_{{\rm eff},i}/H_{{\rm eff},i}^{2}-{\bf s}_{i}$, with the probability $1-\alpha$. The parameter $\alpha$ plays the role of the effective relaxation constant. We mainly use the value $\alpha=0.03$ that provides the overall fastest convergence.

The dipolar part of the effective field takes the longest time to compute. The method uses Fast Fourier Transform (FFT) in the whole sample as one program step. Since the dipolar field is much weaker than the exchange, several cycles of spin alignment can be performed before the dipolar field is updated, which increases the computation speed. The total charge $Q$ of the topological defect has been computed numerically using the lattice-discretized version of Eq. (2).

Computations were performed with Wolfram Mathematica using compilation. Most of the numerical work has been done on the 20-core Dell Precision T7610 Workstation. The FFT for computing the DDI was performed via Mathematica’s function ListConvolve that implicitly uses many processor cores. For this reason, no explicit parallelization was done in our program. However, we have been able to run several independent computations at the same time.

To demonstrate how well our rescaling method given by Eq. (34) works, we have computed the effective radius $R_{\mathrm{eff}}$ of topological defects with $Q=1$, defined as $\pi R_{\mathrm{eff}}^{2}=S_{z}/2$, where $S_{z}$ is the total spin of the bubble defined by $S_{z}=\intop dxdy\left[s_{z}(x,y)+1\right]$. The “original” system is a grid of $316\times 316\times 32$ spins with $D/J=0.01$ and the lattice constant $a$. The rescaled system is a grid of $100\times 100\times 10$ spins with $D^{\prime}/J^{\prime}=0.1$ and the lattice constant $b=\sqrt{10}a=3.16a$. In both cases $\beta=1$. Fig. 6 shows a perfect agreement between $R_{\mathrm{eff}}$ and $R^{\prime}_{\mathrm{eff}}$. (Note that $\delta^{\prime}=\delta$ thus practically $R^{\prime}_{\mathrm{eff}}=R_{\mathrm{eff}}$) in the plot.

Subsequently we performed the energy minimization for biskyrmions with $Q=2$ in a system with $N_{x}\times N_{y}\times N_{z}=1000\times 1000\times 100$, $D/J=0.001$, and $\beta=1$ at different values of the applied field $H$. The separation $d$ was found numerically as the distance between the maxima of $s_{z}$ in the biskyrmion. Then $\lambda$ was extracted with the help of Eq. (12) using the numerically found saddle-point value of $s_{z}$.

In Fig. 7 one can see that there is always a finite separation, $d>0$, in the biskyrmion. For stronger fields, the biskyrmion is smaller and closer to the BP shape, as the energies $E_{ex}$ and $\Delta E$ are close to $8\pi J$. Biskyrmions collapse as the field reaches the stability threshold. With decreasing the magnitude of field the biskyrmion is expanding, gradually transforming into a thin-wall bubble. The ratio $d/\lambda$ changes from 1.5 on the left side of Fig. 7 to 1.1 on its right side. The transition from BP biskyrmions to biskyrmion bubbles can be seen in a significant deviation of energy from $E/(4\pi J)=2$. As the field strength further decreases, the bubble loses its circular shape and transforms into a laminar domain.

Typical spin configurations of biskyrmions are illustrated (zoomed) in Fig. 8 that provides more details as compared to Fig. 1. The in-plane spin components in a BP biskyrmion shown in the upper panel of Fig. 8 decay as $1/r^{2}$ in accordance with Eq. (9) and Eq. (13). On the contrary, in a biskyrmion bubble shown in the lower panel of Fig. 8 the in-plane spin components decay exponentially and are hardly visible away from the bubble. For topological defects with $Q=-2$ (antiskyrmions) no finite separation $d$ was detected in our computations.

Next we performed computations for a similar model with a stronger DDI, $\beta=0.5$. The DDI in excess of the PMA forces the spins into the film’s plane, which suppresses skyrmions. To prevent this from happening, a stronger negative magnetic field has to be applied. The results in Fig. 9 show a larger separation, $d/\lambda\apprge 2$, and the shape close to the BP shape as $E_{ex}/(4\pi J)$ is rather close to 2. On the right side of the figure, the instability of the uniform state with spins down occurs on decreasing the field’s strength.

For the model in which PMA is stronger than DDI, $\beta>1$, no stable skyrmions or biskyrmions were found as they were collapsing even at $H=0$. The topological structures can be stabilized by the negative magnetic field in the range $\beta^{*}<\beta<1$ that depends on the film’s thickness $N_{z}$. Numerical studies show that the range of $\beta$ narrows down for thin films. In the latter, the effect of the DDI is similar to that of the easy-plane PMA, so that one can introduce the effective anisotropy that includes both PMA and DDI, $\widetilde{D}=D(1-1/\beta)$. This effective anisotropy changes its sign at $\beta=1$. This results in the extremely weak distortion of the BP shape of biskyrmions and extremely small controlling fields $H$. However, the model with a single effective anisotropy cannot support stable topological defects, including skyrmions and biskyrmions, at any $\beta\neq 1$. Since in real materials $\beta$ cannot be tuned, very thin non-chiral films with no DMI are not a good medium for the skyrmions. To the contrary, for thicker films, competition of the short-range PMA and long-range ($\sim 1/r$) DDI creates a range of $\beta$ in which skyrmions and biskyrmions can exist.

In particular, for $N_{z}=20$ with $D/J=0.001$, stable biskyrmions were found only for $\beta=1$. The results in Fig. 10 show very weakly distorted BP biskyrmions with practically the same ratio $d/\lambda=1.62$ in the whole range of a very weak field $H$. For $\beta=0.5$ and even for $\beta=0.75$ either the skyrmions collapse if the negative field is too strong or the background spin-down state becomes unstable if the negative field is too weak. For the system of $500\times 500\times 10$ spins, biskyrmions were still found at $\beta=1$ but the range of $H$ was much narrower. One can see that here the exchange energy is much closer to 2 than for $N_{z}=20$, thus the shape of the biskyrmion is much closer to the BP shape. For $H/J<-3.8\times 10^{-5}$ the biskyrmion collapses. For $H/J>-2\times 10^{-5}$ the spin-down background becomes unstable.

Computations on the $500\times 500$ monolayer found no biskyrmions at all even for $\beta=1$. Initial states in the form of biskyrmions evolve into states with separate skyrmions far away from each other that exist in an extremely narrow region of $H$.

The bottom line of this investigation is that biskyrmions should be searched for in sufficiently thick films, such as $N_{z}=100$ for $D/J=0.001$. For other values of the PMA, the scaling of Eq. (34) can be used. For a stronger PMA, biskyrmions can be supported by thinner films.

IV Discussion

We have studied biskyrmions in non-chiral ferromagnetic films of finite thickness with perpendicular magnetic anisotropy (PMA) with account of dipole-dipole interaction (DDI), and discreteness of the atomic lattice. In agreement with experimental findings Yu-2014 ; Zhang2016 we have found that biskyrmions are stable above a certain threshold in the film thickness. In films of insufficient thickness, (e.g., with the number of the atomic layers $N_{z}=10$ for $D/J=0.001$, stable biskyrmions exist only at $\beta=$ PMA/DDI$\,=1$, when the DDI nearly compensates the effect of the PMA. Since this condition is impossible to satisfy in practice, we conclude that stable biskyrmions do not exist in thin non-chiral magnetic films that contain just a few atomic layers. For the monolayer, $N_{z}=1$, no biskyrmions at all were found.

In thicker films (with, e.g., $N_{z}\sim 100$ for $D/J=0.001$ and thus $\delta/a\lesssim N_{z}$ ) stable biskyrmions exist within a finite range of $\beta$ below $\beta=1$. For stronger stabilizing fields, smaller biskyrmions of sizes $\lambda\lesssim d\lesssim\delta$ and the shape close to that provided by the Belavin-Polyakov 2D exchange model, Eqs. (9) - (11), have been observed in our numerical studies. Their energy is close to $8\pi J$. Such BP biskyrmions collapse at a critical value of the stabilizing field. In the opposite limit, when the magnitude of the stabilizing field decreases, the BP biskyrmion transforms into a bigger thin-wall biskyrmion bubble. Regardless of the film thickness and other parameters we did not find any stable biantiskyrmions with $d>0$.

The study presented here has focused on the stability of individual biskyrmions. Biskyrmions observed in experiments Yu-2014 ; Zhang2016 formed distorted triangular lattices. They were obtained from labyrinth domains on increasing the magnetic field in a manner similar to how lattices of $Q=1$ skyrmions have been observed. While our numerical method easily generated lattices of $Q=1$ skyrmions in the films with DMI, numerical generation of stable biskyrmion lattices in centrosymmetric 2D systems remains a challenging problem.

In our computations, stability of a single biskyrmion required an external field while the biskyrmion lattice has been observed even in a zero field Yu-2014 . This means that a sufficiently dense biskyrmion lattice minimizes the sum of the exchange, DDI, and PMA energies at $H=0$, similarly to what happens in the domain state.

For practical applications one has to be able to generate and manipulate individual biskyrmions. It has been demonstrated that skyrmions can be created, annihilated and moved by current-induced spin-orbit torques Yu-NanoLet2016 ; Fert-Nature2017 ; Legrand-Nanolet2017 . Individual $Q=1$ skyrmion bubbles have been generated by pushing elongated magnetic domains through a constriction using an in-plane current Jiang-Sci2015 ; Hoffmann-PhysRep2017 . These methods may not be suited for creating biskyrmions.

It has been shown that small skyrmions can be written and deleted in a controlled fashion with local spin-polarized currents from a scanning tunneling microscope Romming-Sci2013 . It has been also demonstrated that light-induced heat pulses of different duration and energy can write skyrmions in a magnetic film in a broad range of temperatures and magnetic fields Berruto-PRL2018 . These methods can be better suited for creating biskyrmions if experimentalists find the way of using a scanning tunneling microscope with a double tip or heat pulses of the shape resembling biskyrmions.

Recently, it has been experimentally demonstrated and confirmed through micromagnetic computations that stripe domains in a film can be cut into $Q=1$ skyrmions by the magnetic field of the tip of a scanning magnetic force microscope (MFM) Senfu-APL2018 . Writing individual $Q=1$ skyrmions by the MFM tip have been studied theoretically in Ref. AP-2018, . One can also use for that purpose magnetic nanoparticles of the kind used in nanocantilevers for mechanical magnetometry deLoubens .

A simple modification of the above method, tailored to biskyrmions, may be developed by using a double MFM tip consisting of two single tips in close proximity to each other, or a nanocantilever with two magnetic nanoparticles next to each other. We have tested this numerically in the same manner as is described in detail in Ref. AP-2018, . A biskyrmion created this way in the numerical experiment relaxes to the equilibrium size and shape determined by the parameters of the film.

Since biskyrmions carry magnetic moments similar to those of $Q=1$ skyrmions, they can be utilized in a similar way for data storage and information processing. The lack of the rotational symmetry in a biskymion makes its orientation another useful parameter in addition to the magnetic moment. It can open other functionalities in manipulating such information carriers as well. For example, the magnitude of the force exerted on a biskyrmion by the spin polarized current would depend on the direction of the current.

V Acknowledgements

This work has been supported by the grant No. OSR-2016-CRG5-2977 from King Abdullah University of Science and Technology.