Theory of surface-induced multiferroicity in magnetic materials, thin films and multilayers

TL;DR

This paper develops a theoretical framework to understand how electric polarization can emerge at surfaces of magnetic materials, thin films, and multilayers, highlighting two distinct mechanisms involving surface-induced interactions.

Contribution

It introduces a theoretical model distinguishing two pathways for surface-induced multiferroicity based on the presence of Dzyaloshinskii-Moriya interactions in magnetic materials.

Findings

Surface-induced electric polarization can occur via two different mechanisms.

The presence or absence of Dzyaloshinskii-Moriya interactions determines the multiferroic behavior.

Experimental implications of the theory are discussed.

Abstract

We present a theoretical study of the onset of electric polarization close to a surface in magnetic materials and in thin films and multilayers. We consider two different paths that lead to the onset of multiferroic behavior at the boundary in materials that are bulk collinear ferromagnets or antiferromagnets. These two paths are distinguished by the presence or absence of a surface induced Dzyaloshinskii-Moriya interaction which can be taken into account through Lifshitz invariants in the free energy of the system. Experimental consequences are discussed in the light of the developed theory.

Equations

$$\frac{a_0^2}{6}\frac{{\partial }^{2}m(z)}{\partial z^2}+(1-\tau )m(z)-\beta m^3(z)=0$$

$${\xi }^2\frac{{\partial }^{2}f(z)}{\partial z^2}+f(z)-f^3(z)=0$$

$$\delta F=F_{ME}+F_E=\mathbf{p}\cdot (\gamma \nabla ({\mathbf{m}}^2)$$

$$+\gamma \prime [\mathbf{m}(\nabla \cdot \mathbf{m})-(\mathbf{m}\cdot \nabla )\mathbf{m}]+...)+\frac{\mathbf{p}(\mathbf{r}{)}^2}{2{\chi }_E}$$

$$p(z)=-{\chi }_E(\gamma \nabla ({\mathbf{m}}^2)+\gamma \prime \left[\mathbf{m}(\nabla \cdot \mathbf{m})-(\mathbf{m}\cdot \nabla )\mathbf{m}\right]+...)$$

$$p(z)=-m_{\infty }^2\frac{\sqrt{2}}{\xi }{\chi }_E\gamma \left[\tanh (\frac{z+a_0}{\sqrt{2}\xi })[1-{\tanh }^2(\frac{z+a_0}{\sqrt{2}\xi })]\right]$$

$$\mathcal{F}=A\sum _i{\left(\frac{\partial \mathbf{m}}{\partial x_i}\right)}^2+Df(r)\mathcal{L}-K{m}_z^2+{\mathcal{F}}_{ME}+{\mathcal{F}}_E$$

$$\mathcal{F}=A(\frac{d\varphi }{dz}{)}^2+Df(z)(1+\gamma p_z)\frac{d\varphi }{dz}+\frac{p_z^2}{2{\chi }_E}-kcos(n\varphi ).$$

$$2A\frac{d^2\varphi }{d{z}^2}+D\frac{df}{dz}(1+\gamma p_z)$$

$$p_z$$

$$\mathcal{L}+{\mathcal{L}}_{me}=m_z\frac{\partial m_x}{\partial x}-m_x\frac{\partial m_z}{\partial x}+m_z\frac{\partial m_y}{\partial y}-m_y\frac{\partial m_z}{\partial y}+\gamma m_z(p_x\frac{\partial m_y}{\partial z}-p_y\frac{\partial m_x}{\partial z})$$

$$A\left[\frac{{\partial }^2\theta }{\partial z^2}+\frac{{\partial }^2\theta }{\partial {\rho }^2}+\frac{1}{\rho }\frac{\partial \theta }{\partial \rho }-\frac{sin\theta cos\theta }{{\rho }^2}\right]-Df(z)\frac{si{n}^2\theta }{\rho }-Ksin\theta cos\theta =0$$

$$p_x=-{\chi }_{E}Df\gamma \sin \varphi {\cos }^2\theta \frac{\partial \theta }{\partial z}and{p}_y={\chi }_{E}Df\gamma \cos \varphi {\cos }^2\theta \frac{\partial \theta }{\partial z}$$

$$m(\mathbf{l})=B_s\left(-\frac{JS}{k_{B}T}\sum _{\hat{\delta }}m(\mathbf{l}+\hat{\delta })\right)$$

$$\coth \theta =\frac{1}{\theta }+\frac{\theta }{3}-\frac{{\theta }^3}{45}+.......$$

$$B_s(x)=\frac{(s+1)}{3s}\{x-\frac{1}{15s^2}[s(s+1)+\frac{1}{2}]x^3+......\}.$$

$$\sum _{\delta }m(\mathbf{l}+\delta )\simeq 6m(\mathbf{l})+a_0^2{\nabla }^{2}m(\mathbf{l})$$

$$\frac{a_0^2}{6}\frac{{\partial }^{2}m(z)}{\partial z^2}+(1-\tau )m(z)-\beta m^3(z)=0$$

$$A\left(\frac{{\partial }^2\theta }{\partial {\rho }^2}+\frac{{\partial }^2\theta }{\partial z^2}+\frac{1}{\rho }\frac{\partial \theta }{\partial \rho }-\frac{\sin (\theta )\cos (\theta )}{{\rho }^2}\right)+Df(z)\frac{{\sin }^2(\theta )}{\rho }-K\sin (\theta )\cos (\theta )=0\text{mbox}in\Omega$$

$$A\left(\frac{{\partial }^2\theta }{\partial {\rho }^2}+\frac{{\partial }^2\theta }{\partial z^2}+\frac{1}{\rho }\frac{\partial \theta }{\partial \rho }\right)-\frac{1}{2}\left(\frac{A}{{\rho }^2}+K\right)\sin (2\theta )+D\frac{f(z)}{\rho }{\sin }^2(\theta )=0\text{mbox}in\Omega .$$

$$L\theta =A\left(\frac{{\partial }^2\theta }{\partial {\rho }^2}+\frac{{\partial }^2\theta }{\partial z^2}+\frac{1}{\rho }\frac{\partial \theta }{\partial \rho }\right)$$

$$f(\theta )=-\frac{1}{2}\left(\frac{A}{{\rho }^2}+K\right)\sin (2\theta )+D\frac{f(z)}{\rho }{\sin }^2(\theta )$$

$$L\theta +f(\theta )=0\text{mbox}in\Omega .$$

$$L{\theta }_{ij}+f({\theta }_{ij})=0i=2,\dots ,N_r,j=1,\dots ,N_z+1.$$

$$\frac{{\partial }^2\theta }{\partial {\rho }^2}$$

$$\frac{{\partial }^2\theta }{\partial z^2}$$

$$\frac{\partial \theta }{\partial \rho }$$

$$\begin{gathered} \left(\frac{A}{\Delta {\rho }^2}-\frac{A}{2(i-1)\Delta {\rho }^2}\right){\theta }_{i-1,j}+\left(\frac{A}{\Delta {\rho }^2}+\frac{A}{2(i-1)\Delta {\rho }^2}\right){\theta }_{i+1,j}-\left(\frac{2A}{\Delta {\rho }^2}+\frac{2A}{\Delta z^2}\right){\theta }_{ij} \\ +\frac{A}{\Delta z^2}{\theta }_{i,j-1}+\frac{A}{\Delta z^2}{\theta }_{i,j+1}+f({\theta }_{ij})=0i=2,\dots ,N_r,j=1,\dots ,N_z+1. \end{gathered}$$

$$\left(\begin{array}{cccc}{M}_{11}& M_{13}& M_{1b}& M_{1I}\\ M_{31}& M_{33}& M_{3b}& M_{3I}\\ M_{b1}& M_{b3}& M_{bb}& M_{bI}\\ M_{I1}& M_{I3}& M_{Ib}& M_{II}\end{array}\right)\left(\begin{array}{c}{\mathbf{\theta }}_1\\ {\mathbf{\theta }}_3\\ {\mathbf{\theta }}_b\\ {\mathbf{\theta }}_I\end{array}\right)+f(\mathbf{\theta })=0$$

$$\left(\begin{array}{cc}{M}_{11}+M_{13}& M_{1I}\\ M_{I1}+M_{I3}& M_{II}\end{array}\right)\left(\begin{array}{c}{\mathbf{\theta }}_1\\ {\mathbf{\theta }}_I\end{array}\right)+f\left(\left(\begin{array}{c}{\mathbf{\theta }}_1\\ {\mathbf{\theta }}_I\end{array}\right)\right)=-\left(\begin{array}{c}{M}_{1b}\\ M_{Ib}\end{array}\right){\mathbf{\theta }}_b$$

$$\hat{M}\hat{\mathbf{\theta }}+f(\hat{\mathbf{\theta }})=\mathbf{b}.$$

$$J_{\hat{M}}({\hat{\mathbf{\theta }}}^{(k)})\mathbf{\delta }{\hat{\mathbf{\theta }}}^{(k)}=-\left(\hat{M}{\hat{\mathbf{\theta }}}^{(k)}+f({\hat{\mathbf{\theta }}}^{(k)})-\mathbf{b}\right)$$

$${\hat{\mathbf{\theta }}}^{(k+1)}={\hat{\mathbf{\theta }}}^{(k)}+\mathbf{\delta }{\hat{\mathbf{\theta }}}^{(k)}$$

Full text

I Introduction.

Magnetoelectric multiferroic materials continue to attract much interest, due to both scientific as well as technological importance Khomskii ; Wang ; lottermoser . Typically the ferroelectric transition temperature is much higher than the magnetic one and a coupling between the two order parameters is weak. Representative examples are the transition metal perovskites BiFeO3, BiMnO3, belonging to type I class of multiferroics (more details on the classification Khomskii ; cheong ; eerenstein ; Ramesh ; ederer ). When the two ordering temperatures are close or even coincide, such as in TbMnO3 kimura or TbMn2O5 hur strong multiferroic behavior is expected. From the symmetry point of view, the necessity to break both inversion and time-reversal symmetries suggests different possible mechanisms which have been actually realised. In type II multiferroics, magnetism drives the onset of the ferroelectric order parameter, either due to the presence of spin-orbit coupling (SOC) and magnetic frustration e.g. in Ni3V2O6 lawes ; mostovoy ; katsura , or exchange striction e.g. in TbMnO3 and Ca3CoMnO6 sergienko ; Wu or ”phase dislocated” spin density waves e.g. in YMn2O5 chapon ; betouras .

Currently the role of surfaces and interfaces in the properties of materials is the focus of systematic studies Panagopoulos1 ; Gerhard ; Maruyama ; Zutic ; Hu . In technological applications there are important prospects and the theoretical understanding is developing Hellman ; Sander . Experimental techniques have been advanced such that novel phenomena can be detected as a result of the higher precision and resolution. Recent advancements led to the detection of new properties, by distinguishing surface from bulk phenomena or going to the atomic scale UO2 ; Sonntag .

In this work, we study the effects of boundaries in the development of multiferroic behavior in bulk materials, thin films and multilayers. This is a complementary effort to first-principles calculations on the magnetoelectric coupling close to surfaces Tsymbal or monolayers Nakamura of specific materials. A straightforward Ginzburg-Landau (GL) free-energy analysis with appropriate boundary conditions demonstrates that collinear magnetism can generate a ferroelectric polarization near surfaces, even without invoking the mechanism of phase dislocation betouras . In addition, due to the absence of inversion symmetry close to surfaces, a term that promotes the Dzyaloshinskii-Moriya interaction (DMI) can be present, leading to multiferroic behaviour through the formation of spiral magnetic order mostovoy ; katsura . The underlying assumption is that we deal with predominantly magnetic materials with non-zero coupling between magnetic and ferroelectric order parameters. In the following, we analyse separately the two mechanisms.

II Ginzburg-Landau analysis.

We focus on a simple cubic ferromagnet, for simplicity, without frustration. Each spin has an interaction with its six nearest neighbours, according to an isotropic ferromagnetic Heisenberg-type interaction $J$ and the expectation value of the z-component of the magnetisation (spin) at site ${\bf l}$, $m({\bf l})=\langle S_{z}({\bf l})\rangle/S$ is the order parameter of the system. The crystal is assumed to have a (001) surface. In the continuous space approximation, the GL equation of the magnetization $m({\bf l})$ which depends only on $z$, reads Mills (for completeness the details are in Appendix A):

[TABLE]

where $a_{0}$ is the lattice parameter and $\tau=\frac{T}{T_{c}}$ is the reduced temperature ($T_{c}$ is the Curie temperature) and $\beta=\frac{3}{5}[s(s+1)+\frac{1}{2}]/(s+1)^{2}$. In the limit of $(z\to\infty)$, far from the surface, the order parameter takes the bulk value $m_{\infty}=(1-\tau)^{1/2}/\beta^{1/2}$. The free energy expansion, leads to the same equation for both a ferromagnet or an antiferromagnet Mills . Defining $m(z)\equiv m_{\infty}f(z)$ and $\xi^{2}\equiv\frac{1}{6}\frac{a_{0}^{2}}{(1-\tau)}$ then:

[TABLE]

The boundary conditions are $f(\infty)=1$ and $f(0)=a_{0}\frac{\partial f}{\partial z}(0)$. Then the solution of Eq. (1) reads: $m(z)=m_{\infty}\tanh(\frac{z+a_{0}}{\sqrt{2}\xi})$.

We now address the question whether the onset of ferroelectricity is possible due to the existence of this surface which results in the change of the magnetisation. The physical argument is that as the inversion symmetry is broken due to the surface, the onset of multiferroic behavior is possible.

The free energy is supplemented by a term due to the coupling between magnetization and electric polarization p and a term which is the electronic part of the free energy that depends only on p. Since we are interested in systems where the magnetic order is the primary one, the second term is sufficient to be quadratic in p. The free energy reads:

[TABLE]

Note that, in the case of collinear magnetic structure the term proportional to $\gamma$ is non-zero, while the term proportional to $\gamma^{\prime}$ (a Lifshitz invariant for cubic lattice) is zero. Taking the dielectric susceptibility as constant and using $m(z)$ as the magnetic order parameter, then the minimisation with respect to $p(z)$ results in:

[TABLE]

Using the solution of Eq. (1) for the magnetisation, $p(z)$ becomes:

[TABLE]

This function is plotted in Fig. 1, interestingly it is peaked at a distance from the surface $z=\sqrt{2}\xi\tanh^{-1}(1/\sqrt{3})$ with $p_{max}=-\frac{2\sqrt{2}}{3\sqrt{3}}\frac{1}{\xi}\chi_{E}\gamma m_{\infty}^{2}\approx-\frac{1}{\xi}\chi_{E}\gamma m_{\infty}^{2}$, where the negative sign denotes a direction opposite to the direction of $\nabla({\bf m}^{2})$. The distance over which polarization is developed and the location of its peak is also controlled by the temperature through the magnetic correlation length. For a reasonable estimate of the effect we need the range of possible values of $\gamma\chi_{E}$ which in the literature Coupling1 ; Coupling2 ; Coupling3 ; Coupling4 can be found to lie between 10*-23* and 10*-13* sm/A, a typical correlation length Michels $\xi\approx$ 10 nm and a typical value for the bulk magnetization $m_{\infty}\approx$ 10-100 kA/m. Then the range of the possible values of $p_{max}$ is between 10*-10* to 102 $\mu$C/cm2 or to produce a measurable polarisation a value for $\gamma\chi_{E}$ larger than 10*-19* sm/A is enough explanation .

III DMI in magnetic thin films and multilayers.

The lack of inversion symmetry close to surfaces can lead to the induction of DMI. The direction of the d vector in this case can cause the directions of the spins of the nearest neighbours to change in such a way as to break the chiral symmetry close to the surface xia ; crepieux . For the purpose of our investigation of the chiral nature induced by the DMI, the antisymmetric exchange interaction is described by a Lifshitz invariant term which is linear in the spatial derivatives of the magnetisation m(r) of the form $m_{i}\frac{\partial m_{j}}{\partial x_{k}}-m_{j}\frac{\partial m_{i}}{\partial x_{k}}$, where $x_{l}$ denotes a spatial coordinate. These interactions are responsible for breaking the chiral symmetry and stabilizing localized magnetic vortices, with certain chirality which has been observed experimentally in noncentrosymmetric ferromagnetic and antiferromagnetic materials examples . It is possible to observe these effects in centrosymmetric crystals where stresses or applied magnetic fields Bogdanov-2001 or anisotropic frustrated magnetic interactions leonov induce chiral magnetic couplings and vortices/skyrmions. Chiral effects, as a consequence of the DMI energy are not so strong in the bulk, but they can become fundamentally important in magnetic thin films and multilayers or near the surface of a larger crystal where the local symmetry is low. Taking into account experimetal facts, the chiral couplings should also be inhomogeneous Bogdanov-2001 within a magnetic structure with low local symmetry. A phenomenological term for the corresponding chiral energy density is ${\mathcal{F}}_{D}=Df({\bf r}){\mathcal{L}}({\bf m})$ where $D$ is a constant, ${\mathcal{L}}$ is a Lifshitz invariant and the function $f({\bf r})$ is a function describing the inhomogeneous distribution of the magnetic chiral energy. $f({\bf r})$ was interpreted as another field, in addition to the magnetisation Bogdanov-2001 , but it essentially indicates the strength profile of the DMI as a function of the distance from the surface. To demonstrate the physics clearly, we take two functions as the profile function $f(z)$, where $z$ is the distance from the surface: a function exponentially decaying in $z$ and a function $1-tanh(z/\lambda)$, both with maximum at the surface. This behavior has been verified when the magnetisation in a finite-width slab was computed Bogdanov-2001 .

We consider a uniaxial magnetic anisotropy that contributes to the energy density a term ${\mathcal{F}}_{an}=-Km_{z}^{2}$. Then the free energy density reads:

[TABLE]

the first term represents the magnetic exchange interaction with a stiffness constant A (in , the second is the energy of the electric polarization with susceptibility $\chi_{E}$, the third term is due to DMI coupling, the fourth due to anisotropy and ${\mathcal{F}}_{ME}$ is the term of the free energy that couples ${\bf m}$ and ${\bf p}$. In the following all the quantities are dimensionless and for that purpose the magnetization $\vec{m}$ is normalized by its amplitude $M_{0}$, D is in units of stiffness A per m, $chi_{E}$ is in units of vacuum permittivity $\epsilon_{0}$, K is in units of A per m2, polarization p in units of $\sqrt{\epsilon_{0}\mu_{0}}M_{0}$ where $\mu_{0}$ is vacuum permeability, length is in units of lattice constant, $\gamma$ is in units of inverse polarization and the coupling constant in ${\mathcal{F}}_{ME}$ is the product D $\gamma$ (instead of only $\gamma$).

III.1 Anisotropic term with $K<0$

When $K<0$ in Eq. (5) then $\bf{m}$ lies in the $xy$ plane. The Lifshitz invariant ${\mathcal{L}}$ can be taken as $(m_{x}\frac{dm_{y}}{dz}-m_{y}\frac{dm_{x}}{dz})$. The free energy term that couples ${\bf m}$ and ${\bf p}$ then takes the form ${\bf p}\cdot({\bf m}\times\frac{d{\bf m}}{dz})$ which in this particular geometry reads: $D\gamma p_{z}(m_{x}\frac{dm_{y}}{dz}-m_{y}\frac{dm_{x}}{dz})$. It is convenient to work with the angle $\phi$ to describe the vector ${\bf m}$ that lies in the $xy$ plane. The related part of the free energy density becomes:

[TABLE]

The last term describes an in-plane anisotropy with $n$ an even integer, depending on lattice symmetry and/or homogeneous strain Bogdanov-2001 . Minimizing the free energy with respect to ${\bf m}$ and ${\bf p}$, we obtain:

[TABLE]

Inserting Eq. (7) into Eq. (6), we solve numerically for $\phi$ and $p_{z}$ as a function of the distance from the surface. The results are presented in Fig. 2. This physics provides a second mechanism to generate a finite polarization close to the surface which comes from the non-zero value of $\frac{d\phi}{dz}$ as a consequence of the DMI which is maximum at the surface. $f(z)$ is the profile that controls the strength of the DMI. We have checked that both profiles of $f(z)$ lead to the same physics qualitatively and we have used $1-tanh(z/\lambda)$ in Figs 2 and 3.

For an estimate of the effect, the range of values of the product $D\gamma\chi_{E}$ is taken between 10*-23* and 10*-13* sm/A, the intralayer spacing $\Delta z\approx\cal{O}$ (1 nm), typical values for the magnetization Panagopoulos1 ; Boulle ; Takeuchi $M_{0}\approx$ (10-1000) kA/m and typical $\Delta\phi\approx 0.1$ so that $\frac{\Delta\phi}{\Delta z}\approx$ 105 m*-1* while $f$ is of $\cal{O}$(1) close to the surface. Then the range of the possible values of the maximum polarization is $p_{z}\approx$ 10*-10* to 104 $\mu$C/cm2. Again with a reasonable value for ${\mathcal{D}}\gamma\chi_{E}$ of 10*-17* sm/A or larger, the polarization can be detected explanation .

III.2 Anisotropic term with $K>0$

When $K>0$ in the free energy, the relevant Lifshitz invariants may involve gradients along all three directions depending on the respected symmetry. There is a Lifshitz invariant term ${\mathcal{L}}$ purely magnetic as well as a term ${\mathcal{L}}_{me}$ that mixes ${\bf m}$ and ${\bf p}$ in the free energy. In the case of twofold or fourfold symmetry about the z axis we take the Lifshitz invariants to be:

[TABLE]

Using spherical coordinates for the magnetisation ${\bf m}=(sin\theta cos\phi,sin\theta sin\phi,cos\theta)$, cylindrical coordinates for the spatial vector ${\bf r}=(\rho cos\zeta,\rho sin\zeta,z)$ and focusing on the magnetic part of the free energy, the problem has axisymmetric localized solutions $\phi=\zeta$ and $\theta=\theta(\rho,z)$ with $\theta(0)=\pi$ and $\theta(\infty)=0$. The part of the free energy proportional to ${\mathcal{L}}_{me}$ reads: ${\mathcal{F}}_{me}=Df(z)\gamma\cos^{2}\theta\frac{\partial\theta}{\partial z}\left(p_{x}\sin\phi-p_{y}\cos\phi\right)$. Minimizing the free energy with respect to $\theta$ and $p_{x}$, $p_{y}$, we obtain:

[TABLE]

As ${\bf p}$ is the sub-dominant order parameter, to simplify the calculation it is sufficient to neglect its effect on $\theta$. Then the solution for ${\bf m}$ is similar to the purely magnetic problem Bogdanov-2001 . As a result, the magnitude of ${\bf p}$ at any given point is $p=\chi_{E}Df\gamma\cos^{2}\theta|\frac{\partial\theta}{\partial z}|$. In Fig. (3) we present $\theta$ as a function of $\rho$ and $z$ as well as $p_{x}$ and $p_{y}$ (differing by a phase difference) (more details in Appendix B). Similarly to the case (i), using similar typical values and reasonable estimate of $|\frac{\partial\theta}{\partial z}|\approx$ 104 m*-1*, then the range of the possible values of the amplitude of ${\bf p}$ is $p\approx$ 10*-8* to 104 $\mu$C/cm2, which suggests a detectable value explanation .

IV Discussion.

We present a detailed study of the multiferroic behavior in thin films, multilayers and close to the surface of magnetic materials with symmetries (inversion and time reversal) that do not allow the onset of polarization in the bulk. The onset of polarization can be achieved both in the presence or absence of a DMI close to the surface. One mechanism is through the non-zero gradient of the magnetization, while the second is through the change of orientation of the magnetization (and not its amplitude) as a function of the distance from the surface. Our study is essentially the inverse effect of Ref. (valencia, ), where surface-induced magnetisation was detected in the archetypal ferroelectric BaTiO3. In the same spirit, there exist recent intensive efforts to synthesize multiferroic heterostructures and control the interfacial DMI (e.g. Refs. [Vaz, ; Jain, ; hu, ; Ding, ; Hrabec, ; Kuepferling, ]). In thin films DMI originating from strong SOC of interfacial atoms neighboring the magnetic layer can be engineered and controlled. DMI has been engineered in a ferromagnet (FM) interfaced with two different heavy metals, such as Pt/Co/Ir Chen ; Moreau-Luchaire , or in magnetic layers inserted between a heavy metal and an oxide, such as Pt/CoFe/MgO Emori ; Emori2 . In thin films, the deposition conditions would make the couplings vary, given that the strength of the DMI has been recently shown to vary in a model system Pt/Co/Pt, depending e.g. on the temperature variation during deposition Wells . An elegant and reliable method for determining the magnitude of the DMI from static domain measurements even in the presence of hybrid chiral structures was recently demonstrated Legrand while electrical detection of single magnetic skyrmions has been achieved at room temperature in metallic multilayers Maccariello .

We estimate the range of values of the polarization, that can be detected even if the coupling constant that mixes the magnetic and ferroelectric order parameters is several orders of magnitude smaller than the highest ones reported in the literature. The polarization we predict is detectable within the current accuracy of the experimental techniques. Techniques such as polarised neutron reflectivity can be used to detect magnetisation gradients over a length scale of a few nanometers. In case of out-of-plane components, magnetic force microscopy is ideal to detect the magnetisation components. X-ray magnetic circular/linear dichroism can detect the contribution of individual magnetic ions if required. Scanning transmission X-ray microscopy, has been used effectively to study skyrmion dynamics high temporal and spatial resolution Stoll . Recently, the linear magnetoelectric phase in ultrathin MnPS3 was probed by optical second harmonic generation Chu . These techniques in combination with e.g. electrostatic force microscopy or ellipsometry make feasible the detection of both order parameters.

Acknowledgments.

We thank N. Banerjee, A. Bogdanov, P. Borisov, D. Efremov, N. Gidopoulos, T. Hesjedal, D. Khomskii, P. King, P. Radaelli, I. Rousochatzakis and J. van den Brink for useful discussions and communications. The work is supported partly by EPSRC through the grant EP/P003052/1 (JJB) and a scholarship from the Regional Government of Kurdistan-Iraq (ART). JJB also thanks the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the programme ”Mathematical design of new materials”, supported by EPSRC grant no EP/R014604/1, where a part of the work was done.

Appendix A: Derivation of the Ginzburg-Landau equation.

Using molecular field theory, the expectation value of the z-component of the magnetisation (spin) at site ${\bf l}$, $m({\bf l})=\langle S_{z}({\bf l})\rangle/S$ which is the order parameter of the system, is written as:

[TABLE]

where $k_{B}$ is Boltzmann’s constant, $T$ is the temperature, $S$ is the spin of the magnetic ions and $B_{s}(x)$ is the Brillouin function $B_{s}(x)=(1+\frac{1}{2S})\coth[(1+\frac{1}{2S})x]-\frac{1}{2S}\coth(\frac{x}{2S})$. The sum over ${\hat{\bf\delta}}$ ranges over the six nearest neighbours of the spin at site ${\bf l}$. We can then expand $\coth(\theta)$ in powers of $\theta$:

[TABLE]

The Brillouin function then reads:

[TABLE]

We define the reduced temperature $\tau=\frac{T}{T_{c}}$ (for the model we discuss $T_{c}$ is the Curie temperature $T_{c}=2k_{B}J(S+1)/3S$). The order parameter is then defined through the relations $m({\bf l})=<S({\bf l})>/S$ with $m({\bf l})=B_{s}\left(-\frac{\sum{m({\bf l}+{\bf\delta})}}{2\tau k_{B}(s+1)}\right)$. In the continuous space approximation of the lattice, $m({\bf l})$ becomes also function of the continuous ${\bf l}$ and

[TABLE]

where $a_{0}$ is the lattice parameter. Retaining only first order in $\bigtriangledown^{2}m$ terms and using $\beta=\frac{3}{5}[s(s+1)+\frac{1}{2}]/(s+1)^{2}$ and the fact that $m({\bf r})$ depends only on $z$, the GL equation becomes:

[TABLE]

which is Eq. (1) of the main text. The free energy expansion, leads to the same equation for both a ferromagnet or an antiferromagnet.

Appendix B: Details of the calculation of magnetisation for K$>$0.

As explained in the text, the magnetisation in the case K$>$0 is determined by the function $\theta=\theta(\rho,z)$ such that

[TABLE]

where $\Omega=(0,r)\times(-l,l)$. The boundary conditions are

on $\Gamma_{4}=\{(\rho,z)\,:\,\rho=0,\;-l\leq z\leq l\}$, we impose $\theta=\pi$

• -

on $\Gamma_{2}=\{(\rho,z)\,:\,\rho=r,\;-l\leq z\leq l\}$, we impose $\theta=0$

• -

on $\Gamma_{1}=\{(\rho,z)\,:\,0\leq\rho\leq r,\;z=-l\}$ and on $\Gamma_{3}=\{(\rho,z)\,:\,0\leq\rho\leq r,\;z=l\}$, we impose periodic boundary conditions: $\theta(\rho,z=-l)=\theta(\rho,z=l)$.

Equation (15) can be equivalently reformulated as

[TABLE]

Defining the linear operator

[TABLE]

and the nonlinear operator

[TABLE]

we can write the problem in abstract form as: find $\theta$ such that

[TABLE]

Numerical solution using the finite difference method

We introduce the computational grid as the one shown in the following figure:

Let $\Delta\rho$ and $\Delta z$ be the discretisation steps along the $\rho$ and the $z$ axes, respectively, and let $N_{r}$ and $N_{z}$ be the number of intervals in the $\rho$ and $z$ direction. Thus, the point $(i,j)$ has coordinates $((i-1)\Delta\rho,(j-1)\Delta z)$, $i=1,\ldots,N_{r}+1$, $j=1,\ldots,N_{z}+1$.

We also introduce a global numbering of the nodes in such a way that the node $(i,j)$ is associated to the number $k=(j-1)(N_{r}+1)+i$.

We now discretise the equation (19) using finite differences. Considering the nature of the boundary conditions, we will write the equation in all the nodes except those on the boundaries $\Gamma_{2}$ and $\Gamma_{4}$ where we impose Dirichlet boundary conditions. More precisely, we will consider

[TABLE]

Discretisation of the linear term

To discretise the linear term $L$ we introduce the following second-order approximations of the derivatives

[TABLE]

Then, with the help of some algebra, the discretisation of (19) becomes

[TABLE]

Considering that on $\Gamma_{1}$ and $\Gamma_{3}$ we impose periodic boundary conditions, we identify the nodes characterised by indices $(i,1)$ and $(i,N_{z}+1)$ (shown in blue in the figure above). Moreover, we make the assumption that we can identify the nodes $(i,0)$ with the nodes $(i,N_{z})$ (green nodes), and the nodes $(i,N_{z}+1)$ with the nodes $(i,2)$ (red nodes). Thus, in the matrix $M$ associated to (23) we can identify

[TABLE]

and we can reduce the system to

[TABLE]

For simplicity and with obvious choice of notation, we can rewrite (24) as

[TABLE]

Newton’s method for the nonlinear system

To solve the nonlinear system (25), we consider now the Newton’s method: for $k\geq 0$ until convergence, we solve the linear system

[TABLE]

and set

[TABLE]

where $J_{\widehat{M}}(\widehat{\boldsymbol{\theta}}^{(k)})$ is the Jacobian matrix in $\widehat{\boldsymbol{\theta}}^{(k)}$ defined as

[TABLE]

and $f^{\prime}$ is the Gâteaux derivative of the nonlinear operator $f$ defined as

[TABLE]