# Generalized FDTD Scheme for the Simulation of Electromagnetic Scattering in
Moving Structures

Zoé-Lise Deck-Léger1  Amir Bahrami2  Zhiyu Li3,2 and Christophe Caloz2
1Polytechnique Montréal, Canada, 2KU Leuven, Belgium, 3Xi’an Jiaotong
University, China<EMAIL_ADDRESS>
††journal: oe

Electromagnetic scattering in moving structures is a fundamental topic in
physics and engineering. Yet, no general numerical solution to related
problems has been reported to date. We introduce here a generalized FDTD
scheme to remedy this deficiency. That scheme is an extension of the FDTD
standard Yee cell and stencil that includes not only the usual, physical
fields, but also auxiliary, unphysical fields allowing a straightforward
application of moving boundary conditions. The proposed scheme is illustrated
by four examples – a moving interface, a moving slab, a moving crystal and a
moving gradient – with systematic validation against exact solutions.

##  1 Introduction

Moving structures are pervasive in electromagnetics. They may involve both
moving bodies, such stars or vehicles, and moving perturbations, such as fluid
or elastic waves. The category of moving-body structures has a nearly 300-year
history and has lead to the discovery of many fundamental physical phenomena,
such as Bradley aberration [1], Doppler frequency shifting [2], Fizeau
dragging [3], relativity [4], magneto-electric coupling [5], medium
bianisotropy [6] and gravity emulation [7]. The category of moving-
perturbation structures is more recent and particularly amenable to practical
applications; it produces effects such as parametric amplification [8, 9, 10],
nonreciprocal transmission [11, 12], deflected reflection [13, 14], space-time
frequency transitions [15, 16], pseudo Fizeau dragging [17, 18],
spatiotemporal bandgaps [19, 20, 21] and accelerated metamaterial bending
[22].

Surprisingly, no general numerical tool is currently available to simulate
electromagnetic moving structures, whether of the moving-body or moving-
perturbation type. The development of such a tool would be highly desirable to
simulate the plethora of moving structures mentioned above as well as future
ones, particularly the noncanonical – and typically more practical –
structures that do not admit analytical solutions. This issue is well-known.
There have been only two workarounds in the litterature so far [23, 24, 25,
26, 27], both based on the FDTD technique [28, 29], probably selected for its
natural incarnation of both spatial and temporal variations in Maxwell’s
equations. However, one of these approaches is restricted to non-penetrable
objects [23, 24], while the other one implies cumbersome Lorentz frame
transformations [25, 26, 27].

We present here a novel, general and efficient FDTD scheme for simulating
electromagnetic scattering in moving structures. That scheme, contrarily to
that in [23, 24], also applies to penetrable media, allowing hence to handle
gradient structures and metamaterials[17, 30], without requiring problematic
numerical transformations, contrarily to that in [25, 26, 27]. This scheme is
based on a generalized Yee cell and stencil, which include not only the usual,
physical fields, but also auxiliary, unphysical fields that carry the velocity
information, so as to automatically satisfy the moving boundary conditions.

##  2 Types of Systems

The proposed method applies to both moving-matter and moving-perturbation
media [18], the only difference between the two, in terms of modeling, being
that in the former case, motion transforms media that are isotropic at rest
into bianisotropic media [31, 32], whereas in the latter case, motion alters
the values of the permittivity or permeability parameters without affecting
the isotropy of the medium [33, 15, 34]. Figure 1 depicts some of these
structures111In the case of moving-perturbation, these structures are
essentially _traveling-wave_ modulation-type structures, while _standing-wave_
modulation-type structures, such as those studied in [35, 21], which consist
combinations of purely spatial and purely temporal interfaces can be handled
by standard FDTD algorithm. Indeed, the FDTD algorithm naturally enforces
continuity of the tangential $\mathbf{E}$ and $\mathbf{H}$ fields at
stationary interfaces and continuity of $\mathbf{D}$ and $\mathbf{B}$ fields
at temporal interfaces., which will be used for illustration and validation
purposes in Sec. 6: a simple interface [36] [Fig. 1(a)], a slab [37] [Fig.
1(b)], a bilayer crystal [Fig. 1(c)], which can be operated in the Bragg
regime [9] or in the metamaterial (homogeneous) regime [17], and a gradient
[38] [Fig 1(d)], all moving at a velocity $v$ in a direction selected as being
$z$ in a Cartesian coordinate system. The different parts of these structures
can naturally represent arbitrary bianisotropic media, with general tensor
$\overline{\overline{\chi}}=\left[\overline{\overline{\epsilon}},\overline{\overline{\xi}};\overline{\overline{\zeta}},\overline{\overline{\mu}}\right]$,
which may be all different from each other and may be arbitrarily truncated in
the $xy$-plane. Note that all these structures can be seen as a succession of
interfaces, including the gradient considering a decomposition in
subwavelength slabs. For this reason, we shall focus on an interface in the
next section, which deals with boundary conditions.

Figure 1: Types of moving structures, of either the moving-matter or moving-
perturbation type, that can be simulated by the proposed generalized FDTD
scheme. (a) Simple interface. (b) Slab. (c) Bilayer crystal. (d) Gradient.

##  3 Generalized Yee Cell

The boundary conditions at a moving interface, whether associated to moving
modulation or moving matter, read [39, 40]

$\hat{\mathbf{n}}\times(\mathbf{E}_{2}^{*}-\mathbf{E}_{1}^{*})=0,$ (1a)
$\hat{\mathbf{n}}\times(\mathbf{H}_{2}^{*}-\mathbf{H}_{1}^{*})=\mathbf{J}_{\text{s}},$
(1b) $\hat{\mathbf{n}}\cdot(\mathbf{D}_{2}-\mathbf{D}_{1})=\rho_{\text{s}},$
(1c) $\hat{\mathbf{n}}\cdot(\mathbf{B}_{2}-\mathbf{B}_{1})=0,$ (1d)

with

$\mathbf{E}^{*}=\mathbf{E}+\mathbf{v}\times\mathbf{B}\quad\text{and}\quad\mathbf{H}^{*}=\mathbf{H}-\mathbf{v}\times\mathbf{D},$
(2)

where $1$ and $2$ label the media at the two sides the interface [as in Fig.
1(c)], $\mathbf{E}$, $\mathbf{H}$, $\mathbf{D}$, $\mathbf{B}$ are the usual
electromagnetic fields, $\mathbf{J}_{\text{s}}$ and $\rho_{\text{s}}$ are the
usual surface current and charge densities, respectively, $\hat{\mathbf{n}}$
is the unit vector normal to the interface pointing towards medium 1, and
$\mathbf{v}$ is the velocity of the interface. Note that these conditions, as
expected, reduce to the stationary boundary conditions for $\mathbf{v}=0$,
viz., $\hat{\mathbf{n}}\times(\mathbf{E}_{2}-\mathbf{E}_{1})=0$,
$\hat{\mathbf{n}}\times(\mathbf{H}_{2}-\mathbf{H}_{1})=\mathbf{J}_{\text{s}}$,
$\hat{\mathbf{n}}\cdot(\mathbf{D}_{2}-\mathbf{D}_{1})=\rho_{\text{s}}$ and
$\hat{\mathbf{n}}\cdot(\mathbf{B}_{2}-\mathbf{B}_{1})=0$.

The fields in Eq. (2) may decomposed into tangential and normal components
with respect to the interface. The tangential components are

$\mathbf{E}^{*}_{\text{tan}}=\mathbf{E}_{\text{tan}}+\mathbf{v}\times\mathbf{B}-(\hat{\mathbf{n}}\cdot(\mathbf{v}\times\mathbf{B}))\hat{\mathbf{n}}\quad\text{and}\quad\mathbf{H}^{*}_{\text{tan}}=\mathbf{H}_{\text{tan}}-\mathbf{v}\times\mathbf{D}-(\hat{\mathbf{n}}\cdot(\mathbf{v}\times\mathbf{D}))\hat{\mathbf{n}},$
(3a) while the normal components are
$\mathbf{E}^{*}_{\text{norm}}=\mathbf{E}_{\text{norm}}+(\hat{\mathbf{n}}\cdot(\mathbf{v}\times\mathbf{B}))\hat{\mathbf{n}}\quad\text{and}\quad\mathbf{H}^{*}_{\text{norm}}=\mathbf{H}_{\text{norm}}+(\hat{\mathbf{n}}\cdot(\mathbf{v}\times\mathbf{D}))\hat{\mathbf{n}}.$
(3b)

In the case where the interface is perpendicular to the direction of motion,
i.e., $\mathbf{n}\|\mathbf{v}$, the right-most terms in the right hand sides
of Eqs. (3a) and Eqs. (3b) vanish. This is the scenario that is represented in
Fig. 1 and that will be assumed in the sequel of the paper222If we had
$\mathbf{n}\cancel{\|}\mathbf{v}$, the interface would be slanted with respect
to the motion direction. Related problems were treated analytically in [41,
42]. The paper will further assume non-conducting interfaces, and hence
$\mathbf{J}_{\text{s}}=\rho_{\text{s}}=0$ in Eq. (1).

The dynamic boundary conditions (1) differ from the static or stationary
boundary equations only insofar as they involve the auxiliary fields
$\mathbf{E}_{1,2}^{*}$ and $\mathbf{H}_{1,2}^{*}$ instead of the fields
$\mathbf{E}_{1,2}$ and $\mathbf{H}_{1,2}$, whose tangential components include
velocity-dependent extra terms [Eq. (3a)]. This observation suggests the
generalization of the standard Yee cell [37] to a Yee cell with
correspondingly modified $\mathbf{E}_{\text{tan}}$ and
$\mathbf{H}_{\text{tan}}$ fields, without alteration of the conventional
electric and magnetic field staggering arrangement. The result is shown in
Fig. 2, for a pair of generalized Yee cells separated by a moving interface.
The so-defined generalized Yee cell straightforwardly allows to satisfy the
field continuity conditions at a moving interface, i.e., Eqs. (1a) and (1b),
by having the field components $E_{x}^{*}$, $E_{y}^{*}$ and $H_{z}$ uniquely
defined at the permittivity interface of Fig. 2 (or $H_{x}^{*}$, $H_{y}^{*}$
and $E_{z}$ in the case of a corresponding magnetic interface).

Figure 2: Pair of generalized Yee cells, separated by a moving interface [Fig.
1(a)], with the generalization consisting in the substitution of the
tangential physical fields $\mathbf{E}$ and $\mathbf{H}$ by the corresponding
tangential auxiliary fields $\mathbf{E}^{*}$ and $\mathbf{H}^{*}$ in Eq. (2).

This simple generalization of the Yee cell is the core idea of the paper. This
idea might a priori appear to be trivial, but it is not so trivial insofar as
it implies a _mixture of physical and unphysical fields_ , namely, as seen in
Eqs. (1) and (3), the physical fields $\mathbf{D}$, $\mathbf{B}$,
$\mathbf{E}_{\mathrm{norm}}$ and $\mathbf{H}_{\mathrm{norm}}$ and the
unphysical fields $\mathbf{E}_{\mathrm{tan}}^{*}$ and
$\mathbf{H}_{\mathrm{tan}}^{*}$, which will also imply unphysical mixed-field
electromagnetic equations. This strange situation is somewhat reminiscent to
the unphysical anisotropic media appearing in Bérenger’s Perfect Matched
Layers (PMLs) [43] and, as the PMLs, it does not pose any fundamental problem.
The corresponding unphysical Maxwell’s equations constitutive relations, which
will be derived and discretized in Sec. 4 for establishing the generalized
FDTD update equations, are mathematical mappings of the physical Maxwell’s
equations via the simple changes of variable in Eqs. (3), and the physical
fields $\mathbf{E}$ and $\mathbf{H}$ can be found anytime from the auxiliary
fields by simply inverting these relations.

##  4 Generalized Scheme

In this section, we shall derive the generalized FDTD scheme corresponding to
the generalized Yee cell established in Sec. 3. This scheme must, of course,
model the physics of the problem, and hence follow the usual Maxwell-Faraday
and Maxwell-Ampère equations,

$\frac{\partial\mathbf{B}}{\partial
t}=-\nabla\times\mathbf{E}\quad\text{and}\quad\frac{\partial\mathbf{D}}{\partial
t}=\nabla\times\mathbf{H},$ (4)

as well as the general constitutive relations,

$\mathbf{D}=\overline{\overline{\epsilon}}(\mathbf{r}-\mathbf{v}t)\cdot\mathbf{E}+\overline{\overline{\xi}}(\mathbf{r}-\mathbf{v}t)\cdot\mathbf{H}\quad\text{and}\quad\mathbf{B}=\overline{\overline{\zeta}}(\mathbf{r}-\mathbf{v}t)\cdot\mathbf{E}+\overline{\overline{\mu}}(\mathbf{r}-\mathbf{v}t)\cdot\mathbf{H},$
(5)

whose bianisotropic form [44] pertains to moving matter333In this case, for
motion in the $z$ direction ($\mathbf{v}=v\hat{\mathbf{z}}$), the tensors in
Eq. (5) are [44]
$\overline{\overline{\epsilon}}=\epsilon^{\prime}\begin{bmatrix}\alpha&0&0\\\
0&\alpha&0\\\
0&0&1\end{bmatrix},\quad\overline{\overline{\mu}}=\mu^{\prime}\begin{bmatrix}\alpha&0&0\\\
0&\alpha&0\\\
0&0&1\end{bmatrix},\quad\overline{\overline{\xi}}=\begin{bmatrix}0&\chi/c&0\\\
-\chi/c&0&0\\\
0&0&0\end{bmatrix},\quad\overline{\overline{\zeta}}=\overline{\overline{\xi}}^{T},$
(6) with
$\alpha=\frac{1-\beta^{2}}{1-\beta^{2}{n^{\prime}}^{2}},\qquad\chi=\beta\frac{1-{n^{\prime}}^{2}}{1-\beta^{2}{n^{\prime}}^{2}},\qquad\beta=v/c\qquad\text{and}\qquad
n^{\prime}/c=\sqrt{\epsilon^{\prime}\mu^{\prime}},$ where $\epsilon^{\prime}$
and $\mu^{\prime}$ are constant scalars that represent the permittivity and
permeability of the medium at rest ($v=0$), which is assumed to be isotropic,
nondispersive and linear. and which reduce to simple homoisotropic relations
$\left(\overline{\overline{\epsilon}}=\epsilon,~{}\overline{\overline{\mu}}=\mu,~{}\overline{\overline{\xi}}=\overline{\overline{\zeta}}=0\right)$
in the case of moving modulation [33, 15, 34].

In order to apply the generalized Yee cell established in Sec. 3, we must
express the fields $\mathbf{E}$ and $\mathbf{H}$ in Eqs. (4) and (5) in terms
of the fields $\mathbf{E}^{*}$ and $\mathbf{H}^{*}$ using Eq. (2). This
results into the mixed-field Maxwell’s equations

$\frac{\partial\mathbf{B}}{\partial
t}=-\nabla\times\mathbf{E}^{*}+\nabla\times(\mathbf{v}\times\mathbf{B}),$ (7a)
$\frac{\partial\mathbf{D}}{\partial
t}=\nabla\times\mathbf{H}^{*}+\nabla\times(\mathbf{v}\times\mathbf{D}),$ (7b)

and constitutive relations

$\mathbf{D}=\overline{\overline{\epsilon}}(\mathbf{r}-\mathbf{v}t)\cdot\left(\mathbf{E}^{*}-\mathbf{v}\times\mathbf{B}\right)+\overline{\overline{\xi}}(\mathbf{r}-\mathbf{v}t)\cdot\left(\mathbf{H}^{*}+\mathbf{v}\times\mathbf{D}\right),$
(8a)
$\mathbf{B}=\overline{\overline{\zeta}}(\mathbf{r}-\mathbf{v}t)\cdot\left(\mathbf{E}^{*}-\mathbf{v}\times\mathbf{B}\right)+\overline{\overline{\mu}}(\mathbf{r}-\mathbf{v}t)\cdot\left(\mathbf{H}^{*}+\mathbf{v}\times\mathbf{D}\right).$
(8b)

From this point, we shall restrict our developments, for the sake of
simplicity and without any loss of generality, to structures that are i) based
on moving modulation (as opposed to moving matter), ii) one-dimensional (1D)
and iii) illuminated by a plane wave propagating in the direction of the
modulation (i.e., normal to the interface(s)). Using a Cartesian coordinate
system with $z$ corresponding to the direction of the modulation, the non-zero
field components reduce to $E_{x}$ and $H_{y}$ and Eqs. (7) and (8) become
then

$\frac{\partial B_{y}}{\partial t}=-\frac{\partial E_{x}^{*}}{\partial
z}-v\frac{\partial B_{y}}{\partial z},$ (9a) $\frac{\partial D_{x}}{\partial
t}=-\frac{\partial H_{y}^{*}}{\partial z}-v\frac{\partial D_{x}}{\partial z},$
(9b)

and

$D_{x}=\epsilon(z-vt)\left(E^{*}_{x}+vB_{y}\right),$ (10a)
$B_{y}=\mu(z-vt)\left(H^{*}_{y}+vD_{x}\right),$ (10b)

where we have used the fact that the mapping (2) preserves the polarization in
the 1D case444To show this, we insert
$\mathbf{B}=-\int\left(\nabla\times\mathbf{E}\right)\,\mathrm{d}t$, from (4),
into (2), which yields, using the vectorial identity
$\mathbf{A}\times(\nabla\times\mathbf{B})=\nabla(\mathbf{A}\cdot\mathbf{B})-(\mathbf{A}\cdot\nabla)\mathbf{B}-(\mathbf{B}\cdot\nabla)\mathbf{A}-\mathbf{B}\times(\nabla\times\mathbf{A})$,
$\mathbf{E}^{*}=\mathbf{E}-\smallint\left(\mathbf{v}\times(\nabla\times\mathbf{E})\right)\,\mathrm{d}t=\mathbf{E}-\smallint\left(\nabla(\mathbf{v}\cdot\mathbf{E})-(\mathbf{v}\cdot\nabla)\mathbf{E}-(\mathbf{E}\cdot\nabla)\mathbf{v}-\mathbf{E}\times(\nabla\times\mathbf{v})\right)\,\mathrm{d}t.$
The last two last terms in the last expression both vanish in the assumed
regime of constant velocity. Given $\mathbf{E}=E_{x}\hat{\mathbf{x}}$ and
$\mathbf{v}=v\hat{\mathbf{z}}$, the term with $\mathbf{v}\cdot\mathbf{E}$ also
vanishes since $\mathbf{v}\perp\mathbf{E}$. Finally, the term
$(\mathbf{v}\cdot\nabla)\mathbf{E}$ reduces to $\hat{\mathbf{x}}\partial
E_{x}/\partial z$ and is hence parallel to $\mathbf{E}$. Thus,
$\mathbf{E}^{*}\|\mathbf{E}$. It may be similarly shown that
$\mathbf{H}^{*}\parallel\mathbf{H}$., so that only the star field components
$E_{x}^{*}$ and $H_{y}^{*}$, corresponding to $E_{x}$ and $H_{y}$, are
involved,

We can now proceed to the discretization of Eqs. (9) and (10). The static
($v=0$) parts can be discretized in the usual fashion [29], which yields

$B_{y}|^{n}_{k+\frac{1}{2}}=B_{y}|^{n-1}_{k+\frac{1}{2}}-S\left(E^{*}_{x}|^{n-\frac{1}{2}}_{k+1}-E^{*}_{x}|^{n-\frac{1}{2}}_{k}\right)-v\Delta
t\frac{\partial B_{y}}{\partial z},$ (11a)
$D_{x}|^{n+\frac{1}{2}}_{k}=D_{x}|^{n-\frac{1}{2}}_{k}-S\left(H^{*}_{y}|^{n}_{k+\frac{1}{2}}-H^{*}_{y}|^{n}_{k-\frac{1}{2}}\right)-v\Delta
t\frac{\partial D_{x}}{\partial z},$ (11b) and
$E^{*}_{x}|_{k}^{n+\frac{1}{2}}=\frac{D_{x}|_{k}^{n+\frac{1}{2}}}{\epsilon|_{k}^{n+\frac{1}{2}}}-vB_{y},$
(11c)
$H^{*}_{y}|_{k+\frac{1}{2}}^{n}=\frac{B_{y}|_{k+\frac{1}{2}}^{n}}{\mu|_{k+\frac{1}{2}}^{n}}-vD_{x},$
(11d)

where $S=\Delta t/\Delta z$.

In contrast, the dynamic ($v\neq 0$) terms, which have not been discretized
yet in Eqs. (11), require a special treatment for the numerical stability and
minimal dispersion of the overall scheme. We empirically found that the
discretization choices given in Tab. 1 are appropriate in these regards. Note
that these schemes are different for positive velocities, i.e.,
$\mathbf{v}\|\mathbf{\hat{z}}$ or $v>0$ and negative velocities,
$\mathbf{v}\|-\mathbf{\hat{z}}$ or $v<0$.

Positive velocity ($v>0$) |  | Negative velocity ($v<0$) |
---|---|---|---

$\dfrac{\partial B_{y}}{\partial z}=\frac{B_{y}|^{n-1}_{k+\frac{1}{2}}-B_{y}|^{n-1}_{k-\frac{1}{2}}}{\Delta z}$ (12a) |  $\dfrac{\partial B_{y}}{\partial z}=\frac{B_{y}|^{n-1}_{k+\frac{3}{2}}-B_{y}|^{n-1}_{k+\frac{1}{2}}}{\Delta z}$ (12b)
---|---

$\dfrac{\partial D_{x}}{\partial z}=\frac{D_{x}|^{n-1/2}_{k}-D_{x}|^{n-1/2}_{k-1}}{\Delta z}$ (13a) |  $\dfrac{\partial D_{x}}{\partial z}=\frac{D_{x}|^{n-1/2}_{k+1}-D_{x}|^{n-1/2}_{k}}{\Delta z}$ (13b)
---|---

$B_{y}=\frac{B_{y}|^{n}_{k-1/2}+B_{y}|^{n}_{k-3/2}}{2}$ (14a) |  $B_{y}=\frac{B_{y}|^{n}_{k+3/2}+B_{y}|^{n}_{k+1/2}}{2}$ (14b)
---|---

$D_{x}=\frac{D_{x}|_{k+1}^{n-1/2}+D_{x}|_{k}^{n-1/2}}{2}$ (15a) |  $D_{x}=\frac{D_{x}|_{k+1}^{n-1/2}+D_{x}|_{k}^{n-1/2}}{2}$ (15b)
---|---

Table 1: Discretization of the dynamic terms in Eqs. (11).

Figure 3 shows the stencil of the overall proposed generalized FDTD scheme,
which provides a handy pictorial representation of Eqs. (11) with Eqs.
(12)-(15).

Figure 3: FDTD stencil corresponding to the FDTD scheme in Eqs. (11) with Eqs.
(12)-(15).

##  5 Stability

In this section, we demonstrate the stability of the generalized scheme
presented in Sec. 4 using von Neumann’s approach [45]. This approach consists
in inserting plane-wave test fields with the form

$\Psi|_{k}^{n}=\Psi_{0}\xi^{k}\zeta^{n},\quad\text{where}\quad\xi=e^{ik_{z}\Delta
z}\quad\text{and}\quad\zeta=e^{-\alpha\Delta t}$ (14)

into the update equations and determining the condition under which the
$\zeta(\xi)$ solution of the resulting matrix system is smaller than one – but
as close as possible to one to avoid numerical attenuation555In Eq. (14), one
must be carefull not to confuse $k_{z}$, the $z$-component of the wavenumber,
with $k$, the spatial index appearing in the FDTD update equations. No related
ambiguity should exist in the paper insofar as the two quantities are always
distinguishable from each other by the presence or absence of the $z$
subscript..

Inserting the test fields of the above form into Eqs. (11) with Eqs.
(12)(a)-(15)(a) or (12)(b)-(15)(b), setting the determinant of the matrix
associated with the resulting homogeneous system of equations to zero, and
finding the $\zeta$ roots of the corresponding characteristic polynomial
yields (see details in Appendix A)

$\begin{split}\zeta_{\pm}&=\frac{\xi^{-3}}{4n^{2}}\left[Sd\xi^{4}+4\left(b-4S^{2}\right)\xi^{3}+S\left(3\beta
n^{2}+2S\right)\xi^{2}-\beta n^{2}S\xi\right]\\\
&\quad\pm\frac{\xi^{-2}}{4n^{2}}(1-\xi)S\sqrt{\xi^{4}d^{2}+4\left(b-2S^{2}\right)\xi^{3}+4(n^{4}\beta^{2}-6n^{2}S\beta-2S^{2})\xi^{2}-4\beta^{2}n^{2}S\xi}\end{split},$
(15)

where $\beta=v/c$, $\xi=e^{ik_{z}\Delta z}$, $b=n^{2}(4-3\beta S)$ and
$d=2S+\beta n^{2}$. Setting $\zeta=1$ in Eq. (15) provides the stability
threshold of the scheme, with stability achieved for $\zeta<1$ (and
unstability for $\zeta>1$).

The fact that the stability analysis admits a _double solution_ ,
$\zeta_{\pm}$ [Eq. (15)], may a priori appear suprising, given that $\zeta$
solutions are usually unique. However, this can be understood by realizing
that waves that are _co-moving and contra-moving_ with respect to the
modulation see different numerical “landscapes”. Specifically, if $v>0$, then
the scheme uses Eq. (12a)-(15a) [first column in Tab. 1] for the dynamic
parts, which correspond to the dotted segments (with the dashed segment
removed) in Fig. 3; the resulting stencils for the first three equations [Eqs.
(11a)-(11c) with Eqs. (12a)-(14a)] are clearly asymmetric with respect the $z$
direction and therefore, although the stencil for the fourth equation [Eq.
(11d) with Eq. (15a)] is symmetric, the stencil is _overall asymmetric_ with
respect to $z$, so that co-moving and contra-moving waves will see different
stencils. A similar conclusion holds for the case $v<0$, so we expect that,
for a given choice of the modulation direction, $\zeta_{+}$ will correspond to
one of the two (co-moving or contra-moving) regime, while $\zeta_{-}$ will
correspond to the other regime.

The two regimes may be identified by comparing the signs of the spatial and
temporal phases of the test plane-wave in Eq. (14), with the spatial phase and
temporal phase residing in $\xi$ and $\zeta_{\pm}$, respectively. For this
purpose, let us consider a specific numerical position and time, say for
simplicity $k=1$ and $n=1$, and analyze the test plane-wave for the two
possible propagation directions. We have, upon writing
$\alpha_{\pm}=\omega_{\text{i}\pm}\pm i\omega_{\text{r}\pm}$ in $\zeta_{\pm}$
(with $\omega_{\text{i}\pm}>0$ assuming stability),

$\begin{split}\frac{\Psi|_{k=1}^{n=1}}{\Psi_{0}}=\xi\zeta_{\pm}&=e^{ik_{z}\Delta
z-\alpha_{\pm}\Delta t}=e^{ik_{z}\Delta z-(\omega_{\text{i}\pm}\pm
i\omega_{\text{r}\pm})\Delta t}\\\ &=e^{-\omega_{\text{i}\pm}\Delta
t}e^{i(k_{z}\Delta z\mp\omega_{\text{r}\pm}\Delta
t)}=e^{-\omega_{\text{i}\pm}\Delta t}e^{ik_{z}\Delta z}e^{\mp
i\omega_{\text{r}\pm}\Delta t}\\\ &=e^{ik_{z}\Delta
z}e^{-\omega_{\text{i}\pm}\Delta t}\left[\cos(\omega_{\text{r}\pm}\Delta t)\mp
i\sin(\omega_{\text{r}\pm}\Delta t)\right]=\xi\left(\zeta_{\text{r}\pm}\mp
i\zeta_{\text{i}\pm}\right).\end{split}$ (16)

In these relations, the first expression of the second line reveals that the
upper sign, which is associated with $\zeta_{+}$, corresponds to a forward
wave, propagating in the $+z$ direction, while the lower sign, which is
associated with $\zeta_{-}$, corresponds to a backward wave, so that
$\zeta_{+}$ and $\zeta_{-}$ respectively correspond to the co-moving and
contra-moving regimes for $v>0$ (and conversely for $v<0$). Note that the
scheme must be stable for both waves, since typical problems involve both of
them666This is even true for a simple mismatched moving interface, where the
incident and transmitted waves are co-moving and the reflected wave is contra-
moving if the incident wave propagates in the direction of the modulation..
The stability criterion should then be taken as that corresponding to the more
restrictive of the two $\zeta$ solutions.

We shall now show that the $\pm$ signs in Eq. (16) correspond to $\pm$ signs
in Eq. (15). We have to do this numerically because we have not been able to
separate the real and imaginary parts of $\zeta$ in Eq. (15) as we did in Eq.
(16). Inserting the parameters $S=0.5$, $\beta=0.3$, $\epsilon=4$, $\mu=1$ and
$\Delta zk_{z}=2\pi/5$ into Eq. (15) yields $\zeta_{+}=0.925-i0.33$ and
$\zeta_{-}=0.917+i0.23$, which is indeed consistent with Eq. (16), and hence
confirms the consistency of the $\pm$ indices of $\zeta$ in Eqs. (15) and
(16). Incidentally, the magnitudes of the $\zeta$’s, as apparent in the last
equation of Eqs. (14), represent the levels of the wave, which are here
$|\zeta_{+}|=0.98$ and $|\zeta_{-}|=0.95$, which reveals that the co-moving
wave is less attenuated than the contra-moving wave.

Figure 4 plots the attenuation rate $\alpha=\ln(\zeta)/\Delta t$ versus the
meshing density $N_{\lambda}=\lambda_{0}/\Delta z=2\pi/(k_{z}\Delta
z)=i2\pi/\ln(\xi)$ obtained by Eq. (15) and by FDTD simulations, along with
FDTD-simulated space-time field evolutions for two meshing densities, for the
case $v>0$ [corresponding to the stencil with Eqs. (12a)-(15a)]. The FDTD
attenuation results closesly match the attenuation results predicted by the
analytical stability analysis, with distinct levels corresponding, as
expected, to the co-moving and contra-moving regimes. Here, with the choice
$v>0$, the appropriate (more restrictive) criterion corresponds to that for
the co-moving waves ($\zeta_{+}$), whereby the contra-moving ($\zeta_{-}$)
waves will experience an extra level of numerical attenuation compared to what
would be required if they were alone. Symmetrically identical results (not
shown) are obtained for the case $v>0$, corresponding to the stencil with Eqs.
(12b)-(15b).

Figure 4: Attenuation rate versus meshing density obtained by Eq. (15) and by
FDTD simulation (taking field ratios at different times), for $S=0.5$,
$v=0.3c$ ($v>0$, first column in Tab. 1), $\epsilon=1$ and $\mu=1$. The insets
show FDTD-simulated space-time evolution of a $3\lambda$-wide Gaussian pulse
modulated by a harmonic wave of wavelength $\lambda=2\pi/k_{z}$ and period
$T=2\pi/\omega$, over a distance of $20\lambda$ and time of $45T$, for the
mesh densities $N_{\lambda}=10$ and $N_{\lambda}=80$.

##  6 Illustrative and Validating Examples

Figures 5, 6, 7 and 8 provide illustrative and validating examples,
corresponding to Figs. 1(a), (b), (c) and (d), for the problems of
electromagnetic scattering at a moving interface, a moving slab, a moving
cystal and a moving gradient, respectively.

Figure 5: Results for a contra-moving interface [Fig. 1(a)], for the physical
parameters $v=-0.3c$, $\epsilon_{1}=2$, $\epsilon_{2}=4$ and $\mu=1$, and
numerical parameters $S=0.2$ and $\Delta z=\lambda_{0}/150$. The analytical
solutions are given in Appendix B.2. (a) Space-time evolution of the electric
field for a modulated Gaussian pulse
($E_{\text{i}}=e^{-i\omega_{\text{i}}t}e^{-(t/\tau_{\text{i}})^{2}}$),
obtained by combining spatial waveforms at successive times, with input and
output temporal waveforms plotted on the sides. (b) Spectrum of the input and
output fields, obtained by Fourier transforming the waveforms sampled at the
input and output positions in (a). (c) Transmission and reflection
coefficients, obtained by taking the ratios of the Fourier transforms of the
appropriate fields for a short (wide-band) incident pulse. Figure 6: Same as
in Fig. 5 but for a co-moving slab [Fig. 1(b)], with velocity $v=0.3c$,
$\epsilon_{2}=4$ and length
$\ell=\frac{\lambda_{0}}{4n_{2}}\frac{1+n_{2}\beta}{1-n_{1}\beta}$,
corresponding to a space-time quarter-wave slab [30]. Figure 7: Same as in
Fig. 5 but for a co-moving bilayer crystal [Fig. 1(c)], consiting of 5 unit-
cells each made of two space-time quarter-wave slabs $\epsilon_{1}=1$,
$\epsilon_{2}=4$,
$\ell_{1}=\frac{\lambda_{0}}{4n_{1}}\frac{1+n_{1}\beta}{1-n_{2}\beta}$,
$\ell_{2}=\frac{\lambda_{0}}{4n_{2}}\frac{1+n_{2}\beta}{1-n_{1}\beta}$,
Figure 8: Same as in Fig. 5 but for a co-moving gradient [Fig. 1(d)],
corresponding to a linear increase of the permittivity from $\epsilon_{1}$ to
$\epsilon_{2}$ over a length of $\ell=\lambda_{0}/2$.

In all the cases (Figs. 5, 6, 7 and 8) and for all the quantities (incident,
reflected and transmitted electric fields, and reflection and transmission
coefficients), excellent agreement is observed between the numerical results
of the proposed generalized FDTD scheme and the analytical results (given in
Appendix B), with minor observable discrepancies in the scattering
coefficients, which are explained by the limited bandwidth of the test pulse.

The physics of the dynamic scattering observed in Figs. 5, 6, 7 and 8 involves
reflection Doppler frequency shifting, transmission index contrast frequency
shifting, co-moving attenuation, contra-moving amplification, multiple space-
time scattering and space-time stopbands. All these effects have been
described elsewhere [30, 34, 46] and are therefore not discussed here.

##  7 Conclusion

We have presented a generalized FDTD scheme that can simulate electromagnetic
scattering in a great diversity of moving-medium problems. This scheme, using
a combination of physical and auxiliary fields to satisfy moving boundary
conditions in a natural fashion, is both simple and powerful. It fills an
fundamental and important gap in the toolbox of electromagnetic computational
techniques and is hence expected to find wide applications, particularly in
the emergent area of space-time metamaterials.

## Appendices

##  Appendix A Stability Calculation

This section derives the stability condition for the generalized FDTD scheme,
presented in Sec. 4, using von Neumann’s approach [45]. We will therefore
replace everywhere the electromagnetic fields by their plane-wave test
counterpart, $\Psi|_{k}^{n}=\Psi_{0}\zeta^{n}e^{ik_{z}k\Delta z}$, with
$\zeta=e^{-\alpha\Delta t}$. In addition, we will use the convenient parameter
$S=\Delta t/\Delta z$ wherever appropriate.

We shall restrict our derivations to the positive-velocity regime ($v>0$),
whose dynamic part corresponds to the first column of Tab. 1. The negative-
velocity derivations ($v<0$), whose dynamic part corresponds to the second
column of Tab. 1, can easily be obtained upon following similar steps.

Using the plane-wave test field substitution in Eqs. (11a) and (12a), and
dividing the resulting expression by $\zeta^{n}$ and $e^{ik_{z}k\Delta z}$
yields

$B_{y0}e^{ik_{z}\Delta z/2}=B_{y0}e^{ik_{z}\Delta
z/2}\zeta^{-1}-SE^{*}_{x0}\zeta^{-1/2}e^{ik_{z}\Delta
z/2}\left(e^{ik_{z}\Delta z/2}-e^{-ik_{z}\Delta z/2}\right)\\\
-vSB_{y0}\zeta^{-1}\left(e^{ik_{z}\Delta z/2}-e^{-ik_{z}\Delta z/2}\right),$
(17)

which may be alternatively written, upon multiplying by $\zeta
e^{-ik_{z}\Delta z/2}$, forming sines, grouping terms and transferring
everything to the left hand-side of the equality, as

$E^{*}_{x0}\zeta^{1/2}2iS\sin(k_{z}\Delta
z/2)+B_{0}\left(\zeta-1+ve^{-ik_{z}\Delta z/2}2iS\sin(k_{z}\Delta
z/2)\right)=0.$ (18)

Similarly, Eq. (11b) with Eq. (13a) becomes

$D_{x0}\zeta^{1/2}=D_{x0}\zeta^{-1/2}-SH^{*}_{y0}\left(e^{ik_{z}\Delta
z/2}-e^{-ik_{z}\Delta z/2}\right)-vSD_{x0}\zeta^{-1/2}e^{-ik_{z}\Delta
z/2}\left(e^{ik_{z}\Delta z/2}-e^{-ik_{z}\Delta z/2}\right),$ (19)

or, multiplying by $\zeta^{1/2}$, and doing the same next operations as for
the previous equation,

$D_{x0}\left(\zeta-1+2iSve^{-ik_{z}\Delta z/2}\sin(k_{z}\Delta
z/2)\right)+H_{y0}^{*}\zeta^{1/2}2iS\sin(k_{z}\Delta z/2)=0.$ (20)

Next, Eq. (11c) with (14a) becomes, after factoring out $e^{-ik_{z}\Delta z}$
from the last two terms to form a cosine,

$E^{*}_{x0}\zeta^{1/2}=D_{x0}\zeta^{1/2}/\epsilon-\frac{v}{2}B_{y0}e^{-ik_{z}\Delta
z}(e^{ik_{z}\Delta z/2}+e^{-ik_{z}\Delta z/2}),$ (21)

or, multiplying by $\epsilon\zeta^{-1/2}$, forming cosines, etc.,

$D_{x0}-\epsilon E^{*}_{x0}-\epsilon vB_{y0}\zeta^{-1/2}e^{-ik_{z}\Delta
z}\cos(k_{z}\Delta z/2)=0.$ (22)

Finally, Eq. (11d) with (15a) becomes, after factoring out $e^{-ik_{z}\Delta
z/2}$ from the last terms to form a cosine,

$H^{*}_{y0}e^{ik_{z}\Delta z/2}=B_{y0}e^{ik_{z}\Delta
z/2}/\mu-\frac{v}{2}\zeta^{-1/2}e^{ik_{z}\Delta z/2}D_{x0}(e^{ik_{z}\Delta
z/2}+e^{-ik_{z}\Delta z/2}),$ (23)

or, multiplying by $-\mu e^{-ik_{z}\Delta z/2}$, etc.,

$B_{y0}-\mu H^{*}_{y0}-\mu v\zeta^{-1/2}D_{x0}\cos(k_{z}\Delta z/2)=0.$ (24)

Equations (18), (20), (22) and (24) form a system of equations that can be put
in the matrix form

$\begin{bmatrix}A&0&0&B+vC\\\ 0&A&B+vC&0\\\ -\epsilon&0&1&-\epsilon
vDe^{-ik_{z}\Delta z}\\\ 0&-\mu&-\mu vD&1\\\
\end{bmatrix}\begin{bmatrix}E^{*}_{x0}\\\ H^{*}_{y0}\\\ D_{x0}\\\ B_{y0}\\\
\end{bmatrix}=0,$ (25)

where

$A=2i\zeta^{1/2}S\sin(k_{z}\Delta z/2),$ (26a) $B=\zeta-1,$ (26b)
$C=2iSe^{-ik_{z}\Delta z/2}\sin(k_{z}\Delta z/2),$ (26c)
$D=\zeta^{-1/2}\cos(k_{z}\Delta z/2).$ (26d)

Setting the determinant of this homogeneous system to zero to find a
nontrivial solution for $\zeta$ leads to the characteristic polynomial

$\begin{split}\zeta^{2}&+\frac{\zeta}{2n^{2}}\left(4S^{2}-n^{2}(4-3\beta
S)-4S(S+n^{2}\beta)\cos(k_{z}\Delta z)\right)\\\ &+\frac{\zeta
S\beta}{2}\left(\cos(2k_{z}\Delta z)+2i\sin(k_{z}\Delta z)-i\sin(2k_{z}\Delta
z)\right)\\\ &+e^{-ik_{z}\Delta z}[1-S\beta(1-\cos(k_{z}\Delta
z))][S\beta+(S\beta-1)\cos(k_{z}\Delta z)-i\sin(k_{z}\Delta z)]=0.\end{split}$
(27)

This is a polynominal of the second order in $\zeta$ whose solution is given
by Eq. (15). The corresponding polynomial and solutions turn out to be exactly
identical for the negative-velocity regime, although the related derivation
involves distinct equations.

##  Appendix B Analytical Solutions

This section provides analytical solutions that are used for validation of the
illustrative examples given in Sec. 6. The solutions in Secs. B.2 and B.3 are
closed-form solutions for the canonical problems of an interface and a slab,
corresponding to Figs. 1(a) and 1(b), respectively, while the solutions in
Sec. B.4 are transfer-matrix method solutions for the more complex problems of
a bilayer crystal and a gradient, corresponding to Figs. 1(c) 1(d),
respectively.

###  B.1 General Formulas

In all cases, the incident field is a modulated Gaussian pulse, with temporal
and spectral profiles

$\psi_{\text{i}}(t)=e^{-i\omega_{\text{i}}t}e^{-(t/\tau_{\text{i}})^{2}}\quad\text{and}\quad\tilde{\psi}_{\text{i}}(\omega)=\frac{\tau_{\text{i}}}{4\sqrt{\pi}}e^{\tau_{\text{i}}^{2}(\omega-\omega_{\text{i}})/4},$
(28)

where $\tau_{\text{i}}$ is the full-width at half-maximum duration divided by
$2\sqrt{\ln 2}$ and $\omega_{\text{i}}$ is the modulation frequency. The
spectra of the reflected and transmitted pulses are obtained as

$\tilde{\psi}_{\text{t,r}}(\omega)=\\{T(a_{\text{t}}\omega),\Gamma(a_{\text{r}}\omega)\\}\frac{\tau_{\text{i}}}{4\sqrt{\pi}a_{\text{t,r}}}e^{\tau_{\text{i}}^{2}(\omega-
a_{\text{t,r}}\omega_{\text{i}})^{2}/(4a_{\text{t,r}})},$ (29)

where $T$ and $\Gamma$ are the transmission and reflection coefficients, and
$a_{\text{t,r}}$ are the transmission and reflection compression or expansion
factors. This factors are given by

$a_{\text{t}}=\frac{1-n_{\text{in}}v/c}{1-n_{\text{out}}v/c}\quad\text{and}\quad
a_{\text{r}}=\frac{1-n_{\text{in}}v/c}{1+n_{\text{in}}v/c},$ (30)

where $n_{\text{in}}$ is the refractive index of the medium in which the
incident and reflected waves propagate, and $n_{\text{out}}$ is the refractive
index of the medium in which the transmitted wave propagates777The equations
(29) and (30) may be derived by combining results in [30, 46] for
monochromatic waves; they do not seem to have been published explicitely
anywhere so far..

###  B.2 Interface

The scattering coefficients for a modulated interface are found by applying
the continuity conditions (1a) and (1b). They read [33, 30, 46], assuming wave
incidence from medium 1 towards medium 2,

$T_{21}=\frac{2\eta_{2}}{\eta_{1}+\eta_{2}}\frac{1-n_{1}v/c}{1-n_{2}v/c}\quad\text{and}\quad\Gamma_{121}=\frac{\eta_{2}-\eta_{1}}{\eta_{1}+\eta_{2}}\frac{1-n_{1}v/c}{1+n_{1}v/c},$
(31)

while the formulas for the reverse propagation direction are obtained by a
simnple interchange of the subscripts. Note that the scattering coefficients
here independent of the frequency.

###  B.3 Slab

The scattering coefficients for a modulated slab may be found, as for a
stationnary slab, either by the multiple-reflection method or the boundary-
value method [47]. This results in the scattering coefficients [30]

$T=\frac{T_{21}T_{12}e^{-i\omega_{2}^{+}n_{2}\ell}}{1-\Gamma_{212}\Gamma^{-}_{212}e^{-i(\omega_{2}^{+}+\omega_{2}^{-})n_{2}\ell}}\quad\text{and}\quad\Gamma=\Gamma_{121}\frac{1-e^{-i(\omega_{2}^{+}+\omega_{2}^{-})n_{2}\ell}}{1-\Gamma_{212}\Gamma^{-}_{212}e^{i(\omega_{2}^{+}+\omega_{2}^{-})n_{2}\ell}}.$
(32)

These relations imply different frequencies for waves propagating forward
($\omega_{2}^{+}$) and backward ($\omega_{2}^{-}$), namely

$\omega_{2}^{+}=\frac{1-n_{1}v/c}{1-n_{2}v/c}\omega_{\text{i}}\quad\text{and}\quad\omega_{2}^{-}=\frac{1-n_{1}v/c}{1+n_{2}v/c}\omega_{\text{i}},$
(33)

where $T_{21}$ and $\Gamma_{121}$ are provided in Eq. (31), $T_{12}$ and
$\Gamma_{121}$ are obtained by interchanging $1$ and $2$ in Eq. (31), and
$\Gamma_{212}^{-}$ is obtained by interchanging $1$ and $2$ and changing the
sign of $v$ in $\Gamma_{121}$ (31).

###  B.4 Transfer-Matrix Method

The transfer-matrix method [33, 19, 30] allows to compute electromagnetic
scattering for any traveling-wave modulation profile with uniform (space-time)
slabs in the frequency domain888Note that, however, the uniform slab and
frequency domain requirements of the transfer-matrix methods are very
restrictive. The former requirement precludes the simulation of a space-time
structures composed of slabs with _finite temporal durations_ (and hence
invovling space-time corners), which are obviously a basic constraint in
practical scenarios. The latter requirement essentially precludes the
determination of the spatio-temporal features of the scattered waves because
this would imply an extremely complex set of Fourier transformations. In
contrast, the proposed generalized FDTD scheme can straightforwardly simulate
truncated slabs and exact space-time waveforms.. This is even true for a
smooth gradient modulation profile, which can be approximated as a succession
of subwavelength interfaces with constant parameters.

The transfer-matrix method [33, 19, 30] models interfaces by the discontinuity
matrix

$I_{mn}=\frac{1}{t_{nm}^{-}}\begin{bmatrix}t_{mn}t_{nm}^{-}-\gamma_{mnm}\gamma_{nmn}^{-}&\gamma_{mnm}^{-}\\\
-\gamma_{nmn}&1\end{bmatrix},$ (34)

whose coefficients are given by

$t_{mn}=\frac{2\eta_{m}}{\eta_{m}+\eta_{n}}\frac{1-n_{n}v/c}{1-n_{m}v/c}\quad\text{and}\quad\gamma_{nmn}=\frac{\eta_{m}-\eta_{n}}{\eta_{m}+\eta_{n}}\frac{1-n_{n}v/c}{1-n_{n}v/c}$
(35)

and where the $t_{mn}^{-}$ and $\gamma_{nmn}^{-}$ coefficients are found by
inverting the sign of $v$ in (35) and the $t_{nm}$ and $\gamma_{nmn}$
coefficients are found by inverting $m$ and $n$ in (35). On the other hand, it
models (space-time) slabs by the propagation matrix

$P_{m}=\begin{bmatrix}e^{-i\omega_{m}^{+}n_{m}\Delta z}&0\\\
0&e^{i\omega_{m}^{-}n_{m}\Delta z}\end{bmatrix},$ (36)

where

$\omega_{\text{m}}^{+}=\frac{1-n_{1}v/c}{1-n_{\text{m}}v/c}\omega^{+}_{1},\quad\text{and}\quad\omega_{\text{m}}^{-}=\frac{1-n_{1}v/c}{1+n_{m}v/c}\omega_{1}^{+}$
(37)

are the frequencies of the forward- and backward-propagating waves within a
given slab, and $\Delta z$ is the width of this slab (which is constant under
the assumption of uniform discretization).

The transfer matrix of a structure composed of $N$ slabs is then found by
taking the product of alternating discontinuity and propagation matrices,
viz.,

$[M]_{N1}=[I_{N,N-1}][P_{N-1}][I_{N-1,N-2}][P_{N-2}]...[P_{2}][I_{2,1}].$ (38)

Finally, the scattering coefficients can be obtained by converting the
transfer matrix into a scattering matrix [48], which yields

$T=\frac{AD-BC}{D},\quad\Gamma=\frac{-A}{D},$ (39)

where $A$, $B$, $C$ and $D$ are the $(1,1)$, $(1,2)$, $(2,1)$ and $(2,2)$
components of the $[M]_{N1}$ matrix.

## References

  * [1] J. Bradley, “IV. A letter from the Reverend Mr. James Bradley Savilian Professor of Astronomy at Oxford, and FRS to Dr. Edmond Halley Astronom. Reg. &c. giving an account of a new discovered motion of the fix’d stars,” Philos. Trans. Royal Soc. 35, 637–661 (1729).
  * [2] C. Doppler, _Über das farbige Licht der Doppelsterne und einiger anderer Gestirne des Himmels_ (Calve, 1842).
  * [3] H. Fizeau, “Sur les hypothèses relatives à l’éther lumineux, et sur une expérience qui parait démontrer que le mouvement des corps change la vitesse avec laquelle la lumière se propage dans leur intérieur,” C.R. Acad. Sci. 33, 349 (1851).
  * [4] A. Einstein, “Zur elektrodynamik bewegter Körper,” Ann. Physik 322, 891–921 (1905). Reprinted in [49].
  * [5] W. C. Röntgen, “Ueber die durch Bewegung eines im homogenen electrischen Felde befindlichen Dielectricums hervorgerufene electrodynamische Kraft,” Annalen der Physik 271, 264–270 (1888).
  * [6] H. Minkowski, “Die grundgleichungen für die elektromagnetischen vorgänge in bewegten körpern,” Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse 1908, 53–111 (1908).
  * [7] U. Leonhardt and P. Piwnicki, “Optics of nonuniformly moving media,” Phys. Rev. A 60, 4301–4312 (1999).
  * [8] P. K. Tien, “Parametric amplification and frequency mixing in propagating circuits,” J. Appl. Phys. 29, 1347–1357 (1958).
  * [9] E. S. Cassedy and A. A. Oliner, “Dispersion relations in time-space periodic media: Part I-Stable interactions,” Proc. IEEE 51, 1342–1359 (1963).
  * [10] E. S. Cassedy, “Dispersion relations in time-space periodic media: Part II-Unstable interactions,” Proc. IEEE 55, 1154–1168 (1967).
  * [11] Z. Yu and S. Fan, “Complete optical isolation created by indirect interband photonic transitions,” Nat. Photonics 3, 91–94 (2009).
  * [12] N. Chamanara, S. Taravati, Z.-L. Deck-Léger, and C. Caloz, “Optical isolation based on space-time engineered asymmetric photonic band gaps,” Phys. Rev. B 96, 155409:1–12 (2017).
  * [13] Z.-L. Deck-Léger, A. Akbarzadeh, and C. Caloz, “Wave deflection and shifted refocusing in a medium modulated by a superluminal rectangular pulse,” Phys. Rev. B 97, 104305:1–7 (2018).
  * [14] A. M. Shaltout, V. M. Shalaev, and M. L. Brongersma, “Spatiotemporal light control with active metasurfaces,” Science 364, 1–10 (2019).
  * [15] M. Lampe, E. Ott, and J. H. Walker, “Interaction of electromagnetic waves with a moving ionization front,” Phys. Fluids 21, 42–54 (1978).
  * [16] N. Chamanara, Y. Vahabzadeh, and C. Caloz, “Simultaneous control of the spatial and temporal spectra of light with space-time varying metasurfaces,” IEEE Trans. Antennas Propag. 67, 2430–2441 (2019).
  * [17] P. A. Huidobro, E. Galiffi, S. Guenneau, R. V. Craster, and J. B. Pendry, “Fresnel drag in space-time modulated metamaterials,” PNAS 116, 24943–24948 (2019).
  * [18] C. Caloz, Z.-L. Deck-Léger, A. Bahrami, O. C. Vincente, and Z. Li, “Generalized space-time engineered modulation (GSTEM) metamaterials: A global and extended perspective.” IEEE Antennas Propag. Mag. (2022).
  * [19] F. Biancalana, A. Amann, A. V. Uskov, and E. P. O’Reilly, “Dynamics of light propagation in spatiotemporal dielectric structures,” Phys. Rev. E 75, 046607:1–12 (2007).
  * [20] O. Mattei and G. Milton, “Field patterns without blow up,” New J. Phys. 19, 093022:1–16 (2017).
  * [21] Y. Sharabi, A. Dikopoltsev, E. Lustig, Y. Lumer, and M. Segev, “Spatiotemporal photonic crystals,” Optica 9, 585–592 (2022).
  * [22] A. Bahrami, Z.-L. Deck-Léger, and C. Caloz, “Electrodynamics of metamaterials formed by accelerated modulation,” arXiv:2208.11506 [physics.optics] (2022).
  * [23] F. Harfoush, A. Taflove, and G. A. Kriegsmann, “Numerical implementation of relativistic electromagnetic boundary conditions in a laboratory-frame grid,” J. Comput. Phys. 89, 80–94 (1990).
  * [24] F. Harfoush, A. Taflove, and G. A. Kriegsmann, “A numerical technique for analyzing electromagnetic wave scattering from moving surfaces in one and two dimensions,” IEEE Trans. Ant. Propag. 37, 55–63 (1989).
  * [25] H. Iwamatsu and M. Kuroda, “Over set grid generation method coupled with FDTD method while considering the Doppler effect,” IEEJ Trans. Fundam. Mater. 129, 699–703 (2009).
  * [26] K. Zheng, Z. Mu, H. Luo, and G. Wei, “Electromagnetic properties from moving dielectric in high speed with Lorentz-FDTD,” IEEE Ant. Propag. Lett. 15, 934–937 (2016).
  * [27] Y. Zhao and S. Chaimool, “Relativistic finite-difference time-domain analysis of high-speed moving metamaterials,” Sci. Rep. 8, 7686:1–12 (2018).
  * [28] W. F. Ames, _Numerical Methods for Partial Differential Equations_ (Academic Press, 1977).
  * [29] A. Taflove and S. C. Hagness, _Computational Electrodynamics: The Finite-Domain Time-Domain Method_ (Artech House, 2000), 2nd ed.
  * [30] Z.-L. Deck-Léger, N. Chamanara, M. Skorobogatiy, M. Silveirinha, and C. Caloz, “Uniform-velocity spacetime crystals,” Adv. Photonics 1, 056002:1–26 (2019).
  * [31] D. Cheng and J. A. Kong, “Time-harmonic fields in source-free bianisotropic media,” J. Appl. Phys 39, 5792–5796 (1968).
  * [32] B. Bolotovski and S. Stolyarov, “Current status of the electrodynamics of moving media (infinite media),” Sov. Phys. Usp. 17, 875–895 (1975).
  * [33] C. Tsai and B. Auld, “Wave interactions with moving boundaries,” J. Appl. Phys. 38, 2106–2115 (1967).
  * [34] C. Caloz and Z.-L. Deck-Léger, “Spacetime metamaterials, part I: General concepts,” IEEE Trans. Antennas Propag. 68, 1569–1582 (2020).
  * [35] G. Milton and O. Mattei, “Field patterns: A new mathematical object,” Proc. R. Soc. A 473, 20160819:1–26 (2017).
  * [36] Z. Li, X. Ma, A. Bahrami, Z.-L. Deck-Léger, and C. Caloz, “Generalized total internal reflection at dynamic interfaces,” Phys. Rev. B 107, 115129:1–8 (2016).
  * [37] K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1966).
  * [38] A. B. Shvartsburg, “Optics of nonstationary media,” Phys.-Uspekhi 48, 797–823 (2005).
  * [39] W. Pauli, _Theory of Relativity_ (Courier Corporation, 1981).
  * [40] J. V. Bladel, _Relativity and Engineering_ (Springer Science, 1984).
  * [41] T. Shiozawa and K. Hazama, “General solution to the problem of reflection and transmission by a moving dielectric medium,” Radio Sci. 3, 569–575 (1968).
  * [42] D. Censor, “Scattering of a plane wave at a plane interface separating two moving media,” Radio Sci. 4, 1079–1088 (1969).
  * [43] J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114, 185–200 (1993).
  * [44] J. A. Kong, _Electromagnetic Wave Theory_ (EMW Publication, 2008), seventh ed.
  * [45] J. von Neumann and R. D. Richtmyer, “A method for the numerical calculation of hydrodynamic shocks,” J. Appl. Phys 21, 232–237 (1950).
  * [46] C. Caloz and Z.-L. Deck-Léger, “Spacetime metamaterials, part II: Theory and applications,” IEEE Trans. Antennas Propag. 68, 1583–1598 (2020).
  * [47] A. Sommerfeld, _Optics-Lectures on Theoretical Physics, Vol. 4_ (Academic Press, New York, 1957). Sec.7.
  * [48] M. Born and E. Wolf, _Principles of Optics_ (Pergamon, 1980).
  * [49] H. A. Lorentz, A. Einstein, H. Minkowski, H. Weyl, and A. Sommerfeld, _The Principle of Relativity: a Collection of Original Memoirs on the Special and General Theory of Relativity_ (Dover Publications, 1952).


# Evanescent-Field Assisted Photon Collection from Quantum Emitters under a
Solid Immersion Lens

S.G. Bishop  School of Engineering, Cardiff University, Queen’s Building, The
Parade, Cardiff, UK, CF24 3AA  J.K. Cannon  School of Engineering, Cardiff
University, Queen’s Building, The Parade, Cardiff, UK, CF24 3AA  H.B. Yağcı
School of Engineering, Cardiff University, Queen’s Building, The Parade,
Cardiff, UK, CF24 3AA  R.N. Clark  School of Engineering, Cardiff University,
Queen’s Building, The Parade, Cardiff, UK, CF24 3AA  J.P. Hadden  School of
Engineering, Cardiff University, Queen’s Building, The Parade, Cardiff, UK,
CF24 3AA  W.W. Langbein  School of Physics and Astronomy, Cardiff University,
Queen’s Building, The Parade, Cardiff, UK, CF24 3AA  A.J. Bennett  School of
Engineering, Cardiff University, Queen’s Building, The Parade, Cardiff, UK,
CF24 3AA<EMAIL_ADDRESS>
###### Abstract

Solid-state quantum light sources are being intensively investigated for
applications in quantum technology. A key challenge is to extract light from
host materials with high refractive index, where efficiency is limited by
refraction and total internal reflection. Here we show that an index-matched
solid immersion lens can, if placed sufficiently close to the semiconductor,
extract light coupled through the evanescent field at the surface. Using both
numerical simulations and experiments, we investigate how changing the
thickness of the spacer between the semiconductor and first lens impacts the
collection efficiency. Using automatic selection and measurement of 100s of
individually addressable colour centres in several aluminium nitride samples
we demonstrate spacer-thickness dependent photon collection efficiency
enhancement, with a mean enhancement of a factor of $4.2\text{\,}$ and a
highest measured photon detection rate of $743\pm
4\text{\,}\mathrm{k}\mathrm{c}\mathrm{p}\mathrm{s}$.

###### keywords:

Solid Immersion Lens, Photon Collection Efficiency, Aluminum Nitride, Single
Photon, Quantum Optics, Room Temperature, Colour Centres

School of Physics and Astronomy, Cardiff University, Queen’s Building, The
Parade, Cardiff, UK, CF24 3AA

##  1 Introduction

The number of photons collected from a source is of paramount importance for
any quantum optoelectronic technology. This is especially critical for quantum
light sources such as colour centres, quantum dots and single ions where
increasing the photon flux is necessary to allow measurement of correlation
effects. However, in solid state systems these emitters are embedded in high
refractive index materials such as diamond (index 2.4 @
$600\text{\,}\mathrm{nm}$), gallium arsenide (3.4), silicon carbide (2.6) or
indium phosphide (3.5). Refraction at the air-dielectric surface reduces the
solid angle over which light may be collected inside the material. At high
emission angles, light is lost at the interface due to total internal
reflection (TIR), placing an upper limit on the efficiency. Various micro-
processed architectures have been demonstrated to increase collection
efficiency (CE), including etched micro-pillars1, 2, nanowire antennae 3, 4,
photonic crystals5, 6, dielectric resonators 7, 8, and free-space cavities9,
10. All require the photonic structure to be aligned to the emitter in space
and frequency, which does not easily support the study of randomly positioned,
broadband emitters. In contrast, a millimeter-sized solid immersion lens (SIL)
11, 12, 13, a dielectric half-sphere, offers an inherently wide-field CE
enhancement over a broad spectral range. This enhancement is optimised for
buried emitters when the refractive index of the SIL is matched to the index
of the semiconductor 14, 15.

Quantum dots in gallium nitride are unique in offering quantum light emission
in the blue and UV part of the spectrum 16. Due to the recent discovery of
quantum emitters in gallium nitride (GaN)17, 18 and aluminium nitride (AlN)19,
there is renewed interest in nitrides as a platform for quantum optics in the
visible and near infra-red. Polarised quantum light emission at room
temperature, coupled with established processing technology, presents an
exciting platform combining optoelectronic components, integrated photonics,
and quantum emitters20. In addition, the III-Nitrides present one of the only
known room temperature quantum emitters in the telecom band 18.

In this report, we demonstrate how a SIL enhances the CE from colour centres
in AlN 19. We use confocal microscopy with an index matched ZrO2 hemispherical
SIL (h-SIL) ($n_{\mathrm{SIL}}=2.16$ at $600\text{\,}\mathrm{nm}$) on AlN
templates ($n_{\mathrm{AlN}}=2.15$ at $600\text{\,}\mathrm{nm}$). We engineer
the thickness of a layer of PMMA ($n_{\mathrm{PMMA}}=1.49$ at
$600\text{\,}\mathrm{nm}$) to avoid the formation of an air gap between the
semiconductor and SIL. Using finite difference time domain (FDTD) simulations,
we show that engineering the thickness of the polymer is important for
increasing the CE with a high NA lens due to evanescent coupling through the
PMMA layer.

##  2 Methods

###  2.1 An ideal hemispherical SIL

The hemispherical SIL, a truncated sphere bisected across its origin, exploits
the SIL’s geometry to increase photon CE by limiting refraction at the
semiconductor-SIL interface by reducing the index difference. Additionally,
aplanatic imaging is achieved due to the normal incidence of light at the SIL-
air interface as shown in Fig.1(a). The NA of the imaging system is therefore
increased by the refractive index within the semiconductor, reducing the axial
and lateral size of both the excitation and collection point-spread functions
by a factor of the SIL refractive index and the square of the index,
respectively, relative to imaging through a flat interface 21.

Figure 1: Hemispherical solid immersion lens. a) Illustration of the geometry.
b) Calculated (filled symbols) and simulated (empty symbols) collection
efficiency (CE) as a function of the wavelength of the dipole with (green) and
without (blue) a hemispherical SIL. c,d) XZ slices of the log of the electric
field across the dipole without and with the SIL respectively.

The CE may be analytically estimated from Barnes $et$ $al$ 14 considering a
dipole orientated in the plane of an infinite AlN slab19. The blue and green
plots with filled symbols in Fig. 1(b) demonstrate the analytically determined
CE without and with a ZrO2 SIL, respectively. The reflection at the SIL-air
interface is accounted for with Fresnel equations, with a reflectivity of 13%.
To take into account the effect of other reflections, such as that from the
AlN to sapphire interface, which are not included in the analytical analysis,
we perform Finite Difference Time Domain (FDTD) simulations in Lumerical. We
simulate the AlN layer surrounded by a sapphire substrate ($n=1.7$) and a
homogeneous ZrO2 SIL using outgoing boundary conditions (perfectly matched
layers). Anti-symmetric and symmetric boundary conditions are used across the
dipole to reduce the simulated region. The loss due to the reflection at the
SIL-air interface is included separately. The FDTD simulation predicts a
wavelength-varying and enhanced CE with regards to the analytically determined
value, due to the finite reflection from the AlN-sapphire interface, shown as
the blue and green lines in Fig.1(b) with empty symbols. This is evident from
the electric field cross-section plots through the dipole in Fig. 1(c-d). In
this scenario, the expected enhancement of light collection into a NA of 0.9
is increased by a factor of approximately 8 by the SIL, to a mean value of
$35\text{\,}\mathrm{\char 37\relax}$ across the visible spectrum.

###  2.2 Ensemble scanning

To experimentally quantify the enhancement in CE, we use statistical analysis
over an ensemble of individually addressable emitters. To reduce bias in the
selection of the emitters, we developed an automated routine to locate and
characterise emitters. Initially, confocal sample maps were measured revealing
the presence of emitters in the AlN epilayer using an optical dual-axis 4f
scanning system 19. Optical filtering was used to measure only the
$550\text{\,}$ to $650\text{\,}\mathrm{nm}$ spectral range, suppressing
fluorescence from chromium impurities in the sapphire substrate. Due to the
linearly polarised nature of the excitation laser at $532\text{\,}\mathrm{nm}$
and the fact that each emitter has a linearly-polarised excitation dipole with
unknown orientation 19, three scans maps were taken with the laser’s linear
polarisation angle at $0\text{\,}\mathrm{\SIUnitSymbolDegree}$,
$60\text{\,}\mathrm{\SIUnitSymbolDegree}$ and
$120\text{\,}\mathrm{\SIUnitSymbolDegree}$, shown in Fig.2(a-c). We combine
these scan maps in (d), which is the mean of scans (a-c), to avoid pre-
selecting emitters aligned to the laser polarisation. As reported in our
previous work, AlN colour centres in this samples have a zero-phonon line at
$595\text{\,}\mathrm{nm}$ with a broad phonon side-band extending to
$700\text{\,}\mathrm{nm}$ 19.

Figure 2: Automated selection of emitters in confocal scan maps. a,b,c)
Excitation polarisation resolved confocal scans, where the excitation laser
polarisation is rotated in the plane of the sample at an angle of
$0\text{\,}\mathrm{\SIUnitSymbolDegree}$(a),
$60\text{\,}\mathrm{\SIUnitSymbolDegree}$(b) and
$120\text{\,}\mathrm{\SIUnitSymbolDegree}$(c) respectively. d) Combined scan
map where each pixel represents the mean of the 3 polarisation resolved scan
maps. e) Emitter locations determined with the fast peak finding technique. d)
Corresponding filtered emitter locations after Gaussian fitting to each
emitter’s point-spread function.

The resultant scan map was then fed into an algorithm that locates emitters in
the following manner. Initially, the algorithm removes single pixel noise by
using a $3\times 3$ pixel median filter and further smooths the image by
convolving the image with a two-dimensional Gaussian function. A fast pixel-
by-pixel analysis then locates potential emitter locations based on the
intensity of neighbouring pixels. The 126 candidate emitter locations are
marked in Fig.2(e). It can be observed in the image that the algorithm
overestimates the number of emitters and occasionally selects features smaller
than the diffraction limited lateral spot size. Therefore, a 2-dimensional
Gaussian function is fit to each emitter. The Gaussian fit parameters are also
used to rank the emitters based on the Gaussian width being close to the
expected diffraction limited point-spread function of the microscope, symmetry
in X and Y, amplitude and goodness of fit $R^{2}$ value. This Gaussian fit
identifies the centre of the emitters’ point spread function with greater
precision than the scan step-size. Post filtering, 92 emitters are identified
in the scan map Fig.2(f).

Power dependent saturation measurements are then taken on each candidate
emitter sequentially using the following automated process. Firstly, the
scanning system moves to the approximate emitter location. X,Y and Z line
scans are performed to find the maximum along each axis. The linear
polarisation of the $532\text{\,}\mathrm{nm}$ excitation laser is rotated to
maximally couple to the emitter’s absorption dipole. Then a laser-power $P$
dependent intensity $I(P)$ was measured and used to quantify the fitted
saturation intensity at infinite power $I_{\infty}$ as defined by the fit
function $I(P)=I_{\infty}P/(P+P_{sat})$. Measurements with a saturation power
greater than $2\text{\,}\mathrm{mW}$ are removed from each data-set as this is
a good proxy for linear power dependent behaviour, which is not consistent
with a single colour centre (this occurs with less than
$5\text{\,}\mathrm{\char 37\relax}$ frequency). No polarisation optics were
included in the collection path of our microscope. Throughout the
measurements, photo-bleaching was not observed for any emitters but random
telegraph noise in the detected count rate is common.

##  3 Results and Discussion

###  3.1 Collection efficiency enhancement

Using the automated system we set out to determine the efficiency enhancement
for an ensemble of emitters under a SIL. We emphasise that the presence of an
air-gap between the SIL and AlN suppresses such an enhancement. An air gap can
arise from a non-flat back surface of the SIL and/or a non-flat surface of the
semiconductor. To prevent such a gap, we add a thin layer of a polymer, PMMA,
which is soft enough to fill the gap between the SIL and AlN adapting to their
surfaces, in line with other reports 12, 21. Three samples of AlN-on-sapphire
from the same wafer were studied. One with no SIL (SP1), one with
$212\text{\,}\mathrm{nm}$ of PMMA between the SIL and AlN (SP2) and one with
$70\text{\,}\mathrm{nm}$ of PMMA (SP3). The PMMA thickness was determined
using an ellipsometer prior to the addition of the SIL.

The statistical distribution of emitter saturation intensities in each sample
is presented in Fig.3. In sample SP1, without the SIL, 396 emitters were
measured and the resulting distribution is shown in Fig.3(a). The mean of the
distribution is given as
$\mu_{\mathrm{SP1}}$=$68\text{\,}\mathrm{k}\mathrm{c}\mathrm{p}\mathrm{s}$.
Although a similar mean intensity to the NV centre in diamond, the variance in
the saturation intensities differs significantly to what is seen in the NV
centre in diamond where a narrow distribution is often reported 22, which can
be correlated with the NV orientation. We attribute the larger variation in
emitter brightness to differences in the internal energy levels between
emitters which periodically suppress emission on different timescales. For
example, a metastable shelving state 19 evident in the photon statistics of
nearly all emitters, leads to pump power dependent bunching of emission on
timescales of tens to hundreds of nanoseconds. Additionally, some emitters
show telegraph noise which is observed on timescales of seconds, but may also
be present on millisecond or microsecond scales. The relationship between
these instabilities and the time-averaged intensity will be the subject of a
future study. However, it is notable that the predicted photon CE should lead
to a saturated emission intensity at least an order of magnitude greater than
we measure, based on the extrapolated zero-pump power lifetime
$t_{1}=$4.6\text{\,}\mathrm{ns}$$ as determined from fitting the second order
correlation measurements in Fig.3(f).

In SP1, 396 emitters were sampled over an area of
$150\text{\,}\mathrm{\SIUnitSymbolMicro m}$ by
$150\text{\,}\mathrm{\SIUnitSymbolMicro m}$, corresponding to an area of
$57\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}$ per emitter. Due to the
magnification of the SIL, where the displacement in X and Y of our imaging
system leads to a displacement of $X/n_{\mathrm{SIL}}$ and
$Y/n_{\mathrm{SIL}}$ on the sample surface respectively, we ensure we sample a
similar number of emitters per unit area scanned for SP2 and SP3. 63 and 70
emitters are measured for SP2 and SP3 respectively, with the corresponding
distribution in $I_{\infty}$ shown in in Fig.3(b,c).

The mean of the saturated intensity $I_{\infty}$ distribution in each sample
is shown as the dashed vertical lines in Fig.3(a-c). For SP2 and SP3, the mean
$I_{\infty}$ is $211\text{\,}$ and
$285\text{\,}\mathrm{k}\mathrm{c}\mathrm{p}\mathrm{s}$, respectively. This
corresponds to an enhancement factor of $3.1\text{\,}$ and $4.2\text{\,}$,
relative to SP1. Assuming that the mean of the distribution for SP1 represents
a CE of $4.6\text{\,}\mathrm{\char 37\relax}$, as determined from the
analytical analysis, the projected CE for SP2 and SP3 can be estimated to be
$14.3\text{\,}\mathrm{\char 37\relax}$ and $19.3\text{\,}\mathrm{\char
37\relax}$. We note that the highest detected count rate is measured is an
emitter in SP3 with $I_{\infty}=$743\pm
4\text{\,}\mathrm{k}\mathrm{c}\mathrm{p}\mathrm{s}$$. We show data for this
emitter in Fig.3, labelled in (c) with a star. Pump power dependent auto-
correlation data in (f) shows the presence of anti-bunching near time zero, in
combination with longer timescale bunching at high pump powers. A fit to the
$4\text{\,}\mathrm{\SIUnitSymbolMicro m}$ measurement, not shown in the
figure, illustrates the single photon purity, with $g^{(2)}(0)=0.12$. In the
experimental data we have removed an artefact at 8-$30\text{\,}\mathrm{ns}$
caused by optical cross-talk in our interferometer.

Figure 3: Evanescent-field assisted collection by a hemispherical SIL. a-c)
Histogram representing the fitted saturation intensity at infinite power
power, $I_{\infty}$, for a number of emitters in samples SP1, SP2 and SP3. The
polymer thickness is labelled on each plot. For the brightest emitter in SP3
(marked with a red star) we show the power-dependent intensity saturation
curve (d), the intensity as a function of excitation polarisation (e) and
power dependent auto-correlation histograms (f).

###  3.2 Evanescent coupling

To investigate the observed CE dependence on the polymer thickness we return
to FDTD. Using the aforementioned simulation environment the electric field
from a single dipole embedded in the middle of the AlN is captured on the
surface of a box of monitors enclosing the emitter and spanning the dielectric
stack. The near-field electric field on the XY monitor in the ZrO2 layer is
projected into the far field, where the flux is integrated over the solid
angles of the collection optics to determine the collected power. The CE is
defined as the collected power normalised by the power injected into the
simulation via the single dipole.

Figure 4: Simulated CE demonstrating evanescent coupling through a thin
polymer layer. a) CE as a function of the collection free space NA and the
thickness $d$ of the PMMA polymer. b) CE as a function of the polymer
thickness for NA$=$0.9\text{\,}$$.

The CE as a function of the polymer thickness $d$ and NA of the imaging system
is presented in Fig.4. The simulation shows that for
$d>$300\text{\,}\mathrm{nm}$$ the CE is constant and is well below the
expected CE for an ideal system with NA$=$0.9\text{\,}$$, as determined in
Fig.1. This is due to total internal reflection at the AlN-PMMA interface at
angles greater than $44\text{\,}\mathrm{\SIUnitSymbolDegree}$, limiting the
effective collection NA to 0.69.

The greatest gain in CE can be made by combining a high NA lens with a thin
PMMA layer. The CE for a NA of 0.9 as a function of the polymer thickness is
presented in Fig.4(b). We attribute the increased CE for polymer thickness
less than $200\text{\,}\mathrm{nm}$ to evanescent coupling through the polymer
for light incident at angles greater than $44\text{\,}$ degrees. The
evanescent wave in the polymer is converted back to a propagating wave at the
polymer-SIL interface and collected by our imaging system, a phenomenon known
as frustrated total internal reflection. The electric field in the polymer
decays with the distance $z$ away from the interface (in the absence of the
SIL) in the form,

$E(z)=E(s)e^{-2z/d}\quad\mathrm{and,}\quad
p=\frac{\lambda_{\mathrm{Air}}/n_{\mathrm{AlN}}}{2\pi\sqrt{\left(\frac{n_{\mathrm{AlN}}}{n_{\mathrm{PMMA}}}\right)^{2}\sin^{2}(\theta)-1}}$
(1)

where $E(z)$ is the electric field normal to the interface, $E(s)$ is the
electric field at the surface, $\theta$ is the incidence angle of the light,
$\lambda_{0}$ is the wavelength of the light in a vacuum and $p$ is the
penetration depth of the evanescent wave.

Due to the angle dependence of the coupling through the polymer, we consider
the two limiting angles of the system, just beyond the angle of TIR,
$\theta_{\mathrm{TIR}}=$45\text{\,}\mathrm{\SIUnitSymbolDegree}$$, and the NA
of the imaging system,
$\theta_{\mathrm{NA}}=$65\text{\,}\mathrm{\SIUnitSymbolDegree}$$. The
penetration depths as determined using Eq.1 are given as
$p_{\mathrm{TIR}}=$228\text{\,}\mathrm{nm}$$ and
$p_{\mathrm{NA}}=$53\text{\,}\mathrm{nm}$$, respectively. The determined
penetration depths agree well with the simulated thickness dependent
collection efficiency in Fig.4, where for $d<$200\text{\,}\mathrm{nm}$$ the CE
is increased above the vanishing value at infinite polymer thickness. A fit to
the data in Fig.4(b) using an offset exponential is used to interpolate the
expected CE for SP2 and SP3. The determined CE values of $15\text{\,}$ and
$27\text{\,}\mathrm{\char 37\relax}$ agree well with the projected values from
the statistical measurements.

##  4 Conclusion

The photon collection efficiency of quantum emitters in a AlN-on-sapphire
substrate was enhanced by engineering the thickness of a PMMA adhesion layer
under a hemispherical solid immersion lens. An automated emitter finding
routine was used to determine the collection efficiency enhancement of an
ensemble. We demonstrate experimentally and numerically that at high NA and
low PMMA thickness evanescent coupling through the polymer can maximise CE. A
collection efficiency of $19.3\text{\,}\mathrm{\char 37\relax}$ is inferred
for a polymer thickness of $70\text{\,}\mathrm{nm}$, in good agreement with
numerical simulations. The significant increase in collection efficiency from
these room temperature quantum emitters is a step towards efficient off-the-
shelf room temperature quantum light sources.

## Acknowledgements

We acknowledge financial support provided by EPSRC via Grant No. EP/T017813/1
and the European Union’s H2020 Marie Curie ITN project LasIonDef (GA No.
956387). Device processing was carried out in the cleanroom of the ERDF-funded
Institute for Compound Semiconductors (ICS) at Cardiff University.

## References

  * Santori et al. 2002 Santori, C.; Fattal, D.; Vučković, J.; Solomon, G. S.; Yamamoto, Y. Indistinguishable photons from a single-photon device. _Nature_ 2002, _419_ , 594–597
  * Somaschi et al. 2016 Somaschi, N. et al. Near-optimal single-photon sources in the solid state. _Nat. Photonics_ 2016, _10_ , 340–345
  * Babinec et al. 2010 Babinec, T. M.; Hausmann, B. J. M.; Khan, M.; Zhang, Y.; Maze, J. R.; Hemmer, P. R.; Lončar, M. A diamond nanowire single-photon source. _Nat. Nanotechnol._ 2010, _5_ , 195–199
  * Claudon et al. 2010 Claudon, J.; Bleuse, J.; Malik, N. S.; Bazin, M.; Jaffrennou, P.; Gregersen, N.; Sauvan, C.; Lalanne, P.; Gérard, J.-M. A highly efficient single-photon source based on a quantum dot in a photonic nanowire. _Nat. Photonics_ 2010, _4_ , 174–177
  * Noda et al. 2007 Noda, S.; Fujita, M.; Asano, T. Spontaneous-emission control by photonic crystals and nanocavities. _Nat. Photonics_ 2007, _1_ , 449–458
  * Lodahl et al. 2015 Lodahl, P.; Mahmoodian, S.; Stobbe, S. Interfacing single photons and single quantum dots with photonic nanostructures. _Rev. Mod. Phys._ 2015, _87_ , 347–400
  * Kuznetsov et al. 2016 Kuznetsov, A. I.; Miroshnichenko, A. E.; Brongersma, M. L.; Kivshar, Y. S.; Luk’yanchuk, B. Optically resonant dielectric nanostructures. _Science (80-. )._ 2016, _354_
  * Mignuzzi et al. 2019 Mignuzzi, S.; Vezzoli, S.; Horsley, S. A.; Barnes, W. L.; Maier, S. A.; Sapienza, R. Nanoscale Design of the Local Density of Optical States. _Nano Lett._ 2019, _19_ , 1613–1617
  * Tomm et al. 2021 Tomm, N.; Javadi, A.; Antoniadis, N. O.; Najer, D.; Löbl, M. C.; Korsch, A. R.; Schott, R.; Valentin, S. R.; Wieck, A. D.; Ludwig, A.; Warburton, R. J. A bright and fast source of coherent single photons. _Nat. Nanotechnol._ 2021, _16_ , 399–403
  * Dolan et al. 2010 Dolan, P. R.; Hughes, G. M.; Grazioso, F.; Patton, B. R.; Smith, J. M. Femtoliter tunable optical cavity arrays. _Opt. Lett._ 2010, _35_ , 3556
  * Born and Wolf 1999 Born, M.; Wolf, E. _Princ. Opt. Electromagn. Theory Propag. Interf. Diffr. Light_ ; 1999
  * Vamivakas et al. 2007 Vamivakas, A. N.; Atatüre, M.; Dreiser, J.; Yilmaz, S. T.; Badolato, A.; Swan, A. K.; Goldberg, B. B.; Imamoǧlu, A.; Ünlü, M. S. Strong extinction of a far-field laser beam by a single quantum dot. _Nano Lett._ 2007, _7_ , 2892–2896
  * Chen et al. 2018 Chen, Y.; Zopf, M.; Keil, R.; Ding, F.; Schmidt, O. G. Highly-efficient extraction of entangled photons from quantum dots using a broadband optical antenna. _Nat. Commun._ 2018, _9_ , 1–7
  * Barnes et al. 2002 Barnes, W.; Björk, G.; Gérard, J.; Jonsson, P.; Wasey, J.; Worthing, P.; Zwiller, V. Solid-State Single Photon Sources : Light Collection Strategies. _Eur. Phys. J. D - At. Mol. Opt. Plasma Phys._ 2002, _210_ , 197–210
  * Serrels et al. 2008 Serrels, K. A.; Ramsay, E.; Dalgarno, P. A.; Gerardot, B.; O’Connor, J. A.; Hadfield, R. H.; Warburton, R. J.; Reid, D. T. Solid immersion lens applications for nanophotonic devices. _J. Nanophotonics_ 2008, _2_ , 021854
  * Zhu and Oliver 2016 Zhu, T.; Oliver, R. A. Nitride quantum light sources. _EPL (Europhysics Lett._ 2016, _113_ , 38001
  * Berhane et al. 2017 Berhane, A. M.; Jeong, K.-Y.; Bodrog, Z.; Fiedler, S.; Schröder, T.; Triviño, N. V.; Palacios, T.; Gali, A.; Toth, M.; Englund, D.; Aharonovich, I. Bright Room-Temperature Single-Photon Emission from Defects in Gallium Nitride. _Adv. Mater._ 2017, _29_ , 1605092
  * Zhou et al. 2018 Zhou, Y.; Wang, Z.; Rasmita, A.; Kim, S.; Berhane, A.; Bodrog, Z.; Adamo, G.; Gali, A.; Aharonovich, I.; Gao, W.-b. Room temperature solid-state quantum emitters in the telecom range. _Sci. Adv._ 2018, _4_ , eaar3580
  * Bishop et al. 2020 Bishop, S. G.; Hadden, J. P.; Alzahrani, F. D.; Hekmati, R.; Huffaker, D. L.; Langbein, W. W.; Bennett, A. J. Room-Temperature Quantum Emitter in Aluminum Nitride. _ACS Photonics_ 2020, acsphotonics.0c00528
  * Wan et al. 2020 Wan, N. H.; Lu, T. J.; Chen, K. C.; Walsh, M. P.; Trusheim, M. E.; De Santis, L.; Bersin, E. A.; Harris, I. B.; Mouradian, S. L.; Christen, I. R.; Bielejec, E. S.; Englund, D. Large-scale integration of artificial atoms in hybrid photonic circuits. _Nature_ 2020, _583_ , 226–231
  * Bishop et al. 2022 Bishop, S. G.; Hadden, J. P.; Hekmati, R.; Cannon, J. K.; Langbein, W. W.; Bennett, A. J. Enhanced light collection from a gallium nitride color center using a near index-matched solid immersion lens. _Appl. Phys. Lett._ 2022, _120_ , 114001
  * Babinec et al. 2010 Babinec, T. M.; Hausmann, B. J.; Khan, M.; Zhang, Y.; Maze, J. R.; Hemmer, P. R.; Lončar, M. A diamond nanowire single-photon source. _Nat. Nanotechnol._ 2010, _5_ , 195–199