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

Ahmed R. Tarkhany, 1 Marco Discacciati, 2 and Joseph J. Betouras 1, ∗

1 Department of Physics and Centre for the Science of Materials,

Loughborough University, Loughborough LE11 3TU, UK.

2 Department of Mathematical Sciences, Loughborough University, Loughborough LE11 3TU, UK.

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.

## I. INTRODUCTION.

Magnetoelectric multiferroic materials continue to attract much interest, due to both scientific as well as technological importance 1-3 . 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 BiFeO , 3 BiMnO , belonging to type I class of multiferroics 3 (more details on the classification 1,4-7 ). When the two ordering temperatures are close or even coincide, such as in TbMnO 3 8 or TbMn O 2 5 9 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 Ni 3 V O 2 6 10-12 , or exchange striction e.g. in TbMnO 3 and Ca CoMnO 3 6 13,14 or 'phase dislocated' spin density waves e.g. in YMn O 2 5 15,16 .

Currently the role of surfaces and interfaces in the properties of materials is the focus of systematic studies 17-21 . In technological applications there are important prospects and the theoretical understanding is developing 22,23 . 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 24,25 .

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 26 or monolayers 27 of specific materials. A straightforward Ginzburg-Landau (GL) free-energy analysis with appropriate boundary conditions demonstrates that collinear magnetism can generate a ferroelec- tric polarization near surfaces, even without invoking the mechanism of phase dislocation 16 . In addition, due to the absence of inversion symmetry close to surfaces, a term that promotes the DzyaloshinskiiMoriya interaction (DMI) can be present, leading to multiferroic behaviour through the formation of spiral magnetic order 11,12 . 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 zcomponent of the magnetisation (spin) at site l , m ( ) = l 〈 S z ( ) l 〉 /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 ( ) l which depends only on z , reads 28 (for completeness the details are in Appendix A):

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

where a 0 is the lattice parameter and τ = T T c is the reduced temperature ( T c is the Curie temperature) and β = 3 5 [ s s ( + 1) + 1 2 ] / s ( + 1) . 2 In the limit of ( z → ∞ ), far from the surface, the order parameter takes the bulk value m ∞ = (1 -τ ) 1 / 2 /β 1 / 2 . The free energy expansion, leads to the same equation for both a ferromagnet or an antiferromagnet 28 . Defining m z ( ) ≡ m f z ∞ ( ) and ξ 2 ≡ 1 6 (1 a 2 0 -τ ) then:

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

FIG. 1: Polarization as a function of the distance from the surface.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00602v2-figure-5ab1376e-69da-4696-bc4b-3817076a8114.jpg)

The boundary conditions are f ( ∞ ) = 1 and f (0) = a 0 ∂f ∂z (0). Then the solution of Eq. (1) reads: m z ( ) = m ∞ tanh( z + a 0 √ 2 ξ ).

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:

$$\delta F = F _ { M E } + F _ { E } = \mathbf p \cdot ( \gamma \nabla ( \mathbf m ^ { 2 } ) \quad \begin{pmatrix} \text{cnra} \\ \text{in nc} \\ + \gamma ^ { \prime } [ \mathbf m ( \nabla \cdot \mathbf m ) - ( \mathbf m \cdot \nabla ) \mathbf m ] + \dots ) + \frac { \mathbf p ( \mathbf r ) ^ { 2 } } { 2 \chi _ { E } } \quad ( 3 ) \quad \text{edge} \\ \text{app} \quad \end{pmatrix}$$

Note that, in the case of collinear magnetic structure the term proportional to γ is non-zero, while the term proportional to γ ′ (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:

$$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 ] + \dots ) \text{ within a}$$

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

$$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 ] ^ { \text{genous} } \text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$	\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\tt$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$
\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\tiff$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$                                                                                                                                                                                                         \text$                                                                                                                                                                                                      \text                                                                                                                                                                                                       
                                                                                                                                                                                                        }                                                                                                                                                                                                       .                                                                                                                                                                                                      .$$

This function is plotted in Fig. 1, interestingly it is peaked at a distance from the surface z =

√ 2 ξ tanh -1 (1 / √ 3) with p max = -2 √ 2 1 3 √ 3 ξ χ E γm 2 ∞ ≈ -1 ξ χ E γm 2 ∞ , where the negative sign denotes a direction opposite to the direction of ∇ ( 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 γχ E which in the literature 29-32 can be found to lie between 10 -23 and 10 -13 sm/A, a typical correlation length 33 ξ ≈ 10 nm and a typical value for the bulk magnetization m ∞ ≈ 10-100 kA/m. Then the range of the possible values of p max is between 10 -10 to 10 2 µ C/cm 2 or to produce a measurable polarisation a value for γχ E larger than 10 -19 sm/A is enough 34 .

## III. DMI IN MAGNETIC THIN FILMS AND MULTILAYERS.

) within a magnetic structure with low local symme-] 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 35,36 . 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 ∂m j ∂x k -m j ∂m i ∂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 37 . It is possible to observe these effects in centrosymmetric crystals where stresses or applied magnetic fields 38 or anisotropic frustrated magnetic interactions 39 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 38 try. A phenomenological term for the corresponding chiral energy density is F D = Df ( r ) L ( m ) where D is a constant, L is a Lifshitz invariant and the function f ( r ) is a function describing the inhomogeneous distribution of the magnetic chiral energy. f ( r ) was interpreted as another field, in addition to the magnetisation 38 , 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

FIG. 2: magnetisation ( φ ) and polarization as a function of the distance from the surface for in-plane anisotropy with n = 0 (a) and n = 4 (b). The used parameters are A=1, D=1, k=1, γ =1, χ E =1.5. The peak of the polarization is at a distance from the surface in agreement with Fig. 1.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00602v2-figure-697fa030-9511-4bc4-a515-76e06031b55b.jpg)

f ( z ), where z is the distance from the surface: a function exponentially decaying in z and a function 1 -tanh z/λ ( ), both with maximum at the surface. This behavior has been verified when the magnetisation in a finite-width slab was computed 38 .

We consider a uniaxial magnetic anisotropy that contributes to the energy density a term F an = -Km 2 z . Then the free energy density reads:

$$\mathcal { F } = A \sum _ { i } \left ( \frac { \partial \mathbf m } { \partial x _ { i } } \right ) ^ { 2 } + D f ( r ) \mathcal { L } - K m _ { z } ^ { 2 } + \mathcal { F } _ { M E } + \mathcal { F } _ { E } ( 5 )$$

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 χ E , the third term is due to DMI coupling, the fourth due to anisotropy and F ME is the term of the free energy that couples m and p . In the following all the quantities are dimensionless and for that purpose the magnetization glyph[vector] 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 glyph[epsilon1] 0 , K is in units of A per m , polarization p in units of 2 √ glyph[epsilon1] 0 µ M 0 0 where µ 0 is vacuum permeability, length

For an estimate of the effect, the range of values of the product Dγχ E is taken between 10 -23 and 10 -13 sm/A, the intralayer spacing ∆ z ≈ O (1 nm), typical values for the magnetization 17,44,45 M 0 ≈ (10-1000) kA/m and typical ∆ φ ≈ 0 1 so that . ∆ φ ∆ z ≈ 10 5 m -1 while f is of O (1) close to the surface. Then the range of the possible values of the maximum polarization is p z ≈ 10 -10 to 10 4 µ C/cm . Again with a 2 reasonable value for D γχ E of 10 -17 sm/A or larger, the polarization can be detected 34 .

is in units of lattice constant, γ is in units of inverse polarization and the coupling constant in F ME is the product D γ (instead of only γ ).

## A. Anisotropic term with K &lt; 0

When K &lt; 0 in Eq. (5) then m lies in the xy plane. The Lifshitz invariant L can be taken as ( m x dm y dz -m y dm x dz ). The free energy term that couples m and p then takes the form p · ( m × d m dz ) which in this particular geometry reads: Dγp z ( m x dm y dz -m y dm x dz ). It is convenient to work with the angle φ to describe the vector m that lies in the xy plane. The related part of the free energy density becomes:

$$\underset { \text{tection} } { c e \colon a } \quad \mathcal { F } = A ( \frac { d \phi } { d z } ) ^ { 2 } + D f ( z ) ( 1 + \gamma p _ { z } ) \frac { d \phi } { d z } + \frac { p _ { z } ^ { 2 } } { 2 \chi _ { E } } - k \cos ( n \phi ).$$

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

$$\begin{array} { c c c } ^ { n } - \\ & 2 A \frac { d ^ { 2 } \phi } { d z ^ { 2 } } + D \frac { d f } { d z } ( 1 + \gamma p _ { z } ) & + D f \gamma \frac { d p _ { z } } { d z } - n k s i n ( n \phi ) \\ \cdot \mathcal { F } _ { E } ( 5 ) & = 0 \\ \cdot \cdot \end{array}$$

$$p _ { z } \ = - D \gamma \chi _ { E } f \frac { d \phi } { d z } \quad ( 7 )$$

Inserting Eq. (7) into Eq. (6), we solve numerically for φ 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 dφ 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/λ ( ) in Figs 2 and 3.

## B. Anisotropic term with K &gt; 0

When K &gt; 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 L purely magnetic as well as a term L me that mixes m and p in the free energy. In the case of twofold or fourfold symmetry about the z axis we take the Lifshitz invariants to be:

({∐√̂/(√{]√√}

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00602v2-figure-f0c2b4c7-32e9-4873-a839-30604d462854.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00602v2-figure-bd42f581-4545-4780-a3ee-bd6986a65dd2.jpg)

({∐√̂/(√{]√√}

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00602v2-figure-ee47e140-c22f-43a1-b8eb-3f1e71dcc720.jpg)

∖/

∖/

∖/

({∐√̂/(√{]√√}

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00602v2-figure-a4ffc853-0c2a-482c-886b-f5309ecafb9b.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00602v2-figure-7fe8b376-1995-4889-8184-48e5c85d663c.jpg)

∖/

FIG. 3: Numerical results for K = 1. In (a) and (b) θ as a function of z and ρ with boundary conditions as explained in the text. The surface are at z=L and z=-L, while the bulk is at z=0. (c) The component of polarization p y as a function of z and ρ , (d) p y at fixed ρ = 0.05 R and (e) at fixed z= 0.96 L (and L=4 λ ), where it shows oscillatory behavior as a function of ρ and φ for the second layer from the surface. The x-component of the polarization is out of phase by 3 π/ 2 in φ . The rest of the parameters are the same as in Fig. 2.

$$\mathcal { L } + \mathcal { L } _ { m e } = 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 } )$$

Using spherical coordinates for the magnetisation m = ( sinθcosφ, sinθsinφ, cosθ ), cylindrical coordinates for the spatial vector r = ( ρcosζ, ρsinζ, z ) and focusing on the magnetic part of the free energy, the problem has axisymmetric localized solutions φ = ζ

and θ = θ ρ, z ( ) with θ (0) = π and θ ( ∞ ) = 0. The part of the free energy proportional to L me reads: F me = Df z γ ( ) cos 2 θ ∂θ ∂z ( p x sin φ -p y cos φ ). Minimizing the free energy with respect to θ and p x , p y , we obtain:

$$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 { s i n \theta \cos \theta } { \rho ^ { 2 } } \right ] - D f ( z ) \frac { s i n ^ { 2 } \theta } { \rho } - K s i n \theta \cos \theta = 0 \\ p _ { x } & = - \chi _ { E } D f \gamma \sin \phi \cos ^ { 2 } \theta \frac { \partial \theta } { \partial z } \quad a n d \quad p _ { y } = \chi _ { E } D f \gamma \cos \phi \cos ^ { 2 } \theta \frac { \partial \theta } { \partial z }$$

As p is the sub-dominant order parameter, to simplify the calculation it is sufficient to neglect its effect on θ . Then the solution for m is similar to the purely magnetic problem 38 . As a result, the magnitude of p at any given point is p = χ E Dfγ cos 2 θ | ∂θ ∂z | . In Fig. (3) we present θ as a function of ρ 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 | ∂θ ∂z | ≈ 10 4 m -1 , then the range of the possible values of the amplitude of p is p ≈ 10 -8 to 10 4 µ C/cm , 2 which suggests a detectable value 34 .

## 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. (46), where surface-induced magnetisation was detected in the archetypal ferroelectric BaTiO 3 . In the same spirit, there exist recent intensive efforts to synthesize multiferroic heterostructures and control the interfacial DMI (e.g. Refs. [49-54]). 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 55,56 , or in magnetic layers inserted between a heavy metal and an oxide, such as Pt/CoFe/MgO 57,58 . 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 61 . 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 59 while electrical detection of single magnetic skyrmions has been achieved at room temperature in metallic multilayers 60 .

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 62 . Recently, the linear magnetoelectric phase in ultrathin MnPS 3 was probed by optical second harmonic generation 63 . 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 KurdistanIraq (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 l , m ( ) = l 〈 S z ( ) l 〉 /S which is the order parameter of the system, is written as:

$$m ( \mathbf 1 ) = B _ { s } \left ( - \frac { J S } { k _ { B } T } \sum _ { \hat { \delta } } m ( \mathbf 1 + \hat { \delta } ) \right ) \quad ( 1 1 )$$

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 + 1 2 S ) coth[(1 + 1 2 S ) x ] -1 2 S coth( x 2 S ). The sum over ˆ ranges over the δ six nearest neighbours of the spin at site l . We can then expand coth( θ ) in powers of θ :

$$\coth \theta = \frac { 1 } { \theta } + \frac { \theta } { 3 } - \frac { \theta ^ { 3 } } { 4 5 } + \dots \dots \quad \quad ( 1 2 )$$

The Brillouin function then reads:

$$\underset { \text{all} } { \text{order} } \quad B _ { s } ( x ) = \frac { ( s + 1 ) } { 3 s } \{ x - \frac { 1 } { 1 5 s ^ { 2 } } [ s ( s + 1 ) + \frac { 1 } { 2 } ] x ^ { 3 } + \dots \}. \\ \text{The} \quad \text{with} \quad, \quad \cdot \quad, \quad \cdot \quad \cdot \quad T \quad \text{$(13)}$$

We define the reduced temperature τ = T T c (for the model we discuss T c is the Curie temperature T c = 2 k B J S ( +1) 3 / S ). The order parameter is then defined through the relations m ( ) l = &lt; S ( ) l &gt; /S with m ( ) l = B s ( -∑ m ( + ) l δ 2 τk B ( s +1) ) . In the continuous space approximation of the lattice, m ( ) l becomes also function of the continuous l and

$$\sum _ { \delta } m ( 1 + \delta ) \simeq 6 m ( 1 ) + a _ { 0 } ^ { 2 } \nabla ^ { 2 } m ( 1 ) \quad ( 1 4 )$$

where a 0 is the lattice parameter. Retaining only first order in glyph[triangleinv] 2 m terms and using β = 3 5 [ s s ( +1)+ 1 2 ] / s ( +1) 2 and the fact that m ( r ) depends only on z , the GL equation becomes:

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

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 &gt; 0.

As explained in the text, the magnetisation in the case K &gt; 0 is determined by the function θ = θ ρ, z ( ) such that

$$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 ) + D f ( z ) \frac { \sin ^ { 2 } ( \theta ) } { \rho } - K \sin ( \theta ) \cos ( \theta ) = 0 \ \ \text{in } \Omega$$

where Ω = (0 , r ) × -( l, l ). The boundary conditions are

- -on Γ 4 = ( { ρ, z ) : ρ = 0 , - ≤ l z ≤ } l , we impose θ = π
- -on Γ 2 = ( { ρ, z ) : ρ = r, - ≤ l z ≤ } l , we impose θ = 0
- -on Γ 1 = { ( ρ, z ) : 0 ≤ ρ ≤ r, z = - } l and on Γ 3 = { ( ρ, z ) : 0 ≤ ρ ≤ r, z = l } , we impose periodic boundary conditions: θ ρ, z ( = -l ) = θ ρ, z ( = ). l

Equation (15) can be equivalently reformulated as

$$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{ in } \Omega \,.$$

Defining the linear operator

$$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 )$$

and the nonlinear operator

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

we can write the problem in abstract form as: find θ such that

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

## Numerical solution using the finite difference method

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

Let ∆ ρ and ∆ z be the discretisation steps along the ρ and the z axes, respectively, and let N r and N z be the number of intervals in the ρ and z direction. Thus, the point ( i, j ) has coordinates (( i -1)∆ ρ, ( j -1)∆ ), z i = 1 , . . . , N r +1, j = 1 , . . . , 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 Γ 2 and Γ 4 where we impose Dirichlet boundary conditions. More precisely, we will consider

$$L \theta _ { i j } + f ( \theta _ { i j } ) = 0 \quad i = 2, \dots, N _ { r }, \, j = 1, \dots, N _ { z } + 1 \,.$$

## Discretisation of the linear term

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

$$\frac { \partial ^ { 2 } \theta } { \partial \rho ^ { 2 } } \, \approx \, \frac { \theta _ { i - 1, j } - 2 \theta _ { i j } + \theta _ { i + 1, j } } { \Delta \rho ^ { 2 } }$$

$$\frac { \stackrel { \cdot } { \partial ^ { 2 } \theta } } { \partial z ^ { 2 } } \, \approx \, \frac { \theta _ { i, j - 1 } - \stackrel { \cdot } { 2 \theta _ { i j } } + \theta _ { i, j + 1 } } { \Delta z ^ { 2 } } \, \\ \stackrel { \cdot } { \omega } }$$

$$\frac { \overset { \cdot \cdot } { \partial \theta } } { \partial \rho } \, \approx \, \frac { \theta _ { i + 1, j } - \theta _ { i - 1, j } } { 2 \Delta \rho }$$

FIG. 4: Computational grid.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00602v2-figure-a1830a4d-5fbe-4242-8b4c-b7fd2783450a.jpg)

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

$$\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 { 2 A } { \Delta \rho ^ { 2 } } + \frac { 2 A } { \Delta z ^ { 2 } } \right ) \theta _ { i j } \\ + \frac { A } { \Delta z ^ { 2 } } \theta _ { i, j - 1 } + \frac { A } { \Delta z ^ { 2 } } \theta _ { i, j + 1 } + f ( \theta _ { i j } ) = 0 \quad i = 2, \dots, N _ { r }, \ j = 1, \dots, N _ { z } + 1.$$

Considering that on Γ 1 and Γ 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

$$\begin{pmatrix} M _ { 1 1 } & M _ { 1 3 } & M _ { 1 b } & M _ { 1 I } \\ M _ { 3 1 } & M _ { 3 3 } & M _ { 3 b } & M _ { 3 I } \\ M _ { b 1 } & M _ { b 3 } & M _ { b b } & M _ { b I } \\ M _ { I 1 } & M _ { I 3 } & M _ { I b } & M _ { I I } \end{pmatrix} \begin{pmatrix} \theta _ { 1 } \\ \theta _ { 3 } \\ \theta _ { b } \\ \theta _ { I } \end{pmatrix} + f ( \theta ) = 0$$

and we can reduce the system to

$$\begin{pmatrix} M _ { 1 1 } + M _ { 1 3 } & M _ { 1 I } \\ M _ { I 1 } + M _ { I 3 } & M _ { I I } \end{pmatrix} \begin{pmatrix} \theta _ { 1 } \\ \theta _ { I } \end{pmatrix} + f \begin{pmatrix} \theta _ { 1 } \\ \theta _ { I } \end{pmatrix} \begin{pmatrix} = - \begin{pmatrix} M _ { 1 b } \\ M _ { I b } \end{pmatrix} \theta _ { b }$$

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

$$M \ddot { \theta } + f ( \dot { \theta } ) = b \,.$$

## Newton's method for the nonlinear system

To solve the nonlinear system (25), we consider now the Newton's method: for k ≥ 0 until convergence, we solve the linear system

$$J _ { \widehat { M } } ( \widehat { \theta } ^ { ( k ) } ) \, \delta \widehat { \theta } ^ { ( k ) } = - \left ( \widehat { M } \, \widehat { \theta } ^ { ( k ) } + f ( \widehat { \theta } ^ { ( k ) } ) - b \right )$$

and set

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

where J ̂ M ( ̂ θ ( k ) ) is the Jacobian matrix in ̂ θ ( k ) defined as

$$J _ { \widehat { M } } ( \widehat { \theta } ^ { ( k ) } ) = \widehat { M } + f ^ { \prime } ( \widehat { \theta } ^ { ( k ) } )$$

and f ′ is the Gˆteaux derivative of the nonlinear operator a f defined as

$$f ^ { \prime } ( \theta _ { i j } ) = - \left ( \frac { A } { ( ( i - 1 ) \Delta r ) ^ { 2 } } + K \right ) \cos ( 2 \theta _ { i j } ) + D \frac { f ( ( j - 1 ) \Delta z ) } { ( i - 1 ) \Delta r } \sin ( 2 \theta _ { i j } ) \,.$$

$$\begin{smallmatrix}. \end{smallmatrix}$$

- ∗ Electronic address: J.Betouras@lboro.ac.uk
- Panagopoulos, Nature 539 , 510 (2016).
- 1 D. Khomskii, Physics 2 , 20 (2009).
- 2 K. F. Wang, J.-M. Liu, and Z. F. Ren, Adv. in Phys. 58 , 321 (2009).
- 3 T. Lottermoser, T. Lonkai, U. Amann, D. Hohlwein, J. Ihringer, and M. Fiebig, Nature 430 , 541 (2004).
- 4 S.-W Cheong and M. Mostovoy, Nat. Mater. 6 , 13 (2007).
- 5 W. Eerenstein, N. D. Mathur, and J. F. Scott, Nature 442 , 759 (2006).
- 6 R. Ramesh and N. A Spaldin, Nat. Mater. 6 , 21 (2007).
- 7 C. Ederer and N. A. Spaldin, Nat. Mater. 3 , 849 (2004).
- 8 T. Kimura, T. Goto, H. Shintani, K. Ishizaka, T. Arima, and Y. Tokura, Nature 426 , 55 (2003).
- 9 N. Hur, S. Park, P. A. Sharma, J. S. Ahn, S. Guha, and S-W. Cheong, Nature 429 , 392 (2004).
- 10 G. Lawes, A. B. Harris, T. Kimura, N. Rogado, R. J. Cava, A. Aharony, O. Entin-Wohlman, T. Yildirim, M. Kenzelmann, C. Broholm, and A. P. Ramirez, Phys. Rev. Lett. 95 , 087205 (2005).
- 11 H. Katsura, N. Nagaosa, and A. V. Balatsky, Phys. Rev. Lett. 95 , 057205 (2005).
- 12 M. Mostovoy, Phys. Rev. Lett. 96 , 067601 (2006).
- 13 H. Wu T. Burnus, Z. Hu, C. Martin, A. Maignan, J. C. Cezar, A. Tanaka, N. B. Brookes, D. I. Khomskii, and L. H. Tjeng, Phys. Rev. Lett. 102 , 026404 (2009).
- 14 I. A. Sergienko, C. Sen, and E. Dagotto, Phys. Rev. Lett. 97 , 227204 (2006).
- 15 L. C. Chapon, P. G. Radaelli, G. R. Blake, S. Park, and S. W. Cheong, Phys. Rev. Lett. 96 , 097601 (2006).
- 16 J. J. Betouras, G. Giovannetti, and J. van den Brink, Phys. Rev. Lett. 98 , 257602 (2007).
- 17 A. Soumyanarayanan, N. Reyren, A. Fert, and C.
- 18 L. Gerhard, T. K. Yamada, T. Balashov, A. F. Taka´ cs, R. J. H. Wesselink, M. Daene, M. Fechner, S. Ostanin, A. Ernst, I. Mertig and W. Wulfhekel, Nat. Nanotech. 5 , 792-797 (2010).
- 19 T. Maruyama, Y. Shiota, T. Nozaki, K. Ohta, N. Toda, M. Mizuguchi, A. A. Tulapurkar, T. Shinjo, M. Shiraishi, S. Mizukami, Y. Ando and Y. Suzuki, Nat. Nanotech. 4 ,158-161 (2009).
- 20 I. Zˆ uti´, c A. Matos-Abiague, B. Scharf, H. Dery, K. Belashchenko, Materials Today 22 , 85 (2019).
- 21 J.-M. Hu, C.-G. Duan, C.-W. Nan, and L.-Qi. Chen npj Computational Materials 3 , 18 (2017).
- 22 F. Hellman, A. Hoffmann, Y. Tserkovnyak , G. S. D. Beach, E. E. Fullerton, C. Leighton, A. H. MacDonald, D. C. Ralph, D. A. Arena , H. A . Durr, P. Fischer, J. Grollier, J. P. Heremans, T. Jungwirth, A. V. Kimel, B. Koopmans, I. N. Krivorotov, S. J. May, A. K. Petford-Long, J. M. Rondinelli, N. Samarth, I. K. Schuller, A. N. Slavin, M. D. Stiles, O. Tchernyshyov, A. Thiaville, and B. L. Zink, Rev. Mod. Phys. 89 , 025006 (2017).
- 23 D. Sander, S. O. Valenzuela, D. Makarov, C. H. Marrows , E. E. Fullerton, P. Fischer, J. McCord, P. Vavassori, S. Mangin, P. Pirro, B. Hillebrands, A. D. Kent, T. Jungwirth, O. Gutfleisch, C. G. Kim and A. Berger, J. Phys. D: Appl. Phys. 50 363001 (2017).
- 24 S. Langridge, G. M. Watson, D. Gibbs, J. J. Betouras, N. I. Gidopoulos, F. Pollmann, M. W. Long, C. Vettier and G. H. Lander Phys. Rev. Lett. 112 , 167201 (2014).
- 25 A. Sonntag, J. Hermenau, A. Schlenhoff, J. Friedlein, S. Krause, and R. Wiesendanger, Phys. Rev. Lett. 112 , 017204 (2014).
- 26 C. G. Duan, J. P. Velev, R. F. Sabirianov Z. Zhu, J. Chu, S. S. Jaswal, and E. Y. Tsymbal, Phys. Rev.

- Lett. 101 , 137201 (2008).
- 27 K. Nakamura, R. Shimabukuro, Y. Fujiwara, T. Akiyama, T. Ito, and A. J. Freeman, Phys. Rev. Lett. 102 , 187201 (2009).
- 28 D.L. Mills, Phys. Rev. B 3 , 3887 (1971).
- 29 A. S. Logginov, G. A. Meshkov, A. V. Nikolaev, and A. P. Pyatakov JETP ett. 86 , 115 (2007).
- 30 A. Pyatakov, A. S. Sergeev, E. P. Nikolaeva, T. B. Kosykh, A. V. Nikolaev, K. A. Zvezdin, and A. K. Zvezdin, Phys.-Usp. 58 , 981 (2015).
- 31 V. Risingg˚rd, I. Kulagina, and J. Linder, Sci. Rep. a 6 , 31800 (2016).
- 32 A. P. Pyatakov and A. K. Zvezdin, Phys.-Usp. 55 , 557 (2012).
- 33 A. Michels, R. N. Viswanath, J. G. Barker, R. Birringer, and J. Weissmueller, Phys. Rev. Lett. 91 , 267204 (2003).
- 34 The wide range of more than 10 orders of magnitude of the possible strength of the polarization reflects the wide range of reported values of coupling constants. We focus on the necessary value of coupling constant that would provide measurable polarization.
- 35 K. Xia, W. Zhang, M. Lu, and H. Zhai, Phys. Rev. B 55 , 12561 (1997).
- 36 A. Crepieux and C. Lacroix, J. Magn. Magn. Mat. 182 , 341 (1998).
- 37 A. Zheludev, S. Maslov, G. Shirane, Y. Sasago, N. Koide, and K. Uchinokura, Phys. Rev. Lett. 78 , 4857 (1997); B. Lebech, J. Bernhard, and T. Freltoft, J. Phys. Condens. Matter 1 , 6105 (1989).
- 38 A. N. Bogdanov and U. K. Rossler, Phys. Rev. Lett. 87 , 037203 (2001).
- 39 A. O. Leonov and M. Mostovoy, Nat. Commun. 6 , 8275 (2015).
- 40 Due to the physical assumption that the dominant order parameter is the magnetic one, the stiffness energy part of the polarization can be neglected.
- 41 A. N. Bogdanov, U. K. Rossler, and K.-H. Mueller, J. Magn. Magn. Mat. 238 , 155 (2002).
- 42 A. N. Bogdanov and D. A. Yablonsky, Sov. Phys. JETP 68 , 101 (1989)
- 43 A. N. Bogdanov, JETP Lett. 62 , 247 (1995).
- 44 A. Takeuchi, S. Mizushima, and M. Mochizuki, Sci. Rep. 9 , 9528 (2019).
- 45 O. Boulle, J. Vogel, H. Yang, S. Pizzini, D. de Souza Chaves, A. Locatelli, T. Onur Mentes, A. Sala, L. D. Buda-Prejbeanu, O. Klein, M. Belmeguenai, Y. Roussigne, A. Stashkevich, S. Mourad Cherif, L. Aballe, M. Foerster, M. Chshiev, S. Auffret, I. Mihai Miron, and G. Gaudin, Nat. Nanotech. 11 , 449 (2016).
- 46 S. Valencia, A. Crassous, L. Bocher, V. Garcia, X. Moya, R. O. Cherifi, C. Deranlot, K. Bouzehouane, S. Fusil, A. Zobelli, A. Gloter, N. D. Mathur, A. Gaupp, R. Abrudan, F. Radu, A. Barthelemy, and M. Bibes,
- Nat. Mater. 10 , 753 (2011).
- 47 N. Lee, C. Vecchini, Y. J. Choi, L. C. Chapon, A. Bombardi, P. G. Radaelli, and S-W. Cheong, Phys. Rev. Lett. 110 , 137203 (2013).
- 48 S. Artyukhin, K. T. Delaney, N. A. Spaldin, and M. Mostovoy, Nat. Mat. 13 , 42 (2014).
- 49 C. A. F. Vaz, J. Hoffman, Y. Segal, J. W. Reiner, R. D. Grober, Z. Zhang, C. H. Ahn, and F. J. Walker, Phys. Rev. Lett. 104 , 127202 (2010).
- 50 P. Jain, Q. Wang, M. Roldan, A. Glavic, V. Lauter, C. Urban, Z. Bi, T. Ahmed, J. Zhu, M. Varela, Q. X. Jia, and M. R. Fitzsimmons, Sci. Rep. 5 , 9089 (2015).
- 51 J.-M. Hu, L.-Q. Chen, and C.-W. Nan, Adv. Mater. 28 , 15 (2016).
- 52 S. Ding, A. Ross, R. Lebrun, S. Becker, K. Lee, I. Boventer, S. Das, Y. Kurokawa, S. Gupta, J. Yang, G. Jakob, and M. Klaui, Phys. Rev. B 100 , 100406 (2015).
- 53 A. Hrabec, N. A. Porter, A. Wells, M. J. Benitez, G. Burnell, S. McVitie, D. McGrouther, T. A. Moore, and C. H. Marrows, Phys. Rev. B 90 , 020402(R) (2014).
- 54 M. Kuepferling, A. Casiraghi, G. Soares, G. Durin, F. Garcia-Sanchez, L. Chen, C. H. Back, C. H. Marrows, S. Tacchi, G. Carlotti, arXiv:2009.11830 (2020).
- 55 G. Chen, T. Ma, A. T. N'Diaye, H. Kwon, C. Won, Y. Wu, and A. K. Schmid, Nat. Commun. 4 , 2671 (2013).
- 56 C. Moreau-Luchaire, C. Moutafis, N. Reyren, J. Sampaio, C. A. F. Vaz, N. Van Horne, K. Bouzehouane, K. Garcia, C. Deranlot, P. Warnicke, P. Wohlhuter, J.-M. George, M. Weigand, J. Raabe, V. Cros, and A. Fert, Nat. Nanotech. 11 , 444 (2016).
- 57 S. Emori, U. Bauer, S.-M. Ahn, E. Martinez, and G. S. D. Beach, Nat. Mater. 12 , 611 (2013).
- 58 S. Emori, E. Martinez, K.-J. Lee, H.-W. Lee, U. Bauer, S.-M. Ahn, P. Agrawal, D. C. Bono, and G. S. D. Beach, Phys. Rev. B 90 , 184427 (2014).
- 59 W. Legrand, J.-Y. Chauleau, D. Maccariello, N. Reyren, S. Collin, K. Bouzehouane, N. Jaouen, V. Cros, and A. Fert, Sci. Adv. 4 : eaat0415 (2018).
- 60 D. Maccariello, W. Legrand, N. Reyren, K. Garcia, K. Bouzehouane, S. Collin, V. Cros, and A. Fert, Nat. Nanotech. 13 , 233 (2018).
- 61 A. W. J. Wells, P. M. Shepley, C. H. Marrows, and T. A. Moore, Phys. Rev. B 95 , 054428 (2017).
- 62 H. Stoll, M. Noske, M. Weigand, K. Richter, B. Krueger, R. M. Reeve, M. Hanze, C. F. Adolff, F.U. Stein, G. Meier, M. Klaui, and G. Schutz, Front. Phys. 3 , 26 (2015).
- 63 H. Chu, C. Jae Roh, J. O. Island, C. Li, S. Lee, J. Chen, J.-G. Park, A. F. Young, J. Seok Lee, and D. Hsieh, Phys. Rev. Lett. 124 , 027601 (2020).
Ahmed R. Tarkhany Marco Discacciati Joseph J. Betouras
Theory of surface-induced multiferroicity in magnetic materials, thin films and multilayers
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.