Skyrmion Lattices in Electric Fields

TL;DR

This paper investigates how electric fields influence skyrmion lattices in insulating compounds, revealing that electric fields can stabilize or destabilize skyrmions through weak magnetoelectric coupling, which is crucial for data storage applications.

Contribution

It provides a theoretical analysis of electric field effects on skyrmion lattices, including elastic and inelastic responses, and considers the role of cubic anisotropy in stabilization.

Findings

Electric fields can shift skyrmion lattice energy positively or negatively.

Magnetoelectric coupling effects are captured through elastic and inelastic responses.

Electric field manipulation can enable creation and erasure of skyrmion arrays.

Abstract

This paper studies the influence of electric fields on the skyrmion lattice (SkL) in insulating skyrmion compounds with weak magnetoelectric (ME) coupling. The ME coupling mechanism is an interaction between the external electric field $E$ and local magnetization in the sample. Physically, the $E$-field perturbs the spin modulation wave vectors resulting in the distortion of the SkL and the $E$-field induced shift in energy. Due to the relativistic smallness of ME coupling, the important physics is captured already in the elastic ($\propto$$E$) and inelastic ($\propto$$E^2$) responses. In this spirit, the effect of the fourth-order cubic anisotropy responsible for stabilization of the skyrmion phase is taken into account perturbatively. The shift in energy can be either positive or negative, -- depending on the direction of electric field, -- thus stabilizing or destabilizing the SkL phase. Understanding the $E$-field energetics is important from the viewpoint of creation (writing) and destruction (erasing) of skyrmion arrays over the bulk, which is paramount for developing skyrmion-based logical elements and data storage.

Equations

$$\begin{array}{r}\begin{array}{c}\mathbf{S}(\mathbf{r})=\mathbf{m}+\mu \sum _{{\mathbf{Q}}_n}{\mathbf{S}}_{{\mathbf{Q}}_n}{e}^{i{\mathbf{Q}}_n\mathbf{r}+i{\phi }_n}+{\mathbf{S}}_{{\mathbf{Q}}_n}^{\ast }{e}^{-i{\mathbf{Q}}_n\mathbf{r}-i{\phi }_n},\end{array}\end{array}$$

$$\mathcal{H}={\mathcal{H}}_{JDh}+{\mathcal{H}}_A+{\mathcal{H}}_{\alpha E},$$

$${\mathcal{H}}_{JDh}=J(\nabla \mathbf{S}{)}^2+D\mathbf{S}\cdot (\nabla \times \mathbf{S})-\mathbf{h}\cdot \mathbf{S}$$

$${\mathcal{H}}_A=A(S_x^4+S_y^4+S_z^4).$$

$${\mathcal{H}}_{\text{ae}}=\alpha E_x{S}_y{S}_z+\text{cyclic permutations},$$

$$\begin{split}\begin{gathered}\mathcal{R}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{ccc}1&0&1\\ 1&0&{\text{-}}1\\ 0&\sqrt{2}&0\end{array}\right).\end{gathered}\end{split}$$

$$\begin{array}{r}\begin{array}{c}{\mathcal{H}}_A=A\left(\frac{1}{2}{S}_{x^{\prime }}^4+S_{y^{\prime }}^4+3S_{x^{\prime }}^2{S}_{z^{\prime }}^2+\frac{1}{2}{S}_{z^{\prime }}^4\right),\end{array}\end{array}$$

$$\begin{array}{r}\begin{array}{c}{\mathcal{H}}_{\text{ae}}=\frac{\alpha E}{2}\left(S_{x^{\prime }}^2+2\sqrt{2}{S}_{x^{\prime }}{S}_{y^{\prime }}-S_{z^{\prime }}^2\right),\end{array}\end{array}$$

$$\begin{array}{r}\begin{array}{c}{\mathbf{S}}^4(\mathbf{r})={\left[\mathbf{s}(\mathbf{r})+\mathbf{m}\right]}^4=...+4{\mathbf{s}}^2(\mathbf{r})\mathbf{s}(\mathbf{r})\cdot \mathbf{m}+...,\end{array}\end{array}$$

$$\mathbf{m}⟨{\mathbf{s}}^2(\mathbf{r})\mathbf{s}(\mathbf{r})⟩$$

$$\times {\mu }^3\delta ({\mathbf{k}}_1+{\mathbf{k}}_2+{\mathbf{k}}_3),$$

$$W_0=⟨J{\left[\nabla \mathbf{S}(\mathbf{r})\right]}^2+D\mathbf{S}(\mathbf{r})\cdot \left[\nabla \times \mathbf{S}(\mathbf{r})\right]⟩=\sum _{\mathbf{k}}{\mathbf{S}}_{\mathbf{k}}^{†}{\hat{\mathcal{H}}}_0{\mathbf{S}}_{\mathbf{k}},$$

$$\begin{split}\begin{gathered}\hat{\mathcal{H}}_{0}=\left(\begin{array}[]{ccc}Jk^{2}&\text{-}iDk_{z}&iDk_{y}\\ iDk_{z}&Jk^{2}&\text{-}iDk_{x}\\ \text{-}iDk_{y}&iDk_{x}&Jk^{2}\end{array}\right).\end{gathered}\end{split}$$

$$\begin{split}\begin{gathered}\hat{\mathcal{H}}_{0}=\left(\begin{array}[]{ccc}Jk^{2}&0&iDk_{y}\\ 0&Jk^{2}&\text{-}iDk_{x}\\ \text{-}iDk_{y}&iDk_{x}&Jk^{2}\end{array}\right),\end{gathered}\end{split}$$

$$|{\mathbf{S}}_{\mathbf{k}}^{(0)}⟩$$

$$|{\mathbf{S}}_{\mathbf{k}}^{(1)}⟩$$

$$|{\mathbf{S}}_{\mathbf{k}}^{(2)}⟩$$

$${\epsilon }_{\mathbf{k}}^{(0)}=J{k}^2-Dk,{\epsilon }_{\mathbf{k}}^{(1)}=J{k}^2,{\epsilon }_{\mathbf{k}}^{(2)}=J{k}^2+Dk,$$

$$\begin{split}\begin{gathered}|\mathbf{S}^{(0)}_{{\mathbf{k}}}\rangle=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}\text{-}i\hat{k}_{y}\\ i\hat{k}_{x}\\ 1\end{array}\right),\end{gathered}\end{split}$$

$$k_0=\frac{D}{2J},$$

$${\epsilon }_{{\mathbf{k}}_0}^{(0)}=-D{k}_0/2,{\epsilon }_{{\mathbf{k}}_0}^{(1)}=D{k}_0/2,{\epsilon }_{{\mathbf{k}}_0}^{(2)}=3D{k}_0/2.$$

$$\mathbf{S}(\mathbf{r})\simeq \mu \sum _{i=1}^3\left({\mathbf{S}}_{{\mathbf{Q}}_i}{e}^{{\mathbf{Q}}_i\mathbf{r}}+c.c.\right)+m{\hat{e}}_z.$$

$$\begin{array}{r}\begin{array}{c}|\mathbf{S}(\mathbf{r})⟩=\mu \sum _{\mathbf{k}}^{\{\pm {\mathbf{Q}}_1...\pm {\mathbf{Q}}_3\}}|{\mathbf{S}}_{\mathbf{k}}^{(0)}⟩e^{i\mathbf{k}\mathbf{r}+i{\phi }_{\mathbf{k}}}+m|{\mathbf{S}}_0⟩,\end{array}\end{array}$$

$$\begin{array}{r}\begin{array}{c}⟨S_x(\mathbf{r})⟩=0,⟨S_y(\mathbf{r})⟩=0,⟨S_z(\mathbf{r})⟩=m,\end{array}\end{array}$$

$$⟨S_x^2(\mathbf{r})⟩$$

$$⟨S_y^2(\mathbf{r})⟩$$

$$\begin{array}{r}\begin{array}{c}⟨S_z^2(\mathbf{r})⟩=m^2+\sum _{\mathbf{k}}^{\left\{{\mathbf{q}}_1...{\mathbf{q}}_6\right\}}{S}_{\mathbf{k}}^z{S}_{-\mathbf{k}}^z=m^2+3{\mu }^2.\end{array}\end{array}$$

$$⟨S_{\alpha }{S}_{\beta }⟩=0,\alpha \ne \beta ,$$

$$⟨S_{\alpha }{S}_{\beta }⟩=\frac{3}{2}{\mu }^2{\delta }_{\alpha ,\beta }[1+(1+2m^2/3{\mu }^2){\delta }_{\alpha ,z}],$$

$$⟨{\mathbf{S}}^2⟩=m^2+6{\mu }^2=1.$$

$$W_0^{(h)}$$

$$=\sum _{\mathbf{k}}^{\left\{{\mathbf{q}}_1...{\mathbf{q}}_6\right\}}{\mathbf{S}}_{\mathbf{k}}^{†}{\hat{\mathcal{H}}}_0{\mathbf{S}}_{\mathbf{k}}-hm,$$

$$\begin{array}{r}\begin{array}{c}{W}_0^{(h)}=-3D{k}_0{\mu }^2-hm.\end{array}\end{array}$$

Full text

Introduction

Magnetic skyrmions are quasiparticles built from topologically protected vortices of spins in helimagnets. They were theoretically predicted by Bogdanov in 1989 as thermodynamically stable states of chiral magnets,Bogdanov1989 and discovered twenty years after in the form of a hexagonal skyrmion lattice in the helimagnetic conductor MnSi.Binz2006 ; Muhlbauer2009 The magnetic skyrmions gained their name after T. Skyrme, who introduced such field configurations in the context of low-energy physics of mesons and baryons Skyrme1961 to explain the stability of particles through topological protection against continuous field transformations. Currently, skyrmions in chiral magnets are compelling due to their nanoscale size and the mentioned topological protection: They are envisaged as promising information carriers, and the skyrmion racetrack memory was recently proposed, see e.g. Refs. Fert2013 ; Nagaosa2013 ; Racetrack1 ; Racetrack2 ; Racetrack3

As quasiparticles, skyrmions can form a crystalline state, and solid-state concepts, such as symmetry breaking, order parameter, elementary excitations etc., can be applied, therefore opening doors for simple but very efficient theoretical models. The skyrmion crystal was observed both in the reciprocal and real space in the form of a two-dimensional hexagonal lattice, Muhlbauer2009 ; Tonomura2012 as is sketched in Figure 1. The hallmark of the skyrmion lattice (SkL) as seen by SANS (small angle neutron scattering) is the appearance of a six-fold pattern in reciprocal space.Muhlbauer2009 ; Munzer2010 ; Tokunaga2015 ; Adams2012 ; White2014 The thermodynamical stability of skyrmion lattices is improved experimentally by continuous tuning of coupling parameters with hydrostatic pressure Levatic2016 or by perturbing the skyrmion lattice with electric fields. Okamura2016 The interest in the skyrmion lattices and their stability is justified by the solid-state concept that the lattice could be melted into individual skyrmions.

Manipulation and control of skyrmions have become an active topic of skyrmionics. Recent experiments succeeded manipulation of skyrmions with moderate electric fields, electric currents , and thermal gradients. White2012 ; White2014 ; Jonietz2010 ; Yu2012 ; Jiang2015 ; Everschor2012 ; Mochizuki2014 ; Watanabe2016 To avoid the Ohmic heating effects which are undesirable for electronics, the application of moderate electric fields to insulating skyrmion-host compounds (such as ${\text{Cu}_{2}\text{O}\text{Se}\text{O}_{3}}$ Seki2012 ) is potentially more advantageous for the current-driven devices. These observations motivated theoretical proposals for creation (“writing”) skyrmions in insulating helimagnets with the help of external electric field, Mochizuki2015 ; Mochizuki2016 ; Okamura2016 ; Watanabe2016 and subsequent electric-field guiding. Upadhyaya2015 To date, a key experimental challenge in this direction is stabilization and control of the skyrmion lattice by an electric field. Therefore, there is an urge for simple theories of skyrmion lattice response to the E-field. In particular, what is of interest is the shift of skyrmion lattice energy in electric field, with subsequent stabilization of the skyrmion lattice.Okamura2016

In the present paper we discuss the physics of the skyrmion lattice upon application of $E$-field to an insulating skyrmion-host compound, such as ${\text{Cu}_{2}\text{O}\text{Se}\text{O}_{3}}$. In this material, the effect of magnetoelectric coupling arises due to a hybridization mechanism originating from relativistic spin-orbit interaction (see Refs. Seki2012 ; Liu2013 ; Nagaosa2007 ; Belesi2012 ), which gives rise to an electric dipole moment which can be expressed in terms of the local spin variables. The strength of the effect is however relativistically small, which allows us to build an accurate perturbation theory in $E$-fields.

In this paper, we treat the effect of electric field on skyrmion lattices in the two first orders of perturbation theory. The skyrmion lattice in electric fields becomes slightly distorted, and we introduce the elastic and inelastic distortion vectors. The shift of SkL energy comes from expectation values of the electromagnetic coupling and anisotropic contributions.

This paper is organized as follows. In Sec. I, we write down the effective coarse-grained energy functional, and make a rotation of quantization axes, as demanded by the experimental orientation of magnetic field. In Sec. II, we consider the mean-field treatment for the SkL energy density with the magnetic field applied. In Sec. III, we calculate the distortion of the skyrmion lattice by the electric field in an elastic approximation. In Sec. IV, we first calculate the mean-field expectation value of the electromagnetic coupling term, also taking into account the distortion-induced anisotropic contributions. To be consistent in the first two orders of perturbation expansion, we introduce the inelastic distortion vector, and write down all the terms in first and second orders in dimensionless electric field $\ae$. In the concluding section we highlight the main results and discuss the limitations of the model and possible applications.

I Energy density in a coarse-grained model

The Skyrmion Lattice (SkL) is a long-range-order spin configuration which can be visualized as a hexagonal array of vortices (see Fig. 1). Experimentally, the hallmark of the SkL phase is appearance of a six-fold reflection pattern in reciprocal space, as sketched in Fig. 2, each of the wave vectors are rotated by $2\pi/3$ (see e.g. Refs. Muhlbauer2009 ; White2014 for SANS patterns). In this study, we describe the skyrmion lattice by a coarse-grained magnetization $\mathbf{S}({\mathbf{r}})$, which can be built on the three $\mathbf{Q}$-vectors (Fig. 2). With a good accuracy,Muhlbauer2009 ; Binz2008 the SkL phase can be approximated by the multispiral spin structure,

[TABLE]

where $\mathbf{m}\equiv\langle\mathbf{S}({\mathbf{r}})\rangle$ is a uniform magnitization, with (spatial) average defined as $\langle...\rangle=\int\frac{dV}{V}(...)$ throughout the study, and $\mu$ is the weight of the SkL helical modulations. The sum in (1) runs over the “3Q-structure” (Fig. 2), the relative phases $\varphi_{n}$ in (1) are important for minimization of the SkL energy. The expectation of energy density in the coarse-grained model is given by calculating the spatial average $\langle\mathcal{H}\rangle$ with spin function

[TABLE]

where the helimagnetic term

[TABLE]

takes into account Heisenberg interaction ($J$), Dzialoshinskiy-Moriya interaction ($D$) and the Zeeman coupling to the external magnetic field $\mathbf{h}$. This coarse-grained model works when the following hierarchy of energies is respected: (1) the strongest is Heisenberg exchange parameter $J$ which favors the ferromagnetic alignment; (2) the Dzialoshinskiy-Moriya interaction (DMI) is slightly tilting two adjacent spins thus resulting helical modulations. (3) The cubic magnetocrystalline anisotropy $\mathcal{H}_{A}$ is considered the weakest in this hierarchy. Finally, the weak magnetoelectric term $\mathcal{H}_{\alpha E}$ is perturbatively small as estimated further. The spatial modulation of the ordered phase (the wavelength of the helices in the ground state) is given by a wave length of order $\lambda/a\sim J/D\gg 1$, i.e. much larger than the crystal lattice parameter $a$: for example, in ${\text{Cu}_{2}\text{O}\text{Se}\text{O}_{3}}$, $\lambda=630\,\text{\AA}$, $a=8.9\,\text{\AA}$, which gives $\lambda/a\sim 70$. In such case, the magnetization on neighboring lattice sites is varying very slowly and the physics of the system in the ordered phase is appropriately described by a continuous-limit model as assumed in (3).

In this study, we consider the fourth-order anisotropy as it represents the essential physics of the problem by stabilizing the SkL phase.White2014 ; Muhlbauer2009 The symmetry of ${\text{Cu}_{2}\text{O}\text{Se}\text{O}_{3}}$ is described by the $P2_{1}3$ space group, which allows a fourth-order magneto-crystalline anisotropy of the form $A_{1}(S_{x}^{4}+S_{y}^{4}+S_{z}^{4})+A_{2}(S_{x}^{2}S_{y}^{2}+S_{y}^{2}S_{z}^{2}+S_{z}^{2}S_{x}^{2})$, which can be reduced to

[TABLE]

Indeed, proceeding to the unitary parametrization $\mathbf{S}/|\mathbf{S}|=(\sin\theta\cos\psi,\sin\theta\sin\psi,\cos\theta)$, one obtains $S_{x}^{2}S_{y}^{2}+S_{y}^{2}S_{z}^{2}+S_{z}^{2}S_{x}^{2}=-\frac{1}{2}(S_{x}^{4}+S_{y}^{4}+S_{z}^{4})+\frac{1}{2}$, thus $A=A_{1}-A_{2}/2$.

The magneto-electric coupling arises due to the weak $p$-$d$ hybridization mechanism(see Refs. Seki2012 ; Liu2013 ; Nagaosa2007 ; Belesi2012 ), which gives rise to an electric dipole moment $\mathbf{P}=\alpha_{\lambda}(S_{y}S_{z},S_{z}S_{x},S_{x}S_{y})$, i.e. the electric dipole moment coupled to the spin variables $S_{x},S_{y},S_{z}$. Therefore, in external electric fields the ordered phase is perturbed by $-\mathbf{P}\cdot\mathbf{E}$, or

[TABLE]

where $\mathbf{E}=(E_{x},E_{y},E_{z})$ is the external electric field and for simplicity we absorbed the minus sign into $\alpha=-\alpha_{\lambda}$. For ${\text{Cu}_{2}\text{O}\text{Se}\text{O}_{3}}$ the strength of magneto-electric coupling is estimated as $|\alpha|\sim 10^{-33}\,\text{Jm}/\text{V}$, see Ref.Omrani2014 .

The above expressions (3)-(5) are written in the natural frame. Experimentally, one often needs to apply magnetic and electric fields along directions where particular properties of the system are better revealed. In particular, in some $E$-field rotation experiments,White2014 the magnetic field $\mathbf{h}$ is parallel to $[1\,\bar{1}\,0]$, while the electric field is parallel to $[111]$ or $[\bar{1}\bar{1}\bar{1}]$, with $|E_{x}|=|E_{y}|=|E_{z}|\equiv E$ for simplicity of notations (thus the magnitude is $|\mathbf{E}|=\sqrt{3}E$). It is easier to carry out calculations in the rotated spin frame $(x^{\prime},y^{\prime},z^{\prime})$, with $\hat{z}^{\prime}$ set by the direction of magnetic field $\mathbf{h}$, see Fig. 3. In the case of above-mentioned geometry (Fig. 3), the transformation of the rotated spin frame is given by the rotation matrix

[TABLE]

While the helimagnetic term (3) is a Lifshitz invariant, the fourth-order anisotropy (4) is transformed under rotation (6) into

[TABLE]

The $E$-field perturbation (5) also transforms under rotation ${\mathcal{R}}$,

[TABLE]

and the electric field is directed along $\hat{e}_{x^{\prime}}+\hat{e}_{y^{\prime}}$. For simplicity, we drop the prime signs in the subsequent calculations.

We make a remark on the importance of the fourth-order magnetocrystalline anisotropy (4). The appearance of the $3{\mathbf{Q}}$-structure with ${\mathbf{Q}}_{1}+{\mathbf{Q}}_{2}+{\mathbf{Q}}_{3}=0$ with helices phased in a way to form the two-dimensional skyrmion crystalline occur only due to higher-order energy terms represented in the model. In chiral ferromagnets, the anisotropy of at least the fourth order can be considered as the source of the SkL order parameter. Indeed, we can consider the fourth order anisotropy which contains - fully or partially - the term $\mathbf{S}^{4}({\mathbf{r}})=S_{x}^{4}+S_{y}^{4}+S_{z}^{4}+2(S_{x}^{2}S_{y}^{2}+S_{y}^{2}S_{z}^{2}+S_{z}^{2}S_{x}^{2})$. Introducing now new variable without uniform magnetization $\mathbf{s}(\mathbf{r})=\mathbf{S}(\mathbf{r})-\mathbf{m}$, the fourth-order term will contain a cubic term

[TABLE]

where we have not mentioned the other terms (linear, quadratic, quartic). The spatial average of the cubic term in the left-hand side of (9) can be Fourier-transformed as

[TABLE]

where $S_{{\mathbf{k}}}$ and their phases are defined as in Eq.(1), and ${\mathbf{k}}_{1,2,3}$ can be in principle any wave vector. If we define parameters in such a way that $\mu>0$, $(\mathbf{S}_{{\mathbf{k}}_{1}}\cdot\mathbf{S}_{{\mathbf{k}}_{2}})(\mathbf{m}\cdot\mathbf{S}_{{\mathbf{k}}_{3}})>0$, the hexagonal phase will be minimized if only $\varphi_{{\mathbf{k}}_{1}}+\varphi_{{\mathbf{k}}_{2}}+\varphi_{{\mathbf{k}}_{3}}=\pi$ and ${\mathbf{k}}_{1}+{\mathbf{k}}_{2}+{\mathbf{k}}_{3}=0$. The latest condition gives the so called $3\mathbf{Q}$ structure, ${\mathbf{Q}}_{1}+{\mathbf{Q}}_{2}+{\mathbf{Q}}_{3}=0$, with the three wave vectors equirotated by $2\pi/3$ due to symmetry.Muhlbauer2009 ; Binz2006 This situation is shown on Fig. 1.

II Skyrmion lattice in zero electric field: mean-field treatment

In this section we consider the mean-field treatment of the skyrmion lattice. First, we find the single-helix eigenstates, which give rise to a modulated spin structures with wave length $\lambda=2\pi/k_{0}=4\pi J/D$. After that, we construct the skyrmion lattice order parameter, and calculate the mean-field energy of the skyrmion lattice.

We start from considering the interplay between the Heisenberg term and the DMI coupling,

[TABLE]

where in we used the Fourier transform of the spatial average to reciprocal space. Here $\hat{\mathcal{H}}_{0}$ is written in spin representation $\mathbf{S}_{{\mathbf{k}}}=(S_{\mathbf{k}}^{x},S_{\mathbf{k}}^{y},S_{\mathbf{k}}^{z})^{T}$ as a matrix operator

[TABLE]

The propagation vectors skyrmion lattice lie in the plane which is perpendicular to the magnetic field $\mathbf{h}$. Consequently, each of the six helices of the skyrmion lattice is parametrized as ${\mathbf{k}}=(k_{z},k_{y},0)$ in the rotated frame, thus the problem is effectively two-dimensional. The energy matrix (”hamiltonian”) (13) is

[TABLE]

and is diagonalized on the eigenstates

[TABLE]

where $\hat{k}_{x,y,z}=k_{x,y,z}/|\mathbf{k}|$ and we introduced ”bra” and ”ket” notations as shortcuts to write the perturbation formulas of the matrix mechanics Green1965 in a familiar way. Consistent with our previous notations, we denote $\langle\mathbf{S}_{\mathbf{k}}|...|\mathbf{S}_{\mathbf{k}}\rangle=\sum_{\mathbf{k}}\mathbf{S}_{{\mathbf{k}}}^{{\dagger}}...\mathbf{S}_{{\mathbf{k}}}$, which is just a Fourier-transform of the corresponding spatial averaging $\langle...\rangle$ as in Eq.(12).

The spectrum of matrix $\hat{\mathcal{H}}_{0}$ consists of the three equidistant energy solutions, with the energy separation $\pm Dk$, see Fig. 4a,

[TABLE]

For the positive $J$, $D$ the lowest energy solution is the helix

[TABLE]

which corresponds to the eigenvalue $\varepsilon^{(0)}_{{\mathbf{k}}}$. The modulation vector $k_{0}$ is determined by further minimization of the energy, ${\partial\varepsilon^{(0)}_{{\mathbf{k}}}}/{\partial k}=0$. Consequently, the minimum of the dispersion $\varepsilon^{(0)}_{{\mathbf{k}}}=Jk^{2}-Dk$ is satisfied if

[TABLE]

see also Fig. 3a,b. Finally, it is instructive to rewrite one more time the eigenspectrum of $\mathcal{H}_{0}$ as we use it in perturbation formulas further,

[TABLE]

As a remark, it is usually convenient to measure the energy of the system in units of $Dk_{0}$.

Next we consider the superposition of helices into a skyrmion lattice. Neglecting higher-order harmonics, the skyrmion lattice is approximated as uniform magnetization $m\hat{e}_{z}$ along magnetic field and a 3$Q$ multispiral configuration,

[TABLE]

In the ”braket” notation, the real-space skyrmion phase has a shortcut notation

[TABLE]

where $|\mathbf{S}^{(0)}_{{\mathbf{k}}}\rangle$ is given by Eq.(19) and the ferromagnetic component is $|\mathbf{S}_{0}\rangle\equiv(0\,0\,1)^{\intercal}$.

The skyrmion lattice possess some interesting properties. First, using the explicit form (22), the first-order moments are

[TABLE]

so the uniform magnetization (the ”ferromagnetic component”) $\mathbf{m}\equiv\langle\mathbf{S}({\mathbf{r}})\rangle=m\,\hat{e}_{z}$ is the only non-vanishing first-order moment. The second-order in-plane moments are

[TABLE]

while the longitudinal (z) second-order moment is

[TABLE]

One can verify that all the mixed moments vanish,

[TABLE]

so, in general, one has

[TABLE]

where $\alpha,\beta={x,y,z}$. Formula (29) is handy for further expectation value calculations. For a fixed temperature one can use normalization $\langle\mathbf{S}^{2}({\mathbf{r}})\rangle=1$, thus the following constraint holdsWhite2014

[TABLE]

Finally, it is interesting to note that in this approximation the fourth-order moment for the skyrmion lattice is surprisingly $\langle\mathbf{S}^{4}\rangle\neq\langle\mathbf{S}^{2}\rangle^{2}$, i.e. the magnetization field is “soft”.

Now we consider the skyrmion lattice in finite magnetic fields by adding the source term $-\mathbf{h}\cdot\mathbf{S}$,

[TABLE]

which however doesn’t change the mean-field eigenstates (15)-(17) as it contains no spatial derivatives. Thus the expectation value (31) is modified only in terms of elongating the ferromagnetic component $m$, but the topology of the skyrmion order parameter is not affected. The mean-field treatment of the SkL energy (31) yields

[TABLE]

In the mean-field treatment, one may minimize the total energy with respect to all the components of the magnetic moment, including $m$ and $\mu$, if a comparison between two phases is needed. For example, minimizing now Eq.(32) with constraint (30), one gets the mean-field estimates $m_{\text{MF}}=h/Dk_{0}$ and $\mu_{\text{MF}}=\sqrt{(1-m_{\text{MF}}^{2})/6}$. The corresponding energy is therefore given as $W^{(h)}_{0}=-Dk_{0}\left(m^{2}_{\text{MF}}+1\right)/2$.

Finally, we calculate perturbatively the anisotropic contribution around the mean-field. For the anisotropy of form (4) (without spatial derivatives), the wave vectors of helices are not renormalized. Therefore, the contribution of anisotropy to the energy of the skyrmion lattice is calculated by using the same order parameter (23). Decoupling the ferromagnetic contribution $\mathbf{m}$ from Eq.(7), one obtains

[TABLE]

A direct calculation within the unperturbed state Eq.(23) gives

[TABLE]

Note that the term containing $\mu^{3}$, is responsible for the minimization of the $3Q$ structure, as was discussed in Section I. For $A>0$, $m>0$, $\mu>0$, expression (34) is minimized for $\varphi_{1}+\varphi_{2}+\varphi_{3}=\pi$, which gives

[TABLE]

The energy of the skyrmion lattice in the absence of electric field is thus given by (32) and (35) in the mean-field approximation.

III Skyrmion Lattice in Electric Fields: Elastic distortion

In this section we consider the shift in the SkL energy caused by distortion of the skyrmion vectors. We start by re-writing the magneto-electric coupling in the rotated frame (8) in symmetrized form, which in units $Dk_{0}$ is simply

[TABLE]

where we have introduced the dimensionless electric field æ,

[TABLE]

which plays the role of the small parameter of the theory. To illustrate its smallness, we use typical electric fields $E=5\times 10^{6}\,\text{V/m}$, and ${\text{Cu}_{2}\text{O}\text{Se}\text{O}_{3}}$ parameters as $J=4.85\times 10^{-23}\,\text{Jm}/\text{A}$, $k_{0}=D/2J=10^{8}\,\text{m}$, and ME couplingMochizuki2015 is $\alpha\sim 10^{-14}\,\text{J}/\text{m}^{2}\text{V}$. This gives $\ae\sim 0.01$. Thus throughout this study, we build the perturbation theory in orders of $\ae^{1}$ and $\ae^{2}$, which is sufficient for describing both the symmetric and asymmetric responses in $E$-fields.

First, we consider the helix vectors in external electric field. Considering $E$-field as a small perturbation (36) on top of the $JD$-matrix (14), matrix perturbation theory givesGreen1965

[TABLE]

where $\varepsilon^{(n)}_{{\mathbf{k}}}$ are given by Eq.(21) and $|{\mathbf{S}}^{(n)}_{\mathbf{k}}\rangle$, $n=0,1,2$, are eigenstates of $\mathcal{H}_{0}$ as given by Eqs.(14-16). The perturbed helix (38) is by construction normalized on unity up to terms of order ${\mathcal{O}}(\ae^{2})$,

[TABLE]

A direct calculation for the new helix, using formulas (38), (36), (15)-(17), gives

[TABLE]

where ${\mathbf{F}}_{{\mathbf{k}}}={\mathbf{F}}(\hat{k}_{x},\hat{k}_{y})\equiv{\mathbf{F}}_{\mathbf{k}}(\phi)$ is the elastic distortion vector (here for each helix $\hat{k}_{x}=\cos\phi$, $\hat{k}_{y}=\sin\phi$, that is $\phi$ is angle between $\mathbf{k}$ and $\hat{x}$), so that

[TABLE]

with angle-dependent components

[TABLE]

The components $F^{x}_{{\mathbf{k}}}$, $F^{y}_{{\mathbf{k}}}$, $F^{x}_{{\mathbf{k}}}$ of the elastic distortion vector are $\pi$ and $2\pi$ periodic functions, and are illustrated in Fig. 6. Note that the $x$ and $y$ components of both ${\mathbf{S}}^{(0)}_{\mathbf{k}}$ and ${\mathbf{S}}^{(\ae)}_{\mathbf{k}}$ are imaginary, and $z$ component is real for both ${\mathbf{S}}^{(0)}_{\mathbf{k}}$ and ${\mathbf{S}}^{(\ae)}_{\mathbf{k}}$. Physically it means that the field-induced distortion of the skyrmion lattice does not change the orientation of the skyrmion plane. One can verify both analytically and numerically that Eqs.(42)-(43) together with (41) satisfy normalization (39).

IV Energy shift of the Skyrmion Lattice in Electric Fields

In this section, we calculate the shift of the SkL energy in the first two orders in terms of dimensionless electric field $\ae$. This shift is contributed by the expectation value of the magnetoelectric coupling, and the anisotropic contributions due to the skyrmion lattice distortion.

IV.1 Magneto-Electric Response

The first contribution to the energy shift of the skyrmion lattice in electric field comes from taking the expectation value of (36) in the appropriate order in $\ae$. In this study, we consider $\ae$ to be small, thus the second order is sufficient for capturing the essential physics in both antisymmetric (field reversion $E\to-E$ gives the energy shifts of different signs) and symmetric (field reversion $E\to-E$ gives the energy shifts of same signs) cases of interest. Using the explicit expression (40)–(44), one obtains

[TABLE]

The expression (45) is the leading contribution in case if the anisotropy of the system is small. Restoring the dimensional units, the shift of the SkL energy $W$ in electric field $E$ is therefore given by

[TABLE]

where $M_{s}$ is saturation magnetization (in $\text{A}/\text{m}$), $a$ is the lattice constant. Note that the helimagnetic term $H_{JDh}$ is not giving shift contributions to this order.

IV.2 Anisotropy response

Now, we calculate the direct contribution of anisotropy to the energy shift. The main physical mechanism here is the distortion of SkL in electric field, which results to perturbation in the anisotropic energy. To proceed, we explicitly use formulas (40)–(44) to calculate the expectation value of the anisotropic term in the new ground state. In the leading order one therefore obtains

[TABLE]

Here $\phi_{0}$ is the angle between the first skyrmion helix and $\hat{x}$ (note that due to $2\pi/3$ symmetry in SkL rotation, and $6\phi$-arguments in Eq.(47), one can take $\phi_{0}$ as the angle between any of the skyrmion helices and $\hat{x}$). The same angular dependence [third term in Eq.(47)] were reported in study White2014 , where the minimization of SkL energy in electric fields leads to the SkL rotation in real space with respect to $\hat{x}$, if six-order anisotropy is considered, however all the angular-independent anisotropic contributions were there neglected.White2014 We notice that expression (47) is dependent on the relative phases of the helices $\varphi_{i}$. The same situation occurs in (34), when the mean-field energy is dependent on phase. We take again $\cos(\varphi_{1}+\varphi_{2}+\varphi_{3})=-1$, which minimizes the mean-field energy for fixed positive $m$, $\mu$, therefore

[TABLE]

This contribution gives the linear anisotropic response to the external field. However, for the completeness of discussion, we need to take into account further terms $\mathcal{O}(\ae^{2})$, which come both from the elastic distortion of the SkL and the inelastic (quadratic in $\ae$) distortion of the SkL.

IV.2.1 Mixed Inelastic response

To consider the nonlinear anisotropic response, we re-define the perturbed helix eigenstates up to the second order, which are given by a perturbative expansion

[TABLE]

This expansion leads to re-definition of (40) by adding a nonelastic distortion of the skyrmion lattice,

[TABLE]

where ${\mathbf{G}}_{{\mathbf{k}}}={\mathbf{G}}(\hat{k}_{x},\hat{k}_{y})\equiv{\mathbf{G}}_{\mathbf{k}}(\phi)$ is the main inelastic distortion vector,

[TABLE]

The direct calculation of the inelastic distortion vector, by using (49), (36), (15)-(17), gives

[TABLE]

These dependencies are shown in Fig.6. The comparison between the elastic and inelastic distortion vectors is shown on Figure 7. Finally, for the numerical consistence of calculations, one can verify that the new ground state is normalized on unity,

[TABLE]

Therefore, the mixed elastic-inelastic response $\mathcal{O}(\ae^{2})$ is obtained by calculating the expectation value in the new basis (50). The direct calculation gives

[TABLE]

where we have introduced the following dimensionless factors

[TABLE]

The mixed elastic-inelastic contribution (54) gives the second-order correction in $\ae$ which may be important in some particular cases, for example, when the first order correction vanishes. Then, the electric field response does not depend on the field polarity ($\pm E$ is absorbed in $E^{2}$).

Therefore, we have three contributions to the shift in the SkL mean-field energy:

[TABLE]

[TABLE]

[TABLE]

Together, formulas (58)-(60) give the contribution up to the second order in electric field $E$. In case of the weak anisotropy $A\ll D^{2}/J$, the term (60) can be usually neglected. Note also that in the main-order approximation ($\ae^{1}$) the $E$-field induced shift in energy is also the shift in free energy of the skyrmion lattice for a fixed temperature near $T_{C}$.

It is therefore possible to stabilize (if $\Delta W_{\text{total}}<0$) or destabilize (if $\Delta W_{\text{total}}>0$) the skyrmion lattice by choosing the appropriate magnitudes and polarities of electric fields. This coincides with the previous idea of Mochizuki Mochizuki2016 of creating single magnetic skyrmions with external $E$-fields in multiferroic ${\text{Cu}_{2}\text{O}\text{Se}\text{O}_{3}}$. The present study is therefore a bridge towards this idea in the bulk samples, where the skyrmions usually exist in the form of long-range-ordered or partially-disordered skyrmion arrays. The experimental research into writing and erasing skyrmions with electric fields in bulk ${\text{Cu}_{2}\text{O}\text{Se}\text{O}_{3}}$ is on the final stage in our lab and will be published elsewhere.

V Summary

The present model describes the shift of the mean-field energy of the SkL, which is either positive or negative depending on the direction of the electric field. In a particular situation, when the first order terms vanish, the energy shift is of the same sign for both field polarities.

An interesting output of the calculation is the existence of the stationary points if a helix is directed in a proper way (Fig. 7). This feature comes both in the elastic and nonelastic distortions of the skyrmion lattice. However, as the SkL is constructed on the three helices, the mean-field energy of the SkL is still shifted.

Finally, we sketch the limitations of the calculation. First, this study describes the first two perturbative corrections to the mean-field energy of the skyrmion lattice in the multispiral approximation, without comparing the free energies of different possible phases (helical, conical) in the system. Second, the energy functional is taken in the quasiclassical continuous-field limit, which limits the use of the model only to the sufficiently low DMI parameter ($D/J\ll 1$) and excludes the quantum regime (low $T$). Third, the effect of electric field on critical fluctuations on top of the mean-field SkL solution is not considered as it comes as a higher-order contribution in the critical correlation length.

The main physical consequence of the phenomenon under study is the field-induced stabilization of the skyrmion phase in the bulk, which has been so far indirectly observed in ${\text{Cu}_{2}\text{O}\text{Se}\text{O}_{3}}$. Okamura2016 This mechanism if further developed would allow one either to write or erase the skyrmion array over the full sample if it is properly placed in the $H$, $T$ phase diagram, and thus opens further routes for skyrmion-based racetrack logical elements and data storage devices.

Acknowledgments. - The work was supported by the Swiss National Science Foundation, its Sinergia network Mott Physics Beyond the Heisenberg Model (MPBH). The authors would like to thank Achim Rosch and Jiadong Zang for useful discussion.