# Radiation of a short linear antenna above a topologically insulating half-
space

M. Ibarra-Meneses<EMAIL_ADDRESS>Instituto de Ciencias Nucleares,
Universidad Nacional Autónoma de México, 04510 Ciudad de México, México  A.
Martín-Ruiz<EMAIL_ADDRESS>Instituto de Ciencias Nucleares,
Universidad Nacional Autónoma de México, 04510 Ciudad de México, México

###### Abstract

The topological magnetoelectric effect (TME) is a unique macroscopic
manifestation of quantum states of matter possessing topological order and it
is described by axion electrodynamics. In three-dimensional topological
insulators, for instance, the axion coupling is of the order of the fine
structure constant, and hence a perturbative analysis of the field equations
is plenty justified. In this paper we use Green’s function techniques to
obtain time-dependent solutions to the axion field equations in the presence
of a planar domain-wall separating two media with different topological order.
We apply our results to investigate the radiation of a short linear antenna
near the domain-wall.

##  I Introduction

Topological materials are not distinguished by local ordered parameters, but
by topological quantum numbers which could be described and computed using
topological field theories in the low energy limit. The best studied of these
are the topological insulators (TIs), where the corresponding bulk band gap is
similar to an ordinary insulator, but the boundary of the material hosts
gapless states which are protected by the topological numbers of the bulk
geometry Qi and Zhang (2011). The experimental discovery of TIs in HgTe
quantum wells König  _et al._ (2007) and stoichiometric crystals Chen  _et
al._ (2009); Xia  _et al._ (2009); Sato  _et al._ (2010) encourages further
investigations regarding its exotic properties.

Apart from their spectroscopic distinguishing features, TIs have
unconventional electromagnetic response described by axion electrodynamics,
which is defined by the standard Maxwell Lagrangian plus the additional
$\theta$-term (in SI units)
$\mathcal{L}_{\theta}=(\alpha/\pi)\sqrt{\epsilon_{0}/\mu_{0}}\,\theta\,{\bf{E}}\cdot{\bf{B}}$,
where ${\bf{E}}$ and ${\bf{B}}$ are the electromagnetic fields, $\alpha$ is
the fine structure constant, and $\theta$ is an angular variable known in
particle physics as the axion field or axion angle Wilczek (1987). The axion
field was originally proposed to solve the strong-CP problem in QCD Peccei and
Quinn (1977); Weinberg (1978); Wilczek (1978), but its weak interaction
strength with standard model particles and an elegant production mechanism in
the early universe make it a promising candidate to explain dark matter as
well Preskill  _et al._ (1983); Dine and Fischler (1983); Abbott and Sikivie
(1983). In the context of TIs, time-reversal (TR) symmetry indicates that
$\theta=0,\pi$ (mod $2\pi$), and hence the $\theta$-term has no effect on
Maxwell equations in the bulk. The nontrivial topological property, a surface
half-integer quantum Hall effect, manifests only at the boundary between two
manifolds with different topologies, in which case $\theta$ is not globally
constant. Therefore, from the macroscopic point of view, the axion angle
$\theta$ can be viewed as an effective parameter characterizing the media, in
the same footing as its permittivity $\epsilon$ and permeability $\mu$.

The $\theta$-term encodes the most salient feature of the electromagnetic
response of TIs: the topological magnetoelectric effect, where an electric
field can induce a magnetic polarization and a magnetic field can induce an
electric polarization Fiebig (2005). This is why, in the condensed matter
literature, the axion angle is termed the topological magnetoelectric
polarization (TMEP). A number of topological magnetoelectric effects have been
predicted on the basis of this theory, namely, the appearance of charge and
current densities of magnetic monopoles Qi  _et al._ (2009); Karch (2009);
Zang and Nagaosa (2010); Martín-Ruiz  _et al._ (2015, 2019, 2016); Martín-Ruiz
(2018); Meier  _et al._ (2019), the magneto-optical Faraday and Kerr rotation
Chang and Yang (2009); Tse and MacDonald (2011); Crosse (2016); Wu  _et al._
(2016); Filipini and Cambiaso (2023), the Casimir effect Grushin and Cortijo
(2011); Grushin  _et al._ (2011); Rodriguez-Lopez and Grushin (2014); Chen and
Wan (2012); Rodriguez-Lopez (2011); Nie  _et al._ (2013); Martín-Ruiz  _et
al._ (2016), as well as the influence of TIs on atomic systems nearby Song
_et al._ (2014); Fang  _et al._ (2015); Martín-Ruiz and Urrutia (2018);
Martín-Ruiz and Chan-López (2017); Castro-Enriquez  _et al._ (2020); Castro-
Enriquez and Martín-Ruiz (2020).

Time-independent solutions to the axion field equations have been elaborated
by using the electric-magnetic duality of $\mathbb{Z}_{2}$ TIs Karch (2009) as
well as Green’s function techniques Martín-Ruiz  _et al._ (2015, 2016).
However, general solutions to the time-dependent case are still under
investigation Franca  _et al._ (2019); Franca and Urrutia (2022); Franca and
Buhmann (2022). In this work we aim to fill in this gap by computing, in a
perturbative fashion, solutions to the time-dependent axion field equations in
the presence of a planar domain wall separating two media with different
topologies. The validity of a perturbative scheme for the solution of the
field equations is as follows. In the high-energy context the validity is
plenty justified: the axion-to-two- photon coupling $g_{\theta\gamma\gamma}$
which converts photons into axions is expected to be very small Ouellet and
Bogorad (2019). In the condensed matter framework, the smallness of the
$\theta$-term is controlled by the coupling $\alpha\theta$, which is small
since $\alpha\approx 1/137$ and $\theta=\pi$ for a TI with a single Dirac
fermion on its surface. In the hypothetical situation of a TI with a large
number of Dirac fermions on its surface, our perturbative approach ceases to
be valid as $\alpha\theta$ approaches unity. In view that exact solutions to
the static axion field equations have been reported in the literature, we
shall first validate our perturbative scheme by reproducing well-known results
in this case. Next we will safely extend our analysis to the time-dependent
case. Finally, we apply our results to study the radiation of a short linear
antenna near the domain-wall.

The paper is organized as follows. In Sec. II we briefly review the basics of
axion electrodynamics. In Sec. III we solve, in a perturbative fashion, the
static field equations and verify its consistency with previously reported
results for the Green’s function for two ponderable topological insulators
separated by a planar interface. Section IV is devoted to the evaluation of
the time-dependent Green’s function by assuming equal optical properties in
the two media. This configuration allow us to isolate the topological effects.
The details of technical calculations are relegated to the appendices.
Finally, in Sec. V we apply our result to investigate the radiation of a short
linear antenna near a $\theta$ domain-wall.

##  II Topological field theory

The effective field theory governing the electromagnetic response of
topological insulators, independently of the microscopic details, is defined
by the action (in SI units) Qi  _et al._ (2008):

$\displaystyle S=\int
d^{4}x\,\left[\,\frac{1}{2}\left(\epsilon{\bf{E}}^{2}-\frac{1}{\mu}{\bf{B}}^{2}\right)+\frac{\alpha}{\pi}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\theta\,{\bf{E}}\cdot{\bf{B}}\,\right],$
(1)

where $\alpha=e^{2}/(4\pi\epsilon_{0}\hbar c)\simeq 1/137$ is the fine
structure constant, $\epsilon$ and $\mu$ are the permittivity and permeability
of the material, respectively, and $\theta$ is the topological magnetoelectric
polarizability (axion field). Under periodic boundary conditions, the
partition function and all physical quantities are invariant under shifts of
$\theta$ by any multiple of $2\pi$. Since ${\bf{E}}\cdot{\bf{B}}$ is odd under
TR, the only symmetry allowed values of $\theta$ are $0$ or $\pi$ (modulo
$2\pi$). This leads to the $\mathbb{Z}_{2}$ classification of three-
dimensional TR invariant TIs. Clearly, the $\theta$-term in Eq. (1) has no
effect on Maxwell equations in the bulk, since it is a total derivative term.
The nontrivial topological property, a surface half-integer quantum Hall
effect, manifests only when a TR-breaking perturbation is induced on the
surface to gap the surface states, thereby converting the material into a full
insulator. This can be achieved by introducing magnetic dopants to the surface
Zhang  _et al._ (2013) or by the application of an external static magnetic
field Maciejko  _et al._ (2010). In this situation, $\theta$ is quantized in
odd integer values of $\pi$, where the magnitude and sign of the multiple is
controlled by the strength and direction of the TR-breaking perturbation.

The electromagnetic response of a system in the presence of the axionic term
is still described by the ordinary Maxwell equations, but the constituent
relations which define the electric displacement ${\bf{D}}$ and the magnetic
field ${\bf{H}}$ acquire an extra term proportional to $\theta$, i.e.
${\bf{D}}=\epsilon{\bf{E}}+\alpha(\theta/\pi)\,\sqrt{\epsilon_{0}/\mu_{0}}{\bf{B}}$
and
${\bf{H}}=\mu^{-1}{\bf{B}}-\alpha(\theta/\pi)\,\sqrt{\epsilon_{0}/\mu_{0}}\,{\bf{E}}$.
Here, the $\theta$-dependent terms, which arise from the axionic term in Eq.
(1), describe the topological magnetoelectric effect. The frequency components
of the electromagnetic fields (introduced as usual: by making a Fourier
decomposition in time of the sources and fields) obey the inhomogeneous
equations

$\displaystyle\nabla\cdot\left[\epsilon({\bf{r}},\omega)\,{\bf{E}}({\bf{r}},\omega)\right]=\rho({\bf{r}},\omega)-\frac{\alpha}{\pi}\,\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\nabla\theta({\bf{r}})\cdot{\bf{B}}({\bf{r}},\omega),$
(2)
$\displaystyle\nabla\times\left[\mu^{-1}({\bf{r}},\omega)\,{\bf{B}}({\bf{r}},\omega)\right]+i\omega\epsilon({\bf{r}},\omega){\bf{E}}({\bf{r}},\omega)={\bf{J}}({\bf{r}},\omega)+\frac{\alpha}{\pi}\,\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\nabla\theta({\bf{r}})\times{\bf{E}}({\bf{r}},\omega),$
(3)

while the homogeneous equations are untouched due to gauge invariance of the
action (1):

$\displaystyle\nabla\cdot{\bf{B}}({\bf{r}},\omega)=0,$ (4)
$\displaystyle\nabla\times{\bf{E}}({\bf{r}},\omega)=i\omega{\bf{B}}({\bf{r}},\omega).$
(5)

One can further verify that the divergence equations (2) and (4) are subsumed
within the curl equations (3) and (5), which can be verified by taking the
divergence of the latters. We observe that equations (2)-(3) depend only on
the space gradient of the axion field. So, the $\theta$-term does not modify
the field equations in the bulk of topological insulators. However, for a TI
in contact with a conventional insulator, as shown in Fig. 1, the axion field
is piecewise constant across the interface $\Sigma$ and then
$\nabla\theta({\bf{r}})=(\theta_{1}-\theta_{2})\delta(\Sigma)\hat{{\bf{n}}}$,
where $\hat{{\bf{n}}}$ is the outward unit normal to $\Sigma$. These are the
equations we aim to solve in a perturbative fashion along this paper. Our
general approach will be to first Taylor expand the electromagnetic potentials
in powers of the fine structure constant and then solve recursively the
resulting equations. Before embarking to the solution of the time-dependent
field equations, we first tackle the time-independent case and demonstrate
consistency of our results with previously reported ones.

Figure 1: Geometry of two semi-infinite topologically insulating media
separated by the surface $\Sigma$. The media have different magnetoelectric
polarizabilities, permittivities and permeabilities, as depicted.

##  III Solution to the static axion field equations

Let us consider the geometry shown in Fig. 1 with the interface $\Sigma$ being
the plane $z=0$. The half-space $z<0$ is occupied by a topological medium with
dielectric constant $\epsilon_{2}$, magnetic permeability $\mu_{2}$ and axion
angle $\theta_{2}$, whereas the half-space $z>0$ is occupied by a medium with
dielectric constant $\epsilon_{1}$, magnetic permeability $\mu_{1}$ and axion
angle $\theta_{1}$. Then the inhomogeneity in the permittivity, permeability
and the axion field are limited to a finite discontinuity across the planar
surface $z=0$.

To proceed further we write the field equations (2) and (3) in terms of the
electrostatic potential $\phi({\bf{r}})$ and vector potential
${\bf{A}}({\bf{r}})$. In the Coulomb gauge, the static electromagnetic
potentials satisfy

$\displaystyle-\nabla\cdot\left[\epsilon({\bf{r}})\,\nabla\phi({\bf{r}})\right]=\rho({\bf{r}})-\tilde{\alpha}\,\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\delta(z)\,\hat{{\bf{e}}}_{z}\cdot\nabla\times{\bf{A}}({\bf{r}}),$
(6)
$\displaystyle\nabla\times\left[\mu^{-1}({\bf{r}})\,\nabla\times{\bf{A}}({\bf{r}})\right]={\bf{J}}({\bf{r}})-\tilde{\alpha}\,\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\delta(z)\,\hat{{\bf{e}}}_{z}\times\nabla\phi({\bf{r}}),$
(7)

where $\tilde{\alpha}=\alpha(\theta_{1}-\theta_{2})/\pi$. To obtain a
perturbative solution for the electromagnetic potentials in the presence of
arbitrary external sources, we introduce the expansions
$\phi({\bf{r}})=\phi^{(0)}({\bf{r}})+\sum_{n=1}^{\infty}\phi^{(n)}({\bf{r}})$
and
${\bf{A}}({\bf{r}})={\bf{A}}^{(0)}({\bf{r}})+\sum_{n=1}^{\infty}{\bf{A}}^{(n)}({\bf{r}})$,
where $\phi^{(n)}$ and ${\bf{A}}^{(n)}({\bf{r}})$ are understood to vanish as
$\tilde{\alpha}^{n}$. We now substitute the power expansion into equations (6)
and (7), and comparing coefficients of each power of $\tilde{\alpha}$ results
in an infinite series of simultaneous equations. The zeroth-order equations
are simply the Maxwell equations in the absence of the $\theta$-term,

$\displaystyle-\nabla\cdot\left[\epsilon({\bf{r}})\,\nabla\phi^{(0)}({\bf{r}})\right]=\rho({\bf{r}}),\qquad\qquad\nabla\times\left[\mu^{-1}({\bf{r}})\,\nabla\times{\bf{A}}^{(0)}({\bf{r}})\right]={\bf{J}}({\bf{r}}).$
(8)

The $n$-order equations, which are the main concern in this work, are

$\displaystyle\nabla\cdot\left[\epsilon({\bf{r}})\,\nabla\phi^{(n)}({\bf{r}})\right]=\tilde{\alpha}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\delta(z)\,\hat{{\bf{e}}}_{z}\cdot\nabla\times{\bf{A}}^{(n-1)}({\bf{r}}),\quad\nabla\times\left[\mu^{-1}({\bf{r}})\,\nabla\times{\bf{A}}^{(n)}({\bf{r}})\right]=-\tilde{\alpha}\,\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\delta(z)\,\hat{{\bf{e}}}_{z}\times\nabla\phi^{(n-1)}({\bf{r}}).$
(9)

The process of solving these equations begins with the solution to the zeroth-
order equations (8), which involves the well-known electro- and magnetostatic
Green’s functions. We now adopt the following notation. Since translational
invariance is broken in the $z$ direction due to the presence of the
interface, it is natural decompose any vector ${\bf{C}}$ into components
parallel and perpendicular to the normal direction $\hat{{\bf{e}}}_{z}$, i.e.
${\bf{C}}_{z}=C_{z}\hat{{\bf{e}}}_{z}$ and
${\bf{C}}_{\perp}=C_{x}\hat{{\bf{e}}}_{x}+C_{y}\hat{{\bf{e}}}_{y}$,
respectively. So, the general solution for the electrostatic potential is

$\displaystyle\phi^{(0)}({\bf{r}})=\int
G_{\epsilon}({\bf{r}},{\bf{r}}^{\prime})\,\rho({\bf{r}}^{\prime})\,d^{3}{\bf{r}}^{\prime},$
(10)

where the electric Green’s function is given by

$\displaystyle
G_{\epsilon}({\bf{r}},{\bf{r}}^{\prime})=\frac{1}{4\pi\epsilon_{1}}\left[\frac{1}{\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}{}_{\\!\\!\perp}|^{2}+(z-z^{\prime})^{2}}}+\frac{\epsilon_{1}-\epsilon_{2}}{\epsilon_{1}+\epsilon_{2}}\frac{1}{\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}{}_{\\!\\!\perp}|^{2}+(|z|+z^{\prime})^{2}}}\right]$
(11)

in which ${\bf{r}}_{\perp}=x\hat{{\bf{e}}}_{x}+y\hat{{\bf{e}}}_{y}$ and
${\bf{r}}^{\prime}{}_{\\!\\!\perp}=x^{\prime}\hat{{\bf{e}}}_{x}+y^{\prime}\hat{{\bf{e}}}_{y}$
J. Schwinger and Tsai (1998). The solution to the magnetostatic problem is a
bit more subtle. Decomposing the current density and the vector potential into
their components parallel and perpendicular to $\hat{{\bf{e}}}_{z}$, i.e.
${\bf{J}}={\bf{J}}_{\perp}+{\bf{J}}_{z}$ and
${\bf{A}}={\bf{A}}_{\perp}+{\bf{A}}_{z}$, respectively, the general solution
for the vector potential, in the half-space $z>0$, can be written as

$\displaystyle{\bf{A}}^{(0)}_{\xi}({\bf{r}})=\int
G_{\mu}^{(\xi)}({\bf{r}},{\bf{r}}^{\prime})\,{\bf{J}}_{\xi}({\bf{r}}^{\prime})\,d^{3}{\bf{r}}^{\prime},\qquad\xi=\,\perp,z$
(12)

while for the region $z<0$ the solution is simply
${\bf{A}}^{(0)}({\bf{r}})=\int
G_{\mu}^{(\perp)}({\bf{r}},{\bf{r}}^{\prime})\,{\bf{J}}({\bf{r}}^{\prime})\,d^{3}{\bf{r}}^{\prime}$.
The corresponding magnetic Green’s function is given by J. Schwinger and Tsai
(1998)

$\displaystyle G_{\mu}^{(\xi)}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\frac{\mu_{1}}{4\pi}\left[\frac{1}{\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}{}_{\\!\\!\perp}|^{2}+(z-z^{\prime})^{2}}}+\chi\,\frac{\mu_{2}-\mu_{1}}{\mu_{2}+\mu_{1}}\frac{1}{\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}{}_{\\!\\!\perp}|^{2}+(|z|+z^{\prime})^{2}}}\right],$
(13)

where $\chi=+1\,(-1)$ for $\xi=\,\perp(z)$.

With all of the above we are ready to solve equations (9). This is a simple
but not straightforward task, so we relegate the details to Appendix A. The
final results, up to second order in the parameter $\tilde{\alpha}$, are

$\displaystyle\phi({\bf{r}})=\phi^{(0)}({\bf{r}})+\int\left[G_{0i}^{(1)}({\bf{r}},{\bf{r}}^{\prime})\,J_{i}({\bf{r}}^{\prime})+G_{00}^{(2)}({\bf{r}},{\bf{r}}^{\prime})\,\rho({\bf{r}}^{\prime})\right]d^{3}{\bf{r}}^{\prime},$
(14) $\displaystyle
A_{i}({\bf{r}})=A_{i}^{(0)}({\bf{r}})+\int\left[G_{0i}^{(1)}({\bf{r}},{\bf{r}}^{\prime})\,\rho({\bf{r}}^{\prime})+G_{ij}^{(2)}({\bf{r}},{\bf{r}}^{\prime})\,J_{j}({\bf{r}}^{\prime})\right]d^{3}{\bf{r}}^{\prime},$
(15)

where the indices $i,j$ range over the transverse coordinates $x,y$, and

$\displaystyle G_{0i}^{(1)}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\frac{\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}}{2\pi(\epsilon_{1}+\epsilon_{2})(1/\mu_{1}+1/\mu_{2})}\,\frac{\epsilon_{zji}\,({\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime})_{j}}{|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|^{2}}\,\left[1-\frac{|z|+z^{\prime}}{\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|^{2}+(|z|+z^{\prime})^{2}}}\right],$
(16) $\displaystyle G_{00}^{(2)}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\frac{(\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}\,)^{2}}{2\pi(\epsilon_{1}+\epsilon_{2})^{2}(1/\mu_{1}+1/\mu_{2})}\,\frac{1}{\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}{}_{\\!\\!\perp}|^{2}+(|z|+z^{\prime})^{2}}},$
(17) $\displaystyle G_{ij}^{(2)}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\frac{(\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}\,)^{2}}{2\pi(\epsilon_{1}+\epsilon_{2})(1/\mu_{1}+1/\mu_{2})^{2}}\;\left[\frac{\delta_{ij}}{\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}{}_{\\!\\!\perp}|^{2}+(|z|+z^{\prime})^{2}}}+\partial_{i}K_{j}({\bf{r}},{\bf{r}}^{\prime})\right],$
(18)

with the definition

$\displaystyle
K_{j}({\bf{r}},{\bf{r}}^{\prime})=\frac{({\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime})_{j}}{|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|^{2}}\,\left(\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|^{2}+(|z|+z^{\prime})^{2}}-(|z|+z^{\prime})\right).$
(19)

Clearly, the electric (vector) potential is even (odd) under time-reversal
operation, as expected.

### Image magnetic monopole effect

As a consistency check of our results, we shall verify the image magnetic
monopole effect, which has been extensively discussed in the literature by
using different approaches, for example, the method of images Qi  _et al._
(2009), the electric-magnetic duality of TIs Karch (2009) and Green’s function
techniques Martín-Ruiz  _et al._ (2016); Martín-Ruiz (2018). This elusive
effect consists in the appearance of an image magnetic monopole when an
electric charge is placed at a certain distance from the TI’s surface.

We now consider for the charge density a pointlike electric charge $q$ at
$z=z_{0}$. Equations (14) and (15), together with the Green’s functions
(16)-(18), give the electric and vector potentials. For $z>0$, the electric
potential can be interpreted as due to two pointlike electric charges, one of
strength $q$ at $z=z_{0}$, and the other, the image charge, of strength

$\displaystyle
q^{\prime}=q\left[\frac{\epsilon_{1}-\epsilon_{2}}{\epsilon_{1}+\epsilon_{2}}+\frac{2\epsilon_{1}(\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}\,)^{2}}{(\epsilon_{1}+\epsilon_{2})^{2}(1/\mu_{1}+1/\mu_{2})}\right]$
(20)

at $-z_{0}$. For $z<0$ only one pointlike electric charge appears, of strength
$q^{\prime}+q$ located at $z_{0}$. Unremarkably, the above expression reduce
to the standard textbook one in the nontopological limit ($\tilde{\alpha}=0$)
J. Schwinger and Tsai (1998). Due to the magnetoelectricy of TIs, a magnetic
field is induced. The corresponding vector potential is found to be

$\displaystyle{\bf{A}}({\bf{r}})=\frac{g}{4\pi}\,\frac{x\hat{{\bf{e}}}_{y}-y\hat{{\bf{e}}}_{x}}{r_{\perp}^{2}}\,\left[1-\frac{|z|+z_{0}}{\sqrt{r_{\perp}^{2}+(|z|+z_{0})^{2}}}\right],$
(21)

which we identify as the vector potential of a straight vortex line or Dirac
string over the $z$ axis. For $z>0$ ($z<0$) the string is over the negative
(positive) $z$ axis ending at a magnetic monopole at $-z_{0}$ ($z_{0}$). In
other words, the magnetic field for $z>0$ can be interpreted as that of a
magnetic monopole of strength

$\displaystyle
g=\frac{2q\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}}{(\epsilon_{1}+\epsilon_{2})(1/\mu_{1}+1/\mu_{2})}$
(22)

at the image point $-z_{0}$. Similarly, for $z<0$ the magnetic field can be
understood as originated by a monopole of strength $-g$ located at $z_{0}$.
Our results are in a agreement with the ones reported in Refs. Qi  _et al._
(2009); Karch (2009); Martín-Ruiz  _et al._ (2016); Martín-Ruiz (2018). This
can be further verified by power expanding the image electric and magnetic
charges up to order $\tilde{\alpha}^{2}$, thus yielding our results (20) and
(22), respectively.

##  IV Time-dependent solutions to the axion field equations

In this Section we consider the same configuration as in Sec. III, but we
shall investigate now time-dependent solutions to the field equations (2)-(5).
In order to isolate the contribution arising from the surface Hall effect, we
restrict ourselves to the case in which the permittivities are the same, i.e.
$\epsilon_{1}=\epsilon_{2}\equiv\epsilon$, but the topological magnetoelectric
polarizability is a piecewise constant function: $\theta_{1}$ for $z>0$ and
$\theta_{2}$ for $z<0$. We assume the sources located in the region $z>0$.

We now introduce the field-potential relationships in the frequency domain,
${\bf{B}}({\bf{r}},\omega)=\nabla\times{\bf{A}}({\bf{r}},\omega)$ and
${\bf{E}}({\bf{r}},\omega)=-\nabla\phi({\bf{r}},\omega)+i\omega{\bf{A}}({\bf{r}},\omega)$.
In the Lorentz gauge,
$\nabla\cdot{\bf{A}}({\bf{r}},\omega)=i\frac{\omega}{v^{2}}\phi({\bf{r}},\omega)$
with $v=1/\sqrt{\epsilon\mu}$, the electromagnetic potentials satisfy

$\displaystyle-\left(\nabla^{2}+k_{0}^{2}\right)\phi({\bf{r}},\omega)$
$\displaystyle=\frac{\rho({\bf{r}},\omega)}{\epsilon}-\frac{\tilde{\alpha}}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\delta(z)\,\hat{{\bf{e}}}_{z}\cdot\nabla\times{\bf{A}}({\bf{r}},\omega),$
(23)
$\displaystyle-\left(\nabla^{2}+k_{0}^{2}\right){\bf{A}}({\bf{r}},\omega)$
$\displaystyle=\mu\,{\bf{J}}({\bf{r}},\omega)-\tilde{\alpha}\mu\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\delta(z)\,\hat{{\bf{e}}}_{z}\times\left[\nabla\phi({\bf{r}},\omega)-i\omega{\bf{A}}({\bf{r}},\omega)\right],$
(24)

where $k_{0}=\omega/v$. To solve this system of equations we follow the same
perturbative program than before: we introduce the expansions
$\phi({\bf{r}},\omega)=\phi^{(0)}({\bf{r}},\omega)+\sum_{n=1}^{\infty}\phi^{(n)}({\bf{r}},\omega)$
and
${\bf{A}}({\bf{r}},\omega)={\bf{A}}^{(0)}({\bf{r}},\omega)+\sum_{n=1}^{\infty}{\bf{A}}^{(n)}({\bf{r}},\omega)$,
where $\phi^{(n)}({\bf{r}},\omega)$ and ${\bf{A}}^{(n)}({\bf{r}},\omega)$
vanish as $\tilde{\alpha}^{n}$, and solve the system recursively. The zeroth-
order equations are exactly Maxwell equations for $\tilde{\alpha}=0$, i.e.

$\displaystyle-\left(\nabla^{2}+k_{0}^{2}\right)\phi^{(0)}({\bf{r}},\omega)$
$\displaystyle=\frac{\rho({\bf{r}},\omega)}{\epsilon},$ (25)
$\displaystyle-\left(\nabla^{2}+k_{0}^{2}\right){\bf{A}}^{(0)}({\bf{r}},\omega)$
$\displaystyle=\mu\,{\bf{J}}({\bf{r}},\omega),$ (26)

while the $n$-order potentials are determined by

$\displaystyle-\left(\nabla^{2}+k_{0}^{2}\right)\phi^{(n)}({\bf{r}},\omega)$
$\displaystyle=-\frac{\tilde{\alpha}}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\delta(z)\,\hat{{\bf{e}}}_{z}\cdot\nabla\times{\bf{A}}^{(n-1)}({\bf{r}},\omega),$
(27)
$\displaystyle-\left(\nabla^{2}+k_{0}^{2}\right){\bf{A}}^{(n)}({\bf{r}},\omega)$
$\displaystyle=-\tilde{\alpha}\mu\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\delta(z)\,\hat{{\bf{e}}}_{z}\times\left[\nabla\phi^{(n-1)}({\bf{r}},\omega)-i\omega{\bf{A}}^{(n-1)}({\bf{r}},\omega)\right].$
(28)

The first step towards solving the higher-order equations (27)-(28) is to
solve the zeroth-order equations (25)-(26). This requires the well-known
frequency-dependent Green’s function of electromagnetism J. Schwinger and Tsai
(1998):

$\displaystyle
G_{\omega}({\bf{r}},{\bf{r}}^{\prime})=\frac{e^{ik_{0}|{\bf{r}}-{\bf{r}}^{\prime}|}}{4\pi|{\bf{r}}-{\bf{r}}^{\prime}|},$
(29)

such that

$\displaystyle\phi^{(0)}({\bf{r}},\omega)=\frac{1}{\epsilon}\int
G_{\omega}({\bf{r}},{\bf{r}}^{\prime})\,\rho({\bf{r}}^{\prime},\omega)\,d^{3}{\bf{r}}^{\prime},\qquad{\bf{A}}^{(0)}({\bf{r}},\omega)=\mu\int
G_{\omega}({\bf{r}},{\bf{r}}^{\prime})\,{\bf{J}}({\bf{r}}^{\prime},\omega)\,d^{3}{\bf{r}}^{\prime}.$
(30)

With these results we can evaluate now equations (27) and (28). In this case
we solve the system up to first order in the parameter $\tilde{\alpha}$. As
shown in the Appendix B.1, for the vector potential one gets:

$\displaystyle A_{i}({\bf{r}},\omega)=A_{i}^{(0)}({\bf{r}},\omega)+\int
G_{ij}^{(1)}({\bf{r}},{\bf{r}}^{\prime})\,J_{j}({\bf{r}}^{\prime},\omega)\,d^{3}{\bf{r}}^{\prime},$
(31)

and for the scalar potential one spot the identity
$i\frac{\omega}{v^{2}}\phi({\bf{r}},\omega)=\nabla\cdot{\bf{A}}({\bf{r}},\omega)$.
Here,

$\displaystyle G_{ij}^{(1)}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=i\frac{1}{\omega}\left[(\hat{\bf{e}}_{z}\times\nabla)_{i}\partial_{j}-k_{0}^{2}\,\epsilon_{ijz}\right]\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime}),$
(32)

where

$\displaystyle\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})=-\frac{\tilde{\alpha}}{8\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\int_{0}^{\infty}dk_{\perp}\frac{e^{i\sqrt{k_{0}^{2}-k_{\perp}^{2}}(|z|+|z^{\prime}|)}}{k_{0}^{2}-k_{\perp}^{2}}k_{\perp}\,J_{0}(k_{\perp}|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}_{\perp}|).$
(33)

A detailed derivation of this expression can be found in the Appendix B.2. It
is worth to mention that the exact electromagnetic potentials (and hence the
fields) must be calculated using Eq. (33), which unfortunately cannot be
evaluated in an analytical fashion. However, for many practical calculations,
one can obtain an approximated expression valid for
$|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}_{\perp}|\gg 1$ and $|z|+|z^{\prime}|\gg
1$. In the Appendix B.3 we derive, by two different methods, the
aforementioned expression. The approximation relays in the fact that the
integrand in Eq. (33) is a rapidly oscillating function for large
$|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}_{\perp}|$ and $|z|+|z^{\prime}|$, and
hence the leading contribution can be obtained by using, for example, the
stationary phase approximation. The final result is

$\displaystyle\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})\approx
i\frac{\tilde{\alpha}}{8\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\frac{e^{ik_{0}\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}_{\perp}|^{2}+(|z|+|z^{\prime}|)^{2}}}}{k_{0}(|z|+|z^{\prime}|)}.$
(34)

Once we have found the indexed Green’s function, we are ready to tackle
problems involving different time-dependent sources near the domain wall.

##  V Radiation of a short linear antenna near a $\theta$ domain-wall

The problem of a short linear antenna near a dielectric half-space has
attracted great attention for over the years due to its importance to many
practical applications, in particular, in the context of THz applications Jeon
and Grischkowsky (2006), near-field optics Novotny (1997), plasmonics Maier
and Atwater (2005) and nanophotonics Gordon (2008). The exact solution to this
problem was given first by Sommerfeld in 1909 Sommerfeld (1909). Since then
numerous investigations have been carried out in this subject. The main goal
of this Section is to address the problem of the radiation pattern produced by
a short linear antenna placed near the boundary between two $\theta$-half
spaces. To this end, we shall employ the theory developed in the previous
Section.

Figure 2: Short linear antenna (of length $2d$ and frequency $\omega$) at a
distance $z_{0}$ from the planar interface between two topologically
insulating media.

The geometry of the problem is shown in Fig. 2. The left half-space $z<0$ is
occupied by a topological insulator with axion angle $\theta_{1}$, while the
right half-space $z>0$ is occupied by a TI with axion angle $\theta_{2}$. We
consider a simple model for a transmitting antenna oriented parallel to the
interface, say along the the $x$-axis and placed at a distance $z_{0}$ from
the interface, with uniform current flowing along a conductor of length $2d$.
We assume that this current varies sinusoidally with time as

$\displaystyle{\bf{J}}({\bf{r}},t)=I_{0}\sin(\omega_{0}t)\,\delta(y)\delta(z-z_{0})\Theta(d-|x|)\,\hat{{\bf{e}}}_{x},$
(35)

where $I_{0}$ and $\omega_{0}$ denotes the current amplitude and angular
frequency, respectively, the $\delta$’s are Dirac deltas and $\Theta$ is the
Heaviside step function. Because the current is discontinuous at the ends,
charge must be deposited there, being of opposite sign at each end, thus
forming an electric dipole. The exact zeroth-order scalar and vector
potentials can be calculated in a simple fashion and can be found in a variety
of textbooks on electrodynamics J. Schwinger and Tsai (1998). The only nonzero
component of the vector potential (which is in the same direction as the
current) in the time-domain is

$\displaystyle A_{x}^{(0)}({\bf{r}},t)=\frac{\mu
I_{0}}{4\pi}\int_{-d}^{+d}dx^{\prime}\frac{1}{|{\bf{r}}-{\bf{r}}^{\prime}|}\,\sin\left[\omega_{0}\left(t-\frac{|{\bf{r}}-{\bf{r}}^{\prime}|}{c}\right)\right],$
(36)

where
${\bf{r}}^{\prime}(x^{\prime})=x^{\prime}\hat{{\bf{e}}}_{x}+z_{0}\hat{{\bf{e}}}_{z}$.
The result for the scalar potential in the time-domain is

$\displaystyle\phi^{(0)}({\bf{r}},t)=-\frac{I_{0}}{4\pi\epsilon\omega_{0}}\left[\frac{\cos\left[\omega_{0}\left(t-\frac{|{\bf{r}}-{\bf{r}}_{+}|}{c}\right)\right]}{|{\bf{r}}-{\bf{r}}_{+}|}-\frac{\cos\left[\omega_{0}\left(t-\frac{|{\bf{r}}-{\bf{r}}_{-}|}{c}\right)\right]}{|{\bf{r}}-{\bf{r}}_{-}|}\right],$
(37)

where ${\bf{r}}_{\pm}={\bf{r}}^{\prime}(\pm d)$. One can easily verify that
these potentials satisfy the Lorentz gauge condition.

Figure 3: Illustration of cross-sections of the normalized electric field
[panels (a) and (b)] and magnetic field [panels (c) and (d)] generated by a
short linear antenna in vacuum at $z_{0}$.

In order to evince the impact of the axion field, we first evaluate
numerically the electromagnetic potentials and fields at first order in
$\theta$ by using Eq. (31) together with the exact Green’s function (32) with
$\mathcal{G}_{\omega}$ given by Eq. (33). To this end we introduce the
dimensionless quantities $X=x/z_{0}$, $Y=y/z_{0}$ and $Z=z/z_{0}$. Figure 3
shows the electric and magnetic fields in the case $\theta=0$, i.e. in the
absence of the topological magnetoelectricity. Panels (a) and (b) show the
electric field lines (in units of $E_{0}=\frac{\mu Iv}{4\pi z_{0}}$) on the
planes $Y=0$ and $X=d/Z_{0}$ (where the end-points of the antenna are located)
for $\omega t=\pi/4$, respectively. They exhibit the expected behavior:
dipole-like in panel (a) and almost point-like in panel (b).Panels (c) and (d)
show the magnetic field lines (in units of $B_{0}=E_{0}/v$) on the planes
$X=d/z_{0}$ and $Y=0$ for $\omega t=\pi/4$, respectively. The former exhibits
the usual current-encircling behavior dictated by the Ampère law, while the
latter exhibits the expected results, whose filed lines invert their direction
periodically. In order to show the full electromagnetic fields, including the
effects of the magnetoelectric coupling, we have to introduce further
dimensionless quantities besides the coordinates $(X,Y,Z)$, namely,
$T=t/\tau$, $D=d/z_{0}$, $\Lambda=\omega_{0}z_{0}/c$ and $\tau=z_{0}/c$. This
is so because the appearance of an additional length scale in the system: the
antenna-surface distance $z_{0}$. In Fig. 4 we show the total electromagnetic
fields for $\theta=\pi$, which is a typical value for topological insulators.
We observe that the cross-sections $X=0$ exhibit the same behavior as in the
case $\theta=0$, for both, the electric (panel (a)) and magnetic (panel (c))
fields; however, as we move from $X=0$ the field lines deviates from the
nontopological case. Panels (b) and (d) show the $Y=0$ cross-section, in which
case we observe a considerably different behavior, which cannot be interpreted
in a simple manner. Indeed, the charge and current densities sourcing these
fields are far from trivial, but they can be understood in terms of the half-
integer Hall current at the interface. The normalized Hall current
${\bf{j}}_{\mbox{\scriptsize Hall }}$ (in units of $\sigma_{\mbox{\scriptsize
Hall }}E_{0}$, with $\sigma_{\mbox{\scriptsize Hall }}=e^{2}/2h$) is

$\displaystyle{\bf{j}}_{\mbox{\scriptsize Hall
}}({\bf{r}}_{\perp},t)=\hat{{\bf{e}}}_{x}\partial_{Y}F(X,Y)+\hat{{\bf{e}}}_{y}\left[\partial_{X}F(X,Y)-\Lambda^{2}G(X,Y)\right],$
(38)

where

$\displaystyle
F(X,Y)=\frac{\cos\left[\Lambda(T-R_{+})\right]}{R_{+}}-\frac{\cos\left[\Lambda(T-R_{-})\right]}{R_{-}},\qquad
G(X,Y)=\int_{-D}^{+D}dX^{\prime}\frac{\cos\left[\Lambda(T-S(X^{\prime}))\right]}{S(X^{\prime})},$
(39)

with $S(X^{\prime})=\sqrt{(X-X^{\prime})^{2}+Y^{2}+1}$ and
$R_{\pm}=\sqrt{(X\mp D)^{2}+Y^{2}+1}$.

Figure 4: Illustration of cross-sections of the normalized electric field
[panels (a) and (b)] and magnetic field [panels (c) and (d)] generated by a
short linear antenna at a distance $z_{0}$ from a $\theta$-domain wall at
$z=0$.

Now we are interested in the fields that survive at large distances from the
source, i.e. in the radiation zone, which requires approximating the
potentials for $r\gg\lambda$ and $r\gg d$, where $\lambda=2\pi c/\omega$ is
the wavelength. To this end we use the expansions
$|{\bf{r}}-{\bf{r}}^{\prime}|\approx r-{\bf{n}}\cdot{\bf{r}}^{\prime}$ and
$|{\bf{r}}-{\bf{r}}_{\pm}|\approx r-\frac{zz_{0}}{r}\mp\frac{xd}{r}$, where
${\bf{n}}={\bf{r}}/r$. All in all, the zeroth order potentials (up to order
$1/r$) become

$\displaystyle A_{x}^{(0)}({\bf{r}},t)\approx\frac{\mu
I_{0}}{2\pi}\sin\left[\omega_{0}\left(t-\frac{r^{2}-z_{0}z}{vr}\right)\right]\sin\left(\frac{\omega_{0}d}{v}\frac{x}{r}\right)\frac{v}{\omega_{0}x},\qquad\phi^{(0)}({\bf{r}},t)=v\frac{x}{r}\,A_{x}^{(0)}({\bf{r}},t).$
(40)

Using these potentials we find the electromagnetic fields in the radiation
zone:

$\displaystyle{\bf{B}}^{(0)}=\frac{\mu
I_{0}}{2\pi}\,\cos\left[\omega_{0}\left(t-\frac{r^{2}-z_{0}z}{vr}\right)\right]\sin\left(\frac{\omega_{0}d}{v}\frac{x}{r}\right)\,\frac{1}{x}\hat{{\bf{e}}}_{x}\times\hat{{\bf{n}}},\qquad{\bf{E}}^{(0)}=-v\,\hat{{\bf{n}}}\times{\bf{B}}^{(0)}.$
(41)

Following the same program one can derive expressions for the electromagnetic
potentials at first-order. The details are relegated to the Appendix C. Using
Eq. (31) one finds that the nonzero components of the vector potential in the
time-domain at first-order become

$\displaystyle A_{x}^{(1)}({\bf{r}},t)$
$\displaystyle=\frac{\tilde{\alpha}}{8\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\frac{I_{0}}{\omega_{0}k_{0}(|z|+z_{0})}\frac{\partial^{2}}{\partial
x\partial
y}\int_{-d}^{+d}dx^{\prime}\sin\left[\omega_{0}\left(t-\frac{|{\bf{r}}^{\ast}-{\bf{r}}^{\prime}_{i}|}{v}\right)\right],$
(42) $\displaystyle A_{y}^{(1)}({\bf{r}},t)$
$\displaystyle=-\frac{\tilde{\alpha}}{8\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\frac{I_{0}}{\omega_{0}k_{0}(|z|+z_{0})}\left(\frac{\partial^{2}}{\partial
x^{2}}+k_{0}^{2}\right)\int_{-d}^{+d}dx^{\prime}\sin\left[\omega_{0}\left(t-\frac{|{\bf{r}}^{\ast}-{\bf{r}}^{\prime}_{i}|}{v}\right)\right],$
(43)

where
${\bf{r}}^{\prime}_{i}(x^{\prime})=x^{\prime}\hat{{\bf{e}}}_{x}-z_{0}\hat{{\bf{e}}}_{z}$
is the image point of ${\bf{r}}^{\prime}(x^{\prime})$ and
${\bf{r}}^{\ast}=x\hat{{\bf{e}}}_{x}+y\hat{{\bf{e}}}_{y}+|z|\hat{{\bf{e}}}_{z}$.
The first-order scalar potential in the time-domain is

$\displaystyle\phi^{(1)}({\bf{r}},t)=-I_{0}\frac{\tilde{\alpha}}{8\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\frac{y}{|z|+z_{0}}\int_{-d}^{+d}dx^{\prime}\frac{\sin\left[\omega_{0}\left(t-\frac{|{\bf{r}}^{\ast}-{\bf{r}}^{\prime}_{i}|}{v}\right)\right]}{|{\bf{r}}^{\ast}-{\bf{r}}^{\prime}_{i}|}.$
(44)

In the radiation zone these potentials can be further simplified. Indeed, one
realizes than only the component $A_{y}^{(1)}$ of the vector potential
contributes to radiation (to be precise, the term proportional to $k_{0}^{2}$
in Eq. (43)), while the others decay faster than $1/r$. Using the expansion
$|{\bf{r}}^{\ast}-{\bf{r}}^{\prime}_{i}|\approx
r^{\ast}-\frac{xx^{\prime}-z_{0}|z|}{r^{\ast}}$ with $r^{\ast}=r$, we get

$\displaystyle
A_{y}^{(1)}({\bf{r}},t)\approx-\tilde{\alpha}\frac{I_{0}}{4\pi\epsilon\omega_{0}}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}\frac{\mu}{\epsilon}}\frac{r/x}{(|z|+z_{0})}\frac{1}{v}\sin\left[\omega_{0}\left(t-\frac{r^{2}+z_{0}|z|}{v\,r}\right)\right]\sin\left(\frac{\omega_{0}d}{v}\frac{x}{r}\right),\qquad\phi^{(1)}({\bf{r}},t)=v\frac{y}{r}A_{y}^{(1)}({\bf{r}},t).$
(45)

Keeping terms of the order $1/r$ the electromagnetic fields become

$\displaystyle{\bf{B}}^{(1)}=-\tilde{\alpha}\frac{\mu
I_{0}}{4\pi}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}\frac{\mu}{\epsilon}}\cos\left[\omega_{0}\left(t-\frac{r^{2}+z_{0}|z|}{v\,r}\right)\right]\sin\left(\frac{\omega_{0}d}{v}\frac{x}{r}\right)\frac{r/x}{|z|+z_{0}}\hat{{\bf{e}}}_{y}\times\hat{{\bf{n}}},\qquad{\bf{E}}^{(1)}=-v\,\hat{{\bf{n}}}\times{\bf{B}}^{(1)}.$
(46)

As shown in Ref. Martín-Ruiz  _et al._ (2015), energy conservation is
unchanged in the presence of the $\theta$ term, i.e. given the Poynting vector
${\bf{S}}$ and energy density $u$ defined as usual
${\bf{S}}={\bf{E}}\times{\bf{H}}$ and
$u=\frac{1}{2}({\bf{E}}\cdot{\bf{D}}+{\bf{B}}\cdot{\bf{H}})$, respectively,
one realizes that $\partial_{t}u+\nabla\cdot{\bf{S}}=0$. Momentum conservation
is modified by a surface term due to the self-induced charge and current
densities arising there, but it does not concern us in this paper. The angular
distribution of the radiated power takes the usual form
$\frac{dP}{d\Omega}=\frac{v}{\mu}r^{2}|{\bf{B}}|^{2}$, with however a total
magnetic field which acquires a contribution arising from the $\theta$-term.
Keeping the lowest order in the magnetoelectric polarizability
$\tilde{\alpha}$ and averaging (in time) over a complete cycle, the angular
distribution of the radiated power by the linear antenna becomes

$\displaystyle\left<\frac{dP}{d\Omega}\right>\approx\frac{v\mu
I_{0}^{2}}{8\pi^{2}}\sin^{2}\left(\frac{\omega_{0}d}{v}\frac{x}{r}\right)\left[\frac{1-(x/r)^{2}}{(x/r)^{2}}-2(\theta/\pi)\sigma_{\mbox{\scriptsize
Hall
}}\sqrt{\frac{\mu}{\epsilon}}\,\frac{y}{x}\frac{r}{|z|}\cos\left(\frac{z_{0}\omega}{v}\frac{z+|z|}{r}\right)\right].$
(47)

The first term corresponds to the usual radiated power by a linear antenna
embedded in an infinite material media described by the permittivity
$\epsilon$ and permeability $\mu$. The second term, which is proportional to
the surface Hall conductivity, exhibits an interesting divergent behavior near
the interface, as is depicted in Fig. 5. The physical origin of this term is
understood from the surface Hall current and charge densities at the
interface, which generates electromagnetic field which contribution to
radiations in the near zone to the interface.

Figure 5: Angular distribution of the radiated power by a thin linear antenna
near a $\theta$ domain-wall. The position of the interface is evident.

##  VI Summary and discussion

Topological phases are an emerging class of materials which have attracted
much attention in condensed matter physics. The remarkable microscopic
properties encoded in the band structure of topological phases produces
exceptional macroscopic electromagnetic response, which is described by a
topological field theory akin to axion electrodynamics. In general, a linear
magnetoelectric media can be described by adding the term
$\sim\theta_{ij}E_{i}B_{j}$ to the free enthalpy of the system, where
$\theta_{ij}$ is the magnetoelectric tensor. In this frame, axion
electrodynamics is a particular case when the ME tensor is of the form
$\theta_{ij}=\theta\delta_{ij}$.

In this paper we derived, in a perturbative fashion, solutions to the time-
dependent axion electrodynamics field equations for a picewise constant axion
angle. The validity of the perturbative program relies in the fact that the
axion coupling, in the case of topological insulators, is proportional to the
fine structure constant. In the particle physics context this condition is
vastly fulfilled, since it is expected to be very small. The difference is
that in the former case the presence of boundaries is a natural task, while in
the latter the axion is assumed to be distributed over the whole space.

To validate our perturbative scheme, we first tackle the static axion field
equations and derive the components of the Green’s matrix up to second order
in the axion coupling. We show that our results are consistent with the
previously reported exact Green’s matrix when Taylor expanded up to second
order in the coupling. We discuss in passing the image magnetic monopole
effect, which is a salient manifestation of the unconventional magnetoelectric
response of topological insulators. In the time-dependent case, the axion
field equations become more involved and no analytical results have been
reported yet, although there are some first attempts towards this direction.
Here we aim to fill in this gap by applying the same perturbative program to
the axion field equations. For the sake of simplicity, in this paper we choose
that the bulks have the same optical properties, namely,
$\epsilon_{1}=\epsilon_{2}\equiv\epsilon$ and $\mu_{1}=\mu_{2}\equiv\mu$,
which is made to highlight and isolate the effects of the axion coupling. In a
forthcoming paper we will relax this condition and generalize our results to
the case of two topologically insulating media characterized by different
nontrivial optical parameters. This generalization is of course far from
trivial. To illustrate the importance of the axion angle in the time-dependent
case considered here, we investigate the radiation of a short linear antenna
above a topologically insulating domain-wall with different values of
$\theta$.

We first computed numerically the electromagnetic fields due to the linear
antenna by using the exact Green’s function and verify that they exhibit the
expected behavior. We interpret the induced fields in terms of the Hall
current at the interface between the two half-spaces characterized by
different $\theta$-angles. Next, to investigate the effects of the
$\theta$-term in the power radiated, we use the asymptotic expression for the
Green’s function and evaluate the fields in the radiation zone. These field
allow us to compute the time-averaged angular distribution of the radiated
power, which as expected, is dominated by the usual Maxwellian term, but
becomes corrected due to the $\theta$-term. Such correction is quite
interesting, since exhibits a divergent-type behavior near the interface,
which arises from the field produced by the surface charge and current
densities there. This is a distinctive behavior with respect to the usual
radiation pattern.

In summary, a general solution to the static axion field equations has been
previously elaborated in terms of Green’s function techniques; however, the
time-dependent case has received less attention mainly due to its complexity.
In this paper we introduced a perturbative scheme which was successfully
validated with the static theory and then we extended the analysis to the
time-dependent case. We obtained general expressions for the corresponding
matrix of Green’s function which we also approximated in the far-field regime.
We expect our method and results will be of considerable relevance and that
they may constitute the basis for numerous other researches which include more
intricate model of antennas and the radiation of particles near the
$\theta$-interface.

###### Acknowledgements.

M. I.-M. was supported by the CONAHCyT fellowship No. 843426. A. M.-R. and M.
I.-M. has been partially supported by DGAPA-UNAM Project No. IA102722 and by
Project CONAHCyT (México) No. 428214. Useful discussions with L. F. Urrutia,
M. Cambiaso and L. Medel are warmly appreciated.

## Data Availability Statement

No datasets were generated or analyzed during the current study.

## Declarations

Conflict of interest The authors declare no conflict of interest.

##  Appendix A Static Green’s function

###  A.1 First-order Green’s functions

Let us consider the first-order corrections. The solution to Eq. (9) are

$\displaystyle\phi^{(1)}({\bf{r}})=\int
G^{(1)}_{0i}({\bf{r}},{\bf{r}}^{\prime})\,J_{i}({\bf{r}}^{\prime})\,d^{3}{\bf{r}}^{\prime},\qquad
A^{(1)}_{i}({\bf{r}})=\int
G^{(1)}_{i0}({\bf{r}},{\bf{r}}^{\prime})\,\rho({\bf{r}}^{\prime})\,d^{3}{\bf{r}}^{\prime},$
(48)

where

$\displaystyle G^{(1)}_{0i}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=-\tilde{\alpha}\,\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\int
G_{\epsilon}({\bf{r}},{\bf{r}}^{\prime\prime})\,\delta(z^{\prime\prime})\,\epsilon_{zji}\,\partial_{j}^{\prime\prime}\,G_{\mu}^{(\perp)}({\bf{r}}^{\prime\prime},{\bf{r}}\,^{\prime})\,d^{3}{\bf{r}}^{\prime\prime},$
(49) $\displaystyle G^{(1)}_{i0}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=-\tilde{\alpha}\,\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\int
G_{\mu}^{(\perp)}({\bf{r}},{\bf{r}}^{\prime\prime})\,\delta(z^{\prime\prime})\,\epsilon_{zji}\,\partial_{j}^{\prime\prime}\,G_{\epsilon}({\bf{r}}^{\prime\prime},{\bf{r}}\,^{\prime})\,d^{3}{\bf{r}}^{\prime\prime}.$
(50)

Noting that partial derivatives apply only in the transverse coordinates
$x^{\prime\prime}$ and $y^{\prime\prime}$, the properties of the Dirac delta
allow us to compute the integrals with respect to the
$z^{\prime\prime}$-coordinate in a simple fashion. Now, if we use
$\partial_{j}^{\prime\prime}\,G_{\mu}^{(\perp)}({\bf{r}}^{\prime\prime},{\bf{r}}\,^{\prime})=-\partial_{j}^{\prime}\,G_{\mu}^{(\perp)}({\bf{r}}^{\prime\prime},{\bf{r}}\,^{\prime})$
and
$\partial_{j}^{\prime\prime}\,G_{\epsilon}({\bf{r}}^{\prime\prime},{\bf{r}}\,^{\prime})=-\partial_{j}^{\prime}\,G_{\epsilon}({\bf{r}}^{\prime\prime},{\bf{r}}\,^{\prime})$
for $j=x,y$, the double primed derivatives can be get out of the integrals as
a minus primed derivative. So we are left with

$\displaystyle G^{(1)}_{0i}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\tilde{\alpha}\,\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\epsilon_{zji}\,\partial_{j}^{\prime}\,\int_{\Sigma}G_{\mu}^{(\perp)}({\bf{r}},{\bf{r}}^{\prime\prime}_{\perp})\,G_{\epsilon}({\bf{r}}^{\prime\prime}_{\perp},{\bf{r}}\,^{\prime})\,d^{2}{\bf{r}}^{\prime\prime}_{\perp},$
(51) $\displaystyle G^{(1)}_{i0}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\tilde{\alpha}\,\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\epsilon_{zji}\,\partial_{j}^{\prime}\,\int_{\Sigma}G_{\epsilon}({\bf{r}},{\bf{r}}^{\prime\prime}_{\perp})\,G_{\mu}^{(\perp)}({\bf{r}}^{\prime\prime}_{\perp},{\bf{r}}\,^{\prime})\,d^{2}{\bf{r}}^{\prime\prime}_{\perp}.$
(52)

Using the electric (11) and magnetic (13) Green functions we realize that

$\displaystyle
G^{(1)}_{0i}({\bf{r}},{\bf{r}}^{\prime})=G^{(1)}_{i0}({\bf{r}},{\bf{r}}^{\prime})=\frac{\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}}{(\epsilon_{1}+\epsilon_{2})(1/\mu_{1}+1/\mu_{2})}\,\epsilon_{zji}\,\partial_{j}^{\prime}\,\int_{\Sigma}\frac{1}{\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}^{\prime\prime}{}_{\\!\\!\perp}|^{2}+z^{2}}}\,\frac{1}{\sqrt{|{\bf{r}}_{\perp}^{\prime\prime}-{\bf{r}}^{\prime}{}_{\\!\\!\perp}|^{2}+z^{\prime\,2}}}\,\frac{d^{2}{\bf{r}}^{\prime\prime}_{\perp}}{(2\pi)^{2}}.$
(53)

To carry out integration, we use the following integral representation:

$\displaystyle\frac{1}{|{\bf{r}}-{\bf{r}}^{\prime}|}=4\pi\int_{\Sigma_{{\bf{k}}_{\perp}}}\frac{d^{2}{\bf{k}}_{\perp}}{(2\pi)^{2}}\;e^{i{\bf{k}}_{\perp}\cdot({\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime})}\,\frac{e^{-k_{\perp}|z-z^{\prime}|}}{2k_{\perp}},$
(54)

where $\Sigma_{{\bf{k}}_{\perp}}$ is the two-dimensional space
$k_{x}$-$k_{y}$. Using this representation twice in Eq. (53) we get

$\displaystyle G^{(1)}_{0i}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=G^{(1)}_{i0}({\bf{r}},{\bf{r}}^{\prime})=\frac{\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}}{(\epsilon_{1}+\epsilon_{2})(1/\mu_{1}+1/\mu_{2})}\,\epsilon_{zji}$
$\displaystyle\hskip
85.35826pt\times\partial_{j}^{\prime}\,\int_{\Sigma_{{\bf{k}}_{\perp}}}\frac{d^{2}{\bf{k}}_{\perp}}{2\pi}\;\frac{e^{-k_{\perp}|z|}}{k_{\perp}}\,\int_{\Sigma_{{\bf{k}}_{\perp}^{\prime}}}\frac{d^{2}{\bf{k}}_{\perp}^{\prime}}{2\pi}\;\frac{e^{-k_{\perp}^{\prime}z^{\prime}}}{k_{\perp}^{\prime}}\,\int_{\Sigma}\,e^{i{\bf{k}}_{\perp}\cdot({\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime\prime})}\,e^{i{\bf{k}}_{\perp}^{\prime}\cdot({\bf{r}}_{\perp}^{\prime\prime}-{\bf{r}}_{\perp}^{\prime})}\,\frac{d^{2}{\bf{r}}^{\prime\prime}_{\perp}}{(2\pi)^{2}}.$
(55)

Integration with respect to the coordinates ${\bf{r}}^{\prime\prime}_{\perp}$
produces a Dirac delta on the momenta ${\bf{k}}_{\perp}$ and
${\bf{k}}_{\perp}^{\prime}$, i.e.

$\displaystyle
G^{(1)}_{0i}({\bf{r}},{\bf{r}}^{\prime})=G^{(1)}_{i0}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\frac{\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}}{(\epsilon_{1}+\epsilon_{2})(1/\mu_{1}+1/\mu_{2})}\,\epsilon_{zji}$
$\displaystyle\hskip
56.9055pt\times\partial_{j}^{\prime}\,\int_{\Sigma_{{\bf{k}}_{\perp}}}\frac{d^{2}{\bf{k}}_{\perp}}{2\pi}\;\frac{e^{-k_{\perp}|z|}}{k_{\perp}}\,\int_{\Sigma_{{\bf{k}}_{\perp}^{\prime}}}\frac{d^{2}{\bf{k}}_{\perp}^{\prime}}{2\pi}\;\frac{e^{-k_{\perp}^{\prime}z^{\prime}}}{k_{\perp}^{\prime}}\,e^{i{\bf{k}}_{\perp}\cdot{\bf{r}}_{\perp}-i{\bf{k}}_{\perp}^{\prime}\cdot{\bf{r}}_{\perp}^{\prime}}\,\delta({\bf{k}}_{\perp}^{\prime}-{\bf{k}}_{\perp})$
$\displaystyle=\frac{\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}}{(\epsilon_{1}+\epsilon_{2})(1/\mu_{1}+1/\mu_{2})}\,\epsilon_{zji}\;\partial_{j}^{\prime}\,\int_{\Sigma_{{\bf{k}}_{\perp}}}\frac{d^{2}{\bf{k}}_{\perp}}{(2\pi)^{2}}\;e^{i{\bf{k}}_{\perp}\cdot({\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime})}\;\frac{e^{-k_{\perp}(|z|+z^{\prime})}}{k_{\perp}^{2}},$
(56)

where we used the properties of the Dirac delta to perform integration over
${\bf{k}}_{\perp}^{\prime}$. The resulting double integral becomes easier to
perform if we express the area element in polar coordinates, i.e.
$d^{2}{\bf{k}}_{\perp}=k_{\perp}dk_{\perp}d\varphi$, and choose the $k_{x}$
axis in the direction of the vector
${\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}$. Noting that
${\bf{k}}_{\perp}\cdot({\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime})=|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|\,k_{\perp}\cos\varphi$,
we can write

$\displaystyle
G^{(1)}_{0i}({\bf{r}},{\bf{r}}^{\prime})=G^{(1)}_{i0}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\frac{\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}}{2\pi(\epsilon_{1}+\epsilon_{2})(1/\mu_{1}+1/\mu_{2})}\,\epsilon_{zji}\;\partial_{j}^{\prime}\,\int_{0}^{\infty}dk_{\perp}\,\frac{e^{-k_{\perp}(|z|+z^{\prime})}}{k_{\perp}}\left(\frac{1}{2\pi}\int_{0}^{2\pi}d\varphi\,e^{i|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|\,k_{\perp}\cos\varphi}\right).$
(57)

The parentheses in this equation enclose an integral representation of the
Bessel function $J_{0}$. The resulting integral is

$\displaystyle
G^{(1)}_{0i}({\bf{r}},{\bf{r}}^{\prime})=G^{(1)}_{i0}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\frac{\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}}{2\pi(\epsilon_{1}+\epsilon_{2})(1/\mu_{1}+1/\mu_{2})}\,\epsilon_{zji}\;\partial_{j}^{\prime}\,\int_{0}^{\infty}dk_{\perp}\,J_{0}(|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|\,k_{\perp})\,\frac{e^{-k_{\perp}(|z|+z^{\prime})}}{k_{\perp}},$
(58)

and taking the primed derivative

$\displaystyle
G^{(1)}_{0i}({\bf{r}},{\bf{r}}^{\prime})=G^{(1)}_{i0}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\frac{\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}}{2\pi(\epsilon_{1}+\epsilon_{2})(1/\mu_{1}+1/\mu_{2})}\,\frac{\epsilon_{zji}\,({\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime})_{j}}{|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|}\,\int_{0}^{\infty}dk_{\perp}\,J_{1}(|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|\,k_{\perp})\,e^{-k_{\perp}(|z|+z^{\prime})}.$
(59)

This integral is well known Gradshteyn and Ryzhik (2007). The final result is
then

$\displaystyle
G^{(1)}_{0i}({\bf{r}},{\bf{r}}^{\prime})=G^{(1)}_{i0}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\frac{\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}}{2\pi(\epsilon_{1}+\epsilon_{2})(1/\mu_{1}+1/\mu_{2})}\,\frac{\epsilon_{zji}\,({\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime})_{j}}{|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|^{2}}\,\left[1-\frac{|z|+z^{\prime}}{\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|^{2}+(|z|+z^{\prime})^{2}}}\right].$
(60)

###  A.2 Second-order Green’s functions

We now turn to the second-order corrections to the GFs. Integration of Eq. (9)
produces

$\displaystyle\phi^{(2)}({\bf{r}})=\int
G^{(2)}_{00}({\bf{r}},{\bf{r}}^{\prime})\,\rho({\bf{r}}^{\prime})\,d^{3}{\bf{r}}^{\prime},\qquad
A^{(2)}_{i}({\bf{r}})=\int
G^{(2)}_{ij}({\bf{r}},{\bf{r}}^{\prime})\,J_{j}({\bf{r}}^{\prime})\,d^{3}{\bf{r}}^{\prime},$
(61)

where

$\displaystyle G^{(2)}_{00}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=-\tilde{\alpha}\,\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\int
G_{\epsilon}({\bf{r}},{\bf{r}}^{\prime\prime})\,\delta(z^{\prime\prime})\,\epsilon_{zji}\,\partial_{j}^{\prime\prime}\,G_{i0}^{(1)}({\bf{r}}^{\prime\prime},{\bf{r}}\,^{\prime})\,d^{3}{\bf{r}}^{\prime\prime},$
(62) $\displaystyle G^{(2)}_{ij}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=-\tilde{\alpha}\,\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\int
G_{\mu}^{(\perp)}({\bf{r}},{\bf{r}}^{\prime\prime})\,\delta(z^{\prime\prime})\,\epsilon_{zki}\,\partial_{k}^{\prime\prime}\,G_{0j}^{(1)}({\bf{r}}^{\prime\prime},{\bf{r}}\,^{\prime})\,d^{3}{\bf{r}}^{\prime\prime}.$
(63)

Using the first-order Green’s function given by Eq. (60), we verify that
$\partial_{j}^{\prime\prime}\,G_{i0}^{(1)}({\bf{r}}^{\prime\prime},{\bf{r}}\,^{\prime})=-\partial_{j}^{\prime}\,G_{i0}^{(1)}({\bf{r}}^{\prime\prime},{\bf{r}}\,^{\prime})$
for $j=1,2$, and hence the double primed derivatives can be get out of the
integrals as a minus primed derivative, as before. Performing the integration
with respect to $z^{\prime\prime}$ we get

$\displaystyle G^{(2)}_{00}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\tilde{\alpha}\,\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\epsilon_{zji}\,\partial_{j}^{\prime}\,\int_{\Sigma}G_{\epsilon}({\bf{r}},{\bf{r}}^{\prime\prime}_{\perp})\,G_{i0}^{(1)}({\bf{r}}^{\prime\prime}_{\perp},{\bf{r}}\,^{\prime})\,d^{2}{\bf{r}}^{\prime\prime}_{\perp},$
(64) $\displaystyle G^{(2)}_{ij}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\tilde{\alpha}\,\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\,\epsilon_{zki}\,\partial_{k}^{\prime}\int_{\Sigma}G_{\mu}^{(\perp)}({\bf{r}},{\bf{r}}^{\prime\prime}_{\perp})\,G_{0j}^{(1)}({\bf{r}}^{\prime\prime}_{\perp},{\bf{r}}\,^{\prime})\,d^{2}{\bf{r}}^{\prime\prime}_{\perp}.$
(65)

Using the integral representations (51) for the $0i$-components and (52) we
find

$\displaystyle G^{(2)}_{00}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=(\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}\,)^{2}\,\epsilon_{zji}\,\epsilon_{zki}\,\partial_{j}^{\prime}\,\partial_{k}^{\prime}\,\int_{\Sigma}\int_{\Sigma}G_{\epsilon}({\bf{r}},{\bf{r}}^{\prime\prime}_{\perp})\,G_{\epsilon}({\bf{r}}^{\prime\prime}_{\perp},{\bf{r}}^{\prime\prime\prime}_{\perp})\,G_{\mu}^{(\perp)}({\bf{r}}^{\prime\prime\prime}_{\perp},{\bf{r}}\,^{\prime})\,d^{2}{\bf{r}}^{\prime\prime\prime}_{\perp}\,d^{2}{\bf{r}}^{\prime\prime}_{\perp},$
(66) $\displaystyle G^{(2)}_{ij}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=(\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}\,)^{2}\,\epsilon_{zki}\,\epsilon_{zlj}\,\partial_{l}^{\prime}\,\partial_{k}^{\prime}\,\int_{\Sigma}\int_{\Sigma}G_{\mu}^{(\perp)}({\bf{r}},{\bf{r}}^{\prime\prime}_{\perp})\,G_{\epsilon}({\bf{r}}^{\prime\prime}_{\perp},{\bf{r}}^{\prime\prime\prime}_{\perp})\,G_{\mu}^{(\perp)}({\bf{r}}^{\prime\prime\prime}_{\perp},{\bf{r}}\,^{\prime})\,d^{2}{\bf{r}}^{\prime\prime\prime}_{\perp}\,d^{2}{\bf{r}}^{\prime\prime}_{\perp},$
(67)

and substituting the electric (11) and magnetic (13) Green functions we get

$\displaystyle G^{(2)}_{00}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\frac{(\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}\,)^{2}}{2\pi(\epsilon_{1}+\epsilon_{2})^{2}(1/\mu_{1}+1/\mu_{2})}\;\mbox{tr}\,[{\bf{I}}({\bf{r}},{\bf{r}}^{\prime})],$
(68) $\displaystyle G^{(2)}_{ij}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\frac{(\tilde{\alpha}\sqrt{\epsilon_{0}/\mu_{0}}\,)^{2}}{2\pi(\epsilon_{1}+\epsilon_{2})(1/\mu_{1}+1/\mu_{2})^{2}}\;I_{ij}({\bf{r}},{\bf{r}}^{\prime}),$
(69)

where we have defined the symmetric matrix
${\bf{I}}=(I_{ij})\in\mathbb{R}^{2\times 2}$, with components

$\displaystyle
I_{ij}({\bf{r}},{\bf{r}}^{\prime})=(\delta_{ij}\delta_{lm}-\delta_{il}\delta_{jm})\,\partial_{l}^{\prime}\,\partial_{m}^{\prime}\,\int_{\Sigma}\int_{\Sigma}\,\frac{1}{\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}^{\prime\prime}{}_{\\!\\!\perp}|^{2}+z^{2}}}\frac{1}{\sqrt{|{\bf{r}}_{\perp}^{\prime\prime}-{\bf{r}}^{\prime\prime\prime}{}_{\\!\\!\perp}|^{2}+0^{2}}}\frac{1}{\sqrt{|{\bf{r}}_{\perp}^{\prime\prime\prime}-{\bf{r}}^{\prime}{}_{\\!\\!\perp}|^{2}+z^{\prime\,2}}}\,\frac{d^{2}{\bf{r}}^{\prime\prime\prime}_{\perp}}{2\pi}\,\frac{d^{2}{\bf{r}}^{\prime\prime}_{\perp}}{2\pi}$
(70)

To evaluate this expression, we use thrice the integral representation (54),
i.e.

$\displaystyle I_{ij}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=(\delta_{ij}\delta_{lm}-\delta_{il}\delta_{jm})\,\partial_{l}^{\prime}\,\partial_{m}^{\prime}\,\int_{\Sigma}\int_{\Sigma}\,\left(\int_{\Sigma_{{\bf{k}}_{\perp}}}\frac{d^{2}{\bf{k}}_{\perp}}{2\pi}e^{i{\bf{k}}_{\perp}\cdot({\bf{r}}_{\perp}-{\bf{r}}^{\prime\prime}_{\perp})}\frac{e^{-k_{\perp}|z|}}{k_{\perp}}\right)\,\left(\int_{\Sigma_{{\bf{k}}_{\perp}}}\frac{d^{2}{\bf{k}}_{\perp}^{\prime}}{2\pi}e^{i{\bf{k}}_{\perp}^{\prime}\cdot({\bf{r}}_{\perp}^{\prime\prime}-{\bf{r}}^{\prime\prime\prime}_{\perp})}\frac{e^{-k_{\perp}^{\prime}(0)}}{k_{\perp}^{\prime}}\right)$
$\displaystyle\hskip
199.16928pt\left(\int_{\Sigma_{{\bf{k}}_{\perp}}}\frac{d^{2}{\bf{k}}_{\perp}^{\prime\prime}}{2\pi}e^{i{\bf{k}}_{\perp}^{\prime\prime}\cdot({\bf{r}}_{\perp}^{\prime\prime\prime}-{\bf{r}}^{\prime}_{\perp})}\frac{e^{-k_{\perp}^{\prime\prime}z^{\prime}}}{k_{\perp}^{\prime\prime}}\right)\,\frac{d^{2}{\bf{r}}^{\prime\prime\prime}_{\perp}}{2\pi}\,\frac{d^{2}{\bf{r}}^{\prime\prime}_{\perp}}{2\pi}.$
(71)

Integrating with respect to ${\bf{r}}^{\prime\prime}_{\perp}$ and
${\bf{r}}^{\prime\prime\prime}_{\perp}$ produces two Dirac deltas in momentum:

$\displaystyle I_{ij}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=(\delta_{ij}\delta_{lm}-\delta_{il}\delta_{jm})\,\partial_{l}^{\prime}\,\partial_{m}^{\prime}\,\int_{\Sigma_{{\bf{k}}_{\perp}}}\int_{\Sigma_{{\bf{k}}_{\perp}}}\int_{\Sigma_{{\bf{k}}_{\perp}}}\frac{d^{2}{\bf{k}}_{\perp}\,d^{2}{\bf{k}}_{\perp}^{\prime}\,d^{2}{\bf{k}}_{\perp}^{\prime\prime}}{2\pi}\frac{e^{-k_{\perp}|z|}e^{-k_{\perp}^{\prime\prime}z^{\prime}}}{k_{\perp}k_{\perp}^{\prime}k_{\perp}^{\prime\prime}}e^{i{\bf{k}}_{\perp}\cdot{\bf{r}}_{\perp}-i{\bf{k}}_{\perp}^{\prime\prime}\cdot{\bf{r}}_{\perp}^{\prime}}\,\delta({\bf{k}}_{\perp}^{\prime}-{\bf{k}}_{\perp})\,\delta({\bf{k}}_{\perp}^{\prime\prime}-{\bf{k}}_{\perp}^{\prime})$
$\displaystyle=(\delta_{ij}\delta_{lm}-\delta_{il}\delta_{jm})\,\partial_{l}^{\prime}\,\partial_{m}^{\prime}\,\int_{\Sigma_{{\bf{k}}_{\perp}}}\frac{d^{2}{\bf{k}}_{\perp}}{2\pi}\frac{e^{-k_{\perp}(|z|+z^{\prime})}}{k_{\perp}^{3}}e^{i{\bf{k}}_{\perp}\cdot({\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime})}$
(72)

where we have used the properties of the Dirac delta. We now introduce the
polar coordinates $(k_{x},k_{y})=k_{\perp}(\cos\varphi,\sin\varphi)$ and
choose the vector ${\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}$ along the
$k_{x}$ direction. So, the angular integral yields

$\displaystyle I_{ij}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=(\delta_{ij}\delta_{lm}-\delta_{il}\delta_{jm})\,\partial_{l}^{\prime}\,\partial_{m}^{\prime}\,\int_{0}^{\infty}dk_{\perp}\,\frac{e^{-k_{\perp}(|z|+z^{\prime})}}{k_{\perp}^{2}}J_{0}(\,|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|k_{\perp}\,).$
(73)

This integral does not converge as is, but its derivatives does. Taking one of
the derivatives we find

$\displaystyle I_{ij}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=(\delta_{ij}\delta_{lm}-\delta_{il}\delta_{jm})\,\partial_{l}^{\prime}\,\left[\frac{({\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime})_{m}}{|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|\;}\,\int_{0}^{\infty}dk_{\perp}\,\frac{e^{-k_{\perp}(|z|+z^{\prime})}}{k_{\perp}}J_{1}(\,|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|k_{\perp}\,)\right].$
(74)

To perform this integral in a convergent manner, we use the following trick:
$\int_{\Lambda}^{\infty}e^{-k_{\perp}\lambda}d\lambda=e^{-k_{\perp}\Lambda}/k_{\perp}$,
and therefore

$\displaystyle I_{ij}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=(\delta_{ij}\delta_{lm}-\delta_{il}\delta_{jm})\,\partial_{l}^{\prime}\,\left[\frac{({\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime})_{m}}{|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|\;}\,\int_{|z|+z^{\prime}}^{\infty}d\lambda\int_{0}^{\infty}dk_{\perp}\,e^{-k_{\perp}\lambda}\,J_{1}(\,|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|k_{\perp}\,)\right],$
(75)

and this integral is well-known; see, for example, Ref. Gradshteyn and Ryzhik
(2007). The result is

$\displaystyle I_{ij}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=(\delta_{ij}\delta_{lm}-\delta_{il}\delta_{jm})\,\partial_{l}^{\prime}\,\left[\frac{({\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime})_{m}}{|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|^{2}}\,\int_{|z|+z^{\prime}}^{\infty}d\lambda\left(1-\frac{\lambda}{\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|^{2}+\lambda^{2}}}\right)\right],$
(76)

from which we finally obtain

$\displaystyle I_{ij}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=(\delta_{ij}\delta_{lm}-\delta_{il}\delta_{jm})\,\partial_{l}^{\prime}\,\left[\frac{({\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime})_{m}}{|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|^{2}}\,\left(\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|^{2}+(|z|+z^{\prime})^{2}}-(|z|+z^{\prime})\right)\right].$
(77)

From this expression we find

$\displaystyle\mbox{tr}\,[{\bf{I}}({\bf{r}},{\bf{r}}^{\prime})]=\frac{1}{\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}{}_{\\!\\!\perp}|^{2}+(|z|+z^{\prime})^{2}}},$
(78)

and

$\displaystyle I_{ij}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\delta_{ij}\,\mbox{tr}\,[{\bf{I}}({\bf{r}},{\bf{r}}^{\prime})]+\partial_{i}\,K_{j}({\bf{r}},{\bf{r}}^{\prime}),$
(79)

where

$\displaystyle
K_{j}({\bf{r}},{\bf{r}}^{\prime})=\frac{({\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime})_{j}}{|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|^{2}}\,\left(\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}|^{2}+(|z|+z^{\prime})^{2}}-(|z|+z^{\prime})\right).$
(80)

##  Appendix B Frequency-dependent Green’s function

###  B.1 Derivation of the Green’s functions

Substituting the zeroth-order solutions (30) into Eqs. (27) and (28) one finds
that the first order solution are given by

$\displaystyle\phi^{(1)}({\bf{r}},\omega)$
$\displaystyle=-\epsilon_{zij}\partial_{i}\int\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})\,J_{j}({\bf{r}}^{\prime},\omega)\,d^{3}{\bf{r}}^{\prime}$
(81) $\displaystyle A^{(1)}_{i}({\bf{r}},\omega)$
$\displaystyle=-\epsilon_{izj}\partial_{j}\int\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})\rho({\bf{r}}^{\prime})d^{3}{\bf{r}}^{\prime}+i\frac{\omega}{v^{2}}\epsilon_{izj}\int\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})\,J_{j}({\bf{r}}^{\prime},\omega)\,d^{3}{\bf{r}}^{\prime},$
(82)

where we have defined

$\displaystyle\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})=\tilde{\alpha}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\int
G_{\omega}({\bf{r}},{\bf{r}}^{\prime\prime}_{\perp})\,G_{\omega}({\bf{r}}^{\prime\prime}_{\perp},{\bf{r}}^{\prime})\,d^{2}{\bf{r}}^{\prime\prime}_{\perp}.$
(83)

We now perform son manipulations in Eqs. (81) and (82). Taking the divergence
of Eq. (82) we obtain

$\displaystyle\partial_{i}A^{(1)}_{i}({\bf{r}},\omega)$
$\displaystyle=i\frac{\omega}{v^{2}}\epsilon_{izj}\partial_{i}\int\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})\,J_{j}({\bf{r}}^{\prime},\omega)\,d^{3}{\bf{r}}^{\prime}=i\frac{\omega}{v^{2}}\phi^{(1)}({\bf{r}},\omega),$
(84)

where we used that $\epsilon_{zij}\partial_{i}\partial_{j}=0$. So, one we have
computed the vector potential, the scalar potential becomes fully determined.
In order to manipulate Eq. (82) we use the local charge conservation:
$\nabla\cdot{\bf{J}}=i\omega\rho$. Substituting the charge density in Eq. (82)
produces

$\displaystyle A^{(1)}_{i}({\bf{r}},\omega)$
$\displaystyle=-\frac{1}{i\omega}\epsilon_{izj}\partial_{j}\int\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})[\partial^{\prime}_{k}J_{k}({\bf{r}}^{\prime},\omega)]d^{3}{\bf{r}}^{\prime}+i\frac{\omega}{v^{2}}\epsilon_{izj}\int\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})\,J_{j}({\bf{r}}^{\prime},\omega)\,d^{3}{\bf{r}}^{\prime}.$
(85)

Integrating by parts and using the identity
$\partial^{\prime}_{i}\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})=-\partial_{i}\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})$
we obtain

$\displaystyle A^{(1)}_{i}({\bf{r}},\omega)$
$\displaystyle=\frac{1}{i\omega}\epsilon_{zij}\partial_{j}\partial_{k}\int\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})\,J_{k}({\bf{r}}^{\prime},\omega)\,d^{3}{\bf{r}}^{\prime}+i\frac{\omega}{v^{2}}\epsilon_{izk}\int\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})\,J_{k}({\bf{r}}^{\prime},\omega)\,d^{3}{\bf{r}}^{\prime}$
$\displaystyle=\frac{i}{\omega}\left(-\epsilon_{zij}\partial_{j}\partial_{k}+\frac{\omega^{2}}{v^{2}}\epsilon_{izk}\right)\int\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})\,J_{k}({\bf{r}}^{\prime},\omega)\,d^{3}{\bf{r}}^{\prime},$
(86)

which can be written as

$\displaystyle A^{(1)}_{i}({\bf{r}},\omega)$
$\displaystyle=\int\mathcal{G}_{ij}({\bf{r}},{\bf{r}}^{\prime})\,J_{j}({\bf{r}}^{\prime},\omega)\,d^{3}{\bf{r}}^{\prime},$
(87)

where

$\displaystyle\mathcal{G}_{ij}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=\frac{i}{\omega}\left(-\epsilon_{zij}\partial_{j}\partial_{k}+\frac{\omega^{2}}{v^{2}}\epsilon_{izk}\right)\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime}).$
(88)

This expression corresponds to Eq. (32) of the main text.

In the following we derive an exact expression for the function
$\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})$, which cannot be expressed
in term of known functions. We also obtain an analytical expression in the
radiation zone (far field approximation) by using two different methods.

###  B.2 Exact expression for the function
$\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})$

The integral we have to evaluate is given by Eq. (83). Substituting the
frequency-dependent Green’s function (29) the integral becomes:

$\displaystyle\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})=\tilde{\alpha}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\int
d^{2}{\bf{r}}^{\prime\prime}_{\perp}\,\frac{e^{ik_{0}|{\bf{r}}-{\bf{r}}_{\perp}^{\prime\prime}|}}{4\pi\,|{\bf{r}}-{\bf{r}}_{\perp}^{\prime\prime}|}\,\frac{e^{ik_{0}|{\bf{r}}_{\perp}^{\prime\prime}-{\bf{r}}^{\prime}|}}{4\pi\,|{\bf{r}}_{\perp}^{\prime\prime}-{\bf{r}}^{\prime}|}.$
(89)

To evaluate this integral we now use the Sommerfeld identity:

$\displaystyle\frac{e^{ik_{0}|{\bf{r}}|}}{4\pi\,|{\bf{r}}|}=\int_{\Sigma_{{\bf{k}}_{\perp}}}\frac{d^{2}{\bf{k}}_{\perp}}{(2\pi)^{2}}e^{i{\bf{k}}_{\perp}\cdot{\bf{r}}_{\perp}}\frac{ie^{i\sqrt{k_{0}^{2}-k_{\perp}^{2}}|z|}}{2\sqrt{k_{0}^{2}-k_{\perp}^{2}}},$
(90)

where $\Sigma_{{\bf{k}}_{\perp}}$ is the two-dimensional space $k_{x}$-$k_{y}$
and ${\bf{r}}={\bf{r}}_{\perp}+|z|\hat{{\bf{e}}}_{z}$. Using this
representation twice in Eq. (89) we get

$\displaystyle\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})=\tilde{\alpha}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\int_{\Sigma_{{\bf{k}}_{\perp}}}\frac{d^{2}{\bf{k}}_{\perp}}{(2\pi)^{2}}\int_{\Sigma_{{\bf{k}}^{\prime}_{\perp}}}\frac{d^{2}{\bf{k}}^{\prime}_{\perp}}{(2\pi)^{2}}\,e^{i({\bf{k}}_{\perp}\cdot{\bf{r}}_{\perp}-{\bf{k}}^{\prime}_{\perp}\cdot{\bf{r}}_{\perp}^{\prime})}\frac{ie^{i\sqrt{k_{0}^{2}-k_{\perp}^{2}}|z|}}{2\sqrt{k_{0}^{2}-k_{\perp}^{2}}}\,\frac{ie^{i\sqrt{k_{0}^{2}-k_{\perp}^{\prime\,2}}|z^{\prime}|}}{2\sqrt{k_{0}^{2}-k_{\perp}^{\prime\,2}}}\int
d^{2}{\bf{r}}^{\prime\prime}_{\perp}\,e^{i({\bf{k}}^{\prime}_{\perp}-{\bf{k}}_{\perp})\cdot{\bf{r}}_{\perp}^{\prime\prime}}$
(91)

Integration over ${\bf{r}}_{\perp}^{\prime\prime}$ produces a Dirac delta in
momentum, which allow us to evaluate one of the momentum integrals. We obtain:

$\displaystyle\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})=-\frac{\tilde{\alpha}}{4}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\int_{\Sigma_{{\bf{k}}_{\perp}}}\frac{d^{2}{\bf{k}}_{\perp}}{(2\pi)^{2}}e^{i{\bf{k}}_{\perp}\cdot({\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime})}\frac{e^{i\sqrt{k_{0}^{2}-k_{\perp}^{2}}(|z|+|z^{\prime}|)}}{k_{0}^{2}-k_{\perp}^{2}},$
(92)

As before, we now introduce the polar coordinates
$(k_{x},k_{y})=k_{\perp}(\cos\varphi,\sin\varphi)$ and choose the vector
${\bf{r}}_{\perp}-{\bf{r}}_{\perp}^{\prime}$ along the $k_{x}$ direction. So,
the angular integral yields

$\displaystyle\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})=-\frac{\tilde{\alpha}}{8\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\int_{0}^{\infty}dk_{\perp}\frac{e^{i\sqrt{k_{0}^{2}-k_{\perp}^{2}}(|z|+|z^{\prime}|)}}{k_{0}^{2}-k_{\perp}^{2}}k_{\perp}\,J_{0}(k_{\perp}|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}_{\perp}|).$
(93)

This is an exact expression for the required function. It corresponds to Eq.
(33) of the main text.

###  B.3 Far-field approximation

Here we shall approximate the integral (93) in the radiation zone, i.e. for
$\rho,|z|\gg 1$. For crosscheking our result, we compute the integral by two
different methods. One of them used the Sommerfeld integral (Method 1) and the
other is the stationary phase approximation (Method 2).

Method 1. We start by using the identity
$q^{-1}=\int_{0}^{\infty}d\xi\,e^{-q\xi}$ in Eq. (93), for
$q=\sqrt{k_{0}^{2}-k_{\perp}^{2}}$ to obtain

$\displaystyle\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=-\frac{\tilde{\alpha}}{8\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\int_{0}^{\infty}d\xi\,\int_{0}^{\infty}dk_{\perp}\,\frac{k_{\perp}}{\sqrt{k_{0}^{2}-k_{\perp}^{2}}}e^{i\sqrt{k_{0}^{2}-k_{\perp}^{2}}(|z|+|z^{\prime}|+i\xi)}\,J_{0}(k_{\perp}|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}_{\perp}|).$
(94)

Using the Sommerfeld integral (90) this expression can be written in the
simple form:

$\displaystyle\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})=i\frac{\tilde{\alpha}}{8\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\int_{0}^{\infty}d\xi\,\frac{e^{ik_{0}R_{\xi}}}{R_{\xi}},$
(95)

where $R_{\xi}\equiv
R_{\xi}({\bf{r}},{\bf{r}}^{\prime})=\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}_{\perp}|^{2}+(|z|+|z^{\prime}|+i\xi)^{2}}$,
such that
$R_{0}=\sqrt{|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}_{\perp}|^{2}+(|z|+|z^{\prime}|)^{2}}$.
Next we shall approximate this integral in the far-field region. To this end
we use the fact that
$R_{\xi}=R_{0}\sqrt{1+\frac{2i\xi(|z|+|z^{\prime}|)-\xi^{2}}{R_{0}^{2}}}$ to
rewritte the integral (95) as

$\displaystyle\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})=i\frac{\tilde{\alpha}}{8\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\int_{0}^{\infty}d\chi\,\frac{e^{i\zeta\sqrt{1+2i\chi\tilde{z}-\chi^{2}}}}{\sqrt{1+2i\chi\tilde{z}-\chi^{2}}},$
(96)

where $\zeta=k_{0}R_{0}$ and $\tilde{z}=(|z|+|z^{\prime}|)/R_{0}$. Numerical
evaluation of this integral for $\zeta\gg 1$ and $\tilde{z}\ll 1$ reveals that
only small values of $\chi$ contribute to the integral. So we can expand the
root and retain solely the leading order, i.e.

$\displaystyle\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})\approx
i\frac{\tilde{\alpha}}{8\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}e^{i\zeta}\int_{0}^{\infty}d\chi\,\frac{e^{-\left(\zeta\tilde{z}\right)\chi}}{1+i\chi\zeta}=i\frac{\tilde{\alpha}}{8\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\frac{1}{\tilde{z}}\left[\pi+i\,\text{Ei}(i\zeta)\right],$
(97)

whose approximation for $\zeta\gg 1$ is

$\displaystyle\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})\approx
i\frac{\tilde{\alpha}}{8\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\frac{e^{ik_{0}R_{0}}}{k_{0}(|z|+|z^{\prime}|)}.$
(98)

This result corresponds to Eq. (34) of the main text.

Method 2. Since the Bessel function and the complex exponential in the
integrand of Eq. (93) are rapidly oscillating functions in the far-field
region, one can use the stationary phase approximation to evaluate the
integral. To this end we first use the identity
$J_{0}(x)=\tfrac{1}{2}[H^{(1)}_{0}(x)+H^{(2)}_{0}(x)]$ together with the
reflection formula $H^{(1)}_{0}(e^{i\pi}x)=-H^{(2)}_{0}(x)$, where
$H^{(1)}_{0}$ and $H^{(2)}_{0}$ are the Hankel functions. In this case, the
integral (93) can be written as

$\displaystyle\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle=-\frac{\tilde{\alpha}}{16\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\int_{-\infty}^{\infty}dk_{\perp}\,\frac{k_{\perp}}{k_{0}^{2}-k_{\perp}^{2}}e^{ik_{z}(|z|+|z^{\prime}|)}\,H^{(1)}_{0}(k_{\perp}|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}_{\perp}|).$
(99)

We now implement the far-field approximation into the Hankel function. For
$x\gg 1$ we have $H^{(1)}_{0}(x)\approx\sqrt{\frac{2}{\pi
x}}e^{ix-i\frac{\pi}{4}}$, and the integral (99) becomes

$\displaystyle\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})$
$\displaystyle\approx-\frac{\tilde{\alpha}}{16\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\int_{-\infty}^{\infty}dk_{\perp}\,\frac{k_{\perp}}{k_{0}^{2}-k_{\perp}^{2}}\sqrt{\frac{2}{\pi
k_{\perp}|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}_{\perp}|}}\exp\left[i|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}_{\perp}|\left(k_{\perp}+\sqrt{k_{0}^{2}-k_{\perp}^{2}}\frac{|z|+|z^{\prime}|}{|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}_{\perp}|}\right)-i\frac{\pi}{4}\right].$
(100)

We observe that, in the far-field, the exponential factor is a rapidly
oscillating function of $k_{\perp}$ and hence we can apply the stationary
phase approximation to evaluate the dominant contribution. A simple
calculation reveals that the stationary phase point locates at

$\displaystyle
k_{\perp}^{\ast}=\frac{k_{0}}{\sqrt{1+\left(\frac{|z|+|z^{\prime}|}{|{\bf{r}}_{\perp}-{\bf{r}}^{\prime}_{\perp}|}\right)^{2}}}.$
(101)

Finally, we use the well-known asymptotic formula

$\displaystyle\int_{\mathbb{R}}g(x)e^{i\lambda f(x)}dx\approx
g(x_{0})e^{i\lambda
f(x_{0})+i\frac{\pi}{4}\mbox{sgn}(f^{\prime\prime}(x_{0}))}\sqrt{\frac{2\pi}{\lambda|f^{\prime\prime}(x_{0})|}},$
(102)

where $\lambda\gg 1$ and $f^{\prime}(x_{0})=0$, to obtain

$\displaystyle\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})\approx
i\frac{\tilde{\alpha}}{8\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\frac{e^{ik_{0}R_{0}}}{k_{0}(|z|+|z^{\prime}|)},$
(103)

which is again the result of Eq. (98). This validates our expression for the
far-field approximation.

##  Appendix C Evaluation of the first-order potentials for a short linear
antenna

Our starting point here is Eq. (31) for the components of the vector
potential. The first order contribution in the frequency-domain is read off as

$\displaystyle A_{i}^{(1)}({\bf{r}},\omega)$
$\displaystyle=i\frac{1}{\omega}\left[(\hat{\bf{e}}_{z}\times\nabla)_{i}\partial_{j}-k_{0}^{2}\,\epsilon_{ijz}\right]\int\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})\,J_{j}({\bf{r}}^{\prime},\omega)\,d^{3}{\bf{r}}^{\prime},$
(104)

where $\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})$ is the corresponding
Green’s function. The potential in the time-domain is computed by inverse
Fourier transforming the above expression. We obtain:

$\displaystyle A_{i}^{(1)}({\bf{r}},t)$
$\displaystyle=i\frac{1}{2\pi}\left[(\hat{\bf{e}}_{z}\times\nabla)_{i}\partial_{j}+\frac{1}{v^{2}}\frac{\partial^{2}}{\partial
t^{2}}\,\epsilon_{ijz}\right]\int d\omega\,e^{-i\omega
t}\frac{1}{\omega}\int\mathcal{G}_{\omega}({\bf{r}},{\bf{r}}^{\prime})\,J_{j}({\bf{r}}^{\prime},\omega)\,d^{3}{\bf{r}}^{\prime}.$
(105)

The current source (35) in the frequency domain is

$\displaystyle{\bf{J}}({\bf{r}},\omega)$ $\displaystyle=\int e^{i\omega
t}{\bf{J}}({\bf{r}},t)dt=I\,\delta(y)\delta(z-z_{0})\Theta(d-|x|)\,\hat{{\bf{e}}}_{x}\frac{\pi}{i}\left[\delta(\omega+\omega_{0})-\delta(\omega-\omega_{0})\right].$
(106)

Substituting (106) into Eq. (105) and integrating the frequency one obtains

$\displaystyle A_{i}^{(1)}({\bf{r}},t)$
$\displaystyle=-\frac{1}{2\pi}\left[(\hat{\bf{e}}_{z}\times\nabla)_{i}\partial_{x}+\frac{1}{v^{2}}\frac{\partial^{2}}{\partial
t^{2}}\,\epsilon_{ixz}\right]\frac{I\pi}{\omega_{0}}\int_{-d}^{+d}dx^{\prime}\,\left[e^{i\omega_{0}t}\mathcal{G}_{-\omega_{0}}({\bf{r}},{\bf{r}}^{\prime\prime})+e^{-i\omega_{0}t}\mathcal{G}_{\omega_{0}}({\bf{r}},{\bf{r}}^{\prime\prime})\right],$
(107)

where
${\bf{r}}^{\prime\prime}=x^{\prime}\hat{{\bf{e}}}_{x}+z_{0}\hat{{\bf{e}}}_{z}$.
Since we are interested in the far-field domain, we use the Green’s function
(34). In such case the vector potential becomes

$\displaystyle
A_{i}^{(1)}({\bf{r}},t)=-\frac{\tilde{\alpha}}{8\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\frac{I}{\omega_{0}k_{0}(|z|+z_{0})}\left[(\hat{\bf{e}}_{z}\times\nabla)_{i}\partial_{x}-\frac{\omega_{0}^{2}}{v^{2}}\,\epsilon_{ixz}\right]\int_{-d}^{+d}dx^{\prime}\sin\left[\omega_{0}\left(t-\frac{|{\bf{r}}^{\ast}-{\bf{r}}^{\prime}_{i}|}{v}\right)\right],$
(108)

where
${\bf{r}}^{\prime}_{i}(x^{\prime})=x^{\prime}\hat{{\bf{e}}}_{x}-z_{0}\hat{{\bf{e}}}_{z}$
is the image point of ${\bf{r}}^{\prime}(x^{\prime})$ and
${\bf{r}}^{\ast}=x\hat{{\bf{e}}}_{x}+y\hat{{\bf{e}}}_{y}+|z|\hat{{\bf{e}}}_{z}$.
Equations (42) and (43) of the main text are derived from this result by
taking $i=x$ and $i=y$, respectively. Now, to derive the scalar potential, we
appeal to the Lorentz gauge condition. Taking the divergence of the vector
potential we obtain

$\displaystyle\partial_{i}A_{i}^{(1)}({\bf{r}},\omega)$
$\displaystyle=-\frac{\tilde{\alpha}}{8\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\frac{I}{\omega_{0}k_{0}(|z|+z_{0})}\frac{\omega_{0}^{2}}{v^{2}}\,\frac{\partial}{\partial
y}\int_{-d}^{+d}dx^{\prime}\sin\left[\omega_{0}\left(t-\frac{|{\bf{r}}^{\ast}-{\bf{r}}^{\prime}_{i}|}{v}\right)\right]=-\frac{1}{v^{2}}\frac{\partial}{\partial
t}\phi^{(1)}({\bf{r}},t),$ (109)

wherefrom we find

$\displaystyle\phi^{(1)}({\bf{r}},t)=-\frac{\tilde{\alpha}}{8\pi}\frac{\mu}{\epsilon}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\frac{I}{k_{0}(|z|+z_{0})}\,\frac{\partial}{\partial
y}\int_{-d}^{+d}dx^{\prime}\cos\left[\omega_{0}\left(t-\frac{|{\bf{r}}^{\ast}-{\bf{r}}^{\prime}_{0i}|}{v}\right)\right],$
(110)

which, after taking the $y$-derivative, corresponds to Eq. (44) of the main
text.

##  Appendix D Derivation of Eq. (47)

The angular distribution of the radiated power is given by
$\frac{dP}{d\Omega}=\frac{v}{\mu}r^{2}|{\bf{B}}|^{2}$. Keeping linear order in
the magnetic field, i.e.
$|{\bf{B}}|^{2}\approx|{\bf{B}}^{(0)}|^{2}+2{\bf{B}}^{(0)}\cdot{\bf{B}}^{(1)}$,
one finds
$\frac{dP}{d\Omega}\approx\frac{v}{\mu}r^{2}\left[|{\bf{B}}^{(0)}|^{2}+2{\bf{B}}^{(0)}\cdot{\bf{B}}^{(1)}\right]$.
Substituting the zeroth (41) and first (46) order magnetic fields we obtain

$\displaystyle\frac{dP}{d\Omega}$
$\displaystyle\approx\frac{v}{\mu}\left(\frac{\mu
I_{0}}{2\pi}\right)^{2}\cos\left[\omega_{0}\left(t-\frac{r^{2}-z_{0}z}{vr}\right)\right]\sin^{2}\left(\frac{\omega_{0}d}{v}\frac{x}{r}\right)\left\\{\cos\left[\omega_{0}\left(t-\frac{r^{2}-z_{0}z}{vr}\right)\right]\,\frac{1-(x/r)^{2}}{(x/r)^{2}}\right.$
$\displaystyle\hskip
199.16928pt\left.+2\tilde{\alpha}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}\frac{\mu}{\epsilon}}\cos\left[\omega_{0}\left(t-\frac{r^{2}+z_{0}|z|}{v\,r}\right)\right]\,\frac{y/x}{|z|+z_{0}}r\right\\},$
(111)

where we used that
$(\hat{{\bf{e}}}_{x}\times\hat{{\bf{n}}})\cdot(\hat{{\bf{e}}}_{x}\times\hat{{\bf{n}}})=1-(x/r)^{2}$
and
$(\hat{{\bf{e}}}_{x}\times\hat{{\bf{n}}})\cdot(\hat{{\bf{e}}}_{y}\times\hat{{\bf{n}}})=-xy/r^{2}$.

Finally, averaging in time over a complete cycle yields

$\displaystyle\left<\frac{dP}{d\Omega}\right>$
$\displaystyle\approx\frac{v}{2\mu}\left(\frac{\mu
I_{0}}{2\pi}\right)^{2}\sin^{2}\left(\frac{\omega_{0}d}{v}\frac{x}{r}\right)\left\\{\frac{1-(x/r)^{2}}{(x/r)^{2}}+2\tilde{\alpha}\sqrt{\frac{\epsilon_{0}}{\mu_{0}}\frac{\mu}{\epsilon}}\cos\left[\frac{(z+|z|)z_{0}\omega}{rv}\right]\,\frac{y/x}{|z|+z_{0}}r\right\\},$
(112)

where we have used:

$\displaystyle\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^{+T}\cos^{2}\left[\omega_{0}\left(t-\frac{r^{2}-z_{0}z}{vr}\right)\right]dt$
$\displaystyle=\frac{1}{2},$ (113)
$\displaystyle\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^{+T}\cos\left[\omega_{0}\left(t-\frac{r^{2}-z_{0}z}{vr}\right)\right]\cos\left[\omega_{0}\left(t-\frac{r^{2}+z_{0}|z|}{vr}\right)\right]dt$
$\displaystyle=\frac{1}{2}\cos\left[\frac{(z+|z|)z_{0}\omega}{rv}\right].$
(114)

These results stablish Eq. (5).

## References

  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, “Topological insulators and superconductors,” Rev. Mod. Phys. 83, 1057–1110 (2011).
  * König  _et al._ (2007) M. König, S. Wiedmann, C. Brüne, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, “Quantum spin Hall insulator state in HgTe quantum wells,” Science 318, 766–770 (2007).
  * Chen  _et al._ (2009) Y. L. Chen, J. G. Analytis, J.-H. Chu, Z. K. Liu, S.-K. Mo, X. L. Qi, H. J. Zhang, D. H. Lu, X. Dai, Z. Fang, S. C. Zhang, I. R. Fisher, Z. Hussain, and Z.-X. Shen, “Experimental realization of a three-dimensional topological insulator, bi¡sub¿2¡/sub¿te¡sub¿3¡/sub¿,” Science 325, 178–181 (2009).
  * Xia  _et al._ (2009) Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, and M. Z. Hasan, “Observation of a large-gap topological-insulator class with a single Dirac cone on the surface,” Nature Physics 5, 398–402 (2009).
  * Sato  _et al._ (2010) T. Sato, K. Segawa, H. Guo, K. Sugawara, S. Souma, T. Takahashi, and Y. Ando, “Direct evidence for the Dirac-cone topological surface states in the ternary chalcogenide ${\mathrm{tlbise}}_{2}$,” Phys. Rev. Lett. 105, 136802 (2010).
  * Wilczek (1987) F. Wilczek, “Two applications of axion electrodynamics,” Phys. Rev. Lett. 58, 1799–1802 (1987).
  * Peccei and Quinn (1977) R. D. Peccei and Helen R. Quinn, “$\mathrm{CP}$ conservation in the presence of pseudoparticles,” Phys. Rev. Lett. 38, 1440–1443 (1977).
  * Weinberg (1978) S. Weinberg, “A new light boson?” Phys. Rev. Lett. 40, 223–226 (1978).
  * Wilczek (1978) F. Wilczek, “Problem of strong $p$ and $t$ invariance in the presence of instantons,” Phys. Rev. Lett. 40, 279–282 (1978).
  * Preskill  _et al._ (1983) J. Preskill, M. B. Wise, and F. Wilczek, “Cosmology of the invisible axion,” Physics Letters B 120, 127–132 (1983).
  * Dine and Fischler (1983) M. Dine and W. Fischler, “The not-so-harmless axion,” Physics Letters B 120, 137–141 (1983).
  * Abbott and Sikivie (1983) L.F. Abbott and P. Sikivie, “A cosmological bound on the invisible axion,” Physics Letters B 120, 133–136 (1983).
  * Fiebig (2005) M. Fiebig, “Revival of the magnetoelectric effect,” Journal of Physics D: Applied Physics 38, R123 (2005).
  * Qi  _et al._ (2009) X.-L. Qi, R. Li, J. Zang, and S.-C. Zhang, “Inducing a magnetic monopole with topological surface states,” Science 323, 1184–1187 (2009).
  * Karch (2009) A. Karch, “Electric-magnetic duality and topological insulators,” Phys. Rev. Lett. 103, 171601 (2009).
  * Zang and Nagaosa (2010) J. Zang and N. Nagaosa, “Monopole current and unconventional Hall response on a topological insulator,” Phys. Rev. B 81, 245125 (2010).
  * Martín-Ruiz  _et al._ (2015) A. Martín-Ruiz, M. Cambiaso, and L. F. Urrutia, “Green’s function approach to Chern-Simons extended electrodynamics: An effective theory describing topological insulators,” Phys. Rev. D 92, 125015 (2015).
  * Martín-Ruiz  _et al._ (2019) A. Martín-Ruiz, O. Rodríguez-Tzompantzi, J. R. Maze, and L. F. Urrutia, “Magnetoelectric effect of a conducting sphere near a planar topological insulator,” Phys. Rev. A 100, 042124 (2019).
  * Martín-Ruiz  _et al._ (2016) A. Martín-Ruiz, M. Cambiaso, and L. F. Urrutia, “Electromagnetic description of three-dimensional time-reversal invariant ponderable topological insulators,” Phys. Rev. D 94, 085019 (2016).
  * Martín-Ruiz (2018) A. Martín-Ruiz, “Magnetoelectric effect in cylindrical topological insulators,” Phys. Rev. D 98, 056012 (2018).
  * Meier  _et al._ (2019) Q. N. Meier, M. Fechner, T. Nozaki, M. Sahashi, Z. Salman, T. Prokscha, A. Suter, P. Schoenherr, M. Lilienblum, P. Borisov, I. E. Dzyaloshinskii, M. Fiebig, H. Luetkens, and N. A. Spaldin, “Search for the magnetic monopole at a magnetoelectric surface,” Phys. Rev. X 9, 011011 (2019).
  * Chang and Yang (2009) M.-C. Chang and M.-F. Yang, “Optical signature of topological insulators,” Phys. Rev. B 80, 113304 (2009).
  * Tse and MacDonald (2011) W.-K. Tse and A. H. MacDonald, “Magneto-optical Faraday and Kerr effects in topological insulator films and in other layered quantized Hall systems,” Phys. Rev. B 84, 205327 (2011).
  * Crosse (2016) J. A. Crosse, “Optical properties of topological-insulator Bragg gratings: Faraday-rotation enhancement for TM-polarized light at large incidence angles,” Phys. Rev. A 94, 033816 (2016).
  * Wu  _et al._ (2016) L. Wu, M. Salehi, N. Koirala, J. Moon, S. Oh, and N. P. Armitage, “Quantized Faraday and Kerr rotation and axion electrodynamics of a 3D topological insulator,” Science 354, 1124–1127 (2016).
  * Filipini and Cambiaso (2023) Sebastián Filipini and Mauro Cambiaso, “New polarization rotation and exact tem wave solutions in topological insulators,” (2023), arXiv:2308.14673 [cond-mat.mes-hall] .
  * Grushin and Cortijo (2011) A. G. Grushin and A. Cortijo, “Tunable Casimir repulsion with three-dimensional topological insulators,” Phys. Rev. Lett. 106, 020403 (2011).
  * Grushin  _et al._ (2011) A. G. Grushin, P. Rodriguez-Lopez, and A. Cortijo, “Effect of finite temperature and uniaxial anisotropy on the Casimir effect with three-dimensional topological insulators,” Phys. Rev. B 84, 045119 (2011).
  * Rodriguez-Lopez and Grushin (2014) P. Rodriguez-Lopez and A. G. Grushin, “Repulsive Casimir effect with Chern insulators,” Phys. Rev. Lett. 112, 056804 (2014).
  * Chen and Wan (2012) L. Chen and S. Wan, “Critical surface band gap of repulsive Casimir interaction between three-dimensional topological insulators at finite temperature,” Phys. Rev. B 85, 115102 (2012).
  * Rodriguez-Lopez (2011) P. Rodriguez-Lopez, “Casimir repulsion between topological insulators in the diluted regime,” Phys. Rev. B 84, 165409 (2011).
  * Nie  _et al._ (2013) W. Nie, R. Zeng, Y. Lan, and S. Zhu, “Casimir force between topological insulator slabs,” Phys. Rev. B 88, 085421 (2013).
  * Martín-Ruiz  _et al._ (2016) A. Martín-Ruiz, M. Cambiaso, and L. F. Urrutia, “A Green’s function approach to the Casimir effect on topological insulators with planar symmetry,” EPL (Europhysics Letters) 113, 60005 (2016).
  * Song  _et al._ (2014) G. Song, J.-P. Xu, and Y.-P. Yang, “Spontaneous emission of a two-level system near the interface of topological insulators,” EPL (Europhysics Letters) 105, 64001 (2014).
  * Fang  _et al._ (2015) W. Fang, Z.-X. Yang, and G.-X. Li, “Quantum properties of an atom in a cavity constructed by topological insulators,” Journal of Physics B: Atomic, Molecular and Optical Physics 48, 245504 (2015).
  * Martín-Ruiz and Urrutia (2018) A. Martín-Ruiz and L. F. Urrutia, “Interaction of a hydrogenlike ion with a planar topological insulator,” Phys. Rev. A 97, 022502 (2018).
  * Martín-Ruiz and Chan-López (2017) A. Martín-Ruiz and E. Chan-López, “Dynamics of a Rydberg hydrogen atom near a topologically insulating surface,” EPL (Europhysics Letters) 119, 53001 (2017).
  * Castro-Enriquez  _et al._ (2020) L. A. Castro-Enriquez, L. F. Quezada, and A. Martín-Ruiz, “Optical response of a topological-insulator–quantum-dot hybrid interacting with a probe electric field,” Phys. Rev. A 102, 013720 (2020).
  * Castro-Enriquez and Martín-Ruiz (2020) L. A. Castro-Enriquez and A. Martín-Ruiz, “Fano resonances induced by the topological magnetoelectric effect,” EPL (Europhysics Letters) 131, 64003 (2020).
  * Franca  _et al._ (2019) O. J. Franca, L. F. Urrutia, and O. Rodríguez-Tzompantzi, “Reversed electromagnetic Vavilov-Čerenkov radiation in naturally existing magnetoelectric media,” Phys. Rev. D 99, 116020 (2019).
  * Franca and Urrutia (2022) O. J. Franca and L. F. Urrutia, “Radiation from a dipole perpendicular to the interface between two planar semi-infinite magnetoelectric media,” Rev. Mex. Fís. 68, 060701 (2022).
  * Franca and Buhmann (2022) O. J. Franca and S. Y. Buhmann, “Modification of transition radiation by three-dimensional topological insulators,” Phys. Rev. B 105, 155120 (2022).
  * Ouellet and Bogorad (2019) J. Ouellet and Z. Bogorad, “Solutions to axion electrodynamics in various geometries,” Phys. Rev. D 99, 055010 (2019).
  * Qi  _et al._ (2008) X.-L. Qi, T. L. Hughes, and S.-C. Zhang, “Topological field theory of time-reversal invariant insulators,” Phys. Rev. B 78, 195424 (2008).
  * Zhang  _et al._ (2013) J. Zhang, C.-Z. Chang, P. Tang, Z. Zhang, X. Feng, K. Li, L. l. Wang, X. Chen, C. Liu, W. Duan, K. He, Q.-K. Xue, X. Ma, and Y. Wang, “Topology-driven magnetic quantum phase transition in topological insulators,” Science 339, 1582–1586 (2013).
  * Maciejko  _et al._ (2010) J. Maciejko, X.-L. Qi, H. D. Drew, and S.-C. Zhang, “Topological quantization in units of the fine structure constant,” Phys. Rev. Lett. 105, 166803 (2010).
  * J. Schwinger and Tsai (1998) K. Milton J. Schwinger, L. DeRaad and W. Tsai,  _Classical Electrodynamics_ (Perseus Books, Boulder, 1998).
  * Jeon and Grischkowsky (2006) T.-I. Jeon and D. Grischkowsky, “THz Zenneck surface wave (THz surface plasmon) propagation on a metal sheet,” Applied Physics Letters 88, 061113 (2006).
  * Novotny (1997) L. Novotny, “Allowed and forbidden light in near-field optics. i. a single dipolar light source,” J. Opt. Soc. Am. A 14, 91–104 (1997).
  * Maier and Atwater (2005) S. A. Maier and H. A. Atwater, “Plasmonics: Localization and guiding of electromagnetic energy in metal/dielectric structures,” Journal of Applied Physics 98, 011101 (2005).
  * Gordon (2008) R. Gordon, “Surface plasmon nanophotonics: A tutorial,” IEEE Nanotechnology Magazine 2, 12–18 (2008).
  * Sommerfeld (1909) A. Sommerfeld, “Über die ausbreitung der wellen in der drahtlosen telegraphie,” Annalen der Physik 333, 665–736 (1909).
  * Gradshteyn and Ryzhik (2007) I. S. Gradshteyn and I. M. Ryzhik,  _Table of integrals, series, and products_ , seventh ed. (Elsevier/Academic Press, Amsterdam, 2007).