Biskyrmions Lattices in Centrosymmetric Magnetic Films

TL;DR

This paper develops a theoretical framework explaining the formation and evolution of biskyrmion lattices in non-chiral magnetic films, matching experimental observations of their structure and stability.

Contribution

It introduces a mathematical model using elliptic functions to describe biskyrmion lattices and studies their evolution under realistic magnetic interactions.

Findings

Biskyrmion lattices are described by elliptic functions in a 2D exchange model.

The evolution leads to a stable triangular lattice of biskyrmion bubbles with Bloch lines.

The model aligns with experimental observations of non-chiral magnetic films.

Abstract

Abstract Theoretical framework is developed that explains biskyrmion lattices observed in non-chiral magnetic films. We study films of finite thickness containing up to $1000\times1000\times100$ spins. Hexatic biskyrmion lattices in a pure 2D exchange model are naturally described by the Weierstrass $\wp$ and $\zeta$ elliptic functions. Starting with such a lattice as an initial state we investigate how it evolves towards the minimum-energy state in the presence of perpendicular magnetic anisotropy and dipole-dipole interaction. In accordance with experiments, we find that the final state is a triangular lattice of biskyrmion bubbles containing Bloch lines.

Equations

$${\mathcal{H}}_{ex}=\frac{J}{2}\int dxdy\left(\frac{\partial s}{\partial x}\cdot \frac{\partial s}{\partial x}+\frac{\partial s}{\partial y}\cdot \frac{\partial s}{\partial y}\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,$$

$$\frac{\partial \mathbf{s}}{\partial x}=\pm \mathbf{s}\times \frac{\partial \mathbf{s}}{\partial y},\frac{\partial \mathbf{s}}{\partial y}=\mp \mathbf{s}\times \frac{\partial \mathbf{s}}{\partial x}.$$

$$E=4\pi J|Q|.$$

$$\omega (z)=\frac{s_x+i{s}_y}{1-s_z}$$

$$\frac{\partial \omega }{\partial x}=\pm i\frac{\partial \omega }{\partial y},$$

$$s_x+i{s}_y=\frac{2\omega }{|\omega {|}^2+1},s_z=\frac{|\omega {|}^2-1}{|\omega {|}^2+1}.$$

$$\omega (z)=\sum _i\frac{{\lambda }_i{e}^{i{\gamma }_i}}{z-z_i},$$

$$\omega (z)={\left(\sum _i\frac{{\lambda }_i{e}^{i{\gamma }_i}}{z-z_i}\right)}^{-1}.$$

$$\omega (z)=\sum _i\frac{{\lambda }_i{e}^{i{\gamma }_i}}{z^{\ast }-z_i^{\ast }},$$

$$\omega (z)={\left(\frac{\lambda }{z^{\ast }-z_0^{\ast }}\right)}^Q{e}^{i\gamma }.$$

$$\omega (z)=\frac{\lambda e^{i{\gamma }_1}}{z^{\ast }+d/2}+\frac{\lambda e^{i{\gamma }_2}}{z^{\ast }-d/2}.$$

$$\wp (z^{\ast })$$

$$p_1=\frac{a^{\prime }}{2},p_2=\frac{-1+i\sqrt{3}}{4}{a}^{\prime },$$

$$\zeta (z^{\ast })$$

$$\omega (z^{\ast })=\frac{{\lambda }^2}{d}{e}^{i\gamma }\left[\zeta \left(z^{\ast }+\frac{d}{2}{e}^{i\eta }\right)-\zeta \left(z^{\ast }-\frac{d}{2}{e}^{i\eta }\right)\right]$$

$$\mathcal{H}$$

$${\Phi }_{ij,\alpha \beta }\equiv a^3{r}_{ij}^{-5}\left(3r_{ij,\alpha }{r}_{ij,\beta }-{\delta }_{\alpha \beta }{r}_{ij}^2\right),$$

$${\varphi }_{\alpha \beta }\left(n_x,n_y,n_z\right)=\frac{3n_{\alpha }{n}_{\beta }-{\delta }_{\alpha \beta }\left(n_x^2+n_y^2+n_z^2\right)}{{\left(n_x^2+n_y^2+n_z^2\right)}^{5/2}}.$$

$${\bar{\varphi }}_{\alpha \beta }\left(n_x,n_y\right)=\frac{1}{N_z}\sum _{i_z,j_z=1}^{N_z}{\varphi }_{\alpha \beta }\left(n_x,n_y,i_z-j_z\right).$$

$${\bar{\varphi }}_{\alpha \beta }\left(n_x,n_y\right)$$

Full text

I Introduction

Skyrmions were initially introduced in nuclear physics as solutions of the nonlinear $\sigma$-model that can describe atomic nuclei SkyrmePRC58 ; Polyakov-book ; Manton-book . They possess topological charge: $Q=\pm 1,\pm 2$, etc. Different $Q$ arise from different homotopy classes of the mapping of a three-component fixed-length field onto the $xy$ plane. $Q=1$ would correspond to a single nucleon, $Q=2$ could describe a deuteron D1 ; D2 , $Q=4$ would provide a model of an $\alpha$-particle, etc.

In magnetic films skyrmions are defects of the ferromagnetic order. They represent a very active field of research due to their nanoscale size and potential for a dense topologically protected data storage and information processing Nagaosa2013 ; Zhang2015 ; Klaui2016 ; Leonov-NJP2016 ; Hoffmann-PhysRep2017 ; Fert-Nature2017 . Much bigger micron-size magnetic bubbles intensively studied in 1970s MS-bubbles ; ODell possessed similar topological properties. They were cylindrical domains surrounded by narrow domain walls. On the contrary, a typical skyrmion would be small compared to the domain wall thickness, making it conceptually similar to the topological objects studied in nuclear physics BelPolJETP75 ; Lectures .

Research on magnetic skyrmions focuses on their stability and dynamics. Magnetic bubble lattices are in effect domain structures MS-bubbles ; ODell ; Ezawa-PRL2010 ; Makhfudz-PRL2012 stabilized by the perpendicular magnetic anisotropy (PMA), dipole-dipole interaction (DDI), and the external magnetic field. On the contrary, in the absence of other interactions, small skyrmions collapse CCG-PRB2012 due to violation of the scale invariance of the 2D exchange model by the atomic lattice. They can be stabilized by Dzyaloshinskii-Moriya interaction that is present in materials lacking inversion symmetry (DMI) Bogdanov-Nature2006 ; Heinze-Nature2011 ; Boulle-NatNano2016 ; Leonov-NJP2016 , or by quenched randomness EC-DG-PRL2018 ; EC-DG-NJP2018 that has effect similar to the DMI.

While $Q=1$ skyrmions and lattices of such skyrmions have been intensively studied, observation of Q = 2 biskyrmions, see Fig. 1, and biskyrmion lattices has been rare. Unlike skyrmions, that have been mostly observed in chiral films, stable biskyrmions lattices have been initially reported down to zero field in two non-chiral films of sufficient thickness: the La2-2xSr1+2xMn2O7 manganite Yu-2014 and the (Mn1-xNix)65Ga35 half Heusler alloy Zhang2016 . Recent works Loudon2019 ; Yao2019 ; Gobel2019 provided further experimental and numerical evidence of stable biskyrmions in centrosymmetric magnetic materials.

Theoretical works on biskyrmions have been scarce. It has been shown GCZ-EPL2017 that topological defects with $Q>1$ appear naturally in centrosymmetric films due to the presence of Bloch lines in labyrinth domains on increasing magnetic field. 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 the studies of triangular spin lattices Leonov-NatCom2015 and in a model based upon Ginzburg-Landau theory Lin-PRB2016 . Metastable biskyrmion configurations have been observed in Landau-Lifshitz dynamics of a frustrated bilayer film Nowak2017 . However, the study of biskyrmion lattices has been absent so far.

In real systems, the magnetic biskyrmions are more complicated than generic $Q=2$ solutions of the Belavin-Polyakov (BP) 2D exchange model. They are formed by a number of competing interactions on top of the exchange, such as PMA and DDI. The latter rules out any meaningful analytical solution. In this paper we adopt the following approach. First, we prepare a biskyrmion lattice that is a solution of the pure exchange model. Fortunately, it can be described in the complex plane, $z=x+iy$, by standard elliptic Weierstrass functions that have been previously used to build multiskyrmion configurations in nuclear physics Manton-book . We then turn on the PMA and DDI and compute numerically evolution of the system towards the spin configuration that corresponds to the energy minimum.

Our calculations are performed on 3D lattices containing $1000\times 1000\times 100$ spins that describe magnetic films of finite thickness. In accordance with real experiments, we find that the lowest energy achieved in our numerical experiments corresponds to a triangular lattice of biskyrmion bubbles containing Bloch lines.

The paper is organized as follows. Analytical expressions for individual biskyrmion and biskyrmion lattices within the BP model are given in Section II. Biskyrmion lattices in a 2D magnetic film with the PMA and DDI are studied numerically within a discrete spin model in Section III. Our results are summarized in Section IV.

II Biskyrmion lattices in a 2D exchange model

The 2D exchange model is described by the energy

[TABLE]

where ${\bf s}$ is a three-component fixed-length spin field, ${\bf s}^{2}=1$, $J$ is the exchange constant, and summation over spin components $\alpha=x,y,z$ is performed. We choose uniform magnetization, ${\bf s}=(0,0,-1)$, in the negative z-direction at infinity. Spin-field configurations belong to homotopy classes characterized by the topological charge

[TABLE]

that takes quantized values $Q=0,\pm 1,\pm 2,...$. The value of $Q$ shows how many times the spin vector circumscribes the full body angle $4\pi$ as the position vector covers the whole $xy$ plane.

In each homotopy class the energy is minimized by ${\bf s}(x,y)$ satisfying

[TABLE]

which means that the spins are collinear with the exchange field. In 2D this is equivalent to

[TABLE]

One recovers Eq. (3) by differentiating the first of these equations over $y$, the second one over $x$, and subtracting one from the other. With the help of Eqs. (1) and (4) it is easy to obtain the relation between the exchange energy and the topological charge BelPolJETP75 :

[TABLE]

Extremal spin configurations can be obtained by mapping the problem onto the complex plane, $z=x+iy$. In terms of the complex function

[TABLE]

Eqs. (4) reduce to linear equations,

[TABLE]

that are familiar Cauchy-Riemann (CR) conditions of the analyticity of the function. The corresponding spin configuration follow from the relations

[TABLE]

The remarkable property of the model is that any analytic function $\omega(z)$ provides solution for the extremal spin configuration.

In particular, the sum of poles

[TABLE]

describes a collection of antiskyrmions with spins at the poles pointing up against the spin-down background, the parameters $z_{i}$, $\lambda_{i}$, $\gamma_{i}$ defining the position, size, and chirality angle of the $i$-th antiskyrmion. This function satisfies the CR condition, Eq. (7), with the plus sign. In-plane spin components are rotating counterclockwise as the observation point is moving clockwise around the center (the pole) of the antiskyrmion. Each antiskyrmion is contributing -1 to the total topological charge $Q$.

The state with the skyrmions whose central spins are pointing down against the spin-up background can be obtained by rotating all spins of the above collection of antiskyrmions by the angle $\pi$ around the $x$-axis. Conformal invariance of the 2D exchange model allows one to achieve the same effect by taking the reciprocal of the $\omega$-function,

[TABLE]

It describes a collection of skyrmions pointing down against the spin-up background, with each skyrmion corresponding to a zero of $\omega(z)$. Since the rotation of all spins preserves the homotopy class, such skyrmions have topological charge $Q=-1$, with in-plane spin components rotating clockwise when the observation point moves clockwise around the skyrmion’s center. This illustrates a less appreciated fact that the sign of $Q$ is determined by both, the topology of the solution and the boundary condition (direction of spins at infinity).

The expression with a complex-conjugated argument

[TABLE]

that satisfies the CR condition Eq. (7) with the minus sign, describes a collection of skyrmions whose central spins are pointing up against the spin-down background. Here, each skyrmion has topological charge $Q=1$.

The energy of the collection of skyrmions or collection of antiskyrmions is entirely defined by the total topological charge $Q$ through Eq. (5). Q determined by the number of poles and is independent of $z_{i}$, $\lambda_{i}$, and $\gamma_{i}$ due to the symmetry of the 2D exchange model. On the contrary, any combination of skyrmions and antiskyrmions yields a non-analytical $\omega(z)$ that violates the CR conditions and, hence, does not satisfy Eq. (3). In the physical language, this means that such a configuration cannot be static. Skyrmions and antiskyrmions must annihilate while conserving the total topological charge $Q$, with the excess energy going into spin waves.

For an arbitrary $Q$, a skyrmion solution centered at $z=z_{0}$ is given by

[TABLE]

The antiskyrmion solution is obtained by replacing $z^{*}$ and $z_{0}^{*}$ with $z$ and $z_{0}$. A biskyrmion with $Q=2$ and separation $d$ between two skyrmions of chirality $\gamma_{1}$ and $\gamma_{2}$ in the biskyrmion, centered at $z=0$, can be written as

[TABLE]

A generic Bloch-type biskyrmion with opposite chiralities $\gamma_{1}=-\gamma_{2}=\pi/2$ is shown in Fig. 1.

Regular lattices of topological defects (skyrmions or antiskyrmions), can be constructed using functions that are double-periodic in the complex plane $z=x+iy$. The elliptic Weierstrass $\wp$-function provides such a solution for the periodic lattice of $Q=2$ skyrmions,

[TABLE]

Here $m,n$ are integers, and $p_{1,2}$ are half periods of the elliptic function satisfying ${\rm Im}(p_{2}/p_{1})\neq 0$. Consequently $\wp(z^{*})=\wp(z^{*}+2mp_{1})=\wp(z^{*}+2np_{2})$, so that $Q=2$ skyrmions are found at $z^{*}=2mp_{1}+2np_{2}$.

Following experimental observations we will be interested in the triangular biskyrmion lattices. They correspond to the choice

[TABLE]

where $a^{\prime}$ is the skyrmion lattice spacing. The choice for the $\omega$-function is $\omega=\lambda^{2}e^{-i\gamma}\wp(z^{*})$, where $\lambda$ and $\gamma$ represent the size and chirality of skyrmions in the lattice. Two such lattices with $\gamma=0$ and $\gamma=\pi/2$ are shown in Fig. 2.

So far the lattices we have built had zero separation $d$ between the two skyrmions in the biskyrmion, see Eq. (13). To build biskyrmion lattices with a finite $d$, one can use Weierstrass $\zeta$-function defined by $d\zeta(z^{*})/dz^{*}=-\wp(z^{*})$. Its explicit form is given by

[TABLE]

This elliptic function is only quasi-periodic, satisfying $\zeta(z^{*}+2mp_{1}+2np_{2})=\zeta(z^{*})+2m\zeta(p_{1})+2n\zeta(p_{2})$. However, the function

[TABLE]

is double-periodic. It describes the lattice of biskyrmions of opposite chiralities that have a finite separation $d$ between skyrmions in a biskyrmion. The value of $\eta$ selects the orientation of the biskyrmion, with $\eta=0$ and $\eta=\pi/2$ corresponding to the horizontal and vertical orientation respectively. Lattices with vertical orientation of Néel and Bloch biskyrmions are shown in Fig. 3.

It is interesting to notice that standard elliptic functions are better suited for the description of $Q=2$ slyrmion and biskyrmion lattices than $Q=1$ skyrmion lattices.

III Biskyrmion lattices in a 2D magnetic film

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

[TABLE]

Here the exchange coupling is $J$ for the nearest neighbors on a simple cubic lattice and zero otherwise, $D$ is the easy-axis PMA constant, $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 in a 3D coordinate space. 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 due to considerable compensation of the effects of the PMA and DDI.

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. (19), 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. That 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.

Starting with a biskyrmion lattice of the previous section as an initial condition we compute its evolution towards the minimum-energy configuration in a system containing up to $1000\times 1000\times 100$ spins. Our 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 most efficient method uses updates of the dipolar field after all spins are updated rather than after updating each individual spin. Computation of the dipolar field uses the Fast Fourier Transform (FFT) algorithm that yields the dipolar fields in the whole sample as one program step. The total topological charge $Q$ of the lattice is computed numerically using the lattice-discretized version of Eq. (2) in which first derivatives are approximated by the four-point formula.

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.

For each skyrmion state we compute the topological charge $Q$, the exchange energy $E_{ex}$, the total energy $E$, and the difference, $\Delta E$, of the total energy from the energy of the uniformly magnetized state. The difference of the exchange energy from $4\pi JQ$ is the measure of the distortion of BP biskyrmions, while the sign of $\Delta E$ is indicative of whether the formation of the biskyrmion lattice can lower the total energy. Our numerical results for the energies are presented in the units of $J$.

Typical results of the energy minimization at $H=0$ and $\beta=1$ for two values of PMA, starting with the same Weierstrass biskyrmion lattice as the initial condition, are illustrated by Figs.  4 and 5. The final state does not depend on the separation $d$ in the biskyrmions but depends on the lattice period which defines the total topological charge of the system. The latter is approximately conserved and changes only slightly due to the use of the open boundaries. At low PMA, $D/J=0.001$, the system relaxes to a biskyrmion lattice shown in Fig. 4. Its energy, however, is above the energy of the uniformly magnetized state, $\Delta E>0$, indicating that it can only be a metastable state stabilized by the conservation of the topological charge.

At $D/J=0.03$ the initial Weierstrass biskyrmion lattice evolves into a lattice of biskyrmion bubbles shown Fig. 5. The energy of such a lattice is below the energy of the uniformly magnetized state, $\Delta E<0$, indicating that it can be a stable domain structure at $H=0$. It consists of regularly spaced magnetic bubbles with $Q=2$, each having two Bloch lines. Notice that Bloch lines in biskyrmion bubbles forming a triangular lattice at $H=0$ have been observed in experiment Yu-2014 .

Even a lower energy can be achieved by increasing the period of the lattice, which leads to bigger biskyrmion bubbles in the final state, see Fig. 6. This finding confirms conjecture of Ref. Loudon2019, that the absolute minimum of the energy is likely to correspond to the limiting case of biskyrmion bubbles of size much greater than the domain wall thickness.

IV Discussion

We have studied biskyrmion lattices in non-chiral magnetic films of finite thickness. While such lattices have been observed in experiments, their theoretical description presented a considerable challenge. Our approach has been based upon observation that representation of ferromagnetically-coupled spins in the complex plane, $z=x+iy$, that uses a complex function $\omega(z)=(s_{x}+is_{y})/(1-s_{z})$, naturally provides biskyrmion lattices when $\omega(z)$ is written in terms of elliptic functions. Starting with such a function as an initial condition, we study evolution of the spin configuration under the action of perpendicular magnetic anisotropy (PMA) and dipole-dipole interaction (DDI).

Our main result is that in a certain a range of parameters, that is relevant to experimental situations, the energy of a triangular lattice of biskyrmion bubbles is significantly lower than the energy of the uniformly magnetized film in a zero magnetic field. Notice in this connection that for a large, highly non-linear, hysteretic magnetic system that has infinite number of local energy minima one cannot find and compare energies of all metastable non-uniformly magnetized states that emerge due to the exchange interaction, PMA, and DDI (such as various kinds of labyrinth domains, various kinds of cylindrical domains or skyrmions, with and without Bloch lines and Bloch points, etc.). Consequently, it is impossible to say with the absolute certainty that a biskyrmion lattice observed in experiments and simulations is the ground state of the system. Its low energy computed in our model, however, makes it plausible that it represents a kind of the energy-minimizing domain structure that has not been observed until recently.

When biskyrmion lattices correspond to the energy minimum, they are likely to consist of bubbles of size much greater than the domain wall thickness. An additional connection between theory and observation is provided by the fact that biskyrmions arranged in a triangular lattice possess Bloch lines in both, experiments and computations. Further experimental studies of magnetic phases of non-chiral films along the hysteresis loop will shed more light on the place of biskyrmion lattices in the phase diagram.

V Acknowledgements

This work has been supported by the grant No. DE-FG02-93ER45487 funded by the U.S. Department of Energy, Office of Science.