Polarizability of electrically induced magnetic vortex plasma

TL;DR

This paper investigates how electric fields can create and manipulate magnetic vortices in magnetoelectric materials, revealing a new form of topological defect plasma with temperature-dependent polarizability.

Contribution

It introduces the concept of electric field-induced magnetic vortex plasma and analytically and numerically studies its density distribution and polarizability.

Findings

Electric fields can create localized magnetic vortex structures.

The polarizability depends on temperature as lpha  1/T^{1-ta}.

The results align with Euclidean random matrix theory predictions.

Abstract

Electric field control of magnetic structures, particularly topological defects in magnetoelectric materials, draws a great attention in recent years, which has led to experimental success in creation and manipulation by electric field of single magnetic defects, such as domain walls and skyrmions. In this work we explore a scenario of electric field creation of another type of topological defects -- magnetic vortices and antivortices, which are characteristic for materials with easy plane (XY) symmetry. Each magnetic (anti)vortex in magnetoelectric materials (such as type-II multiferroics) possesses a quantized magnetic and an electric charge, where the former is responsible for interaction between vortices and the latter couples the vortices to electric field. This property of magnetic vortices opens a peculiar possibility of creation of magnetic vortex plasma by non-uniform electric fields. We show that the electric field, created by a cantilever tip, produces a "magnetic atom" with a localized spatially ordered spot of vortices ("nucleus" of the atom) surrounded by antivortices ("electronic shells"). We analytically find the vortex density distribution profile and temperature dependence of polarizability of this structure and confirm it numerically. We show that electric polarizability of the "magnetic atom" depends on temperature as $\alpha \sim 1/T^{1-\eta}$ ($\eta>0$), which is consistent with Euclidean random matrix theory prediction.

Equations

$$w=w_e+w_m+w_{me}.$$

$$w_e=\frac{P^2}{2{\chi }_e}-\mathbf{E}\mathbf{P}-\frac{{\mathbf{E}}^2}{8\pi }$$

$$w_m=\frac{\alpha }{2}\left[{\left(\frac{\partial \mathbf{M}}{\partial x}\right)}^2+{\left(\frac{\partial \mathbf{M}}{\partial y}\right)}^2\right]$$

$$w_{me}=\gamma \mathbf{P}(\mathbf{M}(\nabla \mathbf{M})-(\mathbf{M}\nabla )\mathbf{M}).$$

$$w=\frac{P^2}{2{\chi }_e}-\mathbf{E}\mathbf{P}-\frac{{\mathbf{E}}^2}{8\pi }-\gamma \mathbf{P}((\mathbf{M}\nabla )\mathbf{M}-\mathbf{M}(\nabla \mathbf{M}))+$$

$$+\frac{\alpha }{2}\left[{\left(\frac{\partial \mathbf{M}}{\partial x}\right)}^2+{\left(\frac{\partial \mathbf{M}}{\partial y}\right)}^2\right].$$

$$\mathbf{P}=\gamma {\chi }_e((\mathbf{M}\nabla )\mathbf{M}-\mathbf{M}(\nabla \mathbf{M}))+{\chi }_e\mathbf{E}.$$

$${\mathbf{M}}(\mathbf{r})=M_{0}\left(\begin{array}[]{ccc}\cos\phi(\mathbf{r})\\ \sin\phi(\mathbf{r})\\ 0\\ \end{array}\right).$$

$${\mathbf{P}}=\gamma\chi_{e}M_{0}^{2}\left(\begin{array}[]{ccc}-\partial_{y}\phi\\ \partial_{x}\phi\\ 0\\ \end{array}\right)+\chi_{e}{\mathbf{E}}.$$

$$w=\frac{1}{2}\left(\alpha M_0^2-{\chi }_e{\gamma }^2{M}_0^4\right)(\nabla \varphi {)}^2-$$

$$-{\chi }_e\gamma M_0^2\left({\partial }_x\varphi E_y-{\partial }_y\varphi E_x\right)-\frac{1}{2}{\chi }_e{E}^2.$$

$${\rho }_s=(\alpha M_0^2-{\chi }_e{\gamma }^2{M}_0^4)h,$$

$$\varphi =k\arctan \frac{y}{x}+{\varphi }_0,$$

$$q_e=2\pi \gamma {\chi }_e{M}_0^{2}h,$$

$$W=\frac{h}{2}\int \left(\alpha M_0^2(\nabla \varphi {)}^2-2{\chi }_e\gamma M_0^2\left({\partial }_x\varphi E_y-{\partial }_y\varphi E_x\right)\right)dS=$$

$$=\frac{h\alpha M_0^2}{2}\int \left(({\partial }_x\varphi -\frac{{\chi }_e\gamma }{\alpha }{E}_y{)}^2+({\partial }_y\varphi +\frac{{\chi }_e\gamma }{\alpha }{E}_x{)}^2\right)dS$$

$$W_{me}=h{\chi }_e\gamma M_0^2\int \left({\partial }_x\varphi {\partial }_y\phi -{\partial }_y\varphi {\partial }_x\phi \right)dS$$

$$W_{me}=h{\chi }_e\gamma M_0^2(\oint \left({\partial }_x\varphi {\hat{n}}_y-{\partial }_y\varphi {\hat{n}}_x\right)\phi dl-$$

$$-\int \left({\partial }_y{\partial }_x\varphi -{\partial }_x{\partial }_y\varphi \right)\phi dS),$$

$$W_{me}=h{\chi }_e\gamma M_0^2\left(\sum _{cores}2\pi k_i\phi ({\mathbf{r}}_i)-\oint {\text{displaylimits}}_{edge}\phi \nabla \varphi d\mathbf{l}\right)$$

$$W_{me}=-q_e\sum _i{k}_i\phi ({\mathbf{r}}_i).$$

$$q_m^2=\pi {\rho }_s\approx \alpha M_0^{2}h.$$

$$W_m=-2q_m^2\sum _{i<j}{k}_i{k}_j\ln \frac{r_{ij}}{a}.$$

$$W_m=(N_1-N_2{)}^2{q}_m^2\ln \frac{R}{a}-2q_m^2\sum _{n<m}{k}_i{k}_j\ln \frac{r_{ij}}{a}.$$

$$\phi (r)=-\frac{Q_{tip}}{\sqrt{r^2+H^2}}$$

$$W=2q_m^2\ln \frac{r}{a}+Q_{tip}{q}_e\left(\frac{1}{\sqrt{r^2+H^2}}-\frac{1}{H}\right).$$

$$\kappa =\frac{Q_{tip}{q}_e}{2q_m^{2}H},$$

$$x=\frac{r^2}{H^2},$$

$$h=\frac{H}{a}.$$

$$W=\kappa \left(\frac{1}{\sqrt{x+1}}-1\right)+\frac{1}{2}\ln x+\ln h$$

$$(x+1{)}^3={\kappa }^2{x}^2.$$

$$\frac{\kappa }{\sqrt{x_0+1}}<\kappa -\ln h-\frac{1}{2}\ln x_0$$

$$r_0\equiv H\sqrt{x_0}\approx \frac{Q_{tip}{q}_e}{2q_m^2}.$$

Full text

I Introduction

Efficient control of magnetic structures by electric fields might revolutionize the field of magnetic memory devices. Compared to the existing technology that switches magnetic states using electric currents (via created magnetic fields or spin-transfer torques) and thus is subject to joule heating Taniyama:2015 , the replacement of current with an electric field promises to reduce energy dissipation by the factor of 100 [Matsukura:2015, ]. Several mechanisms of the electric-field manipulation of magnetism were proposed such as: electric-field changing the carriers concentration, which mediate magnetic interactions in magnetic semiconductors; changing the magnetic anisotropy in insulator/magnetic metal heterostructures; using magnetoelectric coupling in multiferroic materials Matsukura:2015 . In the present work we will focus on the latter mechanism.

Magnetoelectric effect was first discovered “with the point of the pen” by Dzyaloshinskii Dzyaloshinskii:1959 and then this prediction was experimentally observed by Astrov in $\mathrm{Cr_{2}O_{3}}$ [Astrov:1960, ]. Then, after a fast development in 60s-70s, interest to such materials gradually faded away, because of the weakness of the observed effect. But since the experimental successes of 2000-s – discovery of many multiferroics (materials in which two or more ferroic orders, such as ferromagnetism and ferroelectricity, are present Schmid:1994 ) and “revival” of the field, magnetoelectric materials are still in the focus of both experimental and theoretical research. For recent reviews on multiferroics see Khomskii:2016 ; Fiebig:2016 ; Tokura:2014 ; Zvezdin:2012 ; Spaldin:2010 .

According to the classification of Khomskii Khomskii:2009 , there are two types of multiferroics. Type I consists of multiferroics, where both ferromagnetic and ferroelectric orders coexist, being almost not coupled to each other. While in type-II multiferroics the two orders are coupled and usually a magnetic order creates an electric one. This feature can be used backwards: applying an electric field and inducing polarization, magnetic patterns can be also created and manipulated. While this type-II feature was originally coined to multiferroics, it becomes more spread that this coupling between magnetization and electric polarization can appear in any magnetic insulator, where it is allowed by symmetry Khomskii:2016 . In 2008 Dzyaloshinskii theoretically predicted the possibility of creation and manipulation of domain walls by electric fields in conventional and weak ferromagnets Dzyaloshinskii:2008 and the movement of magnetic domain walls under the action of electric field was directly observed Logginov:2008 ; Pyatakov:2011 .

The use of domain walls and other topological defects in multiferroics Seidel:2016 for magnetic storage is particularly appealing since it provides a mechanism of very dense packing of information, while topological nature protects information from loss under the influence of external perturbations such as heating or mechanical action Parkin:2008 ; Fert:2013 . Among the bright examples of such topological defects are skyrmion lattices Tokura:2012 , as well as individual skyrmions observed in magnetoelectric materials Hsu:2016 , magnetic vortices Mostovoy:2009 ; Pyatakov:2012 and their discrete analogues Cheong:2010 ; Mostovoy:2010 ; Artyukhin:2014 , domain walls Logginov:2008 ; Pyatakov:2011 ; Meier:2015 . Usually such defects and textures are observed in thin-film magnetoelectric materials Ramesh:2007 .

In the present paper we consider theoretically a phenomenological model developed by Mostovoy Mostovoy:2006 for a type-II multiferroics and apply it to a thin-film material with easy plane symmetry. We show both analytically and numerically that in this case a strong enough locally-concentrated electric field may cause formation of a “magnetic atom” that consists of vortex “nucleus” surrounded by antivortex “shells”. This idea allows us to propose an experimental realization for 2D Coulomb plasma in a trap and make predictions about its normal modes. The topic of finite Coulomb plasma clusters, and particularly of their normal modes, draws great attention in different areas of physics (see [Melzer:2001, , Melzer:2012, ] and references therein). However, usually particles in 2D plasma clusters interact between themselves via screened or unscreened 3D Coulomb $1/r$ potential. On the contrary, systems that we consider in this work provide a unique possibility of creation of plasma with constituents interacting via 2D Coulomb $\log r$ potential.

Also presented here research has some reminiscence with quantum electrodynamics. It was proposed that space- or time- dependent electric fields can generate electron-positron pairs from vacuum or fill empty electronic shells in an atom, producing a positron Kleinert:2008 . The idea was first theoretically proposed by Sauter [Sauter:1931, ] and Schwinger [Schwinger:1951, ] and attracts considerable interest in the present days (see, for example, [Kleinert:2013, ]). But the experimental observation is still lacking: such strong electric fields can not be obtained in laboratory conditions yet; now steps are being done in this direction: Extreme Light Infrastructure project ELI , which will presumably allow to reach Schwinger limit Dunne:2009 , is currently under construction.

Another example lies in the physics of heavy atoms. For atomic numbers $Z\gtrsim 1/\alpha$ ($\alpha\approx 1/137$ is the fine structure constant) electric field becomes so strong that the nucleus can tear electrons out of vacuum and fill lower electron levels while emitting positrons Pomeranchuk:1945 ; Zeldovich:1971 . This effect was not yet observed since such heavy atoms are very unstable. However, recently a condensed-matter version of the artificial atom-like structures was realized in graphene Mao:2016 .

The present paper proposes an experimental setup for observation of analogous effects in another part of the realm of condensed matter systems, in multiferroics: we propose that vortex-antivortex pairs can be created from the ground state of a type-II multiferroic thin film by applying a strong non-uniform electric field.

The paper is organized as follows. In Sec.II we introduce our basic phenomenological model of quasi-two-dimensional type-II multiferroic-like material. In Sec.III we explore the consequences of applying strong non-uniform electric field created by a tip of cantilever and find “the first critical electric field”, when the first vortex-antivortex pair is created. In Sec.IV we present our analytical results on the number of induced vortices and vortex density distribution in the continuous model. In Sec.V we deal explicitly with the effects of discreteness of vortex distribution and estimate the principal part of temperature dependence of polarizability of artificial “magnetic atom”. Sec.VI presents our numerical results, which confirm theoretical calculations in Sec. IV and V. In Sec.VII we make some estimates for the real materials and show that the critical electric field may be quite moderate and that even weak magnetic anisotropy protects in-plane spin arrangement. Sec.VIII contains our conclusions and some ideas on experimental techniques that can be applied for observation of the described effects.

II Derivation of the model

In order to describe a 2D material with coupled magnetic and electric subsystems we use the following phenomenological model. We write the total free energy density as a sum of parts arising from electric polarization, magnetization, and magnetoelectric coupling respectively:

[TABLE]

Here

[TABLE]

is the energy of electric polarization LL8 , $\chi_{e}$ is the electric susceptibility in the absence of magnetization; we assume that there is no spontaneous polarization in the absence of magnetization, thus the $\sim P^{4}$ term is absent in the free energy (2). We write the magnetic part of the free energy as

[TABLE]

which is magnetic non-uniformity energy LL8 ; contribution $AM^{2}+BM^{4}$ that keeps $|{\mathbf{M}}|=M_{0}=const$ is implied. Since we are modeling thin-film material, we assume that the film thickness $h$ is much less than the exchange length and neglect the variations of ${\mathbf{M}}$ in the transverse to the film direction. We also note that the same kind of analysis can be done also for an antiferromagnetic model, which also possesses the vortex excitations Yanagisawa:2013 .

The most interesting energy part is the coupling energy Mostovoy:2006 , which characterizes type-II multiferroics

[TABLE]

This Lifshits-type term was introduced based on symmetry grounds; it can be present if spatial inversion symmetry is broken. One of the possible underlying microscopic mechanisms is based on the inverse Dzyaloshinskii-Moriya interaction Katsura:2005 .

Combining three contributions (2), (3), (4) we write the total free energy density of the magnetoelectric material as

[TABLE]

Now we assume that ${\mathbf{M}}$ is a primary order parameter, which can induce ${\mathbf{P}}$. So that we can find polarization by minimizing (II) with respect to $\mathbf{P}$

[TABLE]

Because of the thin film geometry of the sample, magnetization tends to be parallel to the film plane, otherwise the magnetic energy of the demagnetizing field would greatly increase; also material can possess additional easy-plane type of symmetry. Therefore we restrict the magnetization to lie in the plane

[TABLE]

Using the angle variable $\phi$ the polarization can be written as

[TABLE]

Substituting (7) and (8) to (II), we see that electric, magnetic, and magnetoelectric parts of energy combine to

[TABLE]

Assume, for a moment, that ${\mathbf{E}}=0$. In this case expression for energy $W=\int wdV$ is similar to one of the XY model $W_{XY}=\frac{1}{2}\rho_{s}\int(\nabla\phi)^{2}dS$ with an effective spin-wave stiffness

[TABLE]

where $h$ is the film thickness. Hence, in the absence of the electric field, the magnetoelectric coupling just renormalizes the spin-stiffness of the XY model. For typical values of parameters (see Sec.VII and especially text after formula (57) for details): $\alpha M_{0}^{2}$ is greater than $\chi_{e}\gamma^{2}M_{0}^{4}$ by several orders of magnitude; therefore, $\rho_{s}$ remains almost unchanged by magnetoelectric interaction and further we use $\rho_{s}\approx\alpha M_{0}^{2}h$.

Despite this, some interesting features appear. It is well known that in the XY and related models the magnetic vortices play the crucial role KT1 . Vortex placed at the origin has the form

[TABLE]

where $k$ is the winding number of the vortex, $\phi_{0}$ is an arbitrary angle. It turns out that magnetoelectric coupling affects vortices. From (8) it follows that polarization of magnetic vortex is ${\mathbf{P}}=-k\gamma\chi_{e}M_{0}^{2}{\mathbf{r}}/r^{2}$. This leads to the fact that a magnetic vortex acquires an electric charge $kq_{e}$ proportional to the winding number $k$ of the vortexMostovoy:2006 , where

[TABLE]

where $\lambda_{e}$ is a vortex charge per unit film thickness.

Consider now what happens in non-zero electric field. The $\phi$-dependent part of energy can be rewritten with the accuracy $O(\gamma)$ in the coupling constant as

[TABLE]

From this we see that the ground state in the electric field is reached when $\nabla\phi=(E_{y},-E_{x})\chi_{e}\gamma/\alpha$.

For example, a weak constant in-plane electric field modifies the ground state configuration into a cycloidal magnetic structure with a period $\lambda=2\pi\rho_{s}/\chi_{e}\gamma E$ and the wave-vector perpendicular to the electric field Most1 . At $T=0$ this magnetic structure creates polarization ${\mathbf{P}}=\chi_{e}^{2}\gamma^{2}M_{0}^{2}{\mathbf{E}}/\rho_{s}$ and gives contribution to electric susceptibility: $\chi_{cycloid}=\chi_{e}^{2}\gamma^{2}M_{0}^{2}/\rho_{s}$. We also see, that $\lambda\rightarrow\infty$ as $E\rightarrow 0$, and hence, finite-size magnetic orientational defects are not influenced by infinitesimal electric field.

When the electric field is not constant, nor weak, it is convenient to write it with the use of the electric field potential ${\mathbf{E}}=-\nabla\varphi$, so the coupling term (4) can be rewritten as

[TABLE]

Excluding singular points from the integration domain and making cuts in order that it was simply connected, we can integrate by parts and get:

[TABLE]

where the boundary integral over the edge of the sample $dl$ contains a normal to it, $\hat{n}$. Since the integration domain is simply connected and does not contain singularities, the second integral in (II) vanishes.

The first integral can be transformed using the fact that $d{\mathbf{l}}\sim(-\hat{n}_{y},\hat{n}_{x})$, so $\oint\left(\partial_{x}\phi\hat{n}_{y}-\partial_{y}\phi\hat{n}_{x}\right)\varphi\,dl=-\oint\varphi\,\nabla\phi d{\mathbf{l}}$. Integration contour includes not only outer boundary of the sample, but also all singular points of $\nabla\phi$. Since integration around the singularities is clockwise, we take out $-1$ and make it counterclockwise; $\oint\nabla\phi d{\mathbf{l}}=2\pi k_{i}$, where $k_{i}$ is the winding number (vorticity) of $i$-th singular point. Integrals along cuts give [math], since $\partial_{i}\phi$ is single-valued (though $\phi$ itself is not). Therefore, we obtain

[TABLE]

The first term in (16) has a form $\sum{q_{i}\varphi({\mathbf{r}}_{i})}$. This means, that similarly to the case of zero electric field (12), in non-zero electric field vortices still have an electric charge proportional to еру winding number of the vortex $q_{i}=k_{i}q_{e}=k_{i}\cdot 2\pi h\chi_{e}\gamma M_{0}^{2}$. Therefore, vortices carry positive charge $q_{e}$ and antivortices carry negative charge $-q_{e}$.

However, from (16) we see that when the external electric field is applied, boundary effects become important. Consider a single vortex, which has electric charge (12) due to the magnetoelectric interaction. From the total electroneutrality it follows that the boundary of the sample also becomes charged; it acquires charge of the same magnitude but the opposite sign compared to the vortex. In a longitudinal external field, the boundary charge effectively shields the charge of the vortex, and the value of screening depends on the geometry of the sample. In this work we assume for simplicity that our sample is a disk with radius $R$. For such a disk screening of longitudinal constant electric field reduces the effective charge of vortices exactly twice Most1 (see Appendix A for further details). For axially symmetric electric fields with potential $\varphi(r)$ the integral in (16) gives just $q_{edge}\varphi(R)$ with $q_{edge}=h\chi_{e}\gamma M_{0}^{2}\oint\nabla\phi d{\mathbf{l}}$. So for magnetically neutral systems (when $\oint\nabla\phi d{\mathbf{l}}=0$) the edge charge is zero and (16) transforms to

[TABLE]

Coupling energy (16) is additive both in $\phi$ and in $\varphi$. Additivity in $\phi$ tells us that magnetoelectric coupling does not influence the vortex-vortex magnetic interaction, which comes only from non-additive term $(\nabla\phi)^{2}$ in (II). Therefore vortices interact as in XY model via $2D$ Coulomb potential KT1 $W=\pm 2q_{m}^{2}\ln(r/a)$, where $a$ is the short-distance cut-off of order of the lattice spacing and $q_{m}$ is the “magnetic charge”:

[TABLE]

Here and further we put for simplicity vorticities to be $k=\pm 1$, since vortices with higher winding numbers are unstable towards decomposition to these elementary ones.

The magnetic energy of a magnetically neutral system of $N$ vortices and $N$ antivortices can be written asChaikin_Lubensky :

[TABLE]

here $r_{ij}$ is the distance between $i$-th and $j$-th vortex cores.

For a sample with the number $N_{1}$ of vortices not equal to the number $N_{2}$ of antivortices formula (19) can be generalized to

[TABLE]

We would like to emphasize that there is also purely electrostatic Coulomb interaction between the vortex cores, but it gives contribution $\sim\gamma^{2}$, which is negligibly small compared to the interaction $\sim\gamma$ with the external electric field. Therefore, in what follows we adopt an approximation, in which a system of vortices and antivortices interact with each other only via the magnetic subsystem, and in addition, they interact with the external electric fields due to the magnetoelectric coupling.

III Breaking the magnetic vacuum with electric field: vortex - antivortex pairs creation

In this section we explore the basic idea of the paper, that a strong non-uniform electric field can break the “magnetic vacuum” of magnetoelectric material via creation of vortex-antivortex pairs.

III.1 Geometry of the system

Fig.1 shows typical experimental setup that we keep in mind (which was proposed in Refs. [Pyatakov:2012, ],[Pyatakov:2010, ]). An electric field is produced by a tip of a cantilever with a fixed point charge $-Q_{tip}<0$. Because of the magnetoelectric interaction, strong enough electric field makes possible creation of magnetic vortices. In this section we work in the limit of low temperatures, where there is no thermally activated vortex-antivortex (VA) pairs (for $T\ll\pi\rho_{s}=q_{m}^{2}$ their number is exponentially small).

Electric potential produced by the tip of the cantilever with charge $-Q_{tip}$ is

[TABLE]

where $H$ is the distance from the tip to the sample, $r$ is the polar radius of the point of observation. We assume that the sample is thin: its thickness $h\ll H$, so we neglect the variation of the electric field in the transverse to the film direction inside the film.

III.2 Critical tip charge for the creation of the first VA pair

In this subsection we find critical value of electric field that creates the first vortex-antivortex pair at $T=0$ in a sample of big enough size to accommodate the pair as a whole (”big sample”).

Let the electric field produce one VA-pair with vortex situated at the center of the disk and antivortex placed at the distance $r$ from its center (in the considered case of big samples $r\ll R$). Energy of the pair as a function of $r$ consists of the magnetic energy (19) of VA-interaction and magnetoelectric part (17):

[TABLE]

Instead of $Q_{tip}$, $r$, and $H$ it is convenient to introduce dimensionless variables $\kappa$, $x$, and $h$:

[TABLE]

The total dimensionless energy expressed with the parameters $\kappa$, $x$, and $h$ is

[TABLE]

Then, minimizing the energy with respect to the distance between the vortex and the antivortex $r$ by the condition $dW/dx=0$ we obtain the following cubic equation

[TABLE]

This equation always has a negative root (Fig.2). When parameter $\kappa$ increases from [math] to some critical value $\kappa_{0}$ the equation also acquires a positive multiple root. This happens when the discriminant of the cubic equation is equal to [math], from which we find $\kappa_{0}^{2}=3\sqrt{3}/2$ and $x_{min}=2$, which means that for this critical value of charge the most favourable distance for VA-pair is $r_{min}=H\sqrt{2}$. For $\kappa>\kappa_{0}$ the potential energy always possesses a local minimum. Whether it will or will not be an absolute one depends on the dimensionless parameter $h=H/a\gg 1$.

The vortex-free state has zero energy. Then, the local minimum state with one VA pair becomes the global one when its energy becomes negative. This condition $W<0$ is fulfilled when:

[TABLE]

where $x_{0}(\kappa)$ is the greatest root of (25) that corresponds to the local minimum of energy; $x_{0}(\kappa)\geq 2$. Since the left hand side of (26) is positive, the above inequality can be fulfilled only if the right hand side is also positive, therefore, at least $\kappa>\ln h\gg 1$. For such big $\kappa$ equation (25) always has 2 positive roots, with the greater one $x_{0}\approx\kappa^{2}$ (and the lesser one $x\approx 1/\kappa$ corresponding to the local maximum of energy). Returning back to the original variables we find the optimal distance between the vortex and the antivortex that minimizes the energy $W$:

[TABLE]

Substituting $x_{0}\approx\kappa^{2}$ in (26) with the accuracy $O(\kappa^{-1})$ we get

[TABLE]

An approximate solution of the equation $\kappa-\ln h-\ln\kappa=1$ is $\kappa\approx\ln h+1+\ln(\ln h+1)\approx\ln h$. From this we deduce the critical charge $Q_{tip}^{crit}$, above which the first VA-pair forms

[TABLE]

For $Q_{tip}\gtrsim Q_{tip}^{crit}$ vortex-antivortex solution corresponds not only to a local but also to the global minimum of energy.

We also define the critical voltage of the tip as voltage produced at the nearest to the tip point of the sample:

[TABLE]

For the same typical values of parameters as we used in Sec.II, we get $2q_{m}^{2}/q_{e}\simeq 10..100\mbox{V}$, so $\varphi_{0}^{crit}\simeq(10..100\mbox{V})\ln(H/a)$. Higher values of tip voltages can lead to systems with several vortex-antivortex pairs.

III.3 Vortex creation in small samples

Above we considered the case of a big sample, when VA-pair lies far from its boundary. The opposite case of small samples is also possible and some estimates were made in Ref. [Pyatakov:2010, ]. If $r_{0}$ becomes greater than $R$ (or $R<Q_{tip}q_{e}/2q_{m}^{2}$), the antivortex leaves the sample, and because the positively electrically charged vortex remains in its center, the boundary of the sample acquires a negative charge (see Appendix A for details). The energy of the system is

[TABLE]

where $\Delta\varphi=Q_{tip}(1/H-1/\sqrt{R^{2}+H^{2}})$ is the electric potential difference between the boundary and the center of the sample. In this case the single vortex is created when energy $W$ becomes negative, i.e. when voltage

[TABLE]

is applied.

IV Vortex distribution in continuous approximation

In this section we construct an analytical theory in the continuous approximation for vortex distribution. It proves that even such a crude approximation captures many effects and allows us to find the vortex distribution profile, the number of vortices for given $Q_{tip}$, and the ground state energy of the system. These results are numerically verified in Sec.VI.

IV.1 Self-consistent calculation of the vortex density $n(r)$ and the number of vortices $N$

The number of vortices in the system is not fixed (as if it were an external parameter), but rather is determined by the condition of the energy minimum. In this section we calculate self-consistently the number of vortices $N$ created by a charge $-Q_{tip}$ at $T=0$.

Consider the continuous limit (which means that at least $N\gg 1$, so $Q_{tip}\gg Q_{tip}^{crit}$, see (29)). Let $n_{v}({\mathbf{r}})$ and $n_{a}({\mathbf{r}})$ be concentrations of vortices and antivortices respectively. Let $n({\mathbf{r}})=n_{v}({\mathbf{r}})-n_{a}({\mathbf{r}})$. At $T=0$, regions of non-vanishing $n_{v}$ and $n_{a}$ does not intersect, since if they did, the vortices in the intersection area would annihilate with antivortices in order to lower the energy of the system.

Consider a ring between radii $r$ and $r+dr$. In the ground state there is a balance between the electric and magnetic interactions. Magnetic interaction pushes vortices outside while external electric force pulls them inside the ring (and vice versa for antivortices). From (21) external electric force acting upon a given area element with electric charge $q_{e}dN=q_{e}ndS$ is:

[TABLE]

From (19) magnetic force is $2q_{m}^{2}/r$; we can now calculate magnetic force acting at the same area element with magnetic charge $dq_{m}=q_{m}ndS$, using the Gauss’s theorem from electrostatics: an axially-symmetric ring acts on the outer charges with such a force as if all charges of the ring were concentrated in its center and does not act on the inner charges, therefore the magnetic force is

[TABLE]

where $q_{m}^{inside}=q_{m}\int n({\mathbf{r}})d^{2}r=2\pi q_{m}\int_{0}^{r}n(r)rdr$ is a net magnetic charge inside the considered ring.

Equating (31) and (32) we get:

[TABLE]

Multiplying by $r$ and differentiating this with respect to $r$ we obtain the vortex distribution density profile:

[TABLE]

The $r$-dependence of $rn(r)$ that is proportional to the number of vortices at distance $r$ from the disk center is shown in Fig.3. Vortices are concentrated inside a small circle of radius $r_{0}=H\sqrt{2}$, antivortices are smeared over the sample outside $r_{0}$ with their density decaying as $rn(r)\sim-1/r^{2}$.

From (33) we get the total topological charge of vortices and antivortices inside a circle of radius $r$:

[TABLE]

Now we can calculate the total number of vortices (which for big samples is the same as the number of antivortices). The border between vortices and antivortices lies at $r_{0}=H\sqrt{2}$. Using this we get

[TABLE]

For a small system with radius $R$ the number of antivortices can be less than the number of vortices, and there will be a linear charge density at the edge of the sample. This edge charge is

[TABLE]

IV.2 Total ground state energy of the system

Here we calculate the magnetic and electric energy of system of vortices and antivortices in the ground state.

First, we calculate the magnetic energy. We use the continuous version of (19), which can be obtained by the substitution $k_{i}\rightarrow n({\mathbf{r}})dS\equiv n(r)dS$

[TABLE]

Integral over $d^{2}r^{\prime}$ can be calculated using Gauss’s theorem for the potential energy: $\int_{r^{\prime}<r}d^{2}r^{\prime}n(r^{\prime})\ln|{\mathbf{r}}-{\mathbf{r}}^{\prime}|/a=\ln(r/a)\int d^{2}r^{\prime}n(r^{\prime})\equiv\ln(r/a)N(r)$, where $N(r)$ is given by (35), so we obtain

[TABLE]

Integrating this we get:

[TABLE]

Second, we find the energy of the magnetic vortices in external electric field (21):

[TABLE]

Summing (37) and (38) we get the total ground state energy of the system in the continuous limit:

[TABLE]

Since it is negative, and the energy of uniform vortex-free state is taken to be [math], then in the continuous limit formation of vortex structures is always energetically favourable. This estimate, as well as the continuous limit itself, is valid only for $Q_{tip}\gg Q_{tip}^{crit}$.

Concluding this section, we note that despite continuous limit being a crude approximation, it gives us the distribution density of induced vortices and their number, which are numerically verified in Sec.VI. The next section addresses some effects that can not be captured in the continuous approximation.

V Calculation of polarizability: discreteness effects

In this section we find the polarizability of big samples with zero total magnetic vorticity, for which even the most distant from the tip antivortices do not feel the boundaries of the sample. As explained below, continuum approach cannot be applied to polarizability calculation and we have to deal explicitly with the effects of discreteness of vortex distribution over the sample.

We study two limiting cases: (A) the case of very low temperatures $T\ll q_{m}$ when only tip-induced vortices are present; (B) the case of higher temperatures $T\lesssim q^{2}_{m}/2\simeq T_{BKT}$, when close to the Berezinskii-Kosterlitz-Thouless temperature there also exist some amount of thermodynamically activated pairs.

V.1 $T\ll q_{m}^{2}$

Consider a system of $N$ vortices and $N$ antivortices. As a first approximation we can neglect the effect of the finite size of the positive vortex “nuclei” and consider it as a point-like one, neglecting its contribution to the polarizability. Then, remaining $2D$ system of $N$ antivortices has $2N$ degrees of freedom, therefore, $2N$ normal modes.

Consider the system near its local energy minimum. Polarizability at $T=0$ is largely determined by the most soft mode of the system. If the system possesses a zero mode that couples with the electric field, then an arbitrary small perturbation leads to huge distortions, so polarizability of such state is infinite. However, temperature fluctuations can effectively smear the low lying modes, which makes the polarizability finite and strongly temperature dependent.

The problem of finding polarizability can be exactly analytically solved for small fixed values of $N$, but the exact solution seems to be inaccessible for arbitrary large $N$. Here we employ the strategy similar to “energy landscape approach” used for overcooled liquids, glasses and spin glasses, proteins folding, and melting of small clusters Wales:book .

Let $r_{i}=r_{i}^{(0)}+\Delta r_{i}$ ($i=1$,…,$2N$), where $\Delta r_{i}$ are deviations of cartesian coordinates of $N$ antivortices from the local minimum, with the first $N$ components giving us $x$-coordinates of antivortices, and the second $N$ components giving their $y$-coordinates: $r_{1}=x_{1}$, …, $r_{N}=x_{N}$; $r_{N+1}=y_{1}$, …, $r_{2N}=y_{N}$. The total energy $W=W^{(0)}+W^{(1)}$, where $W^{(0)}=\sum_{i,j}V(x_{i},y_{i};x_{j},y_{j})+\sum_{i}V_{tip}(x_{i},y_{i})$ is the sum of magnetic energy of antivortices and electrostatic energy interaction of antivortices with the tip. Small perturbation $W^{(1)}=q_{e}E\sum_{i=1}^{N}\Delta r_{i}$ is the energy of interaction with the infinitesimal in-plane electric field $E$ in the $x$-direction.

We employ the Hessian matrix:

[TABLE]

Near the local minima, energy $W^{(0)}$ can be written as a positive definite quadratic form (vector $\Delta{\mathbf{r}}=(\Delta r_{1},...,\Delta r_{2N})$):

[TABLE]

Symmetric matrix $K_{ij}$ can be diagonalized by some orthogonal transformation $O$: $\tilde{K}=OKO^{-1}$, with coordinates $\Delta r_{i}$ being accordingly transformed to normal coordinates $R_{i}$: $\Delta{\mathbf{r}}=O{\mathbf{R}}$. Then $W^{(0)}$ transforms as

[TABLE]

where $\tilde{K}=diag\{K_{1},...,K_{2N}\}$, $K_{i}$ are eigenvalues of the Hessian matrix, which play an important role as we will see below.

Since $\Delta r_{i}=(O{\mathbf{R}})_{i}=\sum_{j=1}^{2N}O_{ij}R_{j}$, $W^{(1)}$ transforms as

[TABLE]

where

[TABLE]

The partition function can be written in terms of normal coordinates $R_{i}$ as

[TABLE]

Here we introduced a cutoff $L$, which takes into account the fact that for big enough $R_{i}$ potential $W^{(0)}$ is no longer a parabolic one. $L$ corresponds to such displacement of antivortices that shifts the point in the configuration space from the local minimum to the nearest saddle point; we expect $L$ to be of order of distance between nearest antivortices.

Polarization of antivortices with negative charges $-q_{e}$ is

[TABLE]

Polarizability can be found as

[TABLE]

Substituting here the partition function (V.1) we get

[TABLE]

As $T\rightarrow 0$ ($\beta\rightarrow\infty$) Gaussian integrals converge very fast and cut-off $L$ does not play any role: we can extend the limits of integration up to infinity. Then from (48) we find

[TABLE]

From this we see that in the limit $T\rightarrow 0$ polarizability $\alpha$ is determined by the most “soft” modes $K_{i}$.

Now consider the case of non-zero $T$. Performing integrations in (48) we see that cut-off $L$ now enters the expressions for polarizability explicitly:

[TABLE]

where we denoted

[TABLE]

playing the role of smooth cut-off function for the terms in the sum (50) with the lowest eigenvalues $K_{i}$. Function $f(x)$ is sketched in Fig.4.

In the limit of big number of vortices and antivortices $N\gg 1$ the sum (50) can be approximated by an integral:

[TABLE]

where $D(K)$ is the density of $K$-modes.

In order to make some estimates of $\alpha(T)$ and deduce its leading temperature dependence we use the following three approximations.

(i) We estimate $D(K)$ using random matrix theory for an Euclidean ensemble Parisi:1999 . Amir et.al. Amir:2010 analytically calculated the density of eigenvalues for certain 1D dilute systems with short range potential and obtained $D(K)\sim 1/K$; for higher dimensions they predicted that there should be some correction

[TABLE]

where $\eta>0$ depends on the dimensionality of the system and on the form of the particle interaction potential. In Sec.VI.3 we check this assumption by exact numerical diagonalization of the Hessian matrix and find the exponent $\eta$.

(ii) We replace the smooth cut-off function with a sharp Heaviside theta function $f(x)\rightarrow\theta(1-x)$, so the lower limit of integration becomes $T/2L^{2}$. This can be done since $f(x)\sim\frac{2}{3}x^{2}$ for $x\rightarrow 0$ and integral (52) converges at the lower limit as $\int K^{1-\eta}dK$.

(iii) Since coupling to the electric filed is determined only by $x\rightarrow-x$ symmetry of the mode, we assume that all modes couple to the electric field with approximately the same strength, apart from fluctuations on small-$K$ scale:

[TABLE]

For example, if the stationary configuration $r^{(0)}$ is symmetric under reflection $x\rightarrow-x$, then all modes are two-fold degenerate and either symmetric (does not couple to the electric field) or antisymmetric (couple to the electric field with the same weight).

Finally, using the above three approximations, we obtain the following temperature dependence of polarizability:

[TABLE]

Since the integral converges we extended the upper limit of integration up to $\infty$. This estimate holds as long as $T\gg K_{min}L^{2}$. For very low temperatures $T\ll K_{min}L^{2}$, $\alpha\sim 1/K_{min}^{1-\eta}=const(T)$.

The exponent $\eta$ is in a sense a non-universal one and it should depend not only on the dimensionality of the system, but also on the form of the particle interaction and external potentials. We numerically test the validity of the above approximations and the $\alpha\sim 1/T^{1-\eta}$ dependence for our particular case in Sec. VI.3.

V.2 Some estimates for $T\lesssim q_{m}^{2}/2$

Now consider the case $T\lesssim q_{m}^{2}/2$. Then there also exist VA pairs created by thermal fluctuations. The profile of the total vorticity $n({\mathbf{r}})=n_{v}({\mathbf{r}})-n_{a}({\mathbf{r}})$ remains the same as it was found in IV.1 from the force balance considerations (formula (34) and Fig.3). The difference is that now the regions of non-vanishing $n_{v}(\mathbf{r})$ and $n_{a}(\mathbf{r})$ intersect with each other. This describes the process of penetration of antivortices in the vortex region and vice versa because of VA pairs formation.

Now there will be two contributions to polarizability:

(i) contribution of tip-induced vortices, but now their interaction will be renormalized by the Kosterlitz-Thouless dielectric function: $W(r)=\pm 2q_{m}^{2}\ln(r/a)/\epsilon(r)$;

(ii) pure contribution of thermally activated dipole pairs away from the tip.

Below we consider only the renormalized contribution (i) (pure contribution of thermally activated vortices (ii) is considered in Appendix B).

Renormalized interaction energy is $W(r)=\pm 2q_{m}^{2}\ln(r/a)/\epsilon_{KT}(r)$. As long as $\chi_{KT}=(\epsilon_{KT}-1)/2\pi\lesssim 1$, which holds in almost the whole range $T<T_{BKT}\approx q_{m}^{2}/2$ except in a small temperature interval in the vicinity of $T_{BKT}$, there will be no qualitative influence of the thermal VA-pairs onto the field induced vortices.

We can estimate the melting temperature for the induced “magnetic atom” patterns. The pattern starts to melt at a temperature $T_{c}$ when the position fluctuations become comparable to intervortex distance $\Delta r\sim 1/\sqrt{n_{v}}$ ($n_{v}$ is a local vortex concentration). Make an estimate for the vortex situated in the center of the disk. Here only electric potential influences, which is almost parabolic:

[TABLE]

Using (36) and $n_{v}\sim\pi r_{0}^{2}/N=2\pi H^{2}/N$ we get:

[TABLE]

Hence, the magnetic pattern is stable up to $T\lesssim T_{BKT}$. For $T>T_{BKT}$ thermodynamically induced vortex pairs dissociate and break the local magnetic order.

VI Numerical simulation

In this section we describe the procedure and the results of numerical simulation.

VI.1 Numerical procedure

We modeled a system of point-like vortices and antivortices, which interact between themselves (with the magnetic energy (19)) and with the external electric field created by cantilever’s tip (with energy (17), electrostatic potential is given by (21)). The system shows “glassy” behavior, i.e. there are many local minima that are close to the global one. In order to find a state corresponding to such a local energy minimum we performed simulated annealing using Monte Carlo (MC) method Kirkpatrick:1983 . Three types of MC-steps were used: movement of a vortex or an antivortex, creation of VA pair, and destruction of VA-pair. We cooled the system from high temperatures $T_{start}\simeq q_{m}^{2}>T_{BKT}\simeq q_{m}^{2}/2$ (which allowed the system to find the optimal number of vortices) down to as low as $T_{stop}=10^{-6}q_{m}^{2}$, depending on the numerical experiment. Geometrical cooling schedule was used Kirkpatrick:1991 :

[TABLE]

with $\gamma=0.99$ and typically $N_{it}=100000$ number of Monte Carlo iterations at a given temperature.

The simulation was performed for the following ratios between the lattice constant $a$, distance $H$ from the tip to the sample, and radius $R$ of the sample: $a=0.01H$, $R=100H$.

VI.2 Number of vortices and the critical field

As we described in detail in Sec.III and IV, charged tip produces magnetic patterns (“magnetic atoms”) in the system under consideration. Fig.5 shows typical low temperature ($T=0.001q_{m}^{2}$) distribution of magnetization in the system for moderate values of tip voltages ($\varphi=Q_{tip}/H=23q_{m}^{2}/q_{e}$), for which $N=3$ vortices and an $3$ antivortices are created.

Fig.6 shows typical low-temperature configuration of the system for higher values of tip’s voltage ($Q_{tip}q_{e}/H=240q_{m}^{2}$). For the depicted configuration the observed number of vortices is $N^{num}=43$, whereas analytical formula (36) gives $N^{analyt}=240/3^{3/2}\approx 46.2$, this mismatch is explained below. Vortices and antivortices are arranged in a highly regular structure that is reminiscent of the structure of an atom: negatively charged antivortex outer concentric circular “shells” surround positive vortex “nucleus”. Fig.6b depicts only this vortex “nucleus” of the “atom”; arrangement of the vortices, interacting via 2D Coulomb-like $\log r$ magnetic potential, is similar to one of the trapped single species Coulomb plasma, where charges are confined to a plane and interact via 3D Coulomb $1/r$ potential Bedanov:1994 ; Radzvilavicius:2011 . Similar patterns of of magnetic skyrmion systems were observed recently Zhao:2016 .

Fig.7 shows $N(Q_{tip})$ dependence for the ground states of the system with varying $Q_{tip}$. We see that analytical result (36) is in a good agreement with the numerical modeling. Analytical theory, which was derived in a continuous approximation, gives the correct slope, but slightly overestimates the number of vortices.

We also see that the analytical result for the critical charge that creates the first VA-pair (29) gives for the considered parameters ($a=0.01H$) $Q_{tip}^{crit}=2\ln 100q_{m}^{2}H/q_{e}\approx 9.2q_{m}^{2}H/q_{e}$, whereas the numerical simulation gives $Q_{tip}^{crit}\approx 9q_{m}^{2}H/q_{e}$ (inset of Fig.7). This mismatch $\Delta Q=Q_{tip}^{crit}$ between the continuous theory and simulation for $Q_{tip}$ at a given $N$ leads to mismatch in $N$ for given $Q_{tip}$, and the average mismatch should be of order of $\Delta N\approx\Delta Qq_{e}/3^{3/2}q_{m}^{2}H\approx 1.7$, moreover real step-like dependence gives an additional $\pm 1$ deviation from the continuous formula (36), which explains the above mismatch $\Delta N\approx 3.2$ between the analytical formula and the numerical simulation for $Q_{tip}q_{e}/H=240q_{m}^{2}$.

VI.3 Polarizability of big samples

In this section we numerically find the temperature dependence of polarizability of big samples. In order to find polarizability we performed cooling at a fixed value of electric field $E$, then independently repeating this for different values of $E$. From the theoretical analysis of Sec.V we expect the polarizability to be strongly temperature-dependent.

Fig.8 shows the variation of the $x$-th and $y$-th components of the dipole moment with the applied electric field for fixed temperature $T=0.001q_{m}^{2}$. Linear fit for $p_{x}(E)$ dependence gives polarizability $\alpha_{num}=1388.9H^{2}q_{e}^{2}/q_{m}^{2}$.

Fig.9 shows the temperature dependence of polarizability in log-log coordinates, which in the limit of small $\beta$ is fitted by $\ln\alpha\approx 0.74\ln(\beta q_{m}^{2})+\operatorname{const}$. This dependence is consistent with the estimate (55), when we choose $\eta\approx 0.26$, which gives then $\alpha\sim 1/T^{0.74}$.

Using the exact numerical diagonalization we have found an eigenspectrum (Fig.10), from which indeed exponent $\eta\approx 0.25$.

We propose that such non-trivial temperature dependence of polarizability of finite clusters is a generic one. For glassy systems potential energy landscape of the system possesses many local minima that are close to the global one; adding one VA pair to the system completely changes this landscape and its local minima, but preserves its Hessian eigenvalue distribution. Then we argue that temperature fluctuations smear the fine details of this energy landscape, effectively cutting off the low-$K$ (low frequency) normal modes $K_{i}\ll T/L^{2}$. This makes polarizability of the system strongly temperature dependent with a non-trivial dependence $\alpha\sim 1/T^{1-\eta}$. The exponent $\eta$ is non-universal one and depends on the density of Hessian eigenvalues. For example, for single-species 2D plasma confined by quadratic external potential, with particles interacting via 3D Coulomb $1/r$ potential $D(K)\sim 1/K^{0.33}$ [Shweigert:1995, ], so from (53) in this case $\eta\approx 0.67$.

VII Possible experimental observation

In this section we consider the possibilities of experimental observation of the described effects and make estimates for the required electric fields. The most plausible candidates possessing magnetoelectric effect seem to be multiferroic materials. However the most of the known multiferroics are antiferromagnets with the staggered magnetization playing the role of the order parameter. We believe that theoretical results presented in the previous sections are applicable to them at least qualitatively if not quantitatively. The important remark should be now kept in mind: by vortex/antivortex we understand now $\pm 2\pi$ rotation of the staggered magnetization vector. Also in order to ease the comparison with the experimental data we use SI units throughout this section; for transferring from Gaussian to SI units substitutions $\chi_{e}\rightarrow\chi_{e}\epsilon_{0}$, $Q_{tip}\rightarrow Q_{tip}/4\pi\epsilon_{0}$ should be made.

VII.1 Theoretical estimates for real materials

First, we consider the potential materials where the predicted effects could be found. In order to observe creation of “magnetic atoms” by electric field in a thin film we need a fairly strong coupling between the electric and magnetic subsystems is needed. Below we make estimates for type-I and type-II multiferroics.

Type-I multiferroics. We consider prototypical type-I multiferroic – $\mathrm{BiFeO_{3}}$. Its single crystals are rhombohedral at room temperature with the space group $R3c$ and pseudocubic unit cell with $a\approx$4\text{,}\mathrm{\SIUnitSymbolAngstrom} and G-type antiferromagnetic structure [[Lebeugle:2008PRL, ](#bib.bib57)], its spontaneous polarization $P\simeq$1\text{\,}\mathrm{C}\mathrm{/}\mathrm{m}^{2} [Lebeugle:2007APL, ,Neaton:2005, ], which induces cycloidal twisting of the staggered magnetisation with the period $\lambda\simeq 620\div$640\text{,}\mathrm{\SIUnitSymbolAngstrom}$$ [Sosnowska:1982, ,Lebeugle:2008PRL, ], dielectric constant $\epsilon=1+\chi_{e}\simeq 50$ [Lu:2010, ],

In order to estimate constant $\alpha$ characterizing the strength of the exchange interaction (see eq. (3) above), we use the exchange constant for nearest neighbor Fe-Fe interaction $J=$10\text{,}\mathrm{m}\mathrm{e}\mathrm{V} [[Delaire:2012, ](#bib.bib63)]. Hence for the BCC pseudocubic setting we obtain exchange stiffness $A\equiv\alpha M_{0}^{2}/2=Ja^{2}/V=$4\text{\,}\mathrm{p}\mathrm{J}\mathrm{/}\mathrm{m} which is consistent with the estimate deduced from the Neel temperature Sando:2013 ($V\approx a^{3}$ is the volume of the pseudocubic unit cell).

Now using (3), (4), (7) and taking $\phi=kx$, where $k=2\pi/\lambda$, we get

[TABLE]

Minimizing this energy density with respect to $k$ we obtain $k=\gamma P_{y}/\alpha=2\gamma P_{y}M_{0}^{2}/A$ or

[TABLE]

Using estimate (57) we find that $\chi_{e}\epsilon_{0}\gamma^{2}M_{0}^{2}\simeq 3\times 10^{-16}\mbox{\, J/m}$, which is much less than $\alpha M_{0}^{2}\simeq 8\times 10^{-12}\mbox{\, J/m}$ whereas is much less. This indeed justifies our simplification of the formula (10) in Sec.II. Now from (12) using $h=a\approx$4\text{,}\mathrm{\SIUnitSymbolAngstrom}$$we find the vortex electric charge per one layer $q_{e}\approx 9\times 10^{-22}\mbox{\,C}\approx 0.006e$ ($e$ is the charge of electron). From (18) we find the value of “magnetic charge” and using (30) we find the critical tip voltage for creation of first VA pair $\varphi_{0}^{crit}\simeq 20\mbox{\,V}\times\ln(H/a)$. Using (29) with $\ln(H/a)\sim 1$ we also estimate the charge of impurity necessary to create one VA pair in a monoatomic film layer as $Q^{crit}_{impurity}\simeq 10^{-18}\mbox{\,C}\approx 6e$.

Type-II multiferroics. Namely, type-II multiferroics with a noncollinear magnetic order seem to be very natural candidates for observation of “magnetic atoms”. Among them are cycloidal helimagnets such as $\mathrm{TbMnO_{3}}$, $\mathrm{DyMnO_{3}}$, $\mathrm{Eu_{0.75}Y_{0.25}MnO_{3}}$, $\mathrm{Ni_{3}V_{2}O_{8}}$, $\mathrm{MnWO_{4}}$ etc., see Kimura:2007 ; Tokura:2014 for reviews. Although the ground state in cycloidal helimagnets is different from the collinear one, vortex excitations of this ground state can be created anyway.

Here we make estimates for the prototypical type-II multiferroic $\mathrm{TbMnO_{3}}$. At room temperature it has orthorhombically distorted perovskite structure with space group $Pbnm$ and the lattice parameters $a=$5.3\text{,}\mathrm{\SIUnitSymbolAngstrom}, $b=$5.86\text{\,}\mathrm{\SIUnitSymbolAngstrom}, $c=$7.49\text{,}\mathrm{\SIUnitSymbolAngstrom} [[Kenzelmann:2005, ](#bib.bib66)]. Below $T_{N}=41\mbox{\,K}$ magnetic order appears – sinusoidal spin density wave caused by frustrations is formed, with no electric polarization. Below $T=28\mbox{\,K}$ magnetic cycloidal state in $bc$ easy plane with the wavevector $k=\pi(0,0.28,0)$ (and $\lambda=2\pi/k\approx$42\text{\,}\mathrm{\SIUnitSymbolAngstrom}$\approx 7b$) becomes energetically favorable and it induces the electric polarization which grows with lowering the temperature and reaches the value $P=8\times 10^{-4}\,$\text{\,}\mathrm{C}\mathrm{/}\mathrm{m}^{2} at $T=10\mbox{\,K}$; dielectric constants are $\epsilon_{a}\simeq 24$, $\epsilon_{b}\simeq 23$, $\epsilon_{c}\simeq 29$ [Kimura:2003, ], for estimates we get the mean value $\epsilon_{bc}=(\epsilon_{b}+\epsilon_{c})/2=26=1+4\pi\chi_{e}$. Using (6) we calculate the magnetoelectric constant

[TABLE]

As expected its value is greater than for type-I multiferroics (57), though the difference is not that drastic. From (12) we also find the vortex electric charge per one layer $q_{e}\approx 1.8\times 10^{-21}\mbox{\,C}\approx 0.01e$.

In order to estimate the exchange stiffness and “magnetic charge” we use values of microscopic exchange constants. Literature data on these is somewhat controversial, with even different relative strengthes of exchange constants along different unit-cell directions Xiang:2008 ; Nagaosa:2010 ; Milstein:2015 ; Fedorova:2015 . For the estimate we use the exchange in the $c$-direction, which is free of magnetic frustrations, as $J_{c}=$1.5\text{,}\mathrm{m}\mathrm{e}\mathrm{V}, consistent with estimates from the Neel temperature [Heffner:2001 ](#bib.bib72) $k_{B}T_{N}\simeq Jz/3$ ($z$ is the number of nearest neighbors for the magnetic atom); therefore the exchange stiffness $A=J_{c}c^{2}/2V\simeq$0.3\text{\,}\mathrm{p}\mathrm{J}\mathrm{/}\mathrm{m}. Then from (30) we find the critical tip voltage for creation of first VA pair $\varphi_{0}^{crit}\simeq 0.3\mbox{\,V}\times\ln(H/a)$. Using (29) with $\ln(H/a)\sim 1$ we estimate the charge of impurity necessary to create one VA pair in a single layer as $Q^{crit}_{impurity}\simeq 1.6\times 10^{-20}\mbox{\,C}\approx 0.1e$.

Other materials. In principle there might be a wide class of materials where magnetic vortices can be created by electric field, with multiferroics substituting only a narrow subclass of them: potentially, any magnetic insulator can exhibit a magnetoelectric phenomena, when symmetry arguments allow this Khomskii:2016 . One microscopic mechanism for this is the Dzyaloshinskii-Moriya interaction, which is present in many materials, such as weak ferromagnets, helical magnets, etc.

VII.2 Magnetic anisotropy

Here we discuss the limitations of our theoretical picture due to possibility of out-of-plane polarization induced by the tip of the cantilever. The charged tip induces an electric field perpendicular to the plane and consequently a transverse electric polarization, the latter can in turn induce a magnetic spiral with out-of-plane magnetisation component and in-plane wavevector.

In order this effect to be minor there should be some uniform easy-plane anisotropy that also holds magnetization in the plane in the absence of the tip. For example, single-ion anisotropy or one connected with the sample geometry (for ferromagnetic systems it is just the demagnetization energy; but for antiferromagnetic systems some weak surface anisotropy effects are also expected Caminale:2014 ).

First, consider a system with tetragonal or orthorhombic symmetry, where an easy plane exists, and magnetic anisotropy energy is relatively big. Estimate the necessary transverse electric field that can induce a spiral in a “hard” plane (plane containing the hard direction). Here we chose as an example $\mathrm{TbMnO_{3}}$, which has a ground state with a cycloidal spiral in $bc$ plane (thus hard axis is $a$); for estimate we take anisotropy parameter $K_{an}\lesssim 0.1J\simeq 0.1\mbox{\, meV}$ [Xiang:2008, ] per one magnetic Mn atom.

Let the electric field be strong enough to induce a cycloidal magnetic spiral with an out-of-plane magnetization component. Using (3) and (4) we get the following estimate for the energy gained by the spiral state with induced out-of-plane polarization $P_{z}=\epsilon_{0}\chi_{e}E_{z}=\epsilon_{0}\chi_{e}E_{0z}/\epsilon\approx\epsilon_{0}E_{0z}$ ($z$-axis is perpendicular to the plane, $E_{z}$ is the field inside film, $E_{0z}$ is the field in vacuum, $\epsilon\approx 1+\chi_{e}\approx\chi_{e}$):

[TABLE]

Minimizing this with respect to wavevector $k$ we obtain $k\simeq\gamma\epsilon_{0}E_{0z}/\alpha$, therefore (59) transforms to

[TABLE]

which is the energy density that system gains when the spiral state is created. However the anisotropy energy density of the system increases by an amount:

[TABLE]

where factor of $4$ comes because of $4$ magnetic Mn atoms per $\mathrm{TbMnO_{3}}$ unit cell of volume $V=abc$ and factor of $2$ in the denominator because of averaging local anisotropy energies $\sim\cos^{2}(kx)$ over spiral propagation direction $x$. Comparing the gained by the spiral state energy (60) and the increase of the anisotropy energy (61), and taking typical values of parameters for $\mathrm{TbMnO_{3}}$ considered above, we obtain the following estimate for electric field, needed to produce a magnetic cycloid rotating in a “hard” plane

[TABLE]

For comparison, if the tip is placed at $H=$10\text{,}\mathrm{n}\mathrm{m} above the film, in order to create such strong electric field that will influence $M_{z}$ component, the tip’s voltage has to be as large as $\varphi=E_{0z}H\approx$200\text{\,}\mathrm{V} (and even larger for greater $H$). This value is well above the critical voltage $0.3\text{\,}\mathrm{V}$ for VA pair formation estimated for $\mathrm{TbMnO_{3}}$ earlier. Thus, we conclude that almost any reasonable nonzero magnetic easy-plane anisotropy protects the in-plane spin arrangement, in spite of the possibly strong perpendicular to the plane electric field component induced by the tip.

For systems with very weak or no single-ion anisotropy, such as $\mathrm{BiFeO_{3}}$, where only surface effects protect the in-plane magnetization, perpendicular component of the electric field will induce the out-of-plane component of the magnetization.

One way to avoid an out-of-plane magnetisation for such systems is to use two symmetrical charged tips from both sides of the film (Fig.11; still $H\ll h$).

However we claim elsewhere that even for single tip geometry (Fig.1), an out-of-plane component of magnetization near the tip would be localized inside the area of the radius $\sim H$ (distance from the tip to film surface). This will not affect parts of the system outside the latter area, away from the tip. Therefore, the “magnetic atom” nucleus consisting of vortices in $r\sim H$ area may drastically change, but antivortex shells will preserve their qualitative features. This will happen because of the topological nature of vortex excitations. For example, if anisotropy is not strong enough, vortices may have core with an out-of-plane component of ${\mathbf{M}}$ – this will make the core electric charge more smeared, but it will not influence its amount. For a single vortex, the boundaries of the sample are negatively charged with the same charge as if the vortex core was “flat”. The boundaries are not influenced by an out-of-plane component of the electric field of the tip, since it decreases as $\sim 1/r^{3}$ away from the center, while in-plane field decreases only as $\sim 1/r^{2}$. Then by the charge conservation, the positive charge at the vortex core has to remain the same.

VIII Conclusions

In magnetoelectric materials interaction of magnetic and electric subsystems makes possible the creation of magnetic patterns by an applied electric field. In this paper we have presented analytical and numerical analysis of such quasi-two-dimensional type-II multiferroic-like materials, where “atom”-like patterns of magnetic vortices and antivortices can be created by a strong electric field, and studied structure and electric polarizability of such magnetic “atoms”.

When the spot-like electric field, produced by a tip of a cantilever, exceeds some critical value, the first vortex-antivortex pair forms; with increasing electric field more vortices and antivortices are created from the “magnetic vacuum” of the magnetoelectric film. We found analytically and confirmed numerically the critical voltage that is required to create the first vortex-antivortex pair, and found the dependence of the number of vortex pairs on the applied voltage. We have shown that a “magnetic atom” forms: a “nucleus” consists of vortices arranged in a lattice-like structure with a short range order; it is surrounded by antivortices forming outer concentric circular “shells”. We obtained the vortex density profile in the continuous approximation.

In the studied here type-II multiferroic-like magnetoelectric materials vortices carry a positive and antivortices – a negative electric charge, therefore additional weak in-plane electric field pulls them in the opposite directions, thus creating a polarization. In this work we have studied the temperature dependence of electric polarizability of such “magnetic atoms”. Polarizability of finite vortex-antivortex cluster is determined by the eigenvalues of its Hessian matrix. We analyzed properties of these normal modes using results from Euclidean random matrix theory. We speculate that the behavior of polarizability $\alpha\sim 1/T^{1-\eta}$ is generic for disordered (“glassy”) finite classical systems with many local minima near the ground state energy; the exponent $\eta$ depends on the form of the potential and dimensionality of the system. Thermodynamical fluctuations wash out the fine details of the energy landscape, effectively cutting off the low-lying eigenmodes. By the exact numerical diagonalization of Hessian matrix we get the density of its eigenvalues, from which we deduce $\alpha\sim 1/T^{1-\eta}$ ($\eta\approx 0.25$) dependence of polarizability and confirm it directly by Monte Carlo simulation.

The presented theoretical and numerical analysis should find its experimental realization. First, such “magnetic atoms” can be created in thin films of type-II multiferroic materials with easy-plane spin arrangement, for example in cycloidal helimagnets (such as $\mathrm{TbMnO_{3}}$, $\mathrm{Ni_{3}V_{2}O_{8}}$ and othersTokura:2014 ). Although the ground state in cycloidal helimagnets is different from the collinear one, vortex excitations upon this ground state can be created anyway. Also, the materials where magnetic vortices can be created by electric field, are not restricted to only multiferroics: potentially, any magnetic insulator can exhibit a magetoelectric phenomena, when symmetry arguments allow this.

Atomic force microscopy (AFM) techniques might be used in order to create high local electric field. The direct observation of the magnetic vortex patterns formation can be done using high-resolution Lorentz TEM imaging Yu:2010 ; Ratcliff:2016 combined with AFM Erts:2001 ; alternatively combined AFM/STM setups can be used Nakayama:2012 ; Hapala:2015 . A number of indirect methods can be used in order to test the eigenmodes of the induced vortex-antivortex system. These eigenmodes can be experimentally found by measuring the electric polarization in the applied transverse electric field (as studied analytically in Sec.V and numerically in Sec.VI), but also by measuring heat capacity, or microwave absorption spectra experiments Onose:2012 . Dynamical effects associated with the vortex “atom” can be tested using the setup, when the cantilever oscillates along the film plane Yazdanian:2009 or applying a time-dependent voltage, and measuring the energy dissipation.

We finally note that the studied here electric field creation of magnetic vortices and antivortices out of magnetoelectric thin film ground state provides a condensed-matter toy-model for studying the long-standing problem of production of electron-positron pairs out of vacuum applying strong electric fields, which arises in quantum electrodynamics and atomic physics.

Acknowledgements.

The authors are grateful to S.A. Brazovskii for the idea to consider field-induced topological superlattices; to D.I. Khomskii and M.V. Mostovoy for helpful discussions. We acknowledge the financial support of the Ministry of Education and Science of the Russian Federation in the framework of Increase Competitiveness Program of NUST MISiS (No K2-2016-067).

Appendix A Energy of a vortex in the in-plane electric field

In this Appendix we calculate the energy of a charged vortex in the electric field. As it was explained in the main text, the magnetic vortex core acquires an electric charge (12) due to magnetoelectric coupling (II). Since our sample is electrically neutral, its edge acquires negative charge of the same magnitude, which should also be taken into account (see Fig.12a). This edge charge effectively shields the vortex charge, modifying its energy.

Let our sample be disk with radius $R$, which contains a single vortex with vorticity $n$, placed at the distance $X_{0}$ from the center of the disk. Chose a coordinate system as shown in Fig.12b; electric field ${\mathbf{E}}=(E_{x},E_{y})$.

From (8) polarization of single vortex configuration (11) is

[TABLE]

Since we calculate the linear response, $\chi_{e}{\mathbf{E}}$ term was omitted from (8). Then the electric energy of magnetic vortex configuration is $W_{me}=-h\int{\mathbf{P}}{\mathbf{E}}\,d^{2}r$

[TABLE]

Since $l(\theta)=\sqrt{R^{2}-X_{0}^{2}\sin^{2}\theta}-X_{0}\cos\theta\,$ is an even function of $\theta$, then the odd term $E_{y}\sin\theta$ vanishes. Thus,

[TABLE]

Using (12) we obtain $W_{me}=-\frac{1}{2}q_{e}E_{x}X_{0}$ or for general displacement of vortex core ${\mathbf{r}}=(X_{0},Y_{0})$

[TABLE]

We see that for the circular sample geometry this energy is exactly a half of the vortex core energy of charged particle in the in-plane electric field. This effective shielding comes from the boundary charge and its exact value depends on the geometry of the sample. We also note that since polarization (62) is linear in phase angle $\phi$, the same result (65) holds for vortex-antivortex configuration, with ${\mathbf{r}}$ denoting a vector that connects cores of the antivortex and the vortex.

Appendix B Thermally activated pairs

In this Appendix we find the polarizability that comes from thermally activated dipole pairs at $T<T_{BKT}$. Since the chemical potential of a pair is quite big, the concentration of pairs is low and we can employ the approximation of non-interacting pairs. Therefore we can find polarizability of one pair, then multiply it by the number of pairs and we will get polarizability of the VA pairs subsystem.

First, find the polarizability of one VA dipole pair. Energy of a dipole pair is the sum of magnetic interaction energy of the vortex and antivortex (19) and the magnetoelectric energy (65), obtained in the Appendix A.

[TABLE]

Here ${\mathbf{r}}=(x,y)$ is the vector connecting the antivortex core to the vortex core; the electric field is chosen to be parallel to the $x$-axis. The partition function of such a dipole:

[TABLE]

where $S=\pi R^{2}$ is the area of the system, which comes from integration over the dipole center of mass. For zero electric field we get

[TABLE]

Average dipole moment of a pair:

[TABLE]

Formally, this integral diverges at large $r$. However, we should keep in mind that we are calculating linear response (i.e., we take the limit $E\rightarrow 0$), so that $r$-divergence is cut off at the radius of the sample $r=R$. Therefore, if we first take derivative of $p(E)$ with respect to $E$ and then put $E=0$, then we obtain the correct result. Hence, polarizability of one pair is

[TABLE]

integrating which we obtain

[TABLE]

Second, estimate the concentration of VA pairs. Since a weak field does not affect their concentration, we make the estimation for $E=0$. Moreover, due to the fact that the chemical potential of a pair is quite big with respect to $T_{BKT}$: ($\mu\simeq-\pi q_{m}^{2}$, $T_{BKT}\simeq q_{m}^{2}/2$, see Ref.[KT1, ]), concentration of VA pairs is always small. Calculate the grand partition function of non-interacting dipole system:

[TABLE]

Then, the number of thermally activated VA pairs

[TABLE]

Finally, combining (67) and (68) we find the polarizability of the sample that comes from all dipole pairs

[TABLE]

From here we can deduce the electric susceptibility $\chi_{e}=\alpha_{thermal}/S$ and we see that it is consistent with Kosterlitz-Thouless’s susceptibility KT1 with the difference that now we have the electric charge $q_{e}/2$ instead of $q_{m}$ in the numerator. This is natural since the response to the electric field depends on the value of electric charge $q_{e}$ whereas screening comes from the interaction of ”magnetic” charges $q_{m}$.

We see that at low temperatures $T\rightarrow 0$, $\alpha_{thermal}$ exponentially vanishes to [math], whereas diverges as $\alpha_{induced}\sim 1/T^{1-\eta}$, therefore $\alpha_{induced}$ dominates. On the contrary, near $T_{BKT}$ $\alpha_{thermal}$ diverges and dominates over $\alpha_{induced}$.