1 Introduction

The radiation produced by an electric dipole near a planar interface has been well studied over the years and has remained a relevant subject of research for physicists and engineers due its relevance in a wide range of phenomena like practical applications in radio communications ^[1], THz Zenneck wave propagation ^[2], near-field optics ^[3], plasmonics ^[4] and nanophotonics ^[5], just to mention a few examples. In 1909, Sommerfeld ^[6] published a theory for a radiating vertically oriented dipole above a planar and lossy ground which formed the basis for subsequent investigations ^{[7,8,9,10,11]}. Probably, by the fact that the early theory of dipole emission near planar interfaces was written in German, although there was an English version summarized in Sommerfeld’s lectures on theoretical physics ^[12], many aspects of the theory were reinvented and clarified over the years ^{[13,14,15,16]}.

Here we consider planar interfaces constructed with linear homogeneous and isotropic magnetoelectric (ME) media giving rise to the so called magnetoelectric effect (MEE), whereby electric (magnetic) fields are able to induce additional magnetization (polarization) in the material. This effect, which was predicted ^[17] and discovered ^[18] in antiferromagnets, has been widely studied along the years in multiferroic materials and it is codified in an additional parameter of the material: the magnetoelectric polarizability (MEP) ^[19]. The recent discovery of three-dimensional topological insulators (TIs) has boosted the interest in this topic by providing new materials where this effect is predicted to occur ^{[20,21,22,23,24,25]}. Generally speaking, TIs belong to a novel state of matter in which the characterization of their quantum states does not fit into the standard paradigm of condensed matter physics whereby the phases of the material are classified according to order parameters arising from the spontaneous symmetry breaking of the corresponding Hamiltonian according to the phenomenological Ginzburg-Landau theory. A distinguished example of this classification are normal superconductors, where gauge invariance is spontaneously broken. Instead, these states are classified according to topological invariants that arise in the Hilbert space generated by the corresponding Hamiltonians in the reciprocal space of the crystal lattice. They are protected by time-reversal-symmetry (TRS) and admit an insulating bulk together with conducting surface (edge) states. The imposition of this symmetry yields two classes of materials: standard insulators labeled by a zero MEP and TIs characterized by a MEP equal to π. These new materials host a number of exceptional electromagnetic properties. Among them we have: (i) they can carry currents along the edge channels without dissipation, (ii) their MEP is quantized, (iii) the conducting edge states can be interpreted as quasi-particles being massless Dirac fermions. This by itself is an important feature that makes contact with high energy physics and which provides the opportunity to investigate the existence of unseen particles like Majorana fermions, for example. (iv) they are predicted to exhibit the quantized photogalvanic effect in which light can induce a quantized current. For an extensive review of the properties of TIs see for example ^[26,27,28]. All these new features provide additional motivation to reconsider the problem of radiation in magnetolectric media.

The effective theoretical framework to deal with magnetoelectric media is motivated by axion electrodynamics ^[29], which consists of adding the term La=gaγγ a(t,x)FμνF~μν to the Maxwell Lagrangian density Lem, plus a kinetic and a mass terms for the pseudoscalar field a(t, x). The so called axion field a(t, x) was introduced in Ref. ^[30] to propose a solution for the strong CP problem in strong interactions ^[31,32]. In the original formulation ^[29], the coupling constant gaγγ arised from a specific grand unification model of the strong and electroweak interactions. Also, Fμν is the electromagnetic tensor and F~μν=(1/2)ϵμναβFαβ is the dual tensor. The well known relation FμνF~μν=-4 E⋅B allows to rewrite La in terms of the electric E and magnetic B fields. As we will show in the following the coupling La encapsulates the MEE which characterizes the electromagnetic response of the materials we consider in this work. Thus we restrict ourselves to a non-dynamical axion field a(t,x)→ϑ(x), to be called the magnetoelectric polarizability (MEP), which we consider as an additional electromagnetic property of the medium, in the same footing as its permittivity and permeability ^[23,33,34]. Following a standard convention we now consider the interaction term Lϑ=-(α/4π2) ϑ(x)E⋅B, where α is the fine structure constant characteristic of the electromagnetic interaction between fermions and photons in the material, which produces this effective term. We call ϑ-Electrodynamics (ϑ-ED) this restriction of axion electrodynamics and our purpose is to study the radiation produced by a dipole oriented vertically with respect to the interface between two semi-infinite planar magnetoelectric media, having different constant MEPs. Excluding important differences in their microscopic structure, we will refer to the medium as a magnetoelectric, or a ϑ-medium, as long as its macroscopic electromagnetic response can be described in the framework of the ϑ-ED.

The paper is organized as follows. In Sec. 2, we present a review of ϑ-ED which also contains a summary of the calculation of the time-dependent Green’s function (GF) for the 4-potential A^μ in our setup. As the source of the electromagnetic fields, in Sec. 3 we introduce an oscillating vertically oriented point-like electric dipole located at a distance z₀ > 0, on the z-axis perpendicular to the interface between the two media. The convolution of this source with the GF is carried out in the Subsec. 3.1 and yields the corresponding electromagnetic potentials Aμ in terms of closed integrals which are calculated in the Appendix A. The next step is to obtain the far-field approximation of those integrals. This is performed using a modified version of the steepest descent method, which is appropriate to the situation where the integrand is not a smooth function in the vicinity of the stationary phase due to the appearance of poles in the steepest descent path at this point.

This approximation,which heavily relies upon Ref. ^[13] is explained in detail and carefully carried out in the Appendix B. These results are summarized in the Subsec. 3.1.

As a consequence of the presence of the pole we find that the 4-potential acquires a contribution from axially symmetric cylindrical waves (denoted also as surface waves) besides the standard spherically symmetric ones. A detailed analysis on the former kind of waves allows us to introduce what we call the discarding angle θ₀, which permits us to divide the space in two regions: V1 where the cylindrical wave contribution can be neglected and V2 where this contribution has to be considered within a certain range of parameters in what is called the intermediate zone in the literature ^[13]. To characterize the relevance of these cylindrical waves we introduce a rapidly decreasing function measuring their amplitude and realize that for observation distances further away from the intermediate zone in the region V2 they turn out to be very much suppressed with respect to the spherical ones. This situation is quantitatively explained in detail also in the Subsec. 3.1. In Subsec. 3.2, the far-field expressions for E and B are calculated for each region. In Sec. 4 we consider the angular distribution, the total radiated power and the energy transport of the dipolar radiation. In the Subsec. 4.1 we establish the numerical parameters to be used in the subsequent applications. Subsection 4.2 comprises a detailed examination of the angular distribution spectrum dP/dΩ in the region V1. In Subsec. 4.3, we calculate the power radiated into the region V1. Subsection 4.4 is devoted to the energy transport of the radiation in the region V2. Here we also give some numerical estimations considering the media already discussed in Subsecs. 4.2 and 4.3, plus some additional hypothetical choices. Finally, Sec. 5 provides a concluding summary and the conclusions from our results. In the Appendix A we derive the exact expressions for the potential A_μ required to calculate the electromagnetic fields. The far-field approximation is carried out in the Appendix B using a modified steepest descent method. The final Appendix C includes a brief review of the Faddeeva plasma dispersion function which arises in the discussion of the cylindrical waves. Throughout this paper we use Gaussian units with ℏ=c=1, the metric signature is (+,—, —, —,) and ε⁰¹²³ = 1. We follow the conventions of Ref. ^[35].

2 ϑ-Electrodynamics

Let us consider two semi-infinite ϑ-media separated by a planar interface located at z = 0, filling the regions U1,z<0 and U2,z>0 of the space. We take both media to be non-magnetic, i.e. μ₁ = μ₂ = 1 and with the same permittivity ϵ1=ϵ2=ϵ. This condition is motivated by the results of Ref. ^[36], which show that the effects of the MEE are substantially enhanced with respect to the optical contributions when both ϑ-media have the same dielectric constant and permeability. Additionally we assume the parameter ϑ to be piecewise constant so that it takes the values ϑ=ϑ1 in the region U1 and ϑ=ϑ2 in the region U2. This is expressed as

where Θ(z) is the standard Heaviside function with Θ(z)=1,    z≥0   and Θ(z)=0,    z<0. The dynamics is governed by the standard Maxwell equations in a material medium ^[35,37]

which require to specify additional constitutive relations characterizing the medium under consideration. In the case of the magnetoelectric media the constitutive relations are

Here α is the fine-structure constant, ϱ   and j are the external sources given by the charge and current densities respectively. Substituting the constitutive relations (3) into the inhomogeneous equations (2) and using the MEP given in (1) yields our final equations

where u^ is the outward unit normal to the region U1 and

In the case of a TI located in the region U2 (ϑ2=π) in front of a regular insulator (ϑ=0) in region U1, we have θ~=α(2m~+1), with m~ being an integer depending on the details of the TRS breaking at the interface between the two materials.

The homogeneous Maxwell equations still enable us to define the electromagnetic fields E and B in terms the electromagnetic potentials Φ and A as

We observe that Eqs. (4) and (5), together with the constitutive relations (3), can also be derived from the action

which clearly incorporates the modified axion term Lϑ discussed in the Introduction. As usual, the electric and magnetic fields in (8) are written in terms of the potentials according to (7).

The Eqs. (4) and (5) explicitly show that there are no modifications to the dynamics in the bulk (z ≠ 0) with respect to standard electrodynamics. Nevertheless, as it is well known, the solution of a system of differential equations depends crucially upon de boundary conditions. In this way, the new physics induced by Lϑ arises from the interface between the media (z = 0) and will be a consequence of the boundary conditions there. Physically, this is a consequence that TIs behave as normal insulators in the bulk, but possess conducting properties at interfaces, as indicated by the MEE. Even though we are dealing with a continuous dielectric, (ϵ1=ϵ2), the different MEP of both media generate effective transmission and reflection coefficients for electromagnetic waves across the interface. Mathematically, this feature is understood because E⋅B in Lϑ is a total derivative, so the only allowed modifications to the standard Maxwell equations arise when the integration by parts produces ∂αϑ≠0, which precisely define the interface in our problem.

Assuming that the time derivatives of the fields are finite in the vicinity of the surface z = 0, the modified Maxwell equations (4) and (5) imply the following boundary conditions (BCs)

for vanishing external sources on the surface z = 0. The notation is Vz=0-z=0+=Vz=0+-Vz=0-, V|z=0=V(z=0), where z=0± indicates the limits z=0±η, with η→0, respectively. The continuous terms in the right-hand side of the first and second equations in Eq. (9) represent self-induced surface charge and surface current densities, respectively, which clearly demonstrate the MEE localized just at the interface between the two media.

A convenient way to deal with the fields produced by arbitrary sources in electrodynamics, in particular in ϑ-ED, is by using the corresponding GF G   νμ(x,x'), which we briefly revise below ^[38-42]. Before going into the details we comment upon the advantages provided by the use of GF methods over different alternatives in electrodynamics: the knowledge of the GF of a given physical system allows a direct calculation of the corresponding electromagnetic fields for arbitrary sources either analytically or numerically just by direct substitution. This clearly avoids the guesswork required when using the image method, which by the way works only in highly symmetrical cases. Also, it saves a lot of work when one needs to consider different sources in a given system by avoiding to solve the equations for each source. This very useful technique extends to many branches of physics like scattering theory, condensed matter physics and quantum field theory, for example.

In what follows we restrict ourselves to contributions of free sources Jμ=(ϱ,j) located outside the interface, and to systems without BCs imposed on additional surfaces, except for those at infinity. A compact formulation of the problem is given in terms of the potentials (7) expressed in their four dimensional form (Φ,A) together with the GF

which satisfies the equation

in the modified Lorenz gauge ϵ∂Φ/∂t+∇⋅A=0, together with the appropriate BCs. The operator ◊υμ is

The detailed calculation of the GF is reported in Sec. 3 of Ref. ^[42]. Here we only recall the results that are written in terms of G-    νμ, which differs from G    νμ only in the term G    00=G-    00/ϵ. Since the GF is time-translational-invariant it is convenient to introduce the corresponding Fourier transform such that

The final result is presented as the sum of three terms, G-    νμx,x';ω=G-ED  νμx,x';ω+G-θ~  νμx,x';ω+G-θ~2  νμx,x';ω, whose explicitly form are

Here R⊥=(x-x')⊥=x-x',y-y', k⊥=(kx,ky) is the momentum parallel to the interface, kα=(ω,k⊥) and k~0=nω where n=ϵ is the refraction index. We observe that in the static limit (ω = 0), the result (14) reduces to the one reported in Ref. ^[38].

3 Electric Dipole perpendicular to the interface

In this section we determine the electric field of an oscillating point-like electric dipole p=p z^ located at a distance z₀ > 0 on the z-axis and perpendicular to the interface. We restrict ourselves to the far-field approximation (k~0r≫1) starting from the GF given by Eqs. (14).

3.1 The Electromagnetic Potential A^μ

The charge and current density for this dipole are

ϱ(x';ω)=-pδ(x')δ(y')δ'(z'-z0),j(x';ω)=-iωpδ(x')δ(y')δ(z'-z0)z^, (15)

respectively, where δ'(u)=dδ(u)/du. After convoluting the sources (15) with the GF (14) we find the following components of A^μ

A0(x;ω)=-pn2ik~0cosθeik~0(r-z0cosθ)r-1n2θ~2p4n2+θ~2H(x,z0;ω)  , (16)

Aa(x;ω)=-2θ~p4n2+θ~2iεabxbρ∂∂ρI(x,z0;ω)+θ~2p4n2+θ~2iωxaρ∂∂ρJ(x,z0;ω)  , (17)

A3(x;ω)=-iωpeik~0(r-z0cosθ)r  , (18)

where ρ=∥x⊥∥=x2+y2, r=x2+y2+z2, a,b∈1,2, εab=-εba, ε12=+1, and {xa}={x1,x2} with x1=x, x2=y. We also have the functions

H(x,z0;ω)=0∞k⊥3dk⊥k~02-k⊥2J0k⊥ρeik~02-k⊥2(|z|+z0)  , (19)

I(x,z0;ω)=0∞k⊥dk⊥k~02-k⊥2J0k⊥ρeik~02-k⊥2(|z|+z0)  , (20)

J(x,z0;ω)=0∞k⊥dk⊥k~02-k⊥2J0k⊥ρeik~02-k⊥2(|z|+z0)  . (21)

The derivation of the above results can be found in the Appendix A.

The next step is to calculate the integrals (19)-(21) in the far-field approximation. To begin with, we recap the main ingredients of the calculation carried out in full detail for the function H in Appendix B. We employ the steepest descent method ^[44,45,46] and incorporate some modifications based on Refs. ^[13,47]. These modifications are required because the current integrals (19)-(21) have poles coinciding with their stationary point (k⊥)s=k~0sinθ at θ=π/2, as can be seen in Eq. (13) of the Appendix B. This means that the factor of the exponential is not a smooth function around the stationary point now, which will prevent a direct application of the method. The main idea to overcome this difficulty is to subtract and add the conflicting pole, as shown in Eqs. (B.15) and (B.25) of the Appendix B. Thereby, we obtain two integrals: one with the divergence removed in the vicinity of the stationary point and another containing the singularity, which can be directly evaluated. The first integral leads to the ordinary stationary phase contributions and the second one gives contributions that are identified as axially symmetric cylindrical waves ^[13]. The final results for the integrals H,I and J, obtained in full detail in the Appendix B, are

H(x,z0;ω)=k~0eik~0rirsin2θeik~0z0|cosθ||cosθ|-12sinθ-sin2θ+2πik~0rsinθk~02ieik~0rsinθπ2erfc-iik~0r1-sinθ  , (22)

J(x,z0;ω)=eik~0rik~0reik~0z0|cosθ||cosθ|-12sinθ-sin2θ+2πik~0rsinθπeik~0rsinθ2ierfc-iik~0r1-sinθ  , (23)

I(x,z0;ω)=eik~0rireik~0z0|cosθ|  , (24)

where erfc(z) denotes the complementary error function. The contributions in curly brackets arise from the terms including the subtracted pole, they are are well behaved at θ=0,π/2,π and describe the standard spherical waves.

On the other hand, the terms from the erfc function in square brackets arise from the pole itself and describe wave propagation corresponding to the amplitude 1/(rsinθ)12expik~0rsinθ×exp-ik~0t=1z12expik~0z×exp-ik~0t. This is clearly a solution of the Helmholtz equation in cylindrical coordinates in the far-zone, thus giving the name of cylindrical waves to this contribution.

The crucial point is that the amplitude of the cylindrical waves is modulated by the erfc function, with argument

iik~0r1-sinθ=iΛ≡ii s, (25)

which will provide us with a quantitative way of appraising the relevance of the cylindrical waves. To this end, we now discuss the behavior of the function

erfc(iΛ)=1iπeik~0r(1-sinθ)Z(eiπ/4s), (26)

that we rewrite in terms of the the Faddeeva plasma dispersion function Z(Λ) discussed in some detail in the Appendix C. Up to an irrelevant normalization constant, we define the function

F(s)=|Z(eiπ/4s)|2/(2π)  , (27)

which controls the amplitude of the electromagnetic fields of the cylindrical waves and can be readily calculated from Eq. (C.7). As can be seen from the following numerical values F(0) = 0.5, F(1) = 0.112, F(3) = 0.174, … and F(s≫1)→0,F(s) is a rapidly decreasing function of s. In this way, when s≫1 the cylindrical wave contribution can be neglected, whereas for s→0 the function F(s) contributes maximally and the cylindrical waves should be taken into account.

If we recall that ρ=rsinθ, s2 can be rewritten as s2=k~0(r-ρ) and can be interpreted as a measure of how far the observer with coordinates (r,θ,ϕ) is from the interface. In other words, s2 determines how far is the spherical radius r from the cylindrical radius ρ at the observation point. This property enables us to define what we call the discarding angle θ₀, which provides a condition to estimate when we can neglect or not the cylindrical wave contribution, according to the magnitude of F(s). We proceed as follows.

For a given observation distance r₀, such that k~0r0 is large enough to describe the far-field regime, we choose as an arbitrary cutoff point the value s = s₀. At the cutoff we define the discarding angle θ₀ such that

s0=k~0r0(1-sinθ0). (28)

More precisely, we can distinguish the upper hemisphere (UH) θ∈[0,π2] from the lower hemisphere (LH) θ∈[π/2,π] by writing

θ0UH=arcsin1-s02k~0r0,θ0LH=π2+arccos1-s02k~0r0, (29)

respectively. Let us observe that both discarding angles are very close to the interface (θ=π/2) in the far-field regime.

For a given observation distance r₀, the rapidly decreasing behavior of F(s) allows us to adopt the following criterion for estimating the relative weight of the spherical versus the cylindrical waves:

Fors>s0,(0<θ<θ0UHandθ0LH<θ<π)cylindrical waves are neglectedFors<s0,(θ0UH<θ<θ0LH)cylindrical waves are taken into account. (30)

In our case we take s₀ = 1,the function F(s) decreased about five times with respect to its maximum at s = 0. Let us emphasize that s₀ can be arbitrarily chosen much larger than one, which will provide a more stringent cutoff.

In other words, for a fixed r₀ and a chosen s₀ = 1 the discarding angles θ₀ define two regions, as shown in Fig. 1. The region V1, where θ∈0,θ0UH∪θ0LH,π, is such that s > 1, while in its complement, the region V2 where θ∈[θ0UH,θ0LH], we have s > 1.

Figure 1 Diagram showing the regions V1 and V2 determined through the discarding angles θ0UH and θ0LH. The region V1, where only the spherical waves (undulated pink arrows) are taken into account is defined by θ∈0,θ0UH∪θ0LH,π. The region V2 (white hatched region), defined by θ∈θ0UH,θ0LH, contains the intermediate region (s<s0=1) where both cylindrical waves (undulated black waves) and spherical waves have to be taken into account. Nevertheless, going further into the radiation zone for each observation point in the intermediate region one can make s>s0=1, thus making the cylindrical waves unobservable. 

To fix ideas let us consider now the UH and examine what happens when we fix the angle and explore the consequences changing the observation distance to a larger value r>>r0,i.e., we go farther into the radiation zone. Suppose that in the region V1 we consider the angle θ1<θ0UH, where we have s1=k~0r0(1-sinθ1)>1 according to our choices. Then, keeping θ₁ fixed and going to a larger distance r>>r0 would only increase the value of s>=k~0r>(1-sinθ1) such that s>>s1>s0=1. That is to say, all observation points in V1 with r>r0 will have s>s0=1 and the cylidrical waves will not be relevant there.

On the other hand, the region V2 shows a mixed behavior. Again, let us consider an angle θ2>θ0UH where we have s2=k~0r0(1-sinθ2)<1 by construction. Nevertheless, an increase in the observation distance to r> can revert the situation yielding a value s>>1. To this end it is enough to take r>>r0/s2. That is to say, the region V2 contains the intermediate region where the cylindrical waves are relevant, but going further into the radiation zone these waves can be safely neglected. A similar situation occurs in the LH, which we do not discuss in detail here.

From the function erfc-iik~0r(1-sinθ) appearing in the Eqs. (22) and (23) we identify the analogous of the Sommerfeld numerical distance S in the standard dipole radiation, which is determined by rewriting the complementary error function as erfc (-iS)^[13]. Thus we have

S=ik~0r(1-sinθ)=is2, (31)

which varies with the angle θ for a fixed observation distance r. This quantity is closely related to the discarding angle θ₀ according to

|S|=1-sinθ1-sinθ0. (32)

Going back to the calculation of the electromagnetic field in the far-field approximation, we now write the corresponding expressions for the electromagnetic potentials obtained from plugging the results (22-24) into the (16-18). In the case of the region V1 we have

A0(x;ω)=-pn2ik~0cosθeik~0(r-z0cosθ)r-1n2θ~2p4n2+θ~2k~0sin2θ|cosθ|eik~0(r+z0|cosθ|)ir  , (33)

Aa(x;ω)=p4n2+θ~2-2θ~ik~0εabxbr+θ~2iω|cosθ|xareik~0(r+z0|cosθ|)r  , (34)

A3(x;ω)=-iωpeik~0(r-z0cosθ)r, (35)

where we dropped terms of higher order. On the other hand, in the intermediate region V2, (s > 1) the contribution of the cylindrical waves near the interface is apparent. The electromagnetic potential now is

A0(x;ω)=∓pn2ik~0ξeik~0(r∓z0ξ)r-1n2θ~2p4n2+θ~2k~0eik~0ririk~0z0-ξ8-1n2θ~2p4n2+θ~2k~022πik~0reik~0rπ2i+πik~0r2ξ  , (36)

Aa(x;ω)=p4n2+θ~2-2iθ~k~0εabxbρeik~0(r+z0ξ)r+θ~2iωxaρeik~0rrik~0z0-ξ8-θ~2ωk~0xaρ2πik~0reik~0rπ2i+πik~0r2ξ  , (37)

A3(x;ω)=-iωpeik~0(r∓z0ξ)r, (38)

where we again dropped terms of higher order. The minus (plus) sign in the first term of the right-hand side in Eq. (36) corresponds to the UH (LH), respectively. We have also performed a first order power expansion around π/2 in the complementary error function arising from Eqs. (22) and (23) in terms of the variable ξ<1 given by

ξUH≡π/2-θUH,ξLH=θLH+π/2, (39)

for the UH and LH, respectively. Let us observe that for both hemispheres we have ξ>0 and also that (1-sinθ)=ξ2/2. In analogous way it is convenient to introduce the corresponding variable ξ0 related to the discarding angle θ0 just by replacing ξ→ξ0 and θ→θ0 in Eq. (39). Using also Eq. (29) we obtain ξ0=2/(nωr0). In this way |S|=ξ2/ξ02<1 for both hemispheres. The expansion in powers of ξ is only valid in the intermediate zone where we have ξ≪ξ0.