Magnetic Skyrmion Lattice by Fourier Transform Method

TL;DR

This paper introduces a rapid Fourier-based numerical method for studying skyrmion lattices in magnetic materials, accurately predicting phase transitions and discovering new metastable states.

Contribution

A novel Fourier transform method for efficiently analyzing skyrmion and spiral phases in magnetic systems with Dzyaloshinsky-Moriya interaction, including metastable state prediction.

Findings

Reproduced known critical fields for phase transitions with high accuracy.

Predicted new metastable skyrmion lattice states stabilized by coupling with superconducting vortices.

Method adaptable to other micromagnetic systems.

Abstract

We demonstrate a fast numerical method of theoretical studies of skyrmion lattice or spiral order in magnetic materials with Dzyaloshinsky-Moriya interaction. The method is based on the Fourier expansion of the magnetization combined with a minimization of the free energy functional of the magnetic material in Fourier space, yielding the optimal configuration of the system for any given set of parameters. We employ a Lagrange multiplier technique in order to satisfy micromagnetic constraints. We apply this method to a system that exhibits, depending on the parameter choice, ferromagnetic, skyrmion lattice, or spiral (helical) order. Known critical fields corresponding to the helical-skyrmion as well as the skyrmion-ferromagnet phase transitions are reproduced with high precision. Using this numerical method we predict new types of excited (metastable) states of the skyrmion lattice, which may be stabilized by coupling the skyrmion lattice with a superconducting vortex lattice. The method can be readily adapted to other micromagnetic systems.

Equations

$${\hat{H}}_M=-J\sum _{⟨ij⟩}{\vec{S}}_i\cdot {\vec{S}}_j-\sum _{ij}{\vec{D}}_{ij}\cdot \left({\vec{S}}_i\times {\vec{S}}_j\right)+g{\mu }_B\sum _j{\vec{S}}_j\cdot \vec{B},$$

$$F\left[\vec{M}\right]$$

$$\tilde{F}\left[\vec{M}\right]$$

$$\mathcal{F}$$

$$\frac{\partial \mathcal{F}}{\partial M_{\mu }}-{\partial }_{\nu }\frac{\partial \mathcal{F}}{\partial \left({\partial }_{\nu }{M}_{\mu }\right)}=0,$$

$$\left({\partial }_x^2+{\partial }_y^2\right)M_x-4{\partial }_y{M}_z-4\lambda M_x=0,$$

$$\left({\partial }_x^2+{\partial }_y^2\right)M_y+4{\partial }_x{M}_z-4\lambda M_y=0,$$

$$\left({\partial }_x^2+{\partial }_y^2\right)M_z-4\left({\partial }_x{M}_y-{\partial }_y{M}_x\right)-4\lambda M_z=-2\beta .$$

$$M_x(x,y)=\sum _{m,n}{X}_{mn}{e}^{-i{\vec{k}}_{mn}\vec{r}},$$

$$\vec{k}_{mn}=\frac{2\pi}{x_{1}y_{2}}\left(\begin{array}[]{c}my_{2}\\ nx_{1}-mx_{2}\end{array}\right)=\left(\begin{array}[]{c}k_{x}\\ k_{y}\end{array}\right),$$

$$M_x=\sum _{\vec{k}}{X}_{\vec{k}}{e}^{-i\vec{k}\vec{r}},$$

$$\lambda =\sum _{\vec{k}}{\lambda }_{\vec{k}}{e}^{-i\vec{k}\vec{r}}.$$

$$\sum _{{\vec{k}}^{\prime }}\left[\left\{-k^{\prime 2}{X}_{{\vec{k}}^{\prime }}+4i{k}_y^{\prime }{Z}_{{\vec{k}}^{\prime }}\right\}{\delta }_{\vec{k}{\vec{k}}^{\prime }}-4X_{{\vec{k}}^{\prime }}{\lambda }_{\vec{k}-{\vec{k}}^{\prime }}\right]=0,$$

$$\sum _{{\vec{k}}^{\prime }}\left[\left\{-k^{\prime 2}{Y}_{{\vec{k}}^{\prime }}-4i{k}_x^{\prime }{Z}_{{\vec{k}}^{\prime }}\right\}{\delta }_{\vec{k}{\vec{k}}^{\prime }}-4Y_{{\vec{k}}^{\prime }}{\lambda }_{\vec{k}-{\vec{k}}^{\prime }}\right]=0,$$

$$\begin{gathered} \sum _{{\vec{k}}^{\prime }}\left[\left\{-k^{\prime 2}{Z}_{{\vec{k}}^{\prime }}+4i\left(k_x^{\prime }{Y}_{{\vec{k}}^{\prime }}-k_y^{\prime }{X}_{{\vec{k}}^{\prime }}\right)\right\}{\delta }_{\vec{k}{\vec{k}}^{\prime }}-4Z_{{\vec{k}}^{\prime }}{\lambda }_{\vec{k}-{\vec{k}}^{\prime }}\right] \\ =-2\beta {\delta }_{\vec{k},0}. \end{gathered}$$

$$\sum _{{\vec{k}}^{\prime }}{\hat{D}}_{\vec{k}{\vec{k}}^{\prime }}{\vec{M}}_{{\vec{k}}^{\prime }}={\vec{I}}_{\vec{k}},$$

$$\vec{I}_{\vec{k}}=\left(\begin{array}[]{c}0\\ 0\\ -2\beta\delta_{\vec{k},0}\end{array}\right),$$

$$\vec{M}_{\vec{k}}=\left(\begin{array}[]{c}X_{\vec{k}}\\ Y_{\vec{k}}\\ Z_{\vec{k}}\end{array}\right),$$

$${\hat{D}}_{\vec{k}{\vec{k}}^{\prime }}={\hat{K}}_{{\vec{k}}^{\prime }}{\delta }_{\vec{k}{\vec{k}}^{\prime }}-4{\lambda }_{\vec{k}-{\vec{k}}^{\prime }}\hat{I},$$

$$\hat{K}_{\vec{k}^{\prime}}=\left(\begin{array}[]{c c c}k^{\prime 2}&0&4ik_{y}^{\prime}\\ 0&k^{\prime 2}&-4ik_{x}^{\prime}\\ -4ik_{y}^{\prime}&4ik_{x}^{\prime}&k^{\prime 2}\end{array}\right).$$

$${\vec{M}}_{\vec{k}}=\sum _{{\vec{k}}^{\prime }}{\hat{D}}_{\vec{k}{\vec{k}}^{\prime }}^{-1}{\vec{I}}_{{\vec{k}}^{\prime }},$$

$$\sum _{\vec{k}}|{\vec{M}}_{\vec{k}}{|}^2$$

$$\sum _{\vec{k}}{\vec{M}}_{\vec{k}}\cdot {\vec{M}}_{\vec{q}-\vec{k}}$$

$$\tilde{F}\left[\vec{M}\right]$$

$$=\frac{1}{{\gamma }^2}{\tilde{F}}_{ex}\left[{\vec{M}}_{\gamma }\right]+\frac{1}{\gamma }{\tilde{F}}_{DM}\left[{\vec{M}}_{\gamma }\right]+{\tilde{F}}_B\left[{\vec{M}}_{\gamma }\right],$$

$$0$$

$$+\int \left[\frac{1}{{\gamma }^2}\frac{\delta {\tilde{\mathcal{F}}}_{ex}}{\delta {\vec{M}}_{\gamma }}+\frac{1}{\gamma }\frac{\delta {\tilde{\mathcal{F}}}_{DM}}{\delta {\vec{M}}_{\gamma }}+\frac{\delta {\tilde{\mathcal{F}}}_B}{\delta {\vec{M}}_{\gamma }}\right]\frac{\partial {\vec{M}}_{\gamma }}{\partial \gamma }\frac{d{x}^{\prime }d{y}^{\prime }}{A^{\prime }}.$$

$$-2{\tilde{F}}_{ex}\left[\vec{M}\right]-{\tilde{F}}_{DM}\left[\vec{M}\right]+\int \frac{\delta \tilde{\mathcal{F}}}{\delta \vec{M}}(\vec{r}\cdot \nabla )\vec{M}(\vec{r})\frac{dxdy}{A}=0,$$

$$\frac{\delta \tilde{\mathcal{F}}}{\delta \vec{M}}=0$$

$$-2{\tilde{F}}_{ex}[{\vec{M}}_{\mathrm{stat}}]-{\tilde{F}}_{DM}[{\vec{M}}_{\mathrm{stat}}]=0.$$

$$\frac{{\tilde{F}}_{DM}}{{\tilde{F}}_{ex}}=-2,$$

$$\tilde{F}=$$

Full text

I INTRODUCTION

Apart from the more familiar ordered magnetic phases, such as ferro- or antiferromagnetism, in recent years a plethora of more complex magnetically ordered states has been unveiled, both theoretically and experimentally. Among these states, a crucial role is played by magnetic skyrmions,Skyrme (1961) characterized by a whirling pattern of localized Heisenberg spins. Bogdanov and Yablonskii (1989); Bogdanov et al. (1989); Bogdanov and Hubert (1994a, b); Barrett et al. (1995); Brey et al. (1995); Bogdanov and Hubert (1999); Gervais et al. (2005); Rößler et al. (2006); Mühlbauer et al. (2009); Neubauer et al. (2009); Münzer et al. (2010); Yu et al. (2010); Nagaosa and Tokura (2013); Romming et al. (2013, 2015); Leonov et al. (2016); Leonov and Kézsmárki (2017) Magnetic skyrmions can be regarded as localized topological defects of the ferromagnetic state that can only be moved around or deformed in shape and size, but cannot be removed by a smooth deformation of the magnetization pattern due to micromagnetic constraints. They either exist as isolated defects, as, for example, in certain types of topological domain walls found in helical magnets, or they themselves give rise to an extended, energetically favorable phase of helical magnets, in the form of a skyrmion lattice.

Skyrmion lattices were first introduced by Klebanov Klebanov (1985) for neutron crystals, whereas Bogdanov and Yablonskii Bogdanov and Yablonskii (1989) discussed the possibility of thermodynamically stable magnetic vortices (now called skyrmions lattices) arising in anisotropic ferromagnets as an intermediate phase between the uniform and the spiral magnetic order. Specific theoretical predictions for skyrmion lattices in, e.g., $\rm{MnSi}$, $\rm{FeGe}$, and $\rm{Fe_{1-x}Co_{x}Si}$, where made by Bogdanov and Hubert,Bogdanov and Hubert (1994a) and subsequently skyrmion lattices were reported experimentally in $\rm{MnSi}$ Mühlbauer et al. (2009); Neubauer et al. (2009) and $\rm{Fe_{1-x}Co_{x}Si}$.Münzer et al. (2010); Yu et al. (2010) The earliest theoretical predictions Bogdanov and Hubert (1994a) and experimental evidences Nagaosa and Tokura (2013) showed that magnetic skyrmions typically form 2D triangular lattices that are perpendicular to an external magnetic field and translation invariant in the parallel direction. More recent theoretical Lin et al. (2015); Rowland et al. (2016) and experimental Kurumaji et al. (2017); Yu et al. (2018) studies have eventually revealed that an easy-plane magnetic anisotropy may also lead to the stabilization of square skyrmion and meron-antimeron Gross (1978); Affleck (1986); Moon et al. (1995); Brey et al. (1996) (half-skyrmion) lattices.

The formation of magnetic skyrmions may be attributed to various microscopic mechanisms, which often cooperate with each other. We can nevertheless classify these mechanisms into two broad classes, according to the relative size of the skyrmions with respect to the microscopic lattice constant. In the first class of mechanisms, skyrmions arise from the competition between the ferromagnetic exchange coupling and an anisotropic exchange interaction originating from relativistic spin-orbit coupling, namely the so-called Dzyaloshinsky-Moriya (DM) interaction,Dzyaloshinsky (1958); Dzyaloshinskii (1964); Moriya (1960) which breaks either inversion or mirror symmetry. While the ferromagnetic exchange tends to align all spins in the same direction, DM interaction favors a non-collinear alignment, thereby twisting neighboring spins. The size of a spiral or skyrmion pattern emerging from this mechanism is typically two orders of magnitude larger than the crystal lattice constant. Hence, in this case skyrmions are topologically robust against lattice defects and can be very well described in the continuum limit approximation. In the second class of mechanisms, instead, a skyrmion pattern may arise as a result of competing ferro- and antiferromagnetic exchange interactions in frustrated magnets,Okubo et al. (2012); Leonov and Mostovoy (2015); Lin and Hayami (2016); Hu et al. (2017) or in the presence of four-spin interactions.Heinze et al. (2011) The main distinction from the first class of skyrmion patterns is the length scale of the spin modulation, which in this case is of the order of the underlying lattice constant, leading to atomic-scale skyrmion patterns that are typically commensurate with the crystal lattice. Note that this type of skyrmions do not require a broken inversion or mirror symmetry as in the case of DM chiral magnets.

In this work we focus on the large-scale, topologically robust, skyrmion patterns that characterize chiral magnets. Moreover, we consider strictly 2D materials as realized, e.g., in thin films of $\rm{Fe_{1-x}Co_{x}Si}$ or $\rm{FeGe}$.Yu et al. (2010); Huang and Chien (2012) Indeed, in 3D chiral magnets the skyrmion lattice phase is thermodynamically stable only at finite temperatures and can be stabilized down to zero temperature only by a strong easy-axis anisotropyButenko et al. (2010) or the presence of multiple types of spin-orbit couplings.Rowland et al. (2016) Instead, in 2D materials the skyrmion phase is stable at zero temperature over a wide range of external magnetic fields even in the absence of magnetic anisotropies.

The model Hamiltonian for Heisenberg spins in the presence of a DM interaction and an external magnetic field is given by,

[TABLE]

where $J$ is the ferromagnetic exchange constant, $\vec{D}_{ij}$ is the DM coupling vector, $\vec{B}$ is the external field, $\mu_{B}$ is Bohr’s magneton, and $\vec{S}_{i}$ is a Heisenberg spin at site $i$. The first term in Eq. (1) represents the Heisenberg exchange interaction, the second term the anisotropic DM interaction, and the last term the Zeeman energy. The specific form of the DM vector $\vec{D}_{ij}$ depends on the type of relativistic spin-orbit coupling present in the system.Rowland et al. (2016) A Dresselhaus spin-orbit coupling results in a broken inversion symmetry, $\vec{r}\to-\vec{r}$, and a DM vector that is parallel to the relative direction $\hat{r}_{ij}\equiv(\vec{r}_{i}-\vec{r}_{j})/|\vec{r}_{i}-\vec{r}_{j}|$ of neighboring spins, $\vec{D}_{ij}=D_{\parallel}\hat{r}_{ij}$. This induces a rotation of the spin direction within the plane that is perpendicular to the inter-spin direction $\hat{r}_{ij}$, leading to a Bloch-like spiral or skyrmion pattern. Instead, a Rashba spin-orbit coupling results in a broken mirror symmetry with respect to a given plane, e.g., $z\to-z$. In this case the DM vector is perpendicular to both the broken symmetry axis and the inter-spin direction, $\vec{D}_{ij}=D_{\perp}\hat{z}\times\hat{r}_{ij}$, and the spins rotate within the plane that is simultaneously parallel to the inter-spin direction and the broken symmetry axis, leading to Néel-like spirals or skyrmions. In this paper we consider only the Dresselhaus DM interaction. However, we would like to stress that the method presented here is completely general in this respect and can be readily modified, at no additional numerical cost, to include both types of DM interactions and any easy-axis or easy-plane magnetic anisotropy.

In the absence of an external magnetic field (or if the field is sufficiently weak so that the Zeeman energy can be neglected), the interplay of DM and exchange interactions leads to a helical (spiral) spin order.Bak and Jensen (1980) For higher fields the Zeeman interaction becomes relevant and the spiral order gives way to a more complex magnetic state (the skyrmion lattice) until eventually at highest fields a ferromagnetic state is stabilized. Skyrmion magnetic order (see Fig. 1) occurs at intermediate external fields as a compromise between DM interaction and Zeeman energy. In the center of a skyrmion, the spin points in the direction opposite to the magnetic field. When moving radially from the center towards the periphery, the direction of the spin rotates in a similar fashion as in the helical state, until it points in the direction of the magnetic field at the skyrmion boundary. In the regions between the skyrmions the spins stay aligned with the external magnetic field. Hence, by varying the inter-skyrmion distance it is possible to arrive at a net energy gain resulting from both the DM and Zeeman interactions, stabilizing a skyrmion lattice as opposed to a spiral or ferromagnetic order.

The goal of this paper is to develop a fast and accurate numerical procedure for the calculation of skyrmion lattice or helical order by utilizing a method similar to the one used by Brandt for the effective numerical implementation of superconducting vortex lattices.Brandt (1995) We additionally increase the speed of the numerical convergence by invoking a virial theorem in order to find the equilibrium lattice constants for the skyrmion lattice. This general procedure also allows for the study of metastable solutions and can be adapted to combine skyrmion lattices with superconducting vortex lattices. Such a combined system may stabilize new types of skyrmion lattices that otherwise would be metastable.

In Section II we present our analytical and numerical approach to the skyrmion lattice problem, including the results obtained for the phase diagram of a 2D non-centrosymmetric ferromagnet and the field dependence of the inter-skyrmion distance; in Section III we discuss alternative metastable solutions for the skyrmion lattice and the possible ways to stabilize them, focusing in particular on the honeycomb skyrmion lattice.

II SKYRMION LATTICE PHASE

II.1 Euler-Lagrange Equations for a Magnetic System in Fourier Space

In this section we write the free energy functional for a two-dimensional non-centrosymmetric magnetic system in an external field and derive the corresponding Euler-Lagrange equations for the magnetization. Afterwards, we Fourier transform the obtained equations and solve them in Fourier space.

We start with the free energy functional written in the continuum limit. This approximation is legitimate as long as the magnetization varies on a length scale much larger than the crystal lattice spacing $a_{0}$, i.e., when $\left|\nabla\vec{M}\right|a_{0}\ll 1$. In Cartesian coordinates, the free energy per unit area reads

[TABLE]

where $J$ is the ferromagnetic exchange constant, $D$ is the DM interaction, and $\vec{B}$ is the external magnetic field applied (note that, in the continuum limit, the couplings $J$, $D$, and $\vec{B}$ have different dimensions than their corresponding lattice counterparts figuring in Eq. (1); moreover, they incorporate, as a multiplicative factor, the microscopic lattice coordination number). Here, $\vec{M}=\vec{M}(\vec{r})=\left(M_{x}(x,y),M_{y}(x,y),M_{z}(x,y)\right)$, obeying $|\vec{M}|^{2}=1$, and $\vec{B}=\left(0,0,B\right)$, with a constant $B$ corresponding to a uniform external field in the $z$-direction. We use the notation $\kappa=\frac{D}{2J}$ for the ratio between the DM interaction and the ferromagnetic exchange coupling. As $\kappa$ has dimensions of an inverse length, dimensionless coordinates result from $\left(\tilde{x},\tilde{y}\right)=\left(\kappa x,\kappa y\right)=\kappa\vec{r}$. This leads to the dimensionless form of the free energy functional $\tilde{F}=F/(2J\kappa^{2})$,Han et al. (2010)

[TABLE]

where $\vec{\beta}=\vec{B}/(2J\kappa^{2})$ is the dimensionless magnetic field. For simplicity, we shall drop all tildes from now on, remembering that lengths are measured in units of $\kappa^{-1}$ and energies in units of $2J$.

Since the magnetization obeys the micromagnetic constraint $|\vec{M}(x,y)|^{2}=1$ at every point in space, we introduce a Lagrange multiplier field, $\lambda(x,y)$, in order to enforce the aforementioned condition. The total functional density then reads

[TABLE]

We now write the Euler-Lagrange equations for each component of the magnetization, $M_{\mu}$,

[TABLE]

with summation convention applied for $\nu$, with $\nu\in\left\{x,y\right\}$. Explicitly, this leads to

[TABLE]

[TABLE]

[TABLE]

Being inspired by Brandt’s approach for Abrikosov vortices in type-II superconductors, Brandt (1995) we Fourier transform the magnetization components (and the Lagrange multiplier) and solve the problem in Fourier space. As the types of magnetic order promoted by the DM interaction, namely the skyrmion lattice and helical order, are periodic structures, the Fourier approach is preferable. The discrete Fourier transform of $M_{x}$ is written as

[TABLE]

where $X_{mn}=X_{\vec{k}_{mn}}$ is the Fourier coefficient for $M_{x}$, and the discrete wave vectors are given by

[TABLE]

with $m$ and $n$ integer indices and, for a triangular lattice with lattice spacing $a$, $x_{1}=a$, $x_{2}=\frac{a}{2}$, $y_{2}=\frac{a\sqrt{3}}{2}$. In order to shorten the notation, we write

[TABLE]

and by analogy for the remaining components of the magnetization. For the Lagrange multiplier we write

[TABLE]

The Fourier transformed Euler-Lagrange equations, (6)–(8), then read

[TABLE]

[TABLE]

[TABLE]

The above equations can be rewritten in matrix notation as

[TABLE]

where the term on the right hand side is proportional to the external magnetic field,

[TABLE]

$\vec{M}_{\vec{k}}$ is the vector containing the Fourier coefficients of the magnetization components for a given $k$,

[TABLE]

and the matrix $\hat{D}_{\vec{k}\vec{k}^{\prime}}$ is given by

[TABLE]

with $\hat{I}$ being an identity $3\times 3$ matrix, and

[TABLE]

Equation (16) can be solved analytically for $\vec{M}_{\vec{k}}$ via matrix inversion,

[TABLE]

yielding a magnetization functional $\vec{M}_{\vec{k}}[\lambda_{\vec{k}}]$ that depends on the Lagrange multiplier field. This solution, however, does not necessarily obey the micromagnetic constraint $|\vec{M}(x,y)|^{2}=1$. In order to enforce it, we must then solve the constraint equation numerically to find the appropriate Lagrange multiplier. In Fourier space, this amounts to solving the following equations for $\lambda_{\vec{k}}$,

[TABLE]

where we use the $\lambda$-dependence of the magnetization, $\vec{M}_{\vec{k}}=\vec{M}_{\vec{k}}[\lambda_{\vec{k}}]$, from Eq. (21). We find that the Lagrange multiplier $\lambda(x,y)$ takes the largest values in regions where the spins are aligned parallel or anti-parallel to the external field and the lowest values in regions where the spins are perpendicular to it (parallel to the $xy$-plane).

II.2 Skyrmion Lattice

Typical solutions obtained by the procedure described in the previous section are shown in Fig. 2. Keeping in mind that the applied field is directed in positive $z$-direction, a typical skyrmion lattice solution appears with spins aligning along the field direction, $z$, in the region between skyrmions, and in the opposite direction in the centers of skyrmions. Hence, the region between skyrmions is mostly ferromagnetic, while in moving away from the center of a skyrmion along any radial direction, the spin orientation rotates continuously within the plane perpendicular to the given radial direction, acquiring a component along the $x$-$y$ plane, in a similar fashion as found in helical magnetic patterns. In this way, the system is able to gain energy from the DM interaction, which favors the twisting of spins, while retaining a large region of ferromagnetic alignment, promoted by the Zeeman energy.

Shown in Fig. 2 are two self-consistently determined solutions for a triangular skyrmion lattice, one close to the phase transition to a spiral state (top panels), and one close to the phase transition to the homogeneous ferromagnetic state (bottom panels). In the first case the spatial modulation of the magnetization vector along a path going from one skyrmion center to a neighboring one approaches that of a spiral modulation, whereas in the latter case the skyrmions become isolated from each other and the inter-skyrmion regions expand. It is also interesting to observe, in these two extreme cases, the different spatial distribution of the $z$-component of the magnetization within a single skyrmion, whose radius can be conventionally identified in the red dashed contour ($M_{z}=0.95$). In the first case the area of negative magnetization (blue region, enclosed by the thick black contour) is approximately equal to the area of positive magnetization (red region, enclosed between the red dashed contour and the thick black contour). Hence, skyrmions in this case carry, as a whole, essentially no net magnetic moment, as in the case of a helix. On the other hand, in the vicinity of the ferromagnetic transition the area of positive magnetization is nearly three times that of negative magnetization, so that skyrmions do carry a finite magnetic moment along the external field direction. As we shall see later on, the shrinking of the negative magnetization region plays an important role in determining the magnetic field dependence of the inter-skyrmion distance.

Solutions to the Fourier transformed Euler-Lagrange equations, (13)–(15), satisfy the condition of stationarity of the functional within a given lattice geometry. However, there is no guarantee that such a solution is truly stationary with respect to arbitrary variations of the magnetization pattern and, even if this is the case, that the solution represents a global minimum of the functional and not just a local one. Hence, it is necessary to devise an additional procedure to distinguish stable solutions from those metastable or not stable at all.

We recall that the method described above is based on a Fourier expansion of the functional, where quantities such as the magnetization and the Lagrange multiplier are assumed to be periodic on a given lattice structure (triangular, for the time being) with lattice spacing $a$. When seeking a solution to the Euler-Lagrange equations we must therefore remember that we are only exploring the subspace of magnetization profiles with a fixed given periodicity. The classes of solutions that we can obtain in this way include (but are not limited to) the helical magnetic order, the triangular skyrmion lattice, and the ferromagnetic order. Hence, in order to find the globally stable solution, i.e., the ground state of the system, we must first minimize, with respect to the lattice spacing $a$, the free energy corresponding to a given class of solutions, and later on compare the optimal energies of each class of solutions, identifying the lowest one.

A typical profile of the free energy as a function of the lattice spacing is shown, for the skyrmion lattice class of solutions, in Fig. 3 (a). The free energy potential $\tilde{F}(a)$ has a distinct minimum, corresponding to the optimal spacing of a triangular skyrmion lattice. The solutions shown in Fig. 2 are calculated for such free energy minima.

Although the value of the optimal spacing can be, in principle, extracted from this curve alone, one must be aware that its precise determination requires the evaluation of the potential at many points in $a$ around the minimum, due to the quadratic dependence of $\tilde{F}(a)$. On the other hand, there is a far more efficient procedure to extract the optimal spacing, and it is based on the virial theorem. Even more importantly, the virial theorem method allows one to discern whether or not the optimal spacing solution for a given lattice is truly stationary with respect to variations in the lattice geometry.

II.3 Virial theorem

In order to formulate the virial theorem for our problem, we rescale the coordinates inside the spatial integrals via a dimensionless parameter $\gamma$, writing $\vec{r}=\gamma\vec{r}^{\,\prime}$. This leads to the following rescaling of the dimensionless free energy, Eq. (II.1),

[TABLE]

where $\tilde{F}_{ex}$ is the ferromagnetic exchange term, $\tilde{F}_{DM}$ is the DM term, $\tilde{F}_{B}$ is the Zeeman term, and $\vec{M}_{\gamma}(\vec{r}^{\,\prime})=\vec{M}(\gamma\vec{r}^{\,\prime})$. Since the above rescaling is just the result of a change of variables in the integrals, Eq. (24) holds for any value of $\gamma$. Hence, the right-hand side of Eq. (24) must be independent of $\gamma$ and, accordingly, its derivative with respect to $\gamma$ must vanish,

[TABLE]

We can now use the arbitrariness in the value of $\gamma$ and set $\gamma=1$ in the above equation, in order to obtain the following relation between free energy terms,

[TABLE]

where we have used $\left.\frac{\partial\vec{M}_{\gamma}}{\partial\gamma}\right|_{\gamma=1}=(\vec{r}\cdot\nabla)\vec{M}(\vec{r})$. The above equation holds for any magnetization profile $\vec{M}(\vec{r})$, which need not be a minimum of the free energy functional. However, when $\vec{M}(\vec{r})$ represents a stationary solution of the Euler-Lagrange equations, the additional condition

[TABLE]

holds. Therefore, when evaluated on a genuinely stationary solution $\vec{M}_{\rm stat}(\vec{r})$, Eq. (26) reads

[TABLE]

This is the virial theorem Doria et al. (1989) formulated for our problem, Bogdanov and Hubert (1994b); Leonov et al. (2016) which can be rewritten in terms of the following “virial ratio”,

[TABLE]

satisfied by stationary magnetization profiles.

In Fig. 3 (b) we plot the virial ratio (29) against the lattice spacing. As expected, the minimum of the free energy potential coincides precisely with the value of spacing for which we have ${\tilde{F}_{DM}}/{\tilde{F}_{ex}}=-2$, indicating that the solution is a genuinely stable configuration. Furthermore, we clearly see that the virial ratio has a linear dependence on the lattice spacing, and this allows for a much more efficient and precise numerical determination of its optimal value $a_{\rm opt}$, which we can now define through the condition $({\tilde{F}_{DM}}/{\tilde{F}_{ex}})|_{a_{\rm opt}}=-2$, rather than through the minimum of the potential.

II.4 Phase Diagram and Critical Points

The free energy can be expressed in terms of the Fourier components via

[TABLE]

This expression allows for a fast and precise evaluation of the free energy without the need to evaluate numerical derivatives.

In the following, we discuss solutions minimizing the free energy for different external magnetic fields, starting from $\beta=0$. In order to get a better understanding of the transitions, it is also instructive to include the free energy functional of metastable solutions next to the stable ones. This is shown in Fig. 4, where energies of metastable helices, skyrmions, and ferromagnetic state are plotted alongside the stable solutions.

Figure 4 also shows the critical values of $\beta$ for phase transitions. In our units the critical field for the helical-skyrmion phase transition is $\beta_{c_{1}}=0.46$ and the critical field for the skyrmion-ferromagnet phase transition is $\beta_{c_{2}}=1.59$. This is in excellent agreement with the experimentally observed phase diagram of the 2D $\rm{Fe_{1-x}Co_{x}Si}$ chiral magnet,Yu et al. (2010) and with previous theoretical investigations such as Monte-Carlo simulations Yi et al. (2009) and analytical approximate methods.Han et al. (2010) We find that the helical-skyrmion transition at $\beta_{c_{1}}$ is clearly first-order, as underlined by the possibility to extend, as metastable configurations, the skyrmion or spiral order into the neighboring phase (blue and red dotted lines in Fig. 4). The skyrmion-ferromagnet transition at $\beta_{c_{2}}$, instead, appears to be second-order in the continuum limit. This is suggested by the impossibility to find metastable skyrmion lattices above $\beta_{c_{2}}$ and by the diverging behavior of the inter-skyrmion distance in approaching the transition, as discussed hereafter. However, we have to keep in mind that the continuum limit is only an approximation and that in any actual material the spins originate from a crystal lattice. Hence, for real materials the second-order transition could easily turn into a weakly first-order one.

The lattice spacing that provides an energetically optimal configuration varies with the external magnetic field, a behavior that recalls the analogy with Abrikosov vortex lattices in type-II superconductors. However, while the lattice spacing of a superconducting vortex lattice is a monotonic function of the external field, scaling as the inverse square root of the magnetic induction and being determined by the magnetic flux quantization,Abrikosov (1957) the optimal spacing of a skyrmion lattice, shown in Fig. 5, exhibits a non-monotonic dependence upon the external field. Starting from the helical-skyrmion transition, the top panels of Fig. 2 indicate that in the vicinity of $\beta_{c_{1}}$ skyrmions are close packed, so that the lattice spacing is roughly twice the skyrmion radius. Moreover, as already discussed in Section II.2, in this region a given skyrmion is characterized by approximately an equal amount of negative and positive $M_{z}$, resulting in a zero net magnetic moment. If we now increase the magnetic field, the area of negative magnetization around the core of the skyrmion starts to shrink with respect to the positive one, giving the skyrmions a finite net magnetic moment along the field direction, and at the same time reducing the skyrmion radius. Since skyrmions at this stage are still close packed, the slightly smaller skyrmion radius results in a smaller lattice spacing, which reaches a minimum at $\beta=0.82$. Upon further increasing the magnetic field, the skyrmion radius slowly continues to decrease. However, skyrmions cease to be close packed and become more and more isolated, separated from each other by increasingly large ferromagnetic domains, as shown in the bottom panels of Fig. 2. Accordingly, the skyrmion lattice spacing increases rapidly when moving towards the ferromagnetic phase. Exactly at the $\beta_{c_{2}}$ transition, the behavior shown in Fig. 5 seems to indicate a divergent lattice spacing, signaling a continuous vanishing of the skyrmion density. This scenario would imply a second-order phase transition, as mentioned earlier. Unfortunately, with our calculations we are not able to assert this with complete certainty due to the finite truncation of the Fourier series required by the numerical implementation. In fact, Eq. (10) shows that, for an increasingly large lattice spacing $a$, the discrete $\vec{k}_{mn}$ vectors defining the Fourier series become closer and closer to each other, so that to achieve the same spatial resolution one has to include more and more $\vec{k}_{mn}$ vectors. In our numerical implementation, instead, we limited the number of $\vec{k}_{mn}$ vectors to a $25\times 25$ grid.

In addition to the skyrmion lattice spacing, in Figure 5 we present the field dependence of $\lambda_{H}=\frac{2\pi}{k_{H}}$, the period of a helix. At zero magnetic field one can prove analytically that $k_{H}=2$ (in dimensionful units, $k_{H}=2\kappa=D/J$), so that the optimal value for a helix is Dzyaloshinsky (1958); Dzyaloshinskii (1964) $\lambda_{H}=\pi$. From the low-$\beta$ branch of the curve in Fig. 5 we can see that, at $\beta=0$, $\lambda_{H}$ obtained from our calculation is exactly equal to $\pi$, showing that our calculations match the analytical predictions. As long as the field is turned on, the helix gets deformed, rotating faster in the regions where the magnetization is opposite to the field direction, and slower where the magnetization points along the field, so that the system can gain some Zeeman energy from the latter regions. Overall, this leads to a larger period of the helix, which thereby increases with the external field.

III METASTABLE SKYRMION LATTICES

In addition to the stable triangular skyrmion lattice phase illustrated above, we shall now discuss a few examples of alternative metastable skyrmion lattice geometries, and a possible way to stabilize them. If for a stable triangular lattice the virial ratio, and specifically its linear dependence on the lattice spacing, provided a numerically efficient and accurate tool to find the optimal spacing, in the case of metastable solutions it plays an even more important role. In fact, for lattice geometries different from the triangular one, there is no guarantee that the solution obtained by minimizing the free energy with respect to the lattice spacing represents a true metastable solution. Changing the lattice spacing is indeed only one of the possible ways to vary the underlying lattice structure, the others being a change in the ratio and in the angle between the unit vectors of the lattice—namely, a change in the lattice geometry. As discussed in Section II.3, any truly stationary solution of the Euler-Lagrange equations, including metastable skyrmion lattices, should fulfill the virial theorem. Hence, we can use this criterion to distinguish true metastable solutions from those that are local minima of the free energy only within the subset of solutions with a given lattice symmetry.

The first type of metastable skyrmion lattice that we are going to discuss is a square lattice, shown in the top panels of Fig. 6. This pattern can be realized, as an alternative to the triangular lattice, in the whole skyrmions phase between $\beta_{c_{1}}$ and $\beta_{c_{2}}$, and it represents a genuine metastable solution in the sense that it is a local minimum of the free energy functional with respect to arbitrary variations of the magnetization pattern. This is confirmed by the virial ratio, which takes the value of $-2$ at the minimum of the free energy as a function of the lattice spacing. In the phase diagram of Fig. 4, the energy of the square lattice (not shown to avoid an overburdening of the plot) is only slightly higher than the ground state energy, lying between the triangular (blue full curve) and honeycomb (black double-dot-dashed curve) skyrmion lattice energies. In approaching the skyrmion-ferromagnet transition, all skyrmion energies merge together at $\beta_{c_{2}}$. Such a behavior is precisely what one would expect from a second-order phase transition, characterized by a vanishing skyrmion density (or, equivalently, a diverging lattice spacing). Indeed, it is clear that for nearly isolated skyrmions the details of the lattice symmetry become unimportant and the energy of the system is controlled only by the skyrmion density.

The second type of metastable phase that we would like to present is, instead, a peculiar metastable solution that was found while investigating the system, namely a honeycomb skyrmion lattice, depicted in the bottom panels of Fig. 6. As it is well known for honeycomb lattices, the underlaying Bravais lattice for this solution is still a triangular lattice, however the unit cell contains now two skyrmions and the lattice spacing is no longer the distance between the centers of nearest skyrmions, but the distance between the centers of honeycombs. In Fig. 7 we show the free energy of a honeycomb skyrmion lattice as a function of the lattice spacing for a given $\beta$. As in the case of a triangular lattice (shown in the same figure for comparison), the free energy has a clear minimum at a certain optimal spacing. However, in stark contrast with the triangular and square lattice solutions presented earlier, the honeycomb lattice does not fulfill the virial theorem at the minimum of the corresponding free energy. This behavior entails a strong message, telling us that this solution cannot be a true local minimum of the free energy functional. Indeed, the free energy of such a configuration does have a local minimum as a function of the lattice spacing, but only as long as one considers lattice configurations with trigonal symmetry. Instead, if the geometry of the lattice is allowed to change, this solution can relax into a lower energy configuration. At first, one might think that configurations with high symmetries, such as the square or trigonal one, should be protected against deformations of the lattice geometry, so that the corresponding magnetization patterns should always be at least stationary points of the free energy functional. However, this is true only for simple Bravais lattices such as the square or triangular lattices, but no longer holds in the case of Bravais lattices with a basis such as the honeycomb one. In the latter case, for example, one realizes that the honeycomb lattice can be smoothly turned into a square lattice, which has a lower energy, by modifying the angle and ratio between the unit vectors in a certain direction. A modification of the angle and ratio in the opposite direction, instead, leads to a very different configuration with a higher energy. The honeycomb lattice would then relax into the square lattice if the lattice symmetry is left unconstrained. On the other hand, it still represents a valid metastable solution if the space of solutions is constrained, by some external mechanism, to the subset of configurations with trigonal symmetry.

As in the case of a square skyrmion lattice, the free energy of this peculiar metastable solution lies just above the free energy of the stable ground state configuration and, in approaching $\beta_{c_{2}}$, the difference between them becomes almost negligible. Hence, while this type of solution has not been experimentally observed in pure skyrmion materials, there is hope that the honeycomb skyrmion lattice might be stabilized by some sort of magnetic perturbation that establishes an underlying periodic structure with trigonal symmetry. In particular, by looking at the spin texture of this solution, we realize that skyrmions are arranged so as to surround and enclose large ferromagnetic domains. These domains are themselves arranged in a triangular lattice with vertices located at the centers of honeycombs. Since we already know that ferromagnetic domains are stabilized by a large magnetic field, a promising way of stabilizing this honeycomb skyrmion lattice is to make the spins experience an “external” magnetic field that is no longer constant, but modulated periodically so that it is largest at the vertices of a triangular lattice.

A noteworthy example of such a spatial pattern is given by the microscopic magnetic field of a type-II superconductor in the Abrikosov vortex phase, Abrikosov (1957) where the magnetic field is largest at the vortex cores and weaker in the interstitial superconducting region between vortices. Hence, provided that the length and magnetic field scales of the skyrmion material and type-II superconductor are compatible, we argue that a promising way to stabilize the honeycomb skyrmion lattice is to realize a bi- or tri-layer consisting of a non-centrosymmetric ferromagnet and a strong type-II superconductor, electromagnetically coupled, where the Abrikosov vortices would be pinned at the centers of the skyrmion honeycombs. Schematic image of such a system can be found in Fig. 8. A similar heterostructure has been very recently proposed, albeit with a different purpose, in Ref. Dahir et al., 2019, of which we became aware only after submitting our manuscript. There, the authors suggest that the coupling of a chiral ferromagnet with a type-II superconductor, characterized by an intertwined lattice of superconducting vortices and anti-vortices, might stabilize an ordinary triangular skyrmion lattice even in the absence of an external magnetic field, i.e., in the region where helical order is otherwise thermodynamically stable.

IV SUMMARY

In this paper we have introduced a numerically fast and stable approach for modeling magnetic systems by combining a Fourier-transform method with a virial theorem. We have demonstrated the accuracy of the method by applying it to a magnetic skyrmion lattice, reproducing the optimal skyrmion spacing and critical fields. Although we have considered only the case of a DM interaction induced by a Dresselhaus spin-orbit coupling, a Rashba type DM interaction, as well as easy-axis or easy-plane magnetic anisotropies, can be straightforwardly implemented in our method without incurring any additional numerical cost. Hence, we believe that our method should be perfectly suitable to describe also the recently predictedLin et al. (2015); Rowland et al. (2016) and observedYu et al. (2018) skyrmion fractionalization into meron-antimeron lattices.

We then have discussed a novel metastable honeycomb lattice solution, and suggested to combine a skyrmion lattice with a superconducting vortex lattice in order to stabilize this solution. The new phase arranges superconducting vortices in the inter-skyrmion regions, whereas skyrmions are placed in the regions where the local magnetic field is suppressed due to superconductivity. Both of these requirements can be fulfilled in a state hosting a skyrmion honeycomb lattice in combination with a triangular Abrikosov vortex lattice. This system would also lead to interesting potential ways of manipulating skyrmions, which could be moved around by drifting the Abrikosov vortices.Rabinovich et al. (2018); Reichhardt and Reichhardt (2018); Becerra and Milošević (2017) Moreover, the inhomogeneous magnetic pattern of skyrmions could be able to induce odd-frequency spin-triplet correlations in the superconductor already in a bilayer system, without the need of two separate ferromagnetic layers as in the case of superconducting spin valves.Kalenkov et al. (2011); Yokoyama and Linder (2015)

Finally, the proximity between skyrmion textures and either s-wave or p-wave superconductors is a promising tool to generate and manipulate topologically non-trivial excitations such as Majorana bound states and tuneable Weyl points.Hals et al. (2016); Yang et al. (2016); Zhang et al. (2016); Pershoguba et al. (2016); Pöyhönen et al. (2016); Takashima and Fujimoto (2016)

Acknowledgements.

ME acknowledges support by EPSRC grant EP/N017242/1.