Magnetoelectricity in two-dimensional materials

ABSTRACT

Since the initial isolation of few-layer graphene, a plethora of two-dimensional atomic crystals has become available, covering almost all known materials types including metals, semiconductors, superconductors, ferro- and antiferromagnets. These advances have augmented the already existing variety of two-dimensional materials that are routinely realized by quantum confinement in bulk-semiconductor heterostructures. This review focuses on the type of material for which two-dimensional realizations are still being actively sought: magnetoelectrics. We present an overview of current theoretical expectation and experimental progress towards fabricating low-dimensional versions of such materials that can be magnetized by electric charges and polarized electrically by an applied magnetic field – unusual electromagnetic properties that could be the basis for various useful applications. The interplay between spatial confinement and magnetoelectricity is illustrated using the paradigmatic example of magnetic-monopole fields generated by electric charges in or near magnetoelectric media. For the purpose of this discussion, the image-charge method familiar from electrostatics is extended to solve the boundary-value problem for a magnetoelectric medium in the finite-width slab geometry using image dyons, i.e. point objects having both electric and magnetic charges. We discuss salient features of the magnetoelectrically induced fields arising in the thin-width limit.

Graphical abstract

KEYWORDS:

• Magnetoelectric effect
• 2D materials
• magnetic monopoles
• method of images

1. Introduction

Magnetoelectric media [1–3] are characterized by unconventional equilibrium responses to an electric field $\mathcal{E}$ and a magnetic field $\mathcal{B}$. While $\mathcal{E}$ induces only an electric polarization $\mathcal{P}$ in ordinary materials, it also creates a magnetization $\mathcal{M}$ in magnetoelectrics. Similarly atypical is an electric polarization arising due to $\mathcal{B}$ in addition to a magnetization. The full range of linear electromagnetic responses occurring in a magnetoelectric material is embodied in the constitutive relations expressing the electric displacement field $\mathcal{D}\equiv {\epsilon }_0\mathcal{E}+\mathcal{P}$ and the field $\mathcal{H}\equiv {\mu }_0^{-1}\mathcal{B}-\mathcal{M}$ in terms of $\mathcal{E}$ and $\mathcal{B}$, which read

(1) $\mathcal{D}=\underset{\_}{\epsilon }ℰ+\underset{\_}{\alpha }ℬ,\mathscr{H}={\underset{\_}{\mu }}^{-1}ℬ-{\underset{\_}{\alpha }}^Tℰ$(1)

Here $\underset{\_}{\epsilon }$ and $\underset{\_}{\mu }$ are the familiar [4] materials tensors of electric permittivity and magnetic permeability, respectively. The magnetoelectric tensor $\underset{\_}{\alpha }$ is associated with the magnetic-field-induced electric polarization, and its transpose ${\underset{\_}{\alpha }}^{\mathrm{T}}$ governs the electric-field-induced magnetization [1–3]. The interested reader can find a more in-depth discussion of constitutive relations for magnetoelectrics in Appendix A, including a juxtaposition of two commonly adopted conventions.

Seminal works on the magnetoelectric effect were largely focused on the symmetry classification of materials with nonvanishing $\underset{\_}{\alpha }$ [5,6], but ramifications of magnetoelectric responses have also been discussed early on [7]. However, it was the advent of multiferroic materials [8,9] with their large magnetoelectric couplings [10–13] that has greatly boosted current interest in gaining a deeper understanding about magnetoelectricity [14], also with the view to enable useful applications [15–19]. The intriguing connection with the particle-physics concept of axion electrodynamics [20–22] has recently become even more relevant through the realization of topological insulators [23–27] and Weyl semimetals [26–29]. The opportunity to explore sizable magnetoelectric effects in real materials has spurred efforts to describe unconventional electromagnetic phenomena exhibited by magnetoelectrics [30–40]. Among these, simulations of magnetic-monopole fields have been paradigmatic [30–32,37,39].

Nanostructuring generally opens up new avenues for controlling materials properties, including those of multiferroics and magnetoelectrics [41,42]. Recent advances in the design of two-dimensional (2D) atomic crystals [43] and their versatile combination [44–46] have transcended the original platform of low-dimensional systems created in semiconductor heterostructures [47,48]. With magnetic 2D materials already realized in metal [49] and semiconductor quantum wells [50] as well as atomic crystals [51–54], the focus has shifted to the possibility of 2D multiferroics and magnetoelectrics [52–58]. In this context, the purpose of this Review is two-fold. In Sec. 2, we present a survey of 2D materials where magnetoelectricity is expected and/or has already been observed. In the subsequent Sec. 3, we discuss the ramification of a 2D material’s finite size on its magnetoelectric responses, focusing especially on magnetic-monopole fields generated by electric charges placed near or in the medium. Our work is intended to raise awareness of the special features associated with magnetoelectricity in low-dimensional materials, which have been discussed only sporadically until now [59].

2. Survey of 2D magnetoelectric materials

The systematic discussion of magnetoelectric materials properties is usually based on consideration of a medium’s free-energy density $F(\mathcal{E},\mathcal{B})$. Its expansion up to quadratic order in electric and magnetic field components is most generally given by [2,3]

(2) $F(\mathcal{E},\mathcal{B})=F(\mathbf{0},\mathbf{0})-{\mathcal{P}}_{\mathrm{s}}\cdot \mathcal{E}-{\mathcal{M}}_{\mathrm{s}}\cdot \mathcal{B}-\frac{1}{2}\mathcal{E}\cdot ({\underset{\_}{\chi }}_{\mathcal{E}}\mathcal{E})-\frac{1}{2}\mathcal{B}\cdot ({\underset{\_}{\chi }}_{\mathcal{B}}\mathcal{B})\text{break}-\mathcal{E}\cdot (\underset{\_}{\alpha }\mathcal{B})-\dots .$(2)

Here ${\mathcal{P}}_{\mathrm{s}}$ and ${\mathcal{M}}_{\mathrm{s}}$ are the material’s spontaneous electric polarization and magnetization, respectively, while ${\underset{\_}{\chi }}_{\mathcal{E}}$ and ${\underset{\_}{\chi }}_{\mathcal{B}}$ denote the conventional [4] electric and magnetic susceptibility tensors. The linear magnetoelectric effect is embodied in the term containing $\underset{\_}{\alpha }$. Contributions to $F(\mathcal{E},\mathcal{B})$ giving rise to nonlinear electromagnetic (including higher-order magnetoelectric [3,60,61]) responses have been omitted from Eq. (2(2) $F(\mathcal{E},\mathcal{B})=F(\mathbf{0},\mathbf{0})-{\mathcal{P}}_{\mathrm{s}}\cdot \mathcal{E}-{\mathcal{M}}_{\mathrm{s}}\cdot \mathcal{B}-\frac{1}{2}\mathcal{E}\cdot ({\underset{\_}{\chi }}_{\mathcal{E}}\mathcal{E})-\frac{1}{2}\mathcal{B}\cdot ({\underset{\_}{\chi }}_{\mathcal{B}}\mathcal{B})\text{break}-\mathcal{E}\cdot (\underset{\_}{\alpha }\mathcal{B})-\dots .$(2) ), as we are not focusing on these in the following. Similarly, possible interactions between simultaneously present spontaneous electric and magnetic orders, such as ${\mathcal{M}}_{\mathrm{s}}$ being affected when switching ${\mathcal{P}}_{\mathrm{s}}$ in multiferroics, are not being considered in the following unless these present a mechanism for generating a novanishing $\underset{\_}{\alpha }$.

General descriptions of materials properties in terms of tensor quantities can be obtained based on fundamental symmetry considerations [62]. For example, magnetoelectricity can occur only in systems where both spatial-inversion and time-reversal symmetries are broken. This allows us to distinguish the magnetoelectric responses we are interested in from the superficially related phenomena of current-induced magnetization [63] or spin-orbit torque [64], which are based on a relation $\mathcal{M}=\underset{\_}{\eta }\mathcal{J}$ between the magnetization and an electric-current density $\mathcal{J}$ [65–68]. If the current originates from an applied electric field and is related to the latter via Ohm’s law, $\mathcal{J}=\underset{\_}{\sigma }\mathcal{E}$, the resulting dependence of the magnetization on the electric field,

(3) $\mathcal{M}=\underset{\_}{\eta }\mathcal{J}=\underset{\_}{\eta }(\underset{\_}{\sigma }\mathcal{E})\equiv {\underset{\_}{\eta }}^{\mathrm{\prime }}\mathcal{E},$(3)

looks similar to the magnetoelectric effect [69]; see Eq. (A1b). However, the basic conditions for the tensor $\underset{\_}{\eta }$ to exist in a material differ from those required for nonvanishing $\underset{\_}{\alpha }$. Most importantly, current-induced magnetization is possible in materials with unbroken time-reversal symmetry. Furthermore, in contrast to the magnetoelectric effect which is an equilibrium phenomenon, nonequilibrium mechanisms are crucial for generating a magnetization accompanying the electric current [63]. Recent theoretical studies have explored current-induced magnetization in 2D materials, e.g. in twisted bilayer graphene [70] and the interfacial 2D electron gas in SrTiO${ }_3$ [71]. Spin-orbit torques in magnetic 2D materials, such as Fe${ }_3$GeTe${ }_2$ [72] and CrI${ }_3$ [73] have also been considered. An intriguing realization of current-induced magnetization in the 2D transition-metal dichalcogenide MoS${ }_2$ [74] utilized the material’s valley-isospin degree of freedom [75].

We focus here only on magnetoelectric 2D materials, i.e. those having finite $\underset{\_}{\alpha }$. To impose some order on the rapidly expanding zoo of 2D magnetoelectrics, materials are grouped into subsections according to the basic mechanisms underlying their magnetoelectric responses. We start by discussing 2D versions of single-phase magnetoelectrics and multiferroics in Sec. 2.1. Following that, Sec. 2.2 provides an overview of 2D multiferroic heterostructures. Sec. 2.3 is dedicated to magnetic 2D materials whose (ferro- or antiferro-)magnetic order is tunable by electric fields. The unconventional magnetoelectricity of 2D charge carriers in semiconductor quantum wells is considered in Sec. 2.4. We conclude the survey of 2D magnetoelectric materials with Sec. 2.5 exploring the ways in which the valley-degree of freedom enables magnetoelectric responses in 2D atomic crystals, such as few-layer sheets of graphene or transition-metal dichalcogenides.

2.1. Single-phase 2D magnetoelectrics

Before the advent of multiferroics, the typical magnetoelectric bulk material was expected to be an inversion-asymmetric time-reversal-breaking antiferromagnet, much like the first and paradigmatic example Cr${ }_2$O${ }_3$ [5,76]. Although time-reversal symmetry is always broken in ferromagnets but not in all antiferromagnets [77], inversion symmetry is rarely broken in ferromagnets unless the underlying crystal structure is already noncentrosymmetric. In contrast, antiferromagnetic order itself can break inversion symmetry, thus greatly multiplying the possibilities for realizing a magnetoelectric material. By now, many antiferromagnetic compounds are known to be magnetoelectric, but the magnitude of their magnetoelectric couplings ${\alpha }_{ij}$ are usually small [3] – the current record holder is TbPO${ }_4$ with $\left|{\alpha }_{ij}\right|\text{~}0.220\sqrt{{\epsilon }_0/{\mu }_0}$ [78,79].

Despite the fact that the smallness of magnetoelectric couplings is generic for single-phase antiferromagnets, those among them that are van-der-Waals crystals can serve as an attractive starting point for fabricating 2D magnetoelectrics. Pertinent examples are MnPS${ }_3$ [80] and NiI${ }_2$ [81], whose intralayer antiferromagnetic order has been shown to persist in the few-layer limit [82–84]. To date, there has been no direct confirmation of the magnetoelectric effect in these 2D allotropes of bulk magnetoelectrics, but the presence of all necessary ingredients is promising and should motivate further exploration of these materials [85] and similar candidates [86,87].

An interesting alternative to planar antiferromagnets is presented by layered magnets where the magnetization direction alternates between adjacent layers. Bilayers of such materials could constitute single-phase antiferromagnets and show the magnetoelectric effect. A prominent example is CrI${ }_3$ for which the magnetoelectric coupling has been measured to be $\left|{\alpha }_{zz}\right|\text{~}0.034\sqrt{{\epsilon }_0/{\mu }_0}$ [88]. See also related experiments [89,90] and theoretical work [91,92]. It can be expected that similar materials such as WSe${ }_2$ and CrTe${ }_2$ also exhbibit magnetoelectricity [92].

Single-phase materials with large magnetoelectric effects are typically multiferroic, i.e. exhibit coexisting spontaneous electric polarizations and magnetizations. There is, however, a paucity of such materials even in bulk, because the most effective mechanisms for generating ferromagnetisms and ferroelectricity are mutually exclusive [93] due to a classic case of chemical contra-indication [94]: ferroelectrics favor having empty d or f orbitals, whereas magnetism typically arises when these orbitals are partially filled. In fact, no 2D materials that are single-phase ferromagnetic-ferrolectric multiferroics have been realized so far [58], but theoretical studies are pointing to a number of promising candidates: CrN and CrB${ }_2$ [95], MXene bilayers [96], monolayer (CrBr${ }_3$)${ }_2$Li [97], TMPCs-CuMP${ }_2$X${ }_6$ [98,99], monolayer Hf${ }_2$VC${ }_2$F${ }_2$ [100], monolayer VOCl${ }_2$ [101] and VOI${ }_2$ [102] (but see also Ref [103]), halogen-intercalated phosphorene bilayers [104], transition-metal-intercalated MoS${ }_2$ bilayers [105], double-perovskite bilayers [106], as well as monolayer $\alpha$-In${ }_2$Se${ }_3$ [107].

2.2. 2D multiferroic heterostructures

Hybrid systems incorporating ferroelectric and ferromagnetic parts provide a promising avenue for overcoming the intrinsic limitations that make it so hard to have both a spontaneous electric polarization and a magnetization in the same material [8]. In fact, gigantic magnetoelectric couplings of $\left|{\alpha }_{ij}\right|\text{~}4.8\times {10}^3\sqrt{{\epsilon }_0/{\mu }_0}$ have recently been realized in the heterostructure multiferroic FeRh/BTO [108]. The expectation that a similar strategy of combining 2D ferroelectrics with 2D ferromagnets will yield robust 2D multiferroics has fuelled intense theoretical efforts. We discuss below a number of promising materials combinations for which magnetoelectic phenomena in the wider sense have been predicted. However, the available theoretical results do not permit us to infer with certainty that the linear magnetoelectric effect occurs, or estimate the magnitude of magnetoelectric-tensor components $|{\alpha }_{ij}|$.

The first-considered van-der-Waals heterostructure of ferromagnetic Cr${ }_2$Ge${ }_2$Te${ }_6$ combined with ferroelectric In${ }_2$Se${ }_3$ showed intriguing properties [109]. Besides demonstrating the ability to switch the magnetization of the Cr${ }_2$Ge${ }_2$Te${ }_6$ layer via reversal of the electric polarization in the In${ }_2$Se${ }_3$ part, the latter became magnetized via the proximity effect – in effect becoming a single-phase 2D multiferroic.

The combination of magnetic bilayer CrI${ }_3$ with ferroelectric monolayer Sc${ }_2$CO${ }_2$ also yielded interesting magnetoelectric responses [110]. Here the reversal of electric polarization in Sc${ }_2$CO${ }_2$ triggered the transition between antiferromagnetic and ferromagnetic order in bilayer CrI${ }_3$. The same mechanism for magnetoelectricity has also been studied theoretically for a bilayer-CrI${ }_3$ /monolayer-In${ }_2$Se${ }_3$ multiferroic heterostructure [111]. Even greater efficiency for electric control of the antiferromagnetic-to-ferromagnetic transition in bilayer CrI${ }_3$ is expected when it forms a van-der-Waals heterostructure with perovskite-oxide ferroelectrics such as BiFeO${ }_3$ [112].

Control of the magnetization in Fe${ }_n$GeTe${ }_2$ layers via the electric polarization of In${ }_2$Se${ }_3$ was predicted to enable enhanced magnetotransport functionality in a tunnel-junction device [113]. Transition-metal-decorated graphene can also be used as the magnetic part of a multiferroic heterostructure formed with monolayer In${ }_2$Se${ }_3$ [114].

2.3. 2D materials with electrically tunable spontaneous magnetization

Magnetic order generally arises from the exchange interaction between microscopic magnetic moments. In systems where the exchange-interaction strength depends on tunable electronic degrees of freedom, e.g. the charge density in a semiconductor [50], the spontaneous magnetization ${\mathcal{M}}_{\mathrm{s}}$ can become a function of the applied electric field $\mathcal{E}$. Electric control of magnetism based on such mechanisms is attracting a lot of interest [115] but has typically not been discussed using the language and formalism of the magnetoelectric effect. Nevertheless, in the spirit of the expansion (2(2) $F(\mathcal{E},\mathcal{B})=F(\mathbf{0},\mathbf{0})-{\mathcal{P}}_{\mathrm{s}}\cdot \mathcal{E}-{\mathcal{M}}_{\mathrm{s}}\cdot \mathcal{B}-\frac{1}{2}\mathcal{E}\cdot ({\underset{\_}{\chi }}_{\mathcal{E}}\mathcal{E})-\frac{1}{2}\mathcal{B}\cdot ({\underset{\_}{\chi }}_{\mathcal{B}}\mathcal{B})\text{break}-\mathcal{E}\cdot (\underset{\_}{\alpha }\mathcal{B})-\dots .$(2) ) for the field-dependent part of the system’s free-energy density, a term

$-{\mathcal{M}}_{\mathrm{s}}\left(\mathcal{E}\right)\cdot \mathcal{B}\equiv -{\mathcal{M}}_{\mathrm{s}}\left(0\right)\cdot \mathcal{B}-\left(\frac{\mathrm{\partial }{\mathcal{M}}_{\mathrm{s}}}{\mathrm{\partial }\mathcal{E}}{|}_{\mathcal{E}\to 0}\mathcal{E}\right)\cdot \mathcal{B}-\dots$

actually incorporates magnetoelectricity with a tensor $\underset{\_}{\alpha }$ whose elements ${\alpha }_{ij}=[\mathrm{\partial }{\mathcal{M}}_{\mathrm{s}j}\left(\mathcal{E}\right)/\mathrm{\partial }{\mathcal{E}}_i{]}_{\mathcal{E}\to 0}$ quantify the electric-field tunability of the spontaneous magnetization. Typical magnitudes $\left|{\alpha }_{ij}\right|\text{~}4.0\times {10}^{-3}\sqrt{{\epsilon }_0/{\mu }_0}$ have been realized in the dilute magnetic semiconductor GaMnAs [116].

Magnetoelectric effects arising from an electric-field-dependent spontaneous magnetization have been discussed theoretically [117,118] and observed experimentally [119,120] in 2D quantum-well structures made from dilute-magnetic semiconductors. Few-layer Cr${ }_2$Ge${ }_2$Te${ }_6$ has recently been shown to exhibit electric-field-tunable ferromagnetism [121–123]. Using experimental data and results of calculations presented in Ref [122], we estimate $|{\alpha }_{ij}|\sim 2\times {10}^{-6}\sqrt{{\epsilon }_{{\hspace{0.17em}}_0}/{\mu }_{{\hspace{0.17em}}_0}}$ for this 2D magnetoelectric material. It would be interesting to explore whether other currently available electrically tunable 2D magnets, e.g. monolayer Fe${ }_2$GeTe${ }_2$ [124], also exhibit the linear magnetoelectric effect.

2.4. Magnetoelectricity of 2D charge carriers in quantum wells

In bulk conductors, screening and nonequilibrium responses associated with the application of an electric field or the generation of an electric polarization generally preclude any discussion of the magnetoelectric effect.¹ Low-dimensional conductors, however, can exhibit equilibrium responses as long as $\mathcal{E}$ and $\mathcal{P}$ are parallel to spatial directions within which charge-carrier motion is quantized [125–127]. Furthermore, 2D-conductor samples fitted with front and back gates allow changing of $\mathcal{E}$ at constant charge density [128–131], thus making it possible to disentangle electric-polarization responses from charging effects in experiments. Conducting 2D materials are therefore ideally suited for a detailed exploration of magnetoelectricity associated with itinerant charge carriers. In this subsection, we focus on the paradigmatic example of 2D conductors realized in semiconductor-heterostructure quantum wells [47]. A pedagogical introduction to the basic principles of band-gap engineering and physical properties of 2D semiconductor structures is available, e.g. in Ref [48] (especially Chapter 4).

Figure 1. Mechanisms underpinning magnetoelectricity in semiconductor-heterostructure quantum wells. Charge carriers can move freely in the $xy$ plane but are confined in $z$ direction by a symmetric potential arising from the spatially varying conduction-band bottom (drawn in black pen). Equilibrium currents [(a)] and bound-state charge distributions [(b)] are indicated in red. The band bending shown schematically in the right panel of subfigure (a) is due to the applied perpendicular electric field. The qualitatively similar, but typically much smaller in magnitude, band bending arising as a consequence of the magnetic-field-induced electric polarization is not shown in the right panel of subfigure (b).

The microscopic origin of magnetoelectricity in quantum wells are quadrupolar equilibrium currents whose distortion by an applied electric field generates a magnetization [132]. Such equilibrium-current distributions arise from the interplay of spatial-inversion asymmetry with an applied magnetic field [132,133], or with an intrinsic exchange field due to ferromagnetic order [132,133], or with the staggered exchange field of a nontrivial antiferromagnetic order [132].² Furthermore, an applied magnetic field shifts the 2D-bound-state charge distribution, thus inducing an electric polarization perpendicular to the well [132,133], as mandated by the duality of magnetoelectric reponses [1,2]. Figure 1 illustrates schematically the microscopic mechanisms for magnetoelectric responses in semiconductor-heterostructure quantum wells. According to theoretical estimates [132], the magnitude of the magnetoelectric-tensor components in quantum wells can reach values $\left|{\alpha }_{ij}\right|\text{~}{10}^{-4}\sqrt{{\epsilon }_0/{\mu }_0}$, and the $\mathcal{E}$-induced magnetic moment per 2D charge carrier can be as large as $0.6$ Bohr magnetons.

2.5. Magnetoelectricity in 2D materials from topology or valley isospin

Exotic types of magnetoelectricity have been suggested to exist in certain materials, even those in which time-reversal and spatial-inversion symmetries are not broken. One example for these are topological insulators [23–27], which have an isotropic magnetoelectric coupling $\underset{\_}{\alpha }=\alpha \underset{-}{1}$ with a quantized magnitude $\left|\alpha \right|\equiv e^2/\left(4\pi \mathrm{\hslash }\right)\approx \left(1/137\right)\sqrt{{\epsilon }_0/{\mu }_0}$. Among the multitude of 2D materials, a topological magnetoelectric effect was predicted [134–136], and its signatures have recently been observed [137,138], in even-layer MnBi${ }_2$Te${ }_4$.³ While interesting from a fundamental point of view, the inability to adjust the magnitude of the topological magnetoelectric coupling limits its relevance for applications – unless proposals for engineering tunable magnetoelectric couplings via symmetry-breaking finite-size effects in thin MnBi${ }_2$Te${ }_4$ layers [139] come to fruition.

Materials with a valley-isospin degree of freedom [75] lend themselves to an intriguing realization of electrically tunable magnetoelectricity. Concrete proposals for 2D materials have been formulated for bilayer graphene [140,141] and transition-metal-dichalcogenide bilayers [142]. In all these cases, electronic states from different valleys in the 2D-material band structure (usually labelled $\mathbf{K}$ and $\mathbf{K}{ }^{\mathrm{\prime }}$) are connected via time reversal or spatial inversion. An imbalance between charge-carrier densities in the two valleys – a finite valley-isospin density – amounts to a breaking of both symmetries [127] and renders the material magnetoelectric. As various ways for addressing valley isospin have been explored as part of a drive to realize electronic devices based on a new valleytronics paradigm [143], these systems offer the unique opportunity for generating electrically tailored magnetoelectric couplings.

3. Electromagnetism of 2D magnetoelectric media

The survey of 2D magnetoelectric media presented in the previous section outlines the great progress being made in designing, fabricating and characterizing such materials. Once the remaining materials-science challenges are overcome, the focus of further investigations will shift to exploring more broadly how such materials behave electromagnetically. In anticipation of that, this section focuses on a paradigmatic property of magnetoelectric media: the simulation of magnetic-monopole fields discussed extensively for bulk systems [30–32,37]. Below we briefly describe how isotropic magnetoelectric responses can be included as unconventional source terms into the Maxwell equations of ordinary electrodynamics, thus mimicking the effects of axion electrodynamics [20–23]. This forms the basis for solving boundary-value problems involving isotropic magnetoelectrics and enables a generalization of the well-known method of images [4] to include image magnetic monopoles alongside image charges [30]. In the following, we refer to this approach as the method of image dyons and use it to obtain electric and magnetic fields generated by an electric charge placed near or in a finite-width magnetoelectric slab.

A configuration of external sources (i.e. electric charges and currents) determines the electromagnetic fields via the inhomogeneous Maxwell equations [4]

(4) $\mathrm{\nabla }\cdot \mathcal{D}=\rho ,\mathrm{\nabla }\times \mathcal{H}-\frac{\mathrm{\partial }\mathcal{D}}{\mathrm{\partial }t}=\mathcal{J}.$(4)

Here $\rho (\mathbf{r})$ and $\mathcal{J}(\mathbf{r})$ denote the charge and current densities, respectively, as a function of position $\mathbf{r}\equiv (x,y,z)$, and $\mathrm{\nabla }\equiv (\mathrm{\partial }/\mathrm{\partial }x,\mathrm{\partial }/\mathrm{\partial }y,\mathrm{\partial }/\mathrm{\partial }z)$ is the gradient operator. Solving (4(4) $\mathrm{\nabla }\cdot \mathcal{D}=\rho ,\mathrm{\nabla }\times \mathcal{H}-\frac{\mathrm{\partial }\mathcal{D}}{\mathrm{\partial }t}=\mathcal{J}.$(4) ), in conjunction with the set of homogeneous Maxwell equations

(5) $\mathrm{\nabla }\cdot \mathcal{B}=0,\mathrm{\nabla }\times \mathcal{E}+\frac{\mathrm{\partial }\mathcal{B}}{\mathrm{\partial }t}=0$(5)

and the constitutive relations (1(1) $\mathcal{D}=\underset{\_}{\epsilon }ℰ+\underset{\_}{\alpha }ℬ,\mathscr{H}={\underset{\_}{\mu }}^{-1}ℬ-{\underset{\_}{\alpha }}^Tℰ$(1) ) that embody the effect of microscopic degrees of freedom in media, yields the fields throughout space.

As a starting point for exploring the interplay between low dimensionality and magnetoelectricity, we focus in the following on the case of media with isotropic magnetoelectric response whose magnetoelectric tensor has the form $\underset{\_}{\alpha }=\alpha \underset{\_}{1}$. For consistency, isotropy is also assumed for the ordinary polarization and magnetization responses, i.e. $\underset{\_}{\epsilon }=\epsilon \underset{\_}{1}$ and $\underset{\_}{\mu }=\mu \underset{\_}{1}$. To make modifications arising from isotropic magnetoelectricity more explicit, the field-source relations (4(4) $\mathrm{\nabla }\cdot \mathcal{D}=\rho ,\mathrm{\nabla }\times \mathcal{H}-\frac{\mathrm{\partial }\mathcal{D}}{\mathrm{\partial }t}=\mathcal{J}.$(4) ) can be rewritten as

(6) $\mathrm{\nabla }\cdot \left(\epsilon \mathcal{E}\right)=\rho -\left(\mathrm{\nabla }\alpha \right)\cdot \mathcal{B},\mathrm{\nabla }\times \left({\mu }^{-1}\mathcal{B}\right)-\frac{\mathrm{\partial }\left(\epsilon \mathcal{E}\right)}{\mathrm{\partial }t}=J+\left(\mathrm{\nabla }\alpha \right)\times \mathcal{E}+\frac{\mathrm{\partial }\alpha }{\mathrm{\partial }t}\mathcal{B}.$(6)

The particular form of the inhomogeneous Maxwell equations given in (6(6) $\mathrm{\nabla }\cdot \left(\epsilon \mathcal{E}\right)=\rho -\left(\mathrm{\nabla }\alpha \right)\cdot \mathcal{B},\mathrm{\nabla }\times \left({\mu }^{-1}\mathcal{B}\right)-\frac{\mathrm{\partial }\left(\epsilon \mathcal{E}\right)}{\mathrm{\partial }t}=J+\left(\mathrm{\nabla }\alpha \right)\times \mathcal{E}+\frac{\mathrm{\partial }\alpha }{\mathrm{\partial }t}\mathcal{B}.$(6) ) coincides with the one derived in the context of axion electrodynamics [20,21] whose realization in condensed-matter systems is currently attracting great interest [26,27]. It shows that, in situations with piecewise-constant $\alpha$, e.g. when certain parts of space are filled with a homogeneous isotropic magnetoelectric medium, the effects of magnetoelectricity can be represented in terms of charges and currents induced at interfaces.

Various theoretical techniques have been used to solve the set of equations describing electromagnetisms with magnetoelectric media [30–35]. Below we discuss in greater detail the method of introducing image charges and monopoles [30] to satisfy boundary conditions of electromagnetic fields at an interface between vacuum and an isotropic magnetoelectric medium. Generalizing this approach to the case with more than a single interface enables us to investigate how the finite width of a 2D magnetoelectric medium affects its unconventional electromagnetism.

3.1. Image dyons for solving boundary-value problems in magnetoelectrics

The method of image charges is an established technique for solving boundary-value problems in electrostatics [4]. It is based on the possibility to use an appropriately placed fictitious (i.e. image) point charge to represent the polarization-charge distribution generated at an interface by the presence of a real electric point charge. This approach turns out to be most useful when considering high-symmetry configurations, e.g. when the real charges are placed in spaces with planar, cylindrical or spherical interfaces between different media. In particular, the image-charge method has been applied to obtain the electric fields generated by a charge near [144] or inside [145] a dielectric plate.

Figure 2. An electric charge $q_0$ placed near a planar interface (located at $z=0$) between an ordinary dielectric (occupying the region $z>0$) and a magnetoelectric medium (present in region $z<0$) generates image dyons on both sides. The values of their electric (magnetic) charges are related via $q_1^{+}=q_1^{-}$ ($g_1^{+}=-g_1^{-}$).

The method of image charges can be extended to describe situations involving magnetoelectrics by allowing magnetic monopoles [146] alongside each of the electric image charges [30], thus effectively rendering the image objects to be dyons [147,148]. See Figure 2 for an illustration of the previously considered case of an electric charge near a single planar boundary between vacuum and a magnetoelectric medium [30,32,149]. Below we juxtapose the known results for this case with those calculated for our situation of interest, which is a charge placed near or inside a finite-width slab of magnetoelectric material. Analogous to the situation of an electric point charge near or inside a dielectric plate where an infinite number of image charges need to be considered [144,145], an infinite number of image dyons are required for describing the finite-slab magnetoelectric. Figure 3 shows the configurational details pertaining to this case. Before presenting our results, we briefly state the boundary-value problem posed by planar interfaces with magnetoelectric media and sketch the formalism underpinning the image-dyon method used for its solution.

Figure 3. An infinite set of image dyons $d_k^{m\nu }$ is needed to describe electric and magnetic fields generated by an electric point charge $q_0$ placed near [panel (a)] or inside [panel (b)] a finite-width slab of magnetoelectric material. Here $m=\mathrm{U}(\mathrm{L})$ labels dyons arising from mirror-imaging at the upper (lower) boundary, $\nu =+(-)$ when the image dyon is located above (below) the $m$ boundary, and $k\in \mathbb{Z}\mathrm{\setminus }0$ is a counter associated with the dyon location. Dyons with $k<0$ arise only for the case with $q_0$ inside the magnetoelectric medium [panel (b)].

Considering the static limit and assuming any field sources entering the r.h.s. of the inhomogeneous Maxwell equations(4)(4) $\mathrm{\nabla }\cdot \mathcal{D}=\rho ,\mathrm{\nabla }\times \mathcal{H}-\frac{\mathrm{\partial }\mathcal{D}}{\mathrm{\partial }t}=\mathcal{J}.$(4) to be localized away from a planar interface that is parallel to the $xy$ plane and intersects the $z$ axis at $z=z_{\mathrm{b}}$, the electromagnetic fields obey the boundary conditions

(7a) ${\mathcal{D}}_z(R,\theta ,z_{\mathrm{b}}^{+})={\mathcal{D}}_z(R,\theta ,z_{\mathrm{b}}^{-}),{\mathcal{H}}_R(R,\theta ,z_{\mathrm{b}}^{+})={\mathcal{H}}_R(R,\theta ,z_{\mathrm{b}}^{-}),$(7a)

(7b) ${\mathcal{E}}_R(R,\theta ,z_{\mathrm{b}}^{+})={\mathcal{E}}_R(R,\theta ,z_{\mathrm{b}}^{-}),{\mathcal{B}}_z(R,\theta ,z_{\mathrm{b}}^{+})={\mathcal{B}}_z(R,\theta ,z_{\mathrm{b}}^{-}).$(7b)

Here cylindrical coordinates $\mathbf{r}\equiv (R,\theta ,z)$ are utilized for positions, and $z_{\mathrm{b}}^{\pm }={lim}_{\epsilon \to 0}(z_{\mathrm{b}}\pm \epsilon )$. Furthermore, the static limit enables both the electric field and the magnetic field to be expressed in terms of scalar potentials,

(8) $\mathcal{E}=-\mathrm{\nabla }V,\mathcal{B}=-\mathrm{\nabla }U.$(8)

The method of image dyons amounts to writing Ansätze for the scalar potentials arising in a given region of space, due to having a source charge $q_0$ at the location ${\mathbf{r}}_0$, in terms of a superposition of electric-point-charge and magnetic-monopole contributions,

(9) $V\left(\mathbf{r}\right)=\frac{1}{4\pi \epsilon { }^{\mathrm{\prime }}}\left[\frac{q_0}{\left|\mathbf{r}-{\mathbf{r}}_0\right|}+\sum _k\frac{q_k}{\left|\mathbf{r}-{\mathbf{r}}_k\right|}\right],U\left(\mathbf{r}\right)=\frac{\mu { }^{\mathrm{\prime }}}{4\pi }\sum _k\frac{g_k}{\left|\mathbf{r}-{\mathbf{r}}_k\right|}.$(9)

Here $q_k$ ($g_k$) denote the electric (magnetic) charge of the fictitious image dyon $d_k$ located at ${\mathbf{r}}_k$. It is understood that the sum ${\sum }_k$ is restricted to include only those image dyons that are relevant for calculating the fields in a given region, i.e. those that are located outside. Furthermore, ${\epsilon }^{\prime }$ (${\mu }^{\prime }$) is the permittivity (permeability) of the region where the source charge is located. To begin with, the electric and magnetic charges making up each image dyon $d_k\equiv (q_k,g_k)$ as well as its location ${\mathbf{r}}_k$ are unknown. However, the location ${\mathbf{r}}_k$ can be predicted based on the geometry of the system if we require it to be independent of the charges $q_k$ and $g_k$. It turns out that ${\mathbf{r}}_k$ is the location of first- or higher-order mirror images of the source charge with respect to the boundary (or boundaries). Since the slab geometry considered in this work has planes perpendicular to the $z$ axis as its boundaries, and since we assume that the electric source charge is located at ${\mathbf{r}}_0=(0,0,\pm z_0)$ with $z_0>0$, the locations of the image dyons are $r_k=(0,0,z_k)$ with the coordinate $z_k$ determined by $z_0$ and the slab width $w$. In particular, if the magnetoelectric medium is semi-infinite (i.e. there is only one boundary, assumed to be located at $z=0$ without loss of generality), we only have two image dyons $d_1^{+}=(q_1^{+},g_1^{+})$ and $d_1^{-}=(q_1^{-},g_1^{-})$ as shown in Figure 2, and Eq. (9)(9) $V\left(\mathbf{r}\right)=\frac{1}{4\pi \epsilon { }^{\mathrm{\prime }}}\left[\frac{q_0}{\left|\mathbf{r}-{\mathbf{r}}_0\right|}+\sum _k\frac{q_k}{\left|\mathbf{r}-{\mathbf{r}}_k\right|}\right],U\left(\mathbf{r}\right)=\frac{\mu { }^{\mathrm{\prime }}}{4\pi }\sum _k\frac{g_k}{\left|\mathbf{r}-{\mathbf{r}}_k\right|}.$(9) simplifies. Furthermore, we find the relations $q_1^{+}=q_1^{-}$ and $g_1^{+}=-g_1^{-}$. Full details of the calculation and results for image-dyon charges are given in Appendix B.1.

If the magnetoelectric medium has the shape of a finite-width slab, an infinite number of image dyons exists as shown in Figure 3. We can see this by trying to satisfy the field-continuity conditions for each of the two boundaries in turn. Assuming the source charge $q_0$ to be outside the slab, we can start by satisfying the upper boundary conditions while ignoring the lower boundary – this yields the image dyons $d_1^{\mathrm{U}+}=(q_1^{\mathrm{U}+},g_1^{\mathrm{U}+})$ at $z=z_0$ and $d_1^{\mathrm{U}-}=(q_1^{\mathrm{U}-},g_1^{\mathrm{U}-})$ at $z=-z_0$ as in the single-interface situation depicted in Figure 2(a). [Here and in the following, the superscript U(L) is used to indicate that the image dyon results from mirroring at the upper (lower) boundary, which is at $z=0$ ($z=-w$), the $+$ ($-$) label next to it specifies that the charge is located above (below) that boundary, and the integer $k$ relates to the image-dyon location as explained below.] Attempting next to satisfy the boundary conditions at the bottom interface while ignoring the upper boundary, we treat the dyon $(q_0+q_1^{\mathrm{U}+},g_1^{\mathrm{U}+})$ at $z_0$ as the source and obtain additional image dyons $d_1^{\mathrm{L}+}$ at $z=z_0$ and $d_2^{\mathrm{L}-}$ at $z=-(z_0+2w)$. However, their introduction causes boundary conditions at the upper interface to no longer be satisfied. To remedy this, two more image dyons need to be created; $d_2^{\mathrm{U}+}$ at $z=z_0+2w$ and $d_2^{\mathrm{U}-}$ at $z=-(z_0+2w)$, whose properties are determined from assuming dyon $d_2^{\mathrm{L}-}$ as the source. This in turn causes boundary conditions at the lower interface to be violated, which needs to be fixed by further iteration of the image-dyon method. In the end, an infinite number of image dyons are needed to satisfy the boundary conditions at both interfaces asymptotically.

For the case when the source charge $q_0$ is inside the slab, even more image dyons are needed as shown in Figure 3(b). Similar as before, we can start by satisfying the upper boundary conditions while ignoring the lower boundary, and obtain image dyons $d_1^{\mathrm{U}+}=$ at $z=z_0$ and $d_1^{\mathrm{U}-}$ at $z=-z_0$. Next, we satisfy the lower boundary conditions while ignoring the upper boundary and obtain $d_{-1}^{\mathrm{L}+}$ at $z=z_0$, $d_{-2}^{\mathrm{L}-}$ at $z=z_0-2w$, $d_1^{\mathrm{L}+}$ at $z=z_0$, and $d_2^{\mathrm{L}-}$ at $z=-(z_0+2w)$ by treating both $q_0$ and the dyon $d_1^{\mathrm{U}+}$ as sources. The subset of these image dyons having $k>0$ are equivalent to the image dyons found for the case with the source charge outside the slab, and further iteration builds up the same structure as shown in Figure 3(a). In contrast, the subset of image dyons labelled by $k<0$ are entirely new, and their subsequent alternate mirroring creates infinitely many additional dyons.

Generally, the location of an image dyon $d_k^{m\nu }$ is at

(10) ${\mathbf{r}}_k^{\nu }=(0,0,\nu \mathrm{s}\mathrm{g}\mathrm{n}(k)\{z_0+2[k-\mathrm{s}\mathrm{g}\mathrm{n}(k)]w\}).$(10)

Similar to the case where the magnetoelectric medium is semi-infinite, we have $q_k^{m+}=q_k^{m-}$ and $g_k^{m+}=-g_k^{m-}$. The general formulae for calculating the image-dyon charges in the above mentioned iterating process are derived in Appendix B.1.2, and the results are given in Appendix B.2. Note that the relevant image dyons for determining electromagnetic fields in different regions are different. The fields present in the part of space $z>0$ ($z<-w$) are equivalent to those created by the source charge and image dyons $d_k^{m+}$ ($d_k^{m-}$). The field inside the medium is equivalent to that created by the source charge and the image dyons $d_k^{\mathrm{U}+}$ and $d_k^{\mathrm{L}-}$. By superposing the fields generated by the source charge and all relevant image dyons and taking the asymptotic limit, we can obtain the field configurations in different regions.

The general principles of the image-dyon method and its detailed results can also be formally obtained by rewriting Eq. (9)(9) $V\left(\mathbf{r}\right)=\frac{1}{4\pi \epsilon { }^{\mathrm{\prime }}}\left[\frac{q_0}{\left|\mathbf{r}-{\mathbf{r}}_0\right|}+\sum _k\frac{q_k}{\left|\mathbf{r}-{\mathbf{r}}_k\right|}\right],U\left(\mathbf{r}\right)=\frac{\mu { }^{\mathrm{\prime }}}{4\pi }\sum _k\frac{g_k}{\left|\mathbf{r}-{\mathbf{r}}_k\right|}.$(9) using the geometry-adapted Sommerfeld-type identity [150–152]

(11) $\frac{1}{|\mathbf{r}-{\mathbf{r}}_k|}=\frac{2}{\pi }{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}d\gamma d\eta \frac{cos\left[\gamma (x-x_k)\right]cos\left[\eta (y-y_k)\right]}{u}{\mathrm{e}}^{-u|z-z_k|},$(11)

where $u\equiv \sqrt{{\gamma }^2+{\eta }^2}$. Due to cylindrical symmetry with respect to the $z$ axis, we can focus on calculating the scalar potentials in one of the planes including $z$ axis and then extend the results to other planes by rotation. Thus assuming $\mathbf{r}=(R,0,z)$ in cylindrical coordinates, and for ${\mathbf{r}}_0=(0,0,z_0)$, we find

(12a) $V\left(\mathbf{r}\right)=\frac{1}{2{\pi }^2{\epsilon }^{\prime }}\underset{0}{\overset{\mathrm{\infty }}{\int }}\underset{0}{\overset{\mathrm{\infty }}{\int }}d\gamma d\eta \frac{cos\left(\gamma R\right)}{u}\left[q_0{\mathrm{e}}^{-u\left|z-z_0\right|}+C_1{\mathrm{e}}^{uz}+C_2{\mathrm{e}}^{-uz}\right],$(12a)

(12b) $U(\mathbf{r})=\frac{\mu { }^{\mathrm{\prime }}}{2{\pi }^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}d\gamma d\eta \frac{cos(\gamma R)}{u}\left[C_3{\mathrm{e}}^{uz}+C_4{\mathrm{e}}^{-uz}\right].$(12b)

After determining the initially unknown coefficients $C_j$ via the boundary conditions, $C_{1,2}$ ($C_{3,4}$) can be interpreted in terms of electric charges $q_k$ (magnetic charges $g_k$) and their locations ${\mathbf{r}}_k$. Furthermore, $C_1=C_3=0$ ($C_2=C_4=0$) in the region above (below) the slab since we require scalar potentials to vanish at infinity. Thus the image-dyon representation is obtained via Eq. (11)(11) $\frac{1}{|\mathbf{r}-{\mathbf{r}}_k|}=\frac{2}{\pi }{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}d\gamma d\eta \frac{cos\left[\gamma (x-x_k)\right]cos\left[\eta (y-y_k)\right]}{u}{\mathrm{e}}^{-u|z-z_k|},$(11) .

3.2. Fields from a point charge placed in or near a magnetoelectric plate

Without loss of generality, we assume a positively charged point source ($q_{{\hspace{0.17em}}_0}>0$) from now on and only calculate electromagnetic fields in the plane $\theta =0$. As a basic illustration of the magnetoelectric effect, we first present results for an artificial situation where the magnetoelectric has the same dielectric constant $\epsilon ={\epsilon }_0$ and magnetic permeability $\mu ={\mu }_0$ as the surrounding medium; i.e. the magnetoelectric medium only differs by having finite $\alpha$. Fields for configurations with these parameter values and the electric source charge located near (inside) the magnetoelectric are shown in Figure 4 (Figure 5). Results for a more realistic set of magnetoelectric-medium parameters ($\epsilon =10{\epsilon }_0$, $\mu ={\mu }_0$, $\alpha =3\times {10}^{-4}\sqrt{{\epsilon }_0/{\mu }_0}$) are shown in Figure 6. To provide further insight into the obtained field configurations, we indicate the locations of maxima in the distributions of interfacial currents and charges. Plots of the interface-current distributions in finite-width slab geometries are shown in Figure 7.

Figure 4. Field lines of the electric ($\mathcal{E}$) and magnetic ($\mathcal{B}$) fields generated by an electric point charge located outside an isotropic magnetoelectric medium occupying the space $0>z>-w$ and having the same dielectric constant and magnetic permeability as the surrounding ordinary medium ($\epsilon ={\epsilon }_0$, $\mu ={\mu }_0$). The black dot indicates the location of the source charge, and horizontal thick blue lines delineate interfaces between ordinary and magnetoelectric media. Green circles are positioned where the interface-current distribution has maxima and show the current direction. Yellow circles are positioned where the interface-charge distribution has maxima and show its sign. The magnetoelectric is assumed to have $\alpha =\sqrt{{\epsilon }_0/{\mu }_0}$.

Figure 5. Field lines of the electric ($\mathcal{E}$) and magnetic ($\mathcal{B}$) fields generated by an electric point charge located inside an isotropic magnetoelectric medium occupying the space $0>z>-w$ and having the same dielectric constant and magnetic permeability as the surrounding ordinary medium ($\epsilon ={\epsilon }_0$, $\mu ={\mu }_0$). The black dot indicates the location of the source charge, and horizontal thick blue lines delineate interfaces between ordinary and magnetoelectric media. Green circles are positioned where the interface-current distribution has maxima and show the current direction. Yellow circles are positioned where the interface-charge distribution has maxima and show its sign. The magnetoelectric is assumed to have $\alpha =\sqrt{{\epsilon }_0/{\mu }_0}$.

Figure 6. Field lines of the electric ($\mathcal{E}$) and magnetic ($\mathcal{B}$) fields generated by an electric point charge located outside an isotropic magnetoelectric medium occupying the space $0>z>-w$ and having $\epsilon =10{\epsilon }_0$, $\mu ={\mu }_0$, and $\alpha =3\times {10}^{-4}\sqrt{{\epsilon }_0/{\mu }_0}$. The black dot indicates the location of the source charge, and horizontal thick blue lines delineate interfaces between ordinary and magnetoelectric media. Green circles are positioned where the interface-current distribution has maxima and show the current direction. Yellow circles are positioned where the interface-charge distribution has maxima and show its sign.

Figure 7. Magnetoelectrically induced interface currents in the slab geometries considered in this work. The blue (yellow) curves show the variation of the current-density component ${\mathcal{J}}_y$ at the upper (lower) boundary as a function of coordinate $x/z_0$. Pairs $(x/z_0,{\mathcal{J}}_y)$ of values associated with positive-current maxima are indicated, allowing also the negative-current maxima to be identified via symmetry as $(-x/z_0,-{\mathcal{J}}_y)$. Panel (a) [(b), (c), (d), (e), (f)] pertains to the situations shown in Figure 4(b) [Figure 6(b), Figure 5(b), Figure 4(c), Figure 6(c), Figure 5(c)].

We start by discussing the idealized situation with $\epsilon ={\epsilon }_0$ and $\mu ={\mu }_0$, i.e. when the differentiating characteristic of the region inside (outside) the magnetoelectric medium is a finite (vanishing) $\alpha$. When the entire half-space $z<0$ is occupied by the magnetoelectric, there is only one boundary in the system, on which the surface currents and surface charges are induced by the source electric charge, as indicated in Figure 4(a) and Figure 5(a). Since the image charges are effective representations of the surface currents and surface charges, the two magnetic image monopoles $g_1^{+}$ and $g_1^{-}$ located at different sides of the boundary have the same magnitudes but opposite signs while the two electric image charges $q_1^{+}$ and $q_1^{-}$ are identical. Therefore, when $\epsilon ={\epsilon }_0$ and $\mu ={\mu }_0$, the magnetic field above the boundary is non-zero and mirror-symmetric to the one below but with opposite directions. The electric field, specifically ${\mathcal{E}}_z$, is not continuous across the boundary, which is a result of the superposition of the electric field generated by the real source charge and the one generated by the electric image charges, with the latter ones being symmetric with respect to the boundary. When the width $w$ of the magnetoelectric medium is finite and the source charge placed above it, the magnetic fields are non-zero above and inside the slab but vanish below, as shown in Figure 4(b) and Figure 4(c). The vanishing of the magnetic field below the slab is due to the cancellation effect between the magnetic fields generated by the surface currents on the upper and lower boundaries since the directions of the induced currents are opposite at the two boundaries. Furthermore, the magnitude of the current density at the lower surface is weaker than the one at the upper surface since the lower boundary is further away from the source charge and is more spread out, as can be seen from Figure 7. The electric field has the same form inside and below the slab, which arises due to the boundary conditions at $z=-w$ where $\mathcal{B}=0$. As the slab width $w$ decreases, the magnetic field more and more resembles the one generated by a magnetic quadrupole, especially when the source charge is located inside the magnetoelectric medium [see Figure 5(c)]. This is because the distributions of the surface currents in terms of their magnitudes at the two boundaries become more and more similar as the slab becomes thinner (see Figure 7). Therefore, we expect the magnetic field to be the same as a quadrupole field when the slab is effectively 2D since the upper and lower surface current will generate the same but opposite magnetic dipole moments.

Real materials typically have a dielectric constant $\epsilon$ that is much larger than that of vacuum (${\epsilon }_0$) and quite small values of $\alpha$ measured in units of $\sqrt{{\epsilon }_0/{\mu }_0}$. For example, in the paradigmatic magnetoelectric Cr${ }_2$ O${ }_3$, $\epsilon \approx 10{\epsilon }_0$, $\mu \approx {\mu }_0$, and $\alpha =3\times {10}^{-4}\sqrt{{\epsilon }_0/{\mu }_0}$ [22].⁴ We present results for such a materials combination in Figure 6. In the previously discussed case with $\epsilon ={\epsilon }_0$, induced charges and currents at the interfaces arose entirely as a result of magnetoelectricity, i.e. the piecewise-constant spatial variation of $\alpha$ [see Eq. (6)(6) $\mathrm{\nabla }\cdot \left(\epsilon \mathcal{E}\right)=\rho -\left(\mathrm{\nabla }\alpha \right)\cdot \mathcal{B},\mathrm{\nabla }\times \left({\mu }^{-1}\mathcal{B}\right)-\frac{\mathrm{\partial }\left(\epsilon \mathcal{E}\right)}{\mathrm{\partial }t}=J+\left(\mathrm{\nabla }\alpha \right)\times \mathcal{E}+\frac{\mathrm{\partial }\alpha }{\mathrm{\partial }t}\mathcal{B}.$(6) ]. Due to the additional variation in $\epsilon$ between the ordinary and magnetoelectric media, there will be an extra electric response. Due to the magnetoelectric effect, this extra electric response leads to an extra contribution to the interface current at the lower boundary. As a result, the magnetic field no longer vanishes in the region below the slab. It rather resembles the magnetic field above the medium, especially when the slab is thin, as illustrated in Figure 6.

Figure 7 illustrates in greater detail parametric dependencies of the dipolar interface-current distributions that are induced magnetoelectrically by an electric source charge in slab geometries. The presented results are calculated brute-force by taking $\mathrm{\nabla }\times \mathcal{B}$. The upper (lower) row of plots is associated with the thick (thin) slabs, with plots paired in the same column pertaining to the same source-charge location and the same material combination. The line shapes rationalize the observed quadrupolar patterns of magnetic field lines, especially in the thin-slab limit.

4. Conclusions

We have reviewed the currently available variety of theoretically proposed and experimentally realized 2D magnetoelectric materials. Both the already well-established 2D-quantum-well systems in semiconductor heterostructures and the more recently fabricated 2D atomic crystals are viable platforms for the detailed exploration and possible device application of magnetoelectricity. At present, however, theoretical proposals vastly outnumber actual experimental realizations. Thus a concerted effort is needed to pursue the promising theoretical leads and overcome basic materials challenges.

Where experimental data are available, we have quoted or estimated the magnitudes of magnetoelectric couplings in 2D materials. Sizable values have been demonstrated for bilayer CrI${ }_3$, but most other measurements are indicative of much weaker magnetoelectricity than is typically exhibited in single-phase bulk magnetoelectrics. We expect that significantly larger magnetoelectric coupling should be observed in 2D multiferroic heterostructures, but these materials still need to be realized.

In addition to the survey of materials, we also discussed the implications of low dimensionality for paradigmatic magnetoelectric responses, focusing specifically on magnetic-monopole fields generated by electric charges near or in magnetoelectric media. To solve the associated boundary-value problem, we employed a generalization of the classical image-charge method where images carry both electric and magnetic charges, thus forming image dyons. Previously applied only to treat a single interface between an ordinary material and a magnetoelectric, we extended the formalism to describe a magnetoelectric plate. Finite-size effects are shown to importantly affect the magnetoelectric responses. These results provide a useful guide for further experimental exploration of magnetoelectricity in 2D materials, and they will hopefully stimulate broader investigation of this interesting current topic.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Correction Statement

This article has been corrected with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This work was supported by the Marsden Fund Council from New Zealand government funding (Contract No. VUW1713), managed by the Royal Society Te Apārangi.

Notes

1. Despite being sometimes loosely linked with magnetoelectricity, the current-induced magnetization [63] that is possible in gyrotropic conductors and the chiral magnetic effect exhibited by Weyl semimetals [28,29] are fundamentally different from the equilibrium-magnetoelectric responses focused on here.

2. The first scenario where the application of a magnetic field causes the quadrupolar equilibrium currents that are the basis for magnetoelectricity technically constitutes a field-induced (i.e. higher-order) magnetoelectric effect [60,132]. Note also that the mechanism for magnetoelectricity to occur via equilibrium quadrupolar currents in a ferromagnetic-semiconductor quantum well is fundamentally different from the magnetoelectric response arising from the charge-density tunability of the spontaneous magnetization discussed in Sec. 2.3.

3. To be precise, even-layer MnBi${ }_2$Te${ }_4$ is actually antiferromagnetic and therefore not a time-reversal-invariant topological insulator but an axion insulator [27]. This distinction matters because the fate of axion electrodynamics in the 2D limit of nonmagnetic topological-insulator slab geometries is rather intricate [25].

4. Strictly speaking, the magnetoelectric response of Cr${ }_2$O${ }_3$ is not isotropic but uniaxial, i.e. described by a tensor $\underset{\_}{\alpha }\equiv \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\alpha }_{xx},{\alpha }_{yy},{\alpha }_{zz})$ with ${\alpha }_{xx}={\alpha }_{yy}\ne {\alpha }_{zz}$ [155–157]. In our model for a realistic isotropic magnetoelectric, we use the value of $\alpha \equiv (2{\alpha }_{xx}+{\alpha }_{zz})/3$ [22].

Appendix A.

Conventions for magnetoelectric constitutive relations

Depending on circumstances, different forms of constitutive relations for the responses to electromagnetic fields in a magnetoelectric medium are preferred [79]. If the electric field $\mathcal{E}$ and the magnetic field $\mathcal{B}$ are adopted as controllable quantities, the form

(A1a) $\mathcal{P}={\underset{\_}{\chi }}_{\mathcal{E}}\mathcal{E}+\underset{\_}{\alpha }\mathcal{B},$(A1a)

(A1b) $\mathcal{M}={\underset{\_}{\chi }}_{\mathcal{B}}\mathcal{B}+{\underset{\_}{\alpha }}^{\mathrm{T}}\mathcal{E}$(A1b)

is most suitable. Here $\mathcal{P}\equiv -\mathrm{\partial }F/\mathrm{\partial }\mathcal{E}$ and $\mathcal{M}\equiv -\mathrm{\partial }F/\mathrm{\partial }\mathcal{B}$ are the electric polarization and the magnetization, respectively, given in terms of derivatives of the magnetoelectric material’s free energy per unit volume $F$. See Equation (2)(2) $F(\mathcal{E},\mathcal{B})=F(\mathbf{0},\mathbf{0})-{\mathcal{P}}_{\mathrm{s}}\cdot \mathcal{E}-{\mathcal{M}}_{\mathrm{s}}\cdot \mathcal{B}-\frac{1}{2}\mathcal{E}\cdot ({\underset{\_}{\chi }}_{\mathcal{E}}\mathcal{E})-\frac{1}{2}\mathcal{B}\cdot ({\underset{\_}{\chi }}_{\mathcal{B}}\mathcal{B})\text{break}-\mathcal{E}\cdot (\underset{\_}{\alpha }\mathcal{B})-\dots .$(2) . Alternatively, in situations when $\mathcal{H}\equiv {\mu }_0^{-1}\mathcal{B}-\mathcal{M}$ is fixed instead of $\mathcal{B}$, it is more practical to use the relations

(A2a) $\mathcal{P}={\underset{\_}{\chi }}_{\mathcal{E}}^{\mathrm{\prime }}\mathcal{E}+{\underset{\_}{\alpha }}^{\mathrm{\prime }}\mathcal{H},$(A2a)

(A2b) $\mathcal{M}={\underset{\_}{\chi }}_{\mathcal{H}}\mathcal{H}+{\mu }_0^{-1}{\underset{\_}{\alpha }}^{{ }^{\mathrm{\prime }}\mathrm{T}}\mathcal{H}.$(A2b)

The derived fields $\mathcal{D}\equiv {\epsilon }_0\mathcal{E}+\mathcal{P}$ and $\mathcal{H}$ can then be expressed in terms of the fundamental fields $\mathcal{E}$ and $\mathcal{B}$ equally well using either one of the two related definitions $\underset{\_}{\alpha }$ and ${\underset{\_}{\alpha }}^{\prime }$ of the magnetoelectric tensor;

(A3a) $\mathcal{D}=\underset{-}{\epsilon }\mathcal{E}+\underset{-}{\alpha }\mathcal{B}\equiv \underset{-}{\epsilon }\mathcal{E}+\underset{\_}{{\alpha }^{\prime }}{\underset{-}{\mu }}^{-1}\mathcal{B},$(A3a)

(A3b) $\mathcal{H}={\mu }^{-1}\mathcal{B}-{\alpha }^T\mathcal{E}\equiv {\mu }^{-1}\left(\mathcal{E}-{{\alpha }^{\prime }}^T\mathcal{E}\right),$(A3b)

with $\underset{-}{\epsilon }={\epsilon }_0\underset{-}{1}+{\underset{-}{\chi }}_{\mathcal{E}}$ and $\underset{-}{\mu }={\mu }_0\left(\underset{-}{1}+{\underset{-}{\chi }}_{\mathcal{H}}\right)\equiv ({\mu }_0^{-1}\underset{-}{1}-{\underset{-}{\chi }}_{\mathcal{B}}{)}^{-1}$. Similarly, the constitutive relations that give $\mathcal{D}$ and $\mathcal{B}$ in terms of $\mathcal{E}$ and $\mathcal{H}$ read (with ${\underset{-}{\epsilon }}^{\prime }={\epsilon }_0\underset{-}{1}+{\underset{-}{\chi }}_{\mathcal{E}}^{\prime }$)

(A4a) $\mathcal{D}=\underset{\_}{{\epsilon }^{\prime }}\mathcal{E}+\underset{\_}{{\alpha }^{\prime }}\mathcal{H}\equiv \underset{\_}{{\epsilon }^{\prime }}\mathcal{E}+\underset{-}{\alpha }\underset{-}{\mu }\mathcal{H},$(A4a)

(A4b) $\mathcal{B}=\underset{\_}{\mu }\mathcal{H}+{\underset{\_}{\alpha }}^{\prime \mathrm{T}}\mathcal{E}\equiv \underset{\_}{\mu }(\mathcal{H}+{\underset{\_}{\alpha }}^{\mathrm{T}}\mathcal{E}).$(A4b)

The susceptibilities ${\underset{\_}{\chi }}_{\mathcal{E}}$, ${\underset{\_}{\chi }}_{\mathcal{B}}$, ${\underset{\_}{\chi }}_{\mathcal{H}}$, and by extension $\underset{-}{\epsilon }$ and $\underset{\_}{\mu }$, are symmetric tensors [62], but $\underset{\_}{\alpha }$ and ${\underset{-}{\alpha }}^{\prime }$ are not necessarily symmetric [1,2]. The necessity to distinguish ${\underset{-}{\chi }}_{\mathcal{E}}^{\prime }$ from ${\underset{\_}{\chi }}_{\mathcal{E}}$ arises in a magnetoelectric medium [153] because the condition $\mathcal{H}=0$ at finite $\mathcal{E}$ requires finite $\mathcal{B}\equiv {\underset{\_}{\alpha }}^{\prime \mathrm{T}}\mathcal{E}$, and the latter makes a contribution to $\mathcal{P}$. More specifically, ${\underset{-}{\chi }}_{\mathcal{E}}^{\prime }={\underset{-}{\chi }}_{\mathcal{E}}+\underset{-}{\alpha }\underset{-}{\mu }{\underset{-}{\alpha }}^{\mathrm{T}}\equiv {\underset{-}{\chi }}_{\mathcal{E}}+{\underset{-}{\alpha }}^{\prime }{\underset{-}{\mu }}^{-1}\underset{-}{\alpha }{\prime }^{\mathrm{T}}$. Hence, in principle, the meaning of the dielectric tensor in a magnetoelectric medium needs to be carefully defined based on the physical situation and/or details of the measurement protocol. In contrast, there is no ambiguity associated with $\underset{\_}{\mu }$ or the magnetoelectric tensor $\underset{\_}{\alpha }\equiv {\underset{\_}{\alpha }}^{\prime }{\underset{\_}{\mu }}^{-1}$. Practically, the fundamental limit [154] $|(\underset{\_}{\alpha }{)}_{ij}|\le \sqrt{{({\underset{\_}{\chi }}_{\mathcal{E}})}_{ii}{({\underset{\_}{\chi }}_{\mathcal{B}})}_{jj}}$ on the magnitude of magnetoelectric-tensor elements also constrains the difference between ${\underset{\mathrm{\_}}{\chi }}_{{\hspace{0.17em}}_{\mathcal{E}}}^{\prime }$ and ${\underset{\_}{\chi }}_{\mathcal{E}}$. For example, in a uniaxial magnetoelectric medium where $(\underset{\_}{\mu }{)}_{ij}={\mu }_{ii}{\delta }_{ij}$, $(\underset{\_}{\alpha }{)}_{ij}={\alpha }_{ii}{\delta }_{ij}$ etc., the inequality

(A5) $\frac{|({\underset{\mathrm{\_}}{\chi }}_{\mathcal{E}}^{\prime }-{\underset{\mathrm{\_}}{\chi }}_{\mathcal{E}}{)}_{ii}|}{{({\underset{\mathrm{\_}}{\chi }}_{\mathcal{E}})}_{ii}}\le \frac{{(\underset{\mathrm{\_}}{\mu })}_{ii}}{{\mu }_0}-1$(A5)

holds whose r.h.s. is typically small. Moreover, the magnitudes of magnetoelectric-tensor elements in real materials are generally well below their maximum value mandated by thermodynamic stability [154].

Appendix B.

Formalism for determining image dyons

B.1. Case of a semi-infinite magnetoelectric medium

B.1.1. Electric point charge as the source

To discuss the case where an electric charge is outside a semi-infinite magnetoelectric medium, we assume the region $z>0$ ($z<0$) to be vacuum (a magnetoelectric) as illustrated in Figure 2(a). We introduce specific Ansätze for the scalar potentials;

(B1a) $V(R,\theta ,z>0)=\frac{1}{4\pi {\epsilon }_0}\left[\frac{q_0}{\sqrt{{\left(z-z_0\right)}^2+R^2}}+\frac{q_1^{-}}{\sqrt{{\left(z+z_0\right)}^2+R^2}}\right],$(B1a)

(B1b) $V(R,\theta ,z<0)=\frac{1}{4\pi {\epsilon }_0}\left[\frac{q_0+q_1^{+}}{\sqrt{{\left(z-z_0\right)}^2+R^2}}\right],$(B1b)

(B1c) $U(R,\theta ,z>0)=\frac{{\mu }_0}{4\pi }\left[\frac{g_1^{-}}{\sqrt{{(z+z_0)}^2+R^2}}\right],$(B1c)

(B1d) $U(R,\theta ,z<0)=\frac{{\mu }_0}{4\pi }\left[\frac{g_1^{+}}{\sqrt{{(z-z_0)}^2+r^2}}\right].$(B1d)

Applying boundary conditions as specified by Eq. (7(7a) ${\mathcal{D}}_z(R,\theta ,z_{\mathrm{b}}^{+})={\mathcal{D}}_z(R,\theta ,z_{\mathrm{b}}^{-}),{\mathcal{H}}_R(R,\theta ,z_{\mathrm{b}}^{+})={\mathcal{H}}_R(R,\theta ,z_{\mathrm{b}}^{-}),$(7a) ), we find the image charges and monopoles as

(B2a) $q_1^{+}=q_1^{-}=-q_0\frac{\left(\epsilon -{\epsilon }_0\right)\left(1/\mu +1/{\mu }_0\right)+{\alpha }^2}{\left(\epsilon +{\epsilon }_0\right)\left(1/\mu +1/{\mu }_0\right)+{\alpha }^2},$(B2a)

(B2b) $g_1^{+}=-g_1^{-}=\frac{q_0}{{\mu }_0}\frac{2\alpha }{\left(\epsilon +{\epsilon }_0\right)\left(1/\mu +1/{\mu }_0\right)+{\alpha }^2}.$(B2b)

Similarly, when the source charge is inside the magnetoelectric medium and located at $(0,0,-z_0)$, where $z_0>0$ [see Figure 2(b)], we use the Ansätze

(B3a) $V(R,\theta ,z>0)=\frac{1}{4\pi \epsilon }\left[\frac{q_0+q_1^{-}}{\sqrt{{\left(z+z_0\right)}^2+R^2}}\right],$(B3a)

(B3b) $V(R,\theta ,z<0)=\frac{1}{4\pi \epsilon }\left[\frac{q_0}{\sqrt{{\left(z+z_0\right)}^2+R^2}}+\frac{q_1^{+}}{\sqrt{{\left(z-z_0\right)}^2+R^2}}\right],$(B3b)

(B3c) $U(R,\theta ,z>0)=\frac{\mu }{4\pi }\left[\frac{g_1^{-}}{\sqrt{{(z+z_0)}^2+R^2}}\right],$(B3c)

(B3d) $U(R,\theta ,z<0)=\frac{\mu }{4\pi }\left[\frac{g_1^{+}}{\sqrt{{(z-z_0)}^2+R^2}}\right].$(B3d)

The image charges satisfying the boundary conditions are found as

(B4a) $q_1^{+}=q_1^{-}=-q_0\frac{\left(-\epsilon +{\epsilon }_0\right)\left(1/\mu +1/{\mu }_0\right)+{\alpha }^2}{\left(\epsilon +{\epsilon }_0\right)\left(1/\mu +1/{\mu }_0\right)+{\alpha }^2},$(B4a)

(B4b) $g_1^{+}=-g_1^{-}=\frac{q_0}{{\mu }_0}\frac{2\alpha }{\left(\epsilon +{\epsilon }_0\right)\left(1/\mu +1/{\mu }_0\right)+{\alpha }^2}.$(B4b)

Note that the image charges from Eq. (B4) are the same as the ones in Eq. (B2) if we simultaneously swap $\epsilon$ with ${\epsilon }_0$ and $\mu$ with ${\mu }_0$. Also, in the limit $\alpha =0$, we recover the known results for the situation where an electric point charge is near/inside a semi-infinite dielectric medium.

B.1.2. Dyon as the source

To assist the calculation for the cases where the magnetoelectric medium is a slab with finite width, it is useful to obtain the results when the source is a dyon, i.e. a point particle located at ${\mathbf{r}}_{\mathrm{d}}=(0,0,\pm z_{\mathrm{d}})$ having both magnetic and electric charges $g_{\mathrm{d}}$ and $q_{\mathrm{d}}$, respectively. Considering this is needed as part of our iterative process for satisfying boundary conditions at interfaces with magnetoelectrics.

If we started with the upper boundary, the dyon is always on the side of the magnetoelectric medium. Assuming that the region $z<0$ is where the dyon is located, the Ansätze are

(B5a) $V(R,\theta ,z>0)=\frac{1}{4\pi \epsilon { }^{\mathrm{\prime }}}\left[\frac{q_{\mathrm{d}}+q_1^{-}}{\sqrt{{\left(z+z_{\mathrm{d}}\right)}^2+R^2}}\right],$(B5a)

(B5b) $V(R,\theta ,z<0)=\frac{1}{4\pi \epsilon { }^{\mathrm{\prime }}}\left[\frac{q_{\mathrm{d}}}{\sqrt{{\left(z+z_{\mathrm{d}}\right)}^2+R^2}}+\frac{q_1^{+}}{\sqrt{{\left(z-z_{\mathrm{d}}\right)}^2+R^2}}\right],$(B5b)

(B5c) $U(R,\theta ,z>0)=\frac{\mu { }^{\mathrm{\prime }}}{4\pi }\left[\frac{g_{\mathrm{d}}+g_1^{-}}{\sqrt{{(z+z_{\mathrm{d}})}^2+R^2}}\right],$(B5c)

(B5d) $U(R,\theta ,z<0)=\frac{\mu { }^{\mathrm{\prime }}}{4\pi }\left[\frac{g_{\mathrm{d}}}{\sqrt{{(z+z_{\mathrm{d}})}^2+R^2}}+\frac{g_1^{+}}{\sqrt{{(z-z_{\mathrm{d}})}^2+R^2}}\right].$(B5d)

Depending on whether the original electric source charge $q_0$ is outside or inside the magnetoelectric slab, we choose ${\epsilon }_0$ (${\mu }_0$) or $\epsilon$ ($\mu$) for $\epsilon { }^{\mathrm{\prime }}$ ($\mu { }^{\mathrm{\prime }}$), respectively, to keep with our overall convention.

The resulting image charges for the case when the original source charge $q_0$ is outside the magnetoelectric slab are

(B6a) $q_1^{+}=q_1^{-}=\frac{-q_{\mathrm{d}}\left[\left(-\epsilon +{\epsilon }_0\right)\left(1/\mu +1/{\mu }_0\right)+{\alpha }^2\right]+g_{\mathrm{d}}\left(2\alpha {\epsilon }_0{\mu }_0/\mu \right)}{\left(\epsilon +{\epsilon }_0\right)\left(1/\mu +1/{\mu }_0\right)+{\alpha }^2},$(B6a)

(B6b) $g_1^{+}=-g_1^{-}=\frac{q_{\mathrm{d}}\left(2\alpha \epsilon /{\mu }_0{\epsilon }_0\right)+g_{\mathrm{d}}\left[\left(\epsilon +{\epsilon }_0\right)\left(1/\mu -1/{\mu }_0\right)+{\alpha }^2\right]}{\left[\left(\epsilon +{\epsilon }_0\right)\left(1/\mu +1/{\mu }_0\right)+{\alpha }^2\right]}.$(B6b)

Alternatively, when the original source charge $q_0$ is inside the magnetoelectric slab, the image charges are

(B7a) $q_1^{+}=q_1^{-}=\frac{-q_{\mathrm{d}}\left[\left(-\epsilon +{\epsilon }_0\right)\left(1/\mu +1/{\mu }_0\right)+{\alpha }^2\right]+g_{\mathrm{d}}\left(2\alpha \epsilon \right)}{\left(\epsilon +{\epsilon }_0\right)\left(1/\mu +1/{\mu }_0\right)+{\alpha }^2}.$(B7a)

(B7b) $g_1^{+}=-g_1^{-}=\frac{q_{\mathrm{d}}\left(2\alpha /\mu \right)+g_{\mathrm{d}}\left[\left(\epsilon +{\epsilon }_0\right)\left(1/\mu -1/{\mu }_0\right)+{\alpha }^2\right]}{\left[\left(\epsilon +{\epsilon }_0\right)\left(1/\mu +1/{\mu }_0\right)+{\alpha }^2\right]}.$(B7b)

B.2. Case of a finite-width magnetoelectric slab

B.2.1. Electric point charge located outside the slab

When $z_0>0$, the Ansätze are

(B8a) $V(R,\theta ,z>0)=\frac{1}{2{\pi }^2{\epsilon }_0}\underset{0}{\overset{\mathrm{\infty }}{\int }}\underset{0}{\overset{\mathrm{\infty }}{\int }}d\gamma d\eta \frac{cos\left(\gamma R\right)}{u}\left(q_0{\mathrm{e}}^{-u\left|z-z_0\right|}+A{\mathrm{e}}^{-uz}\right),$(B8a)

(B8b) $V(R,\theta ,0>z>-w)=\frac{1}{2{\pi }^2{\epsilon }_0}\underset{0}{\overset{\mathrm{\infty }}{\int }}\underset{0}{\overset{\mathrm{\infty }}{\int }}d\gamma d\eta \frac{cos\left(\gamma R\right)}{u}\left(B{\mathrm{e}}^{uz}+C{\mathrm{e}}^{-uz}\right),$(B8b)

(B8c) $V(R,\theta ,z<-w)=\frac{1}{2{\pi }^2{\epsilon }_0}\underset{0}{\overset{\mathrm{\infty }}{\int }}\underset{0}{\overset{\mathrm{\infty }}{\int }}d\gamma d\eta \frac{cos\left(\gamma R\right)}{u}D{\mathrm{e}}^{uz},$(B8c)

(B8d) $U(R,\theta ,z>0)=\frac{{\mu }_0}{2{\pi }^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}d\gamma d\eta \frac{cos(\gamma R)}{u}E{\mathrm{e}}^{-uz},$(B8d)

(B8e) $U(R,\theta ,0>z>-w)=\frac{{\mu }_0}{2{\pi }^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}d\gamma d\eta \frac{cos(\gamma R)}{u}\left(F{\mathrm{e}}^{uz}+G{\mathrm{e}}^{-uz}\right),$(B8e)

(B8f) $U(R,\theta ,z<-w)=\frac{{\mu }_0}{2{\pi }^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}d\gamma d\eta \frac{cos(\gamma R)}{u}H{\mathrm{e}}^{uz}.$(B8f)

Here we absorbed the source-charge terms into $B$, $C$ and $D$ in the Ansätze for $V(R,\theta ,0>z>-w)$ and $V(R,\theta ,z<-w)$. Below we list the solutions when ${\epsilon }_0=\epsilon$ and ${\mu }_0=\mu$ as an example, i.e. when the only difference between the regions inside and outside the slab is the finite magnetoelectricity coupling inside.

(B9) $\begin{array}{cc}A& =q_0\frac{-{\mathrm{e}}^{-u{z}_0}\sigma +{\mathrm{e}}^{-u(z_0+2w)}\sigma }{1-{\mathrm{e}}^{-2uw}\sigma },\\ B& =q_0\frac{4{\epsilon }_0}{4{\epsilon }_0+{\alpha }^2{\mu }_0}\frac{{\mathrm{e}}^{-u{z}_0}}{1-{\mathrm{e}}^{-2uw}\sigma },C=0,\\ D& =B,\\ E& =q_0\frac{2\alpha }{4{\epsilon }_0+{\alpha }^2{\mu }_0}\frac{-{\mathrm{e}}^{-u{z}_0}+{\mathrm{e}}^{-u(z_0+2w)}}{1-{\mathrm{e}}^{-2uw}\sigma },\\ F& =q_0\frac{2\alpha }{4{\epsilon }_0+{\alpha }^2{\mu }_0}\frac{{\mathrm{e}}^{-u{z}_0}}{1-{\mathrm{e}}^{-2uw}\sigma },G=q_0\frac{2\alpha }{4{\epsilon }_0+{\alpha }^2{\mu }_0}\frac{{\mathrm{e}}^{-u(z_0+2w)}}{1-{\mathrm{e}}^{-2uw}\sigma },\\ H& =0,\end{array}$(B9)

with $\sigma =\frac{{\alpha }^2{\mu }_0}{4{\epsilon }_0+{\alpha }^2{\mu }_0}$. We can interpret the above results in terms of image charges shown in Table B1 by applying the following relation:

(B10) $\frac{1}{1-{\mathrm{e}}^{-2uw}\sigma }=\sum _{k=1}^{\mathrm{\infty }}{\mathrm{e}}^{-2kuw}{\sigma }^{k-1}.$(B10)

From expressions given in Table B1, we can see that, when ${\epsilon }_0=\epsilon$ and ${\mu }_0=\mu$, $V(R,\theta ,0>z>-w)$ and $V(R,\theta ,z<-w)$ share the same expression, and $U(R,\theta ,z<-w)=0$. Thus, the electric field is continuous and $\mathcal{B}=0$ in the region where $z<0$, as shown in Figure 4.

Table B1. The potentials in terms of image charges when the point charge is above the magnetoelectric slab. Note that $d_1^{\mathrm{L}-}$ does not exist by construction

	Potential	Image charge when ${\epsilon }_0=\epsilon$ and ${\mu }_0=\mu$
$V(R,\theta ,z>0)$	$\frac{1}{4\pi {\epsilon }_0}\left[\frac{q_0}{\left|\mathbf{r}-{\mathbf{r}}_0\right|}+\sum _{m=\mathrm{U},\mathrm{L}}\sum _{k=1}^{\mathrm{\infty }}\frac{q_k^{m-}}{\left|\mathbf{r}-{\mathbf{r}}_k^{-}\right|}\right]$	$q_k^{\mathrm{L}\pm }=0$$q_k^{\mathrm{U}\pm }=\{\begin{array}{cc}-q_0\sigma & k=1\\ q_0\left(1-\sigma \right){\sigma }^{\left(k-1\right)}& k\ge 2\end{array}$
$V(R,\theta ,0>z>-w)$	$\frac{1}{4\pi {\epsilon }_0}\left[\frac{q_0}{\left|\mathbf{r}-{\mathbf{r}}_0\right|}+\sum _{k=1}^{\mathrm{\infty }}\left(\frac{q_k^{\mathrm{U}+}}{\left|\mathbf{r}-{\mathbf{r}}_k^{+}\right|}+\frac{q_k^{\mathrm{L}-}}{\left|\mathbf{r}-{\mathbf{r}}_k^{-}\right|}\right)\right]$
$V(R,\theta ,z<-w)$	$\frac{1}{4\pi {\epsilon }_0}\left[\frac{q_0}{\left|\mathbf{r}-{\mathbf{r}}_0\right|}+\sum _{m=\mathrm{U},\mathrm{L}}\sum _{k=1}^{\mathrm{\infty }}\frac{q_k^{m+}}{\left|\mathbf{r}-{\mathbf{r}}_k^{+}\right|}\right]$
$U(R,\theta ,z>0)$	$\frac{{\mu }_0}{4\pi }\sum _{m=\mathrm{U},\mathrm{L}}\sum _{k=1}^{\mathrm{\infty }}\frac{g_k^{m-}}{|\mathbf{r}-{\mathbf{r}}_k^{-}|}$	$\begin{array}{cc}\hspace{0.17em}& g_k^{\mathrm{U}+}=-g_k^{\mathrm{U}-}\\ =& -g_k^{\mathrm{L}+}=g_{k+1}^{\mathrm{L}-}\\ =& q_0\frac{2\alpha }{4{\epsilon }_0+{\alpha }^2{\mu }_0}{\sigma }^{k-1}\end{array}$
$U(R,\theta ,0>z>-w)$	$\frac{{\mu }_0}{4\pi }\sum _{k=1}^{\mathrm{\infty }}\left(\frac{g_k^{\mathrm{U}+}}{|\mathbf{r}-{\mathbf{r}}_k^{+}|}+\frac{g_k^{\mathrm{L}-}}{|\mathbf{r}-{\mathbf{r}}_k^{-}|}\right)$
$U(R,\theta ,z<-w)$	$\frac{{\mu }_0}{4\pi }\sum _{m=\mathrm{U},\mathrm{L}}\sum _{k=1}^{\mathrm{\infty }}\frac{g_k^{m+}}{|\mathbf{r}-{\mathbf{r}}_k^{+}|}$

Table B2. The potentials in terms of image charges when the point charge is in the magnetoelectric slab. Note that $d_{-1}^{m-}$, $d_{-1}^{\mathrm{U}+}$ and $d_1^{\mathrm{L}-}$ do not exist by construction

	Potential	Image charge when ${\epsilon }_0=\epsilon$ and ${\mu }_0=\mu$
$V(R,\theta ,z>0)$	$\frac{1}{4\pi \epsilon }\left[\frac{q_0}{\left|\mathbf{r}-{\mathbf{r}}_0\right|}+\sum _{m=\mathrm{U},\mathrm{L}}\sum _{k=1}^{\mathrm{\infty }}\frac{q_{\pm k}^{m-}}{\left|\mathbf{r}-{\mathbf{r}}_{\pm k}^{-}\right|}\right]$	$\begin{array}{rl}& q_k^{\mathrm{U}+}=-q_{-(k+1)}^{\mathrm{U}+}\\ & =q_{-k}^{\mathrm{L}+}=-q_k^{\mathrm{L}+}=-q_0{\sigma }^k;\\ & q_{-(k+1)}^{\mathrm{L}-}=-q_{k+1}^{\mathrm{L}-}\\ & =q_k^{\mathrm{U}-}=-q_{-(k+1)}^{\mathrm{U}-}=-q_0{\sigma }^{k-1}\end{array}$
$V(R,\theta ,0>z>-w)$	$\frac{1}{4\pi \epsilon }\left[\frac{q_0}{\left|\mathbf{r}-{\mathbf{r}}_0\right|}+\sum _{k=1}^{\mathrm{\infty }}\left(\frac{q_{\pm k}^{\mathrm{U}+}}{\left|\mathbf{r}-{\mathbf{r}}_{\pm k}^{+}\right|}+\frac{q_{\pm k}^{\mathrm{L}-}}{\left|\mathbf{r}-{\mathbf{r}}_{\pm k}^{-}\right|}\right)\right]$
$V(R,\theta ,z<-w)$	$\frac{1}{4\pi \epsilon }\left[\frac{q_0}{\left|\mathbf{r}-{\mathbf{r}}_0\right|}+\sum _{m=\mathrm{U},\mathrm{L}}\sum _{k=1}^{\mathrm{\infty }}\frac{q_{\pm k}^{m+}}{\left|\mathbf{r}-{\mathbf{r}}_{\pm k}^{+}\right|}\right]$
$U(R,\theta ,z>0)$	$\frac{\mu }{4\pi }\sum _{m=\mathrm{U},\mathrm{L}}\sum _{k=1}^{\mathrm{\infty }}\frac{g_{\pm k}^{m-}}{|\mathbf{r}-{\mathbf{r}}_{\pm k}^{-}|}$	$\begin{array}{c}{g}_{-k}^{\mathrm{U}+}=g_k^{\mathrm{L}+}\\ =g_k^{\mathrm{L}-}=g_{-k}^{\mathrm{U}-}=0;\\ g_k^{\mathrm{U}+}=-g_{-k}^{\mathrm{L}+}\\ =g_{-\left(k+1\right)}^{\mathrm{L}-}=-g_k^{\mathrm{U}-}=q_0\frac{2\alpha }{4{\epsilon }_0+{\alpha }^2{\mu }_0}{\sigma }^{k-1}\end{array}$
$U(R,\theta ,0>z>-w)$	$\frac{\mu }{4\pi }\sum _{k=1}^{\mathrm{\infty }}\left(\frac{g_{\pm k}^{\mathrm{U}+}}{|\mathbf{r}-{\mathbf{r}}_{\pm k}^{+}|}+\frac{g_{\pm k}^{\mathrm{L}-}}{|\mathbf{r}-{\mathbf{r}}_{\pm k}^{-}|}\right)$
$U(R,\theta ,z<-w)$	$\frac{\mu }{4\pi }\sum _{m=\mathrm{U},\mathrm{L}}\sum _{k=1}^{\mathrm{\infty }}\frac{g_{\pm k}^{m+}}{|\mathbf{r}-{\mathbf{r}}_{\pm k}^{+}|}$

B.2.2. Electric point charge located inside the slab

When $z_0>0$ and ${\mathbf{r}}_0=(0,0,-z_0)$, the Ansätze are

(B11a) $V(R,\theta ,z>0)=\frac{1}{2{\pi }^2\epsilon }\underset{0}{\overset{\mathrm{\infty }}{\int }}\underset{0}{\overset{\mathrm{\infty }}{\int }}d\gamma d\eta \frac{cos\left(\gamma R\right)}{u}A{\mathrm{e}}^{-uz},$(B11a)

(B11b) $V(R,\theta ,0>z>-w)=\frac{1}{2{\pi }^2\epsilon }\underset{0}{\overset{\mathrm{\infty }}{\int }}\underset{0}{\overset{\mathrm{\infty }}{\int }}d\gamma d\eta \frac{cos\left(\gamma R\right)}{u}\left(q_0{\mathrm{e}}^{-u\left|z+z_0\right|}+B{\mathrm{e}}^{uz}+C{\mathrm{e}}^{-uz}\right),$(B11b)

(B11c) $V(R,\theta ,z<-w)=\frac{1}{2{\pi }^2\epsilon }\underset{0}{\overset{\mathrm{\infty }}{\int }}\underset{0}{\overset{\mathrm{\infty }}{\int }}d\gamma d\eta \frac{cos\left(\gamma R\right)}{u}D{\mathrm{e}}^{uz},$(B11c)

(B11d) $U(R,\theta ,z>0)=\frac{\mu }{2{\pi }^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}d\gamma d\eta \frac{cos(\gamma R)}{u}E{\mathrm{e}}^{-uz},$(B11d)

(B11e) $U(R,\theta ,0>z>-w)=\frac{\mu }{2{\pi }^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}d\gamma d\eta \frac{cos(\gamma R)}{u}\left(F{\mathrm{e}}^{uz}+G{\mathrm{e}}^{-uz}\right),$(B11e)

(B11f) $U(R,\theta ,z<-w)=\frac{\mu }{2{\pi }^2}{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}d\gamma d\eta \frac{cos(\gamma R)}{u}H{\mathrm{e}}^{uz}.$(B11f)

Here we absorbed the source-charge terms into $A$ and $D$ in the Ansätze for $V(R,\theta ,z>0)$ and $V(R,\theta ,z<-w)$. Again, we list the solutions when ${\epsilon }_0=\epsilon$ and ${\mu }_0=\mu$ as an example:

(B12) $\begin{array}{cc}A& =q_0\frac{4{\epsilon }_0}{4{\epsilon }_0+{\alpha }^2{\mu }_0}\frac{{\mathrm{e}}^{-u{z}_0}}{1-{\mathrm{e}}^{-2uw}\sigma },\\ B& =q_0\frac{-{\mathrm{e}}^{-u{z}_0}\sigma +{\mathrm{e}}^{u(z_0-2w)}\sigma }{1-{\mathrm{e}}^{-2uw}\sigma },C=q_0\frac{{\mathrm{e}}^{-u(z_0+2w)}\sigma -{\mathrm{e}}^{-u(-z_0+2w)}\sigma }{1-{\mathrm{e}}^{-2uw}\sigma },\\ D& =q_0\frac{4{\epsilon }_0}{4{\epsilon }_0+{\alpha }^2{\mu }_0}\frac{{\mathrm{e}}^{u{z}_0}}{1-{\mathrm{e}}^{-2uw}\sigma },\\ E& =q_0\frac{2\alpha }{4{\epsilon }_0+{\alpha }^2{\mu }_0}\frac{-{\mathrm{e}}^{-u{z}_0}+{\mathrm{e}}^{-u(-z_0+2w)}}{1-{\mathrm{e}}^{-2uw}\sigma },\\ F& =q_0\frac{2\alpha }{4{\epsilon }_0+{\alpha }^2{\mu }_0}\frac{{\mathrm{e}}^{-u{z}_0}}{1-{\mathrm{e}}^{-2uw}\sigma },G=q_0\frac{2\alpha }{4{\epsilon }_0+{\alpha }^2{\mu }_0}\frac{{\mathrm{e}}^{-u(-z_0+2w)}}{1-{\mathrm{e}}^{-2uw}\sigma },\\ H& =q_0\frac{2\alpha }{4{\epsilon }_0+{\alpha }^2{\mu }_0}\frac{{\mathrm{e}}^{-u{z}_0}-{\mathrm{e}}^{u{z}_0}}{1-{\mathrm{e}}^{-2uw}\sigma }.\\ \hspace{0.17em}& \hspace{0.17em}\end{array}$(B12)

The corresponding image charges we get are shown in Table B2.

We can derive from Table B2 that, when ${\epsilon }_0=\epsilon$ and ${\mu }_0=\mu$, the expressions for the potentials can be simplified since ${\sum }_{m=\mathrm{U},\mathrm{L}}{q}_k^{m+}={\sum }_{m=\mathrm{U},\mathrm{L}}{q}_{-k}^{m-}=0$, ${\sum }_{m=\mathrm{U},\mathrm{L}}{q}_{-k}^{m+}={\sum }_{m=\mathrm{U},\mathrm{L}}{q}_k^{m-}=q_0(1-\sigma ){\sigma }^{k-1}$, $\sum _{m=\mathrm{U},\mathrm{L}}{g}_k^{m+}=\sum _{m=\mathrm{U},\mathrm{L}}{g}_{-k}^{m-}=q_0\frac{2\alpha }{4{\epsilon }_0+{\alpha }^2{\mu }_0}{\sigma }^{k-1}$, and $\sum _{m=\mathrm{U},\mathrm{L}}{g}_{-k}^{m+}=\sum _{m=\mathrm{U},\mathrm{L}}{g}_{k+1}^{m-}=-q_0\frac{2\alpha }{4{\epsilon }_0+{\alpha }^2{\mu }_0}{\sigma }^k$, which helps us to verify the relation between the image charges in Table B2 and the expressions of the constants in Eqs. (B12).