# Long-range interactions of quantum emitters coupled through structured baths

Lewis Ruks<EMAIL_ADDRESS>Thomas Busch  Quantum Systems Unit, Okinawa
Institute of Science and Technology Graduate University,
Onna, Okinawa 904-0495, Japan

###### Abstract

Quantum emitters coupled to structured baths are the ideal testbed for
investigating atypical transport stemming from exotic wave propagation
dynamics. In this work, we characterize wave propagation in generic Hermitian
local periodic baths in $d$ dimensions, and investigate the combined effects
of anisotropy and quasi-breaking of periodicity. We derive an asymptotic
decomposition of the Green’s function into resonant travelling waves obeying
power-law decay, and rapidly decaying evanescent waves. The far-field
transport is characterized in terms of the transverse effective mass tensor,
and beyond caustics a universal scaling for the exponential decay length of
ghost waves is obtained. Finally applying an artifical gauge field on the
lattice, we find emission resonant with the corresponding open orbits is
robustly quasi-one dimensional with topological protection. Our results shed
fresh light on wave propagation in discrete structured media.

††preprint: APS/123-QED

Introduction —The paradigmatic setting underpinning quantum optics is that of
quantum emitters (QEs) coupled to a (photonic) bath. Since the discovery of
the Purcell effect led to the general realization that modifying the bosonic
mode structure through a cavity effects the emission profile in the weak-
coupling regime Purcell (1995), structured baths have received a great deal of
attention from both the classical Akahane  _et al._ (2003); Khelif  _et al._
(2004); Krauss (2007) and quantum Rao and Hughes (2007); Yoshie  _et al._
(2004); Nakamura  _et al._ (2004) communities. Spurred on by recent
technological advances in coupling at the micro- and nano-scale Kotsifaki and
Chormaic (2019); Meng  _et al._ (2018); Yu  _et al._ (2019); Ziegler  _et al._
(2019); Yu  _et al._ (2017); Solano  _et al._ (2017), many theoretical
proposals for realising exotic phenomena by coupling QEs to structured baths
have been discussed particularly in the last decade Perczel  _et al._ (2020);
Hood (2017); García-Elcano  _et al._ (2020); Zhou  _et al._ (2017); Sheremet
_et al._ (2021); González-Tudela and Cirac (2017); Ying  _et al._ (2019);
González-Tudela  _et al._ (2015), offering potential applications in quantum
information Shahmoon and Kurizki (2013); Van Loo  _et al._ (2013a); Paulisch
_et al._ (2016) and simulation of many-body physics Scully and Svidzinsky
(2009); Douglas  _et al._ (2015); De Bernardis  _et al._ (2021); Argüello-
Luengo  _et al._ (2019); Kim  _et al._ (2021). Central to all dynamics is the
dispersion relation and dimensionality of the underlying bath, which
determines the emission profile of a weakly coupled QE in space.

Generically, resonant emission into the band of a periodic bath is understood
to consist of a travelling part and a rapidly decaying evanescent part Ying
_et al._ (2019), whose natures are determined by the geometry of the bath
dispersion relation at the resonant frequency Zhou  _et al._ (2017); Yuan  _et
al._ (2019); González-Tudela and Cirac (2018); García-Elcano  _et al._ (2020).
In the one-dimensional case, or for isotropic higher-dimensional cases, the
Green’s function, which determines the dynamics, may be decomposed exactly
within a tight-binding formulation into contributions from travelling waves
(real poles) and evanescent waves (complex poles). Excepting specific analyses
Ying  _et al._ (2019); Kim  _et al._ (2021); Yu  _et al._ (2017); González-
Tudela and Cirac (2017), an exact characterisation of the travelling part in
generic higher-dimensional environments is currently unknown. In this work, we
present a comprehensive picture of the asymptotic dynamics of linear wave
propagation in local periodic baths and illustrate fundamental results using
examples in the square tight-binding lattice. We use differential geometry to
prove an asymptotic decomposition of wave propagation in a generic local
$d$-dimensional ($d$D) bath as a sum of resonant contributions that follow a
power-law decay, and evanescent waves that decay faster than any polynomial.
Our analysis yields an approximation of the Green’s function as an integral
over the resonant level set (Fermi surface) that converges up to exponential
errors with increasing source-receiver separation. Through the geometry of
resonant level sets, we show how excitations of diverging transverse mass
produce caustics beyond which interactions are evanescently decaying, yet
oscillating, ‘ghost waves’ Kuś  _et al._ (1993); Narimanov (2019); Newman  _et
al._ (2018). A universal scaling of ghost wave decay with respect to the off-
caustic angle can be derived using saddle point analysis, and deep into the
ghost wave regime we find that spin-spin exchanges effected through the bath
may be purely incoherent. We probe the associated open orbits Danner  _et al._
(1994) in quasimomentum space via an effective magnetic field on the lattice.
Quasi-1D emission with non-zero transport and macroscopic bath dependence may
be realised, in contrast to isotropic magnetic orbits De Bernardis  _et al._
(2021). Simulating out-of-equilibrium dynamics, we find the emission to
periodically and robustly refocus down to the single site level in the plane
transverse to propagation. Consequently, QEs located at the refocusing points
strongly interact to form high-fidelity bound states in the continuum Hsu  _et
al._ (2016); Chang  _et al._ (2012); Sinha  _et al._ (2020) using only a few
QEs. The topological protection of open orbits enables preservation of
transport in the presence of local obstructions and moderate global disorder.

Formalism —We consider a local and periodic $d$D tight-binding Hamiltonian
$\hat{H}$, which is linear in the bosonic operators $\hat{a}_{i}^{(\dagger)}$,
so that (with $\hbar=1$)
$\hat{H}=\sum_{ij}J_{ij}\hat{a}_{i}^{\dagger}\hat{a}_{j}$ with
$[\hat{a}_{i}^{\dagger},\hat{a}_{j}]=\delta_{ij}$ and each site is coupled to
only finitely many neighbours. Propagation in periodically structured baths
can be analysed through the Green’s function $\hat{{G}}(z)=(z-\hat{H})^{-1}$
with off-diagonal elements given by
${G}(\mathbf{r},\mathbf{r}^{\prime},z)=\langle\mathbf{r}|\hat{G}|\mathbf{r}^{\prime}\rangle$
and expressed as an integral over the Brillouin zone (BZ) of the momentum
space kernel when $\hat{H}$ is hermitian

${G}(\bm{\rho},\Delta)=\mathcal{A}\sum_{\nu}\int_{\text{BZ}}\frac{d^{d}\mathbf{k}}{(2\pi)^{d}}\frac{\psi_{\nu}(\mathbf{k},\mathbf{r})\psi_{\nu}^{*}(\mathbf{k},\mathbf{r^{\prime}})}{\Delta-\omega_{\nu}(\mathbf{k})+i0}.$
(1)

Here $\mathcal{A}$ is the volume of a unit cell of the structure,
$\omega_{\nu}(\mathbf{k})$ describes the dispersion of band $\nu$, and
$\psi_{\nu}(\mathbf{k},\mathbf{r})=\langle\mathbf{r}|\psi_{\nu}(\mathbf{k})\rangle$
is the wavefunction in position space for the diagonal Bloch wave
$|\psi_{\nu}(\mathbf{k})\rangle$. Translational invariance gives dependence
only on the separation vector $\bm{\rho}=\mathbf{r}-\mathbf{r}^{\prime}.$ In
the following, we exploit linearity to treat only a single band, and drop the
band subscript $\nu$. Away from a Van Hove singularity, the Sokhotski-Plemelj
theorem permits a decomposition of $G$ into (in)coherent terms mediated
through off-resonant (resonant) excitations as
$G={\Omega}-\frac{i}{2}{\Gamma}$ with

$\displaystyle{\Gamma}$
$\displaystyle=-2\mathfrak{I}[G]=\mathcal{A}\int_{S}\frac{d^{d-1}\mathbf{k}}{(2\pi)^{d-1}}\frac{\psi(\mathbf{k},\mathbf{r})\psi^{*}(\mathbf{k},\mathbf{r^{\prime}})}{{v}(\mathbf{k})},$
(2) $\displaystyle{\Omega}$
$\displaystyle=\mathfrak{R}[G]=\mathcal{A}\mathcal{P}\int_{\text{BZ}}\frac{d^{d}\mathbf{k}}{(2\pi)^{d}}\frac{\psi(\mathbf{k},\mathbf{r})\psi^{*}(\mathbf{k},\mathbf{r^{\prime}})}{\Delta-\omega(\mathbf{k})},$
(3)

where $S=S(\Delta)=\\{\mathbf{k}:\omega(\mathbf{k})=\Delta\\}$ is the resonant
set at $\Delta$ and $\nabla\omega(\mathbf{k})=\mathbf{v}(\mathbf{k})$ is the
group velocity with $v=|\mathbf{v}|$. Under our assumptions (see Appendix
A.1), $S=S(\Delta)$ is a (union of) smooth $(d-1)D$ dimension manifold(s) Zhao
and Wang (2013). In artifical light-matter systems $S$ is highly controllable
due to the tunability of $\Delta$ through the QEs Zhou  _et al._ (2017);
González-Tudela  _et al._ (2015); Hood (2017); Perczel  _et al._ (2020);
Douglas  _et al._ (2015); Argüello-Luengo  _et al._ (2019), and can be probed
using the same QEs. We consider QEs of resonant frequency $\Delta,$ each
coupled to a single site of the bath $\mathbf{r}_{i}$ with strength
$|g|\ll|J_{ij}|$, the weakly coupled regimeGonzález-Tudela and Cirac (2018);
García-Elcano  _et al._ (2020). In the thermodynamic bath limit one makes the
Born-Markov approximation to obtain the effective spin-spin coupling elements
between spins $i$ and $j$ of the standard master equation as Gardiner  _et
al._ (2004); García-Elcano  _et al._ (2020)

$H_{ij}=g^{2}G(\bm{\rho}_{ij},\Delta).$ (4)

The emitter dynamics are therefore entirely determined through $G,$ and
chiefly through $\omega(\mathbf{k})$ and $\Delta.$ Whilst we present our
results in the context of QEs coupled to tight-binding lattices, we note that
qualitative phenomena are observable Markos and Soukoulis (2008) and
systematically reproducible Matlack  _et al._ (2018) in continuum media,
including ultracold gases Köhl  _et al._ (2005); Argüello-Luengo  _et al._
(2019), photonic crystals Lidorikis  _et al._ (1998); Bayindir  _et al._
(2000); Cerjan  _et al._ (2020), and acoustic meta-crystals Xiao  _et al._
(2015); Xue  _et al._ (2019). In particular we emphasize that our
investigation applies generally to waves in linear and periodic media.

In the current study we investigate transport of emission in a one-band tight-
binding anisotropic square lattice. The nearest neighbour coupling strength is
given by $J_{x(y)}$ for coupling in the $x(y)$ direction along
$\mathbf{a}_{x}=a\hat{\mathbf{x}}$ $(\mathbf{a}_{y}=a\hat{\mathbf{y}})$ as in
Figure 2(b). In the Bloch basis we have
$\omega(\mathbf{k})=2(J_{x}\cos(k_{x}a)+J_{y}\cos(k_{y}a)).$ We take
$-\Delta=J=J_{x}=J_{y}$ well away from the Van Hove singularity at
$\mathbf{k}=\mathbf{0}$ or at the edges of the BZ (Figure 1(a)).

Figure 1: (a) Dispersion relation of the isotropic square tight-binding
lattice in the first BZ. White diamonds indicate Van Hove singularities
$v(\mathbf{k})=0$, whilst the translucent red region gives a tube of radius
$\epsilon$ about $S$, bounded by dashed lines. (b) Normalized coherent part
$J\Omega$ calculated exactly (Appendix A.3) together with the approximation
(A) given through the real part of equation (5) together with the classical
far-field (F) prediction given through equation (7) as
$\bm{\rho}=n\hat{\bm{\rho}}$ is increased along the square diagonal,
$\hat{\bm{\rho}}=\mathbf{a}_{x}+\mathbf{a}_{y}.$ Inset shows relative error
$E=|\Omega-\Omega_{i}|/\Omega,$ ($i=A,F$) on a log scale.

Asymptotic decomposition —The off-resonant interactions (3) are given through
a strongly singular continuum sum that has not received an analytical
treatment to our knowledge, so that their general role in long-range
interactions is unclear beyond one dimensionHood (2017); Sánchez-Burillo  _et
al._ (2020). However, a treatment as a rapidly oscillating integral is
permitted as $\bm{\rho}$ is increased commensurate with the lattice in a fixed
direction. Standard methods of complex analysis used in 1D Katznelson (2004);
Hood (2017); Rao and Hughes (2007) or employed with Wannier functions Brouder
_et al._ (2007); Katznelson (2004) can not be readily applied to the strongly
singular manifold in (3) when $d>1$. However, when $S$ is smooth, results of
differential geometry can be employed. Using a partition of unity, we may
localize the contributions near (away from) $S,$ where the integrand is
strongly singular (analytic), as in Figure 1(a). We apply Weyl’s tube
formulaAlvarez (2017); Gray (2012); Khudaverdian  _et al._ (2007); Choe and
Koiso (2015)to asymptotically obtainsup  the dominant contributions from a
tubular region of $S$. Up to evanescent corrections eva , we obtain an
approximation for $\Omega$, and thus ${G}$, as a regular sum of waves resonant
as $S$ in the limit of large $\rho=|\bm{\rho}|$:

$G\sim-i\mathcal{A}\int_{S}\frac{d^{d-1}\mathbf{k}}{(2\pi)^{d-1}}\Pi_{\rho}(\hat{\mathbf{v}}\cdot\hat{\bm{\rho}})\frac{\psi(\mathbf{k},\mathbf{r})\psi^{*}(\mathbf{k},\mathbf{r^{\prime}})}{{v}(\mathbf{k})}.$
(5)

Here $\hat{\mathbf{v}}$ ($\hat{\bm{\rho}}$) is the normalized group velocity
(separation) direction $\mathbf{v}/v$ ($\bm{\rho}/\rho$) at $\mathbf{k}\in S$,
and $\Pi_{\rho}$ is an amplitude function, defined in sup , that smoothly
varies between 0 and 1 as the angle $\hat{\mathbf{v}}\cdot\hat{\bm{\rho}}$
varies. Crucially, $\Pi_{\rho}(t)\to H(t)$ as $\rho\to\infty$, where $H$ is
the Heaviside function. With the decomposition
$\Pi_{\rho}=\frac{1+\Theta_{\rho}}{2},$ containing (un)scattered part 1/2
($\Theta_{\rho}/2$) corresponding to $\Gamma$ ($\Omega$), equation (5)
generalizes the dominant plane wave contribution $e^{ikx}$ from scattering in
1D Hood (2017); Le Kien and Rauschenbeutel (2017); Pichler  _et al._ (2015);
Rao and Hughes (2007), and places a fundamental restriction on the form of
travelling waves at infinity. For comparison, the classical approximation
obtained sup  by a stationary phase argument (see equation (7), Stamnes
(2002); Zakharevich (2017); Camley and Maradudin (1983); Chigrin (2004)), is
valid up to the leading inverse power of $\rho.$ Through Equation (5) power-
law decays of all orders are retained for $\Omega$, which can be computed on
the same footing as $\Gamma$ through a regular surface integral whilst
converging to evanescent order. In Fig. 1(b), for separation vector
$\bm{\rho}=n(\mathbf{a}_{x}+\mathbf{a}_{y})$ along the square diagonal, the
evanescent versus power-law decay of the error from approximants (5) and (7)
can be seen in the inset. Equation $\eqref{eq:medrange}$ offers particular
utility when strongly subradiant Zhang and Mølmer (2019); Van Loo  _et al._
(2013b); Asenjo-Garcia  _et al._ (2017); Kornovan  _et al._ (2019) spin-spin
dynamics are present in a large system, where additional corrections present
in (3) and beyond those obtained by heuristic arguments Hood (2017); Le Kien
and Rauschenbeutel (2017); Pichler  _et al._ (2015); Eldredge  _et al._ (2016)
have the potential to observably modify dynamics Kornovan  _et al._ (2019).

Diverging transverses mass extending long range excitations.—The standard
dichotomy of emitters weakly coupled to baths in the Markovian approximation
is between emission into travelling waves for resonant frequency within the
band Ying  _et al._ (2019); González-Tudela and Cirac (2017); Kim  _et al._
(2021), and evanescent waves forming bound states outside of the band Bykov
(1975); John and Wang (1990); Douglas  _et al._ (2015); González-Tudela  _et
al._ (2015). However, in $dD$ with $d>2,$ both regimes of decay may occur at a
single frequency, offering greater tunability in existing setups. The
availability of travelling waves is quantified through the density of states.
In a similar manner, one may quantify all the all the resonant contributions
to $G$ with group velocity flux through a local solid angle. In sup , we map a
sum over $S$ in $\mathbf{k}$-space to a sum over the unit sphere in
$\hat{\mathbf{v}}$-space, and arrive at a cross section with respect to
direction of propagation:

$\sigma(\hat{\mathbf{v}},\Delta)\sim\frac{1}{v(\mathbf{k})|K(\mathbf{k})|}=|\det[\mathbf{m}^{T}(\mathbf{k})]|v^{d-2}(\mathbf{k})$
(6)

where $K$ is the (intrinsic) Gaussian curvature of $S(z)$, and an implicit sum
is taken over all $\mathbf{k}$ in $S$ with propagation direction
$\hat{\mathbf{v}}.$ The transverse effective mass tensor $\mathbf{m}^{T}$
restricts the usual effective mass tensor
$[\mathbf{m}^{-1}(\mathbf{k})]_{ij}=\frac{\partial^{2}\omega}{\partial
k_{i}\partial k_{j}}$ to the tangent space of $S$ at $\mathbf{k}$, so that
only the transverse response to transverse impulses determines long-range
propagation. Physically, a weakly coupled emitter produces a long-time bath
population distribution localized around $S(\Delta^{\prime})$ with population
distribution $\sigma(\hat{\mathbf{v}},\Delta^{\prime})$, (see sup ) to lowest
order. Here $\Delta^{\prime}$ is the bath-modified emitter frequency.
(Normalizing and) integrating over the cross section additionally gives the
normalized density of states at $\Delta$ (conservation law for resonant
scattering Chigrin (2004); Yang  _et al._ (2019); Zhou  _et al._ (2017)), so
that (6) also admits an interpretation as density of resonant waves with
respect to propagation direction. Through this analysis we obtain a
correspondence between intrinsic curvature on $S$, transverse effective mass
$\mathbf{m}^{T}$ of excitations, and and scattering cross section. Naturally,
the influence of the effective mass is also present in the Green’s function
elements. Applicationsup  of the classical method of stationary phase to (5),
results in the further approximation

${G}\sim
i\theta\mathcal{A}\frac{\psi(\mathbf{k}_{0},\mathbf{r})\psi^{*}(\mathbf{k}_{0},\mathbf{r^{\prime}})}{v(\mathbf{k}_{0})\sqrt{|K(\mathbf{k}_{0})|}}\bigg{(}\frac{1}{2\pi\rho}\bigg{)}^{(d-1)/2}$
(7)

plus terms of order $O(\rho^{-(1+d)/2})$ in the limit, where a implicit
discrete and finite sum is over the wavevectors $\mathbf{k}_{0}$ that satisfy
the resonance condition $\mathbf{k}_{0}\in S,$ with
$\hat{\mathbf{v}}(\mathbf{k}_{0})$ parallel to ${\bm{\rho}}$ and
${\mathbf{v}}(\mathbf{k}_{0})\cdot{\bm{\rho}}>0$. The phase factor $\theta$ is
given in Appendix B.2. Through (6), a divergence of transverse mass (vanishing
of curvature) gives a singularity in (7), invalidating analysis and allowing
one to move past the generic $\rho^{(1-d)/2}$ scaling along the caustic
directionsChigrin (2004); Camley and Maradudin (1983); sup  corresponding to
group velocities of massive excitations. Beyond the caustic, $\sigma$ drops to
zero, as in Figure 2(b). In this region, oscillating and evanescently decaying
waves are the dominant means of propagation, whose properties can be inferred
from local massive excitations when close to the caustic.

Figure 2: Wave propagation in the anisotropic square lattice. (a) Dispersion
relation for $J_{y}=2J_{x}=2J,$ with $\Delta=-J$. The white dot denotes the
massive excitation $m^{T}(\mathbf{k}_{\infty})\to\infty.$ (b) Cross sectional
density of states given by Equation (6). When more than one wavevector
contributes, the sum is taken. (c) Logarithm of $|JG|$ for varying
$\tan\theta=\rho_{y}/\rho_{x}$. Inset shows $J\Omega$ and $J\Gamma$ for
$\theta=0.$ (d) Inverse decay length $\kappa$ for $\Gamma$ beyond the critical
angle. Diamonds give the exact value through calculation of the Green’s
function, whilst the solid line gives the approximate value obtained by
fitting a decay of the form (8).

Ghost waves.—Beyond the caustic direction, the Green’s function elements decay
evanescently with oscillations — these so-called ghost waves Kuś  _et al._
(1993); Narimanov (2019) typically accompany hyperbolic dispersion in
structured media Cortes and Jacob (2017). Adjusting $\bm{\rho}$ beyond the
caustic allows one to modulate the exponential decay between QEs akin to
standard dynamics near a band edgeHood (2017); González-Tudela  _et al._
(2015); Douglas  _et al._ (2015). Specifically, in Figure 2(c), the
exponential decay length of $G$ is tuned between zero and a maximum value
determined by bath microscopics as $\bm{\rho}$ is varied. Whilst dynamics deep
into the ghost wave regime are dependent upon specifics of the bath, a
universal scaling behaviour near the caustic can generally be obtained for
$\Gamma$. When the separation vector is perturbed beyond a global caustic, the
stationary point of $\mathbf{k}\cdot\bm{\rho}$ on $S$ bifurcates into the
complex plane. Using saddle point analysis sup , one may then arrive at a
generic decay in 2D:

$\Gamma\sim\frac{e^{-\kappa\rho}}{\rho^{1/2}},\ \ \ \ \ \ \kappa\propto
a|\theta-\theta_{c}|^{3/2}.$ (8)

as observed in Figure 2(d), which is consistent with the freedom of complex
$\mathbf{k}$ in (7). In Figure 2(c), the decay of $\Omega$ is at least
dominated by $\Gamma,$ so that $G\sim\Gamma$ within the ghost wave regime.
Remarkably, when $\bm{\rho}\propto\mathbf{x}$, $\Omega$ is many orders of
magnitude smaller than $\Gamma,$ (inset) allowing one to construct nearest-
neighbor spin models with purely incoherent interactions. This behaviour is
not limited to the current setup sup . Furthermore, in coupled systems of QEs
and photonic crystals both detuning and angle can be freely varied Douglas
_et al._ (2015); González-Tudela  _et al._ (2015), giving one dynamical
control over the decay length of interactions in the ghost wave regime in the
same manner as the $1/\sqrt{\Delta}$ scaling of decay length beyond a band
edgeDouglas  _et al._ (2015); González-Tudela  _et al._ (2015); Hood (2017);
Newman  _et al._ (2018).

Figure 3: Emission from quantum emitters in the effective magnetic field, for
a bath 111 (1201) sites wide in the $x(y)$ direction. System parameters as in
Figure 2, with $\alpha=0.01.$ In (a-c), $g=0.1J,$ and in (d-e), $g=0.025J.$
(a) Bath population $\langle\hat{a}^{\dagger}\hat{a}\rangle$ at $t=6\tau$ for
an emitter prepared in excited state $|e\rangle$, and located at the central
white diamond (0,0). (b) Bath population at $t=\tau$. Green lines give
semiclassical trajectories initially located at (0,0) with $\mathbf{k}$
sampled from $S$. Arrows correspond to photonic currentDe Bernardis  _et al._
(2021). (c) Bath population along the ${x}$ axis when ${y}=0,l.$ (d) Total
emitter population for one, two, and three emitters, each separated by $3l.$
The two-atom (three atom) subradiant state
$\frac{|eg\rangle-|ge\rangle}{\sqrt{2}}$
$\Big{(}\frac{|egg\rangle-2|geg\rangle+|gge\rangle}{\sqrt{6}}\Big{)}$ is
prepared. (e) Bath population at $t=10\tau$ for two emitters separated a
distance $3l$ and prepared as in (c).

Quasi-1D transport.—We now restrict our attention to 2D, although an immediate
generalization applies in 3D. Away from Van Hove singularities, we have the
invariant for a generic resonant set:

$\int_{S}\frac{d\mathbf{k}}{m^{T}(\mathbf{k})v(\mathbf{k})}=2\pi n.$ (9)

For a single closed orbit, $n=1.$ This quantity does not change with
variations in parameters until a Van Hove singularity with zero group velocity
is encountered and $S$ self-intersects. In the case of Figure 2(a), we find
$n=0$, corresponding to an open orbit. To probe the orbit, we subject the
lattice to an effective magnetic field. With the absence of Berry curvature,
the semi-classical equations of motion for a wavepacket are Gosselin  _et al._
(2006); Fang  _et al._ (2012):

$\displaystyle\dot{\mathbf{r}}=\mathbf{v}(\mathbf{k}),$ (10)
$\displaystyle\dot{\mathbf{k}}=-\mathbf{v}(\mathbf{k})\times{\mathbf{B}}.$
(11)

We take the effective magnetic field
$\mathbf{B}=\frac{\Phi}{a^{2}}\hat{\mathbf{z}}=B\hat{\mathbf{z}}$ pointing out
of the lattice, corresponding to the physical effect of a gauge field
$\mathbf{A}(\mathbf{r})=\frac{\mathbf{r}\times\hat{\mathbf{z}}}{2a^{2}}$
effecting phase accumulation $J_{ij}\to J_{ij}e^{i\phi_{ij}}$, with
$\phi_{ij}=\int_{\mathbf{r}_{i}}^{\mathbf{r}_{j}}\mathbf{A}(\mathbf{r})\cdot
d\mathbf{r}$ and $\Phi=\sum_{\square}\phi_{ij}$. Conservation of
$\omega(\mathbf{k})$ through (11) dictates that the wave packet traces along
$S$ in $\mathbf{k}$ space, so that the geometry of $S$ will be imprinted into
evolution in $\mathbf{r}$ space. In particular, circular Landau orbits are
obtained for closed orbits of isotropic regionsDe Bernardis  _et al._ (2021);
Page (1930) but the zero winding number (9) of open orbits results in an
average transport of wave packets over time. We may characterise the spatial
(temporal) period $\mathbf{l}\ (\tau)$ of orbit trajectories through the
dimensionless strength of the magnetic field, $\alpha=\frac{\Phi}{2\pi}$,
together with fundamental bath quantities:

$\displaystyle\mathbf{l}=\int
d\mathbf{r}=\frac{1}{B}\int_{S}\frac{d\mathbf{k}\mathbf{v}(\mathbf{k})}{v(\mathbf{k})}=\frac{\hat{\mathbf{a}}_{y}}{\alpha},$
(12) $\displaystyle\tau=\int
dt=\frac{1}{B}\int_{S}\frac{d\mathbf{k}}{v(\mathbf{k})}=\frac{\Gamma}{2\alpha}.$
(13)

The transverse extent of the orbit may be calculated similarly. Provided that
lattice constant is fixed, $\mathbf{l}$ is invariant under microscopic
variations of the lattice until a Van Hove singularity is encountered, so that
such quasi-1D orbits are expected to be present generally when band
transitions and Berry curvature can be neglectedGosselin  _et al._ (2006).
Phenomena due to non-zero $l=|\mathbf{l}|$ are manifest when a QE is coupled
to the bath, as emission is dominated by wavepackets localized on $S$ (as in
Equation (5)) that evolve as Equations (10-13). In Figure 3, the population
plots (a) and (b) reveal quasi-1D transport. Notably, periodicity of the
travelling orbits combined with common initial position at the QE demands
periodic refocussing of emission. In Figure 3(c), the population in a cross-
sectional slice of the bath is given at the QE position, and after a single
period. Excepting the QE position where local evanescent fields contribute,
refocusing at subsequent sites occurs almost down to the single site level,
and independently of bath length scale. Exploiting this single site resolution
can in particular be used effect cavities within the lattice. The Markovian
decay of a single emitter is strongly modified when more than one emitter is
present, due to strong coupling at the refocussing points. Separating emitters
by integer multiples of $l$ in the $y$ direction, an effective cavity can be
formed from QEs prepared in a dark stateSinha  _et al._ (2020); Chang  _et
al._ (2012). A snapshot of bath population in 3(e) is shown at
$t=5\Gamma/\alpha$ for two QEs, and in Figure 3(d), the collective decay of
one, two, and three QEs is shown, prepared in dark states in each case.
Significant emitter and localized bath population between emitters is then
sustained orders of magnitude beyond the single-emitter lifetime due to
destructive interference forming bound states in the continuum Chang  _et al._
(2012); Sinha  _et al._ (2020); Hsu  _et al._ (2016). However, the cavities
modes inherently comprise a continuum of $\mathbf{k}$ modes in 2D (lying on
$S$), whilst population nodes do not occur at the edges of the cavity, in
contrast to the recent established correspondence between coupled QEs and site
vacancies Leonforte  _et al._ (2021). These atypical cavities will be
investigated in a separate work.
Topological robustness.—We finally investigate the robustness of the quasi-1D
transport in the system. Elastic scattering of wavepackets preserves energy,
so that waves scattered from local perturbations are also restricted to
quasi-1D. In Figure 4(a) we observe bidirectional scattering, together with
periodic refocussing of the beam after impinging on an obstruction, similar to
the ‘healing’ observed in Bessel and Airy beamsMcGloin and Dholakia (2005);
Broky  _et al._ (2008), albeit effected by a different mechanism. In addition
to local obstructions, quasi-1D propagation is also robust to global lattice
perturbations through invariant (9). Introducing on-site energy disorder $XJ,$
where $X$ is a random variable sampled uniformly in $[-\chi,\chi],$ for the
disorder parameter $\chi$ and averaging over 500 realizations, the average
bath populations are given in Figure 4(b-f). The mean bath population retains
its behaviour in Fig 4(b), with only slight background noise. In (c-f) a
cross-section is taken at ${y}=2l$, with the solid line showing the average
and the shadow denoting a standard deviation width. The periodic refocussing
is seen to be qualitatively preserved, with bath population maintaining its
peak by many orders of magnitude at the central site. The fluctuations in the
weakly populated sites are additionally many orders of magnitude smaller than
the central site population, which is well preserved with relatively small
fluctuations up to moderate disorder $\chi=0.5.$ In Fig 4(d), the subsequent
evanescent decay rate is additionally shared across all disorders. Even up to
moderate frequency disorders probing the vicinity Van Hove singularities,
invariant (9) protects quasi-1D emission in an manifestation of the Fermi
surface topology Yuan  _et al._ (2019); Zhao and Wang (2013); He  _et al._
(2019) in condensed matter physics.

Figure 4: Effects of disorder. Parameters are as in Fig. 3 with $g=0.025J$.
(a) Bath population at $t=3\tau$ in the presence of an obstruction (white
line) at $y=l/2,x\in[-10a,10a]$ obtained by decoupling all obstruction sites
from the lattice. (b) Averaged bath population at time $t=3\tau$ for disorder
$\chi=0.5.$ (c-f) Average log-populations for the cross-section $y=2l$ in the
presence of varying disorder parameter $\chi.$ The solid lines give the
average whilst the shadow denotes a standard deviation. In (f) the population
in the presence of no disorder corresponding to Fig. 3(a) is shown for
comparison.

Phenomena induced by Fermi surface topology concern the behaviour of
$\mathbf{k}$ on $S$ rather than wavefunctions on the entire band, and Van Hove
singularities in contrast to degeneracies. Whilst studies are generally
limited to 2D, it would be interesting to explore the topological effects
afforded by energy landscape in higher dimensions, possibly in conjunction
with standard wavefunction induced topological effects.

Conclusion.—To summarize, we have studied the propagation of wave dynamics in
structured media, exploring the consequences of (quasi-)breaking of isotropy
and periodicity. We have proved a decomposition of linear wave propagation in
arbitrary periodic local lattices as a sum of a travelling part and evanescent
part. Long-range wave transport can be dominantly characterised by effective
mass of resonant excitations through the geometry of the level set. Inducing
open orbits in the dispersion relation can effect ghost waves as the dominant
transport mechanism, for which we found a universal long range scaling near
caustics. Finally, breaking periodicity with an effective magnetic field
resulted in topologically protected quasi-1D transport in direct contrast to
the usual Landau orbits, offering robust bidirectional transport and atypical
bound states in the continuum. This work offers general insights into the
regimes of wave propagation in generic systems, with potential applications in
emerging metamaterials and quantum optics. Of particular interest would more
detailed studies including transport around Van Hove singularities and
degeneracies, whilst studying interacting non-linear many-body dynamics may
prove particularly fruitful.

###### Acknowledgements.

This work was supported by the Okinawa Institute of Science and Technology.
Simulations in this work were performed using the open source Julia packages
DifferentialEquations.jlRackauckas and Nie (2017), QuadGK.jlJohnson , and
MDBM.jlBachrathy and Stépán (2012). The authors acknowledge discussions with
F. le Kien. We are additionally grateful for the scientific cluster and
resources provided by the Scientific Computing and Data Analysis section of
the Research Support Division at OIST.

##  Appendix A Derivation of the long-range behaviour of the Green’s function

###  A.1 Proving that degeneracies contribute only evanescent decay

When site couplings are local (i.e. a site only couples to finitely many
neighbours), the dispersion relation is analytic on the real axis Des
Cloizeaux (1964) away from degeneracies. Away from Van Hove singularities
$\nabla\omega_{\nu}(\mathbf{k})\neq 0$, the implicit function theorem applies,
and the resonant set $S_{\nu}$ is a smooth manifold of co-dimension $(d-1),$
allowing one to apply the full power of differential geometry. However, at a
degeneracy, it seems that the integrand of

$I=\sum_{\nu}\int_{BZ}\frac{d^{d}\mathbf{k}}{(2\pi)^{d}}\frac{\psi_{\nu}(\mathbf{k},\mathbf{r})\psi^{*}_{\nu}(\mathbf{k},\mathbf{r}^{\prime})}{\Delta-\omega_{\nu}(\mathbf{k})}$
(14)

becomes non-smooth Kato (2013), contributing a power-law decay. Assuming the
degeneracy does not coincide with the resonant set $S$, we show that summing
over all bands removes any non-analyticity in individual bands. This allows us
to consider only a single band in the main text when the resonant set is away
from degeneracies. The derivation follows in a manner similar to Sánchez-
Burillo  _et al._ (2020). Recall the unitary transformation between
eigenfunctions $\psi_{i}$ and the Bloch wave corresponding to an individual
sublattice:

$|\psi_{\nu}(\mathbf{k})\rangle=\sum_{j}U_{\nu
j}(\mathbf{k})|\mathbf{k}_{j}\rangle,$ (15)

where
$|\mathbf{k}\rangle_{i}=\frac{1}{\sqrt{N}}\sum_{r_{i}}e^{i\mathbf{k}\cdot\mathbf{r}_{i}}|\mathbf{r}_{i}\rangle$
with $|\mathbf{r}_{i}\rangle$ the localized excitation at any lattice site
$\mathbf{r}_{i}$ in lattice $i$. In the thermodynamic limit we conversely have
$|\mathbf{r}_{i}\rangle=\int_{BZ}\frac{d^{d}\mathbf{k}}{(2\pi)^{d}}e^{-i\mathbf{k}\cdot\mathbf{r}_{i}}|\mathbf{k}_{i}\rangle,$
and may rewrite the inner product wavefunctions:

$\psi_{\nu}(\mathbf{k},\mathbf{r}_{i})=\langle\mathbf{r}_{i}|\psi_{\nu}(\mathbf{k})\rangle=\sum_{j}\int_{BZ}\frac{d^{d}\mathbf{k}^{\prime}}{(2\pi)^{d}}e^{i\mathbf{k}\cdot\mathbf{r}_{i}}U_{\nu
j}(\mathbf{k})\langle\mathbf{k}^{\prime}_{i}|\mathbf{k}_{j}\rangle=e^{i\mathbf{k}\cdot\mathbf{r}_{i}}U_{\nu
i}(\mathbf{k}).$ (16)

Assuming $\mathbf{r}$ and $\mathbf{r}^{\prime}$ belong to lattices $i$ and
$i^{\prime}$ respectively, we rewrite the integrand

$I=\sum_{\nu}\int_{BZ}\frac{d^{d}\mathbf{k}}{(2\pi)^{d}}\frac{e^{i\mathbf{k}\cdot\bm{\rho}}U_{i^{\prime}\nu}^{*}(\mathbf{k})U_{\nu
i}(\mathbf{k})}{\Delta-\omega_{\nu}(\mathbf{k})},$ (17)

or more succinctly in matrix notation

$I=\int_{BZ}\frac{d^{d}\mathbf{k}}{(2\pi)^{d}}e^{i\mathbf{k}\cdot\bm{\rho}}\big{[}\mathbf{U}^{\dagger}(\mathbf{k})(\Delta-\mathbf{D}(\mathbf{k}))^{-1}\mathbf{U}(\mathbf{k})\big{]}_{i^{\prime}i}.$
(18)

We note the diagonalization
$\mathbf{h}(\mathbf{k})=\mathbf{U}^{\dagger}(\mathbf{k})\mathbf{D}(\mathbf{k})\mathbf{U}(\mathbf{k})$
and recognize the inner element as an inverse of the finite dimensional
resolvant

$I=\int_{BZ}\frac{d^{d}\mathbf{k}}{(2\pi)^{d}}e^{i\mathbf{k}\cdot\bm{\rho}}(\Delta-\mathbf{h}(\mathbf{k}))^{-1}_{i^{\prime}i}=\int_{BZ}\frac{d^{d}\mathbf{k}}{(2\pi)^{d}},$
(19)

where $\mathbf{h}(\mathbf{k})$ is now the finite dimensional matrix with
analytic elements describing coupling between each of the Bloch sublattices.
By standard resolvant analysis Kato (2013), the resolvant is analytic in both,
$\mathbf{k}$ and $\Delta,$ whenever $\Delta$ is not equal to any of the
eigenvalues of $\mathbf{h}(\mathbf{k}),$ i.e. not equal to any of the energy
bands at $\mathbf{k}.$ Thus, when the resonant level set energy $\Delta$ does
not lie at a degeneracy, the integrand is analytic with respect to integration
variable $\mathbf{k},$ and localization of the full Green’s function integral
around a degeneracy but away from the resonant level produces an analytic
integrand with can be treated with the asymptotic analysis in the following
sections.

###  A.2  Long range behaviour of the coherent part of the Green function in
$d$ dimensions

Consider the Green’s function describing propagation through a single mode
$\nu$ (whose subscript we omit in the following) of a periodic continuum bath:

${G}(\bm{\rho},\Delta)=\mathcal{A}\int_{\text{BZ}}\frac{d^{d}\mathbf{k}}{(2\pi)^{d}}\frac{\psi(\mathbf{k},\mathbf{r})\psi^{*}(\mathbf{k},\mathbf{r^{\prime}})}{\Delta-\omega(\mathbf{k})+i0},$
(20)

The integral is understood through the limiting absorption principle, taking a
small imaginary part in the denominator to zero. Whilst the wave function
$\psi$ is in this case a scalar, we note that similar expressions arise for
$\psi$ a tensor, i.e. when polarization is available in an electromagnetic
system. With time reversal symmetry of the bath assumed, we extract the
strongly singular continuum integral corresponding to coherent wave
propagation:

$\displaystyle{\Omega}=\mathcal{A}\mathcal{P}\int_{\text{BZ}}\frac{d^{d}\mathbf{k}}{(2\pi)^{d}}\frac{\psi(\mathbf{k},\mathbf{r})\psi^{*}(\mathbf{k},\mathbf{r^{\prime}})}{\Delta-\omega(\mathbf{k})},$
(21)

and proceed with an asymptotic large
$\bm{\rho}=\mathbf{r}-\mathbf{r}^{\prime}$ treatment, where the direction of
$\bm{\rho}$ is kept fixed. Having shown analyticity of the Kernel around
degeneracies in Appendix A.1, we now address the (near) resonant contributions
from $S$. Under the assumption of locality in Appendix A.1, we may take $S$ to
be a smooth manifold of codimension one, and note that $\omega(\mathbf{k})$ is
analytic for real $\mathbf{k}$ away from degeneracies. With this we may use a
partition of unity (or localization argument). That is, we may choose smooth
bump functions Johnson (2015) with support on arbitrarily small regions near
$S$ and away from $S$ respectively. Away from $S$ the kernel is smooth, and so
evanescent decay of the integrand is obtained.

Let us first factorize the Bloch waves:

${\Omega}=\mathcal{A}\mathcal{P}\int_{\text{BZ}}\frac{d^{d}\mathbf{k}}{(2\pi)^{d}}\frac{\Phi(\mathbf{k})e^{i\bm{\rho}\cdot\mathbf{k}}}{\Delta-\omega(\mathbf{k})},$
(22)

where we increase the separation with $\mathbf{r},\mathbf{r}^{\prime}$ fixed
modulo the lattice, so that
$\Phi(\mathbf{k})=\psi(\mathbf{k},\mathbf{r})\psi^{*}(\mathbf{k},\mathbf{r^{\prime}})e^{-i\bm{\rho}\cdot\mathbf{k}},$
$\Phi$ is unchanging with $\bm{\rho}$ and can be considered purely as a
function of $\mathbf{k}.$ We may parametrise the vicinity of the resonant set
$S=S(\Delta)=\\{\mathbf{k}:\omega(\mathbf{k})=\Delta\\}$ using canonical
(Fermi) coordinates in a tubular neighbourhood of $S$:

$\mathbf{k}=(t,\bm{\xi})=\bm{\xi}+t\hat{\mathbf{n}}({\bm{\xi}}),$ (23)

where $\bm{\xi}$ parameterizes the level set $S$, and $\hat{\mathbf{n}}$ is
the unit normal, which we take facing outwards. As the level set is a smooth
surface of codimension $1$, it is orientable and such a normal always exists.
Note that as group velocity is assumed non-vanishing on the surface, the
velocity vector always points inside or outside. We assume without loss of
generality that the group velocity points outside, although a similar result
is obtained otherwise. The above parametrisation is smooth when $t$ is smaller
than any of the radii of curvature anywhere on $S$ Alvarez (2017), and we
choose $\epsilon$ such that $|t|<\epsilon$ satisfies this condition.

We now use the partition of unity to isolate the contribution from $S$. For
this we consider the integral:

$R=\mathcal{A}\mathcal{P}\int_{\delta
S}\frac{d^{d}\mathbf{k}}{(2\pi)^{d}}\phi_{\delta
S}(\mathbf{k})\frac{\Phi(\mathbf{k})e^{i\bm{\rho}\cdot\mathbf{k}}}{\Delta-\omega(\mathbf{k})},$
(24)

where $\delta S$ is a thin tubular region surrounding $S$, with
$|t|<\epsilon$. This is smoothly facilitated by the bump function
$\phi_{\delta S}\in C^{\infty}$ forming one of our partition of unity. The
bump function can be given entirely through dependence on $t$:

$\phi_{\delta S}(\mathbf{k})=\phi_{\epsilon}(t),$ (25)

where $\phi_{\epsilon}(t)$ is the usual one-dimensional bump function,
normalized to unity at zero:

$\phi_{\epsilon}(x)=\begin{cases}\exp(\frac{x^{2}}{x^{2}-\epsilon^{2}})&|x|<\epsilon\\\
0&\text{Otherwise}\end{cases}$ (26)

We now change coordinates systems, giving the volume element obtained in the
derivation ofWeyl’s tube formula Alvarez (2017); Gray (2012); Khudaverdian
_et al._ (2007); Choe and Koiso (2015):

$d^{d}\mathbf{k}=\det(\mathbf{I}+t\mathbf{K})d\bm{\xi}dt,$ (27)

where $\mathbf{K}$ is the second fundamental form of the surface $S$. In the
case of an implicit surface we have Verpoort (2008)

$\mathbf{K}=\mathbf{K}(\bm{\xi})=-\frac{\mathbf{H}_{T}(\bm{\xi})}{v(\bm{\xi})},$
(28)

for the Hessian

$[\mathbf{H}(\mathbf{k})]_{ij}=\frac{\partial^{2}\omega(\mathbf{k})}{\partial
k_{i}\partial k_{j}},$ (29)

which we restrict to a linear map $\mathbf{H}_{T}(\bm{\xi})$ acting on the
tangent space of $S$ at $\bm{\xi}.$ We obtain

$R=\mathcal{A}\mathcal{P}\int_{S}\frac{e^{i\bm{\xi}\cdot\mathbf{R}}d^{d-1}\bm{\xi}}{(2\pi)^{d-1}}\int_{-\epsilon}^{\epsilon}\frac{dt}{2\pi}\phi_{\epsilon}(t)\frac{\Phi(t,\bm{\xi})e^{it\bm{\rho}\cdot\hat{\mathbf{n}}}}{\Delta-\omega(t,\bm{\xi})}\det(\mathbf{I}+t\mathbf{K}),$
(30)

We note that the only non-smoothness is at $t=0$ and rewrite

$R=\mathcal{A}\int_{S}\frac{e^{i\bm{\xi}\cdot\mathbf{R}}d^{d-1}\bm{\xi}}{(2\pi)^{d-1}}\mathcal{P}\int_{-\epsilon}^{\epsilon}\frac{dt}{2\pi}\frac{L(t,\bm{\xi})e^{it\bm{\rho}\cdot\hat{\mathbf{n}}}\phi_{\epsilon}(t)}{t}\det(\mathbf{I}+t\mathbf{K}).$
(31)

Here, $L(t,\bm{\xi})$ defines a smooth function within $\delta S$ (local
expansion in small $t$ reveals that smoothness can be extended to $t=0$).
Additionally, the determinant may be expanded, upon which we need only the
lowest order term (the higher order terms cancel $t$s and produce evanescent
decay). We can separate out only the contribution at $t=0$, and discard small
contributions to obtain

$\Omega\sim
R\sim\mathcal{A}\int_{S}\frac{e^{i\bm{\xi}\cdot\mathbf{R}}d^{d-1}\bm{\xi}}{(2\pi)^{d-1}}\mathcal{P}\int_{-\epsilon}^{\epsilon}\frac{dt}{2\pi}\frac{L(0,\bm{\xi})e^{it\bm{\rho}\cdot\hat{\mathbf{n}}}\phi_{\epsilon}(t)}{t}+O(\rho^{-\infty}),$
(32)

as $\rho\to\infty$. We use $O(\rho^{-\infty})$ to denote decay faster than any
polynomial which in the main text we term ‘evanescent decay’. The remainder we
subtracted formed a smooth function in $\delta S$, and for $|t|<\epsilon/2$
the surface element is also smooth, so that the remainder will decay
evanescently by the principal of non-stationary phaseZakharevich (2017).

Thus, replacing $L(0,\bm{\xi})$ we may write (ignoring evanescent corrections)

$\Omega=\frac{i\mathcal{A}}{2}\int_{S}\Phi_{\rho}(\hat{\bm{\rho}}\cdot\hat{\mathbf{n}}(\mathbf{k}))\frac{\psi(\mathbf{k},\mathbf{r})\psi^{*}(\mathbf{k},\mathbf{r^{\prime}})d^{d-1}\mathbf{k}}{v(\mathbf{k})(2\pi)^{d-1}}+O(\rho^{-\infty}),$
(33)

where we have relabelled $S\ni\bm{\xi}$ with $\mathbf{k}\in S$ and where we
have the phase function

$\Theta_{\rho}(x)=\int_{-\infty}^{\infty}\frac{dt\sin(tx)\phi_{\rho\epsilon}(t)}{\pi
t}.$ (34)

Figure 5: The Fourier transform $\Theta_{\rho}(x)$ for different values of
$\rho$. Due to the identical scaling of $\Theta_{\rho}$ with $\epsilon$ (in
Equation (34)), we fix $\epsilon=1$ here. $\Theta_{\rho}$ smoothly
interpolates $-1$ and $1$ over the real line, with $\Theta_{\rho}(0)=0$.
$\Theta_{\rho}(x)$ pointwise converges to $\text{Sign}(x)$ with evanescent
corrections as $\rho\epsilon\to\infty.$ We see that the sign changes according
to the sign of $\hat{\mathbf{v}}$ projected onto $\hat{\bm{\rho}}$, or the
angle between the two. Inset shows the variation of $\Theta_{\rho}$ as one
traverses an isotropic level set $\mathbf{k}(t)=(\cos(2\pi t),\sin(2\pi t)).$

$\Theta_{\rho}$ acts as a sigmoid function, varying smoothly from -1 at
$-\infty$ to 1 at $\infty$ (see Figure 5), and as the integral of the Fourier
transform of the bump function Johnson (2015), has the asymptotic behaviour

$\Theta_{\rho}(x)=\text{Sign}(x)+O\big{(}\text{Erf}(|\rho\epsilon
x|^{3/4})\big{)}=\text{Sign}(x)+O\big{(}|\rho\epsilon
x|^{-1/4}\exp(-|\rho\epsilon x|^{1/2})\big{)}$ (35)

with $\epsilon>0$ and $x\to\pm\infty.$ We thus have the full Green’s function
approximant of the main text, converging to $G$ with evanescent errors as
$\bm{\rho}\to\infty:$

$G\sim\frac{i\mathcal{A}}{2}\int_{S}(1+\Theta_{\rho}(\hat{\bm{\rho}}\cdot\hat{\mathbf{n}}))\frac{\psi(\mathbf{k},\mathbf{r})\psi^{*}(\mathbf{k},\mathbf{r^{\prime}})d^{d-1}\mathbf{k}}{v(\mathbf{k})(2\pi)^{d-1}},$
(36)

where the amplitude of the main text is given by

$\Pi_{\rho}(\hat{\bm{\rho}}\cdot\hat{\mathbf{n}})=\frac{(1+\Theta_{\rho}(\hat{\bm{\rho}}\cdot\hat{\mathbf{n}}))}{2}$
(37)

We now briefly discuss when the above approximation can be expected to be
accurate. Most notably, the tubular radius $\epsilon$ should be large enough
to contain a significant fraction of the band, whilst the dispersion relation
and Bloch wavefunctions should not vary significantly throughout the band. If
this is not the case, our approximation will perform poorly at short and
intermediate ranges due to insufficient capturing of the entire band dynamics
away from $S$. Physically, this will occur when strong resonances exist away
from $S$, when $S$ has regions of large curvature, or when $S$ lies close to a
Van Hove singularity or degeneracy. In any of these cases $\Omega$ can be
expected to dominate at shorter distances García-Elcano  _et al._ (2020);
González-Tudela and Cirac (2018); Perczel and Lukin (2020); González-Tudela
and Cirac (2017).

###  A.3  Calculation of the exact Green’s function

The strongly singular manifold present in (1) of the main text presents
significant calculational difficultly, so that to obtain exact results either
recurrence relations Katsura  _et al._ (1971) or a small imaginary part
$0i\to\epsilon i$ are typically used. However, the former suffers from
instabilities for large $\bm{\rho},$ whilst the latter suffers an increase in
computational cost due to the divergence of the integrand near $S$. To
navigate around these difficulties we apply a higher-dimensional variant of
integration by parts that gives us a regular integrand for $\Omega$ (with only
divergent derivative). The integrand then permits standard quadrature
techniques, whilst $\Gamma$ is obtained through the regular surface integral
(2). In one dimension repeated integration by parts on periodic functions
$f,\omega$ gives

$\mathcal{P}\int\frac{dkf(k)}{\Delta-\omega(k)}=-\mathcal{P}\int
dk\log|\Delta-\omega(k)|\bigg{(}\frac{f(k)}{v(k)}\bigg{)}^{\prime}=\mathcal{P}\int
dk(\Delta-\omega(k))\big{[}\log|\Delta-\omega(k)|-1\big{]}\Bigg{(}\frac{1}{v(k)}\bigg{(}\frac{f^{\prime}(k)}{v(k)}\bigg{)}^{\prime}\Bigg{)}^{\prime},$
(38)

when the group velocity $v$ is non-vanishing on the support of $f.$ We note
that endpoint contributions cancel via periodicity, and the log divergences at
either side of singularities can also be shown to cancel. If $v$ vanishes on
the support of $f,$ we can decompose
$f(k)=f(k)\phi_{\epsilon}(\Delta-\omega(k))+f(k)\big{(}1-\phi_{\epsilon}(\Delta-\omega(k))\big{)}$
for the bump function $\phi_{\epsilon}$, such that $\epsilon$ is small enough
that the support of the first term does not contain any points of vanishing
group velocity (i.e. as satisfied by the tubular region). The first term is
then amenable to integration by parts via (38), and the second term is a
regular integrand. In higher dimensions we accordingly make use of the
divergence theorem $\int_{S}\nabla\cdot\mathbf{F}d\mathbf{k}=\int_{\partial
S}\mathbf{F}\cdot\hat{\mathbf{n}}dS,$ noting again that periodic boundary
conditions and divergences of opposite signs cancel on boundaries to see

$\mathcal{P}\int\frac{d\mathbf{k}f(\mathbf{k})}{\Delta-\omega(\mathbf{k})}=-\mathcal{P}\int
d\mathbf{k}\log|\Delta-\omega(\mathbf{k})|\nabla\cdot\bigg{(}\frac{\mathbf{v}(\mathbf{k})f(\mathbf{k})}{v^{2}(\mathbf{k})}\bigg{)},$
(39)

and applying integration by parts again for the final result. The bump
function may be used here to regularize if necessary.

##  Appendix B Characteristics of anisotropic emission

###  B.1 Calculating the directional density of states

We here characterise emission properties in real space through the
distribution of waves on $S$ with respect to the propagation direction. When
$S$ is non-convex (when $S$ has a point of zero curvature) there are multiple
wavevectors contributing to a given direction. In the following we assume that
$S$ is convex, whereby the result is generally obtained by partitioning into
regions where the sign of the curvature is unchanging. For bath quantities
considered as a distribution in $\mathbf{k}$-space, we may instead consider
them as a distribution over the propagation direction
$\hat{\mathbf{v}}(\mathbf{k})=\mathbf{v}(\mathbf{k})/v(\mathbf{k})$, which
maps onto the $(d-1)$-sphere. The mapping $\mathbf{k}\to\hat{\mathbf{v}}$ is
formally known as the Gauss map, whose Jacobian is given as the Gaussian
curvature. The change of coordinates is then

$d^{d-1}\mathbf{k}=\frac{d\Omega}{|K(\mathbf{k})|},$ (40)

for the intrinsic Gaussian curvature of $S$ at $\mathbf{k}.$ We may then
decompose the formula for density of states:

$D(\Delta)=\int_{S^{d-1}}\frac{d\Omega}{(2\pi)^{d}v(\mathbf{k})|K(\mathbf{k})|},$
(41)

so that the number of states with propagation direction through an
infinitesimal cross section $d\Omega$ is obtained by introducing the unit cell
area (the argument $\Delta$ may be dropped when the resonant set $S(\Delta)$
is understood):

$D(\hat{\mathbf{v}},\Delta)=\frac{\mathcal{A}}{(2\pi)^{d}v(\mathbf{k})|K(\mathbf{k})|}.$
(42)

This can further be treated via the introduction of
$\mathbf{m}^{-1}(\mathbf{k})$, the inverse of the effective mass tensor of
excitations in the periodic media, given by

$[\mathbf{m}^{-1}]_{ij}=\frac{\partial^{2}\omega(\mathbf{k})}{\partial
k_{i}\partial k_{j}}.$ (43)

In particular, examination of a perturbation
$\hat{\mathbf{n}}(\mathbf{k})\to\hat{\mathbf{n}}(\mathbf{k}+d\mathbf{k})$
reveals that impulses directed in (perpendicular to) $S$ result in changes in
propagation direction remaining in (perpendicular to) $S$. We have the
explicit form of the Jacobian

$\mathbf{J}_{\hat{\mathbf{n}}}(\mathbf{k})=\frac{\mathbf{H}(\mathbf{k})[\mathbf{I}-\hat{\mathbf{n}}\hat{\mathbf{n}}]}{v(\mathbf{k})},$
(44)

which equals

$\mathbf{J}_{\hat{\mathbf{n}}}(\mathbf{k})=\frac{\mathbf{H}(\mathbf{k})}{v(\mathbf{k})}=\frac{1}{\mathbf{m}^{T}(\mathbf{k})v(\mathbf{k})},$
(45)

when restricted to $S$. The effective mass tensor is now restricted
transversely to $\mathbf{m}^{T}(\mathbf{k}),$ a function acting in the tangent
space of $\mathbf{k}$ where the inverse is generally well defined.

The Gaussian curvature for the implicitly defined resonant set $S$ is then
given through

$\displaystyle
K(\mathbf{k})=-\det[\mathbf{J}_{\hat{\mathbf{n}}}(\mathbf{k})],$ (46)

leading to the state density with respect to the cross section:

$D(\hat{\mathbf{v}},\Delta)=\mathcal{A}\bigg{|}\frac{\det\big{[}v(\mathbf{k})\mathbf{m}^{T}(\mathbf{k})\big{]}}{(2\pi)^{d}v(\mathbf{k})}\bigg{|},$
(47)

which in the following analysis we find proportional to the scattering cross
section in $\hat{\mathbf{v}}$ space.

We begin with the long time population of bath modes as $t\to\infty$González-
Tudela and Cirac (2017); Yang  _et al._ (2019):

$C_{\mathbf{k}}\to-\frac{ge^{-i\omega(\mathbf{k})t}}{\omega(\mathbf{k})-\Delta^{\prime}+i\gamma/2}\sim-\frac{ge^{-i\omega(\mathbf{k})t}}{\omega(\mathbf{k})-\Delta^{\prime}+i\gamma/2},$
(48)

where the frequency shift has been absorbed into $\Delta\to\Delta^{\prime}$,
and $\gamma=g^{2}\Gamma$. The bath population is maximized when
$\omega(\mathbf{k})=\Delta^{\prime}$. In particular, the $\mathbf{k}$-space
cross section is obtained by Yang  _et al._ (2019); Zhou  _et al._ (2017)

$\sigma(\mathbf{k},\Delta)=\frac{\gamma|C_{\mathbf{k}}|^{2}}{v(\mathbf{k})}=\frac{1}{v(\mathbf{k})}\frac{g^{2}\gamma}{(\omega(\mathbf{k})-\Delta^{\prime})^{2}+\frac{\gamma^{2}}{4}}$
(49)

The long-range interactions are composed of the waves resonant with the
emitter for $|g|\ll|J|$, and so the emission pattern in the far-field is
dominated by the angular distribution induced on $S(\Delta^{\prime})\sim
S(\Delta)$, obtained through the coordinate transform:

$\sigma(\hat{\mathbf{v}},\Delta)=\frac{4g^{2}}{v(\mathbf{k})|K(\mathbf{k})|\gamma},$
(50)

whose scaling with $v,K$ is shared by the cross sectional density of states
(47). The bath population $|C_{\mathbf{k}}|^{2}$ scales as
$1/|K(\mathbf{k})|$.

###  B.2 Stationary phase approximation on the resonant level set

To evaluate the asymptotics

$G\sim\frac{i\mathcal{A}}{2}\int_{S}(1+\Theta_{\rho}(\hat{\bm{\rho}}\cdot\hat{\mathbf{n}}))\frac{\psi(\mathbf{k},\mathbf{r})\psi^{*}(\mathbf{k},\mathbf{r^{\prime}})d^{d-1}\mathbf{k}}{v(\mathbf{k})(2\pi)^{d-1}},$
(51)

in the long-range geometric optics limit, we make use of the stationary phase
argument. Whilst the phase function $\Theta_{\rho}$ is also rapidly
oscillating with $\rho,$ integration by parts shows that for our choice of
$\epsilon$ small enough, the oscillations due to the exponent are dominant,
and the leading contribution to the integral come from points where the
exponential argument is at an extremum on $S$. The integral rapidly converges
to a typical exponential integral, and we use the general stationary phase
formula on a manifoldZakharevich (2017):

$\int_{S}d^{d-1}\mathbf{k}f(\mathbf{k})e^{i\lambda\theta(\mathbf{k})}=\Big{(}\frac{2\pi}{\lambda}\Big{)}^{(d-1)/2}\sum_{\mathbf{k}_{0}}f(\mathbf{k}_{0})\frac{e^{\text{Sgn}(H(\mathbf{k}_{0}))\frac{\pi
i}{4}}e^{i\lambda\theta(\mathbf{k}_{0})}}{\sqrt{|\det[H(\mathbf{k}_{0})]|}}+O(\lambda^{(d-3)/2}),$
(52)

as $\lambda\to\infty$, where the sum is over (assumed) discrete critical
points $\mathbf{k}_{0}$ of $\theta$ on $S.$ The Hessian $H$ is also computed
as a Hessian on the manifold $S.$ When applied to (36), we have simply the
phase function $\theta(\mathbf{k})=\mathbf{k}\cdot\hat{\bm{\rho}},$ the height
function with respect to $\hat{\bm{\rho}}.$ According to Izumiya  _et al._
(2016), the Hessian of the height function restricted to $S$ and at a critical
point is simply given as the second fundamental form,
$H(\mathbf{k}_{0})=-\mathbf{K}(\mathbf{k}_{0})$, so that the denominator is
simply the square root of $|K(\mathbf{k}_{0})|$. Plugging in (36) gives the
final long-range approximation:

$G_{R}=i\mathcal{A}\frac{\psi(\mathbf{k}_{0},\mathbf{r})\psi^{*}(\mathbf{k}_{0},\mathbf{r^{\prime}})e^{\frac{i\text{Sgn}[K(\mathbf{k}_{0})]}{4}}}{v(\mathbf{k}_{0})\sqrt{|K(\mathbf{k}_{0})|}}\bigg{(}\frac{1}{2\pi\rho}\bigg{)}^{(d-1)/2}+O(\rho)^{(d-3)/2},$
(53)

We can arrive at the same conclusion by locally parametrising the surface and
using the classical stationary phase formula in Euclidean space:

$\int_{\mathbf{R}^{n}}d\sigma
p(\mathbf{t})e^{iq(\mathbf{t})X}=\sum_{\mathbf{t}_{0}}\frac{p(\mathbf{t}_{0})e^{i(q(\mathbf{t}_{0})X+\text{Sgn}(\mathbf{H}_{q}(\mathbf{t}_{0}))\pi/4)}}{|\text{Det}(\mathbf{H}_{q}(\mathbf{t}_{0}))|^{1/2}}\bigg{(}\frac{2\pi}{X}\bigg{)}^{n/2}+O(X^{-n/2}).$
(54)

###  B.3 Asymptotics beyond the critical angle

We approximate the exponential decay rate just of $\Gamma$ close to but beyond
the critical angle using asymptotic analysis. Near a Van Hove singularity or
when the dispersion varies rapidly, significant contributions to $\Omega$ may
come from regions away from the resonant level set. However, as in the main
text, it is expected that at infinity $\Omega$ and $\Gamma$ become equal up to
phase. We perform calculations in two dimensions for simplicity, but note that
concepts may be extended to higher dimensions. We begin with

${\Gamma}=\mathcal{A}\int_{S}\frac{d\mathbf{k}}{(2\pi)}\frac{\psi(\mathbf{k},\mathbf{r})\psi^{*}(\mathbf{k},\mathbf{r^{\prime}})}{{v}(\mathbf{k})}=\mathcal{A}\int_{S}\frac{d\mathbf{k}}{(2\pi)}\frac{\Phi(\mathbf{k})e^{i\mathbf{k}\cdot\bm{\rho}}}{{v}(\mathbf{k})}.$
(55)

We suppose that $\bm{\rho}_{\infty}$ lies exactly at a simple caustic
corresponding to the wavevector $\mathbf{k}_{\infty}.$ By simple we mean that
$K(\mathbf{k}_{\infty})\neq 0.$ If we then consider a unit length direction
$\hat{\bm{\rho}}=(\hat{\bm{\rho}}_{\infty}+\hat{\bm{\epsilon}})$ just beyond
the caustic with $\epsilon=|\hat{\bm{\epsilon}}|\ll 1$, we have

$\displaystyle\hat{\bm{\epsilon}}\cdot\mathbf{v}(\mathbf{k}_{\infty})=0,$ (56)
$\displaystyle{m}^{T}(\mathbf{k}_{\infty})\epsilon=0,$ (57)
$\displaystyle\bm{\rho}_{\infty}\parallel\mathbf{v}(\mathbf{k}_{\infty}).$
(58)

That is, $\bm{\rho}_{\infty}$ is parallel to the caustic, and
$\hat{\bm{\epsilon}}$ travels orthogonally beyond the caustic. As in the main
text, the caustic is assumed global such that there is no travelling wave with
real vector $\mathbf{k}$ such that
$\hat{\bm{\rho}}\parallel\mathbf{v}(\mathbf{k}).$ Take a periodic
parametrisation $\mathbf{k}(t)$ of $S$ with
$\mathbf{k}_{\infty}=\mathbf{k}(t_{\infty})$, and consider the extension for
$t$ in the complex plane, with small imaginary part. We seek the stationary
point

$\displaystyle(\hat{\bm{\rho}}+\hat{\bm{\epsilon}})\cdot\mathbf{k}^{\prime}(t)=0.$
(59)

For small $\epsilon,$ we expect the stationary point to comprise a
correspondingly small correcting complex part $t^{\prime}$ in $t$ and thus
$\mathbf{k}.$ We expand

$0=\hat{\bm{\rho}}\cdot\mathbf{k}^{\prime}(t)=\hat{\bm{\rho}}\cdot\mathbf{k}^{\prime}(t_{\infty}+t^{\prime})=(\hat{\bm{\rho}}_{\infty}+\hat{\bm{\epsilon}})\cdot\left[\mathbf{k}^{\prime}(t_{\infty})+\mathbf{k}^{\prime\prime}(t_{\infty})t^{\prime}+\frac{t^{\prime
2}}{2}\mathbf{k}^{\prime\prime\prime}(t_{\infty})\right].$ (60)

To lowest order in ${\epsilon}$ and $t^{\prime},$ we obtain

$0=\hat{\bm{\epsilon}}\cdot\mathbf{k}^{\prime}(t_{\infty})+\frac{t^{\prime
2}}{2}\hat{\bm{\rho}}_{\infty}\cdot\mathbf{k}^{\prime\prime\prime}(t_{\infty}).$
(61)

Under assumption of simple caustic, the third derivative is non-zero and the
above equation is balanced for purely imaginary correction $t^{\prime}$:
$i\sqrt{\epsilon}\sim t^{\prime}$. The normalized argument of the exponent at
the new stationary point has a real and imaginary part to leading order

$\displaystyle\mathfrak{R}\big{[}i\hat{\bm{\rho}}\cdot\mathbf{k}(t_{\infty}+t^{\prime})\big{]}\sim-\epsilon^{3/2},$
(62)
$\displaystyle\mathfrak{I}\big{[}i\hat{\bm{\rho}}\cdot\mathbf{k}(t_{\infty}+t^{\prime})\big{]}\sim\hat{\bm{\rho}}\cdot\mathbf{k}(t_{\infty}).$
(63)

Now that the saddle point is obtained, we proceed to deform the contour of
integration. As $\mathbf{k}$ is periodic with respect to $t$, the same period
is respected when extended into the complex plane, in a small strip around the
$t$ axis. We may then shift the contour into the complex plane where the
boundary contributions cancel due to periodicity when the endpoints have the
same imaginary part. We deform to a contour passing through the new stationary
point with the real part of the exponent constant and scaling as
$-\epsilon^{3/2}.$ Stationary phase then applies to the rapidly varying
imaginary part, producing the usual inverse square root contribution in two
dimensions. The difference is now that a small real part is contained in the
exponential, and we obtain the long-range behaviour:

$\displaystyle\Gamma\sim\frac{\exp\big{[}-\epsilon^{3/2}\rho+i\bm{\rho}\cdot\mathbf{k}_{\infty}\big{]}}{\sqrt{\rho}}.$
(64)

In two dimensions in particular, a small rotation of $\bm{\rho}_{\infty}$ by
an angle $\epsilon$ corresponds to a perpendicular shift by $\bm{\epsilon}$ to
leading order, so that the $\epsilon^{3/2}$ scaling of decay rate is also
observed with variance of angle. The above approximation provides the fit in
the main text. In higher dimensions, one may heuristically extend the above
line of reasoning and consider the higher dimensional analogue of steepest
descent to arrive at a scaling in the generic case:

$\displaystyle\Gamma\sim\frac{\exp\big{[}-\epsilon^{3/2}\rho+i\bm{\rho}\cdot\mathbf{k}_{\infty}\big{]}}{\rho^{(d-1)/2}}.$
(65)

However, the geometry of saddle points is far more complex in higher
dimensions, so that when the saddle point is degenerate, or multiple
$\bm{\epsilon}$ are available, the scaling of both $\rho$ and $\epsilon$ can
change, so that the above is limited to a heuristic argument when $d>2$. We
additionally note that the above analysis is limited to small perturbations,
as inherent branch points of the parametrisation will dominate if the
perturbations are too large.

###  B.4 Supplementary numerics

####  B.4.1 Green’s function elements in the honeycomb lattice

Figure 6: Basic analysis of the honeycomb lattice. (a) Dispersion relation,
displaying approximate isotropy at $\omega_{0}=2J.$ (b) Coherent part of
interactions when emitters are separated along a lattice vector
$\bm{\rho}=n\mathbf{a}_{1}.$ The approximate obtained in this work (5)
together with the classical approximation obtained by (7). Relative error is
shown in the inset.

Whilst the square lattice offers the clearest insight into dynamics, we here
present calculations in the honeycomb lattice to further verify our
theoretical results, and show that phenomena observed are not strictly limited
to the square lattice of the main investigation. In particular, we show that
degeneracies do not typically contribute to long-range interactions, and show
how dominantly incoherent interactions can be present in the honeycomb lattice
also.

The honeycomb tight-binding lattice (see Figure 7(c)) Hamiltonian exhibits
generic Dirac-point degeneracies between two bands. We consider QEs of
resonant frequency $\Delta$, coupled to sites in the same sub-lattice.
Assuming a nearest neighbor coupling of magnitude $J_{1}=J$ and next-nearest-
neighbor coupling $J_{2}=TJ$, we obtain the off-diagonal term:

$G=\frac{\mathcal{A}}{2}\sum_{\pm}\int_{\text{BZ}}\frac{d^{2}\mathbf{k}}{(2\pi)^{2}}\frac{e^{i\mathbf{k}\cdot\bm{\rho}}}{\Delta-\omega_{\pm}(\mathbf{k})},$
(66)

where $\omega_{\pm}(\mathbf{k})$ gives the eigenenergy in the Bloch wave basis
for the upper (lower) band, and can be found for $T=0$ in (say) González-
Tudela and Cirac (2018). We take $\Delta=2J,$ where the Born-Markov
approximation for two spins is accurateGonzález-Tudela and Cirac (2018). The
results with $T=0$ can be seen in Fig 6(b), whereby the calculation of the
exact expression (3) is compared with the real part of (5). As in the main
text, we find convergence up to evanescent contributions to the true value of
$G,$ so that the degeneracy does not contribute. We note that when near a
degeneracy, the interactions may become dominantly coherentGonzález-Tudela and
Cirac (2018), or even perfectly coherentGarcía-Elcano  _et al._ (2020) due to
vanishing density of states. In this case, the radius of the tubular region
would necessarily shrink to zero as we tune closer to the emitter, so that
dominant dynamics require a perturbative expansion localized at the
degeneracy, and beyond the Markov approximation.

Figure 7: Wave transport in the anisotropic honeycomb lattice, obtained by
effecting next-nearest-neighbour interactions with hopping strength $0.25J$
(denoted by the blue lines in (c), whilst nearest neighbour interactions are
denoted by orange.) (a) Dispersion relation of the lower branch of the
anisotropic honeycomb lattice. the resonant wavevectors corresponding to the
directions investigated are given by red dots, whilst the (higher order)
massive excitation is given by white dot $k_{\infty}(k_{\infty}^{HO})$. (b)
Decay of the Green’s function elements along the two directions corresponding
the divergences of the transverse effective mass. (c) The cross sectional
density of states. (d) The decay of the relative error between our
approximation and the exact result. (e) Ghost wave decay for the $x$ direction
(with maximal decay), and for the lattice vector $\mathbf{a}_{1}.$

The investigation corresponding to Figure 2 of the main text is given in
Figure 7, where the lattice is now taken with $T=0.25$. In this case an
emitter with $\Delta=-1.5J$ couples to higher-order massive particles,
displayed by the vanishing of both curvature and its derivative of $S$ as in
Figure 7(a). This results in a higher-order divergence in the usual
directional density of states in (c) for the corresponding vertical direction,
and a $n^{-1/4}$ power-law decay in (b), where the error in (d) is
exponentially decaying as predicted. In (e), the ghost wave decay is confirmed
also beyond the global caustics of the honeycomb lattice. Notably, along the
$x$ direction of maximum decay, the coherent interactions are vanishingly
small compared to incoherent interactions. Near the caustic, this is no longer
the case.

## References

  * Purcell (1995) E. M. Purcell, in  _Confined Electrons and Photons_ (Springer, 1995) pp. 839–839.
  * Akahane  _et al._ (2003) Y. Akahane, T. Asano, B.-S. Song, and S. Noda, nature 425, 944 (2003).
  * Khelif  _et al._ (2004) A. Khelif, A. Choujaa, S. Benchabane, B. Djafari-Rouhani, and V. Laude, Applied physics letters 84, 4400 (2004).
  * Krauss (2007) T. F. Krauss, Journal of Physics D: Applied Physics 40, 2666 (2007).
  * Rao and Hughes (2007) V. M. Rao and S. Hughes, Physical Review B 75, 205437 (2007).
  * Yoshie  _et al._ (2004) T. Yoshie, A. Scherer, J. Hendrickson, G. Khitrova, H. Gibbs, G. Rupper, C. Ell, O. Shchekin, and D. Deppe, Nature 432, 200 (2004).
  * Nakamura  _et al._ (2004) H. Nakamura, Y. Sugimoto, K. Kanamoto, N. Ikeda, Y. Tanaka, Y. Nakamura, S. Ohkouchi, Y. Watanabe, K. Inoue, H. Ishikawa, _et al._ , Optics Express 12, 6606 (2004).
  * Kotsifaki and Chormaic (2019) D. G. Kotsifaki and S. N. Chormaic, Nanophotonics 8, 1227 (2019).
  * Meng  _et al._ (2018) Y. Meng, A. Dareau, P. Schneeweiss, and A. Rauschenbeutel, Physical Review X 8, 031054 (2018).
  * Yu  _et al._ (2019) S.-P. Yu, J. A. Muniz, C.-L. Hung, and H. Kimble, Proceedings of the National Academy of Sciences 116, 12743 (2019).
  * Ziegler  _et al._ (2019) J. Ziegler, R. Klaiss, A. Blaikie, D. Miller, V. R. Horowitz, and B. J. Alemán, Nano letters 19, 2121 (2019).
  * Yu  _et al._ (2017) H. Yu, G.-B. Liu, J. Tang, X. Xu, and W. Yao, Science advances 3, e1701696 (2017).
  * Solano  _et al._ (2017) P. Solano, P. Barberis-Blostein, F. K. Fatemi, L. A. Orozco, and S. L. Rolston, Nature communications 8, 1 (2017).
  * Perczel  _et al._ (2020) J. Perczel, J. Borregaard, D. E. Chang, S. F. Yelin, and M. D. Lukin, Physical Review Letters 124, 083603 (2020).
  * Hood (2017) J. D. Hood,  _Atom-light interactions in a photonic crystal waveguide_ , Ph.D. thesis, California Institute of Technology (2017).
  * García-Elcano  _et al._ (2020) I. García-Elcano, A. González-Tudela, and J. Bravo-Abad, Physical Review Letters 125, 163602 (2020).
  * Zhou  _et al._ (2017) M. Zhou, L. Ying, L. Lu, L. Shi, J. Zi, and Z. Yu, Nature communications 8, 1 (2017).
  * Sheremet  _et al._ (2021) A. S. Sheremet, M. I. Petrov, I. V. Iorsh, A. V. Poshakinskiy, and A. N. Poddubny, arXiv preprint arXiv:2103.06824 (2021).
  * González-Tudela and Cirac (2017) A. González-Tudela and J. I. Cirac, Physical Review Letters 119, 143602 (2017).
  * Ying  _et al._ (2019) L. Ying, M. Zhou, M. Mattei, B. Liu, P. Campagnola, R. H. Goldsmith, and Z. Yu, Physical review letters 123, 173901 (2019).
  * González-Tudela  _et al._ (2015) A. González-Tudela, C.-L. Hung, D. E. Chang, J. I. Cirac, and H. Kimble, Nature Photonics 9, 320 (2015).
  * Shahmoon and Kurizki (2013) E. Shahmoon and G. Kurizki, Physical Review A 87, 033831 (2013).
  * Van Loo  _et al._ (2013a) A. F. Van Loo, A. Fedorov, K. Lalumiere, B. C. Sanders, A. Blais, and A. Wallraff, Science 342, 1494 (2013a).
  * Paulisch  _et al._ (2016) V. Paulisch, H. Kimble, and A. González-Tudela, New Journal of Physics 18, 043041 (2016).
  * Scully and Svidzinsky (2009) M. O. Scully and A. A. Svidzinsky, Science 325, 1510 (2009).
  * Douglas  _et al._ (2015) J. S. Douglas, H. Habibian, C.-L. Hung, A. V. Gorshkov, H. J. Kimble, and D. E. Chang, Nature Photonics 9, 326 (2015).
  * De Bernardis  _et al._ (2021) D. De Bernardis, Z.-P. Cian, I. Carusotto, M. Hafezi, and P. Rabl, Physical Review Letters 126, 103603 (2021).
  * Argüello-Luengo  _et al._ (2019) J. Argüello-Luengo, A. González-Tudela, T. Shi, P. Zoller, and J. I. Cirac, Nature 574, 215 (2019).
  * Kim  _et al._ (2021) E. Kim, X. Zhang, V. S. Ferreira, J. Banker, J. K. Iverson, A. Sipahigil, M. Bello, A. González-Tudela, M. Mirhosseini, and O. Painter, Physical Review X 11, 011015 (2021).
  * Yuan  _et al._ (2019) N. F. Yuan, H. Isobe, and L. Fu, Nature communications 10, 1 (2019).
  * González-Tudela and Cirac (2018) A. González-Tudela and J. I. Cirac, Physical Review A 97, 043831 (2018).
  * Kuś  _et al._ (1993) M. Kuś, F. Haake, and D. Delande, Physical review letters 71, 2167 (1993).
  * Narimanov (2019) E. E. Narimanov, Advanced Photonics 1, 046003 (2019).
  * Newman  _et al._ (2018) W. D. Newman, C. L. Cortes, A. Afshar, K. Cadien, A. Meldrum, R. Fedosejevs, and Z. Jacob, Science advances 4, eaar5278 (2018).
  * Danner  _et al._ (1994) G. Danner, W. Kang, and P. Chaikin, Physical review letters 72, 3714 (1994).
  * Hsu  _et al._ (2016) C. W. Hsu, B. Zhen, A. D. Stone, J. D. Joannopoulos, and M. Soljačić, Nature Reviews Materials 1, 1 (2016).
  * Chang  _et al._ (2012) D. E. Chang, L. Jiang, A. Gorshkov, and H. Kimble, New Journal of Physics 14, 063003 (2012).
  * Sinha  _et al._ (2020) K. Sinha, A. González-Tudela, Y. Lu, and P. Solano, arXiv preprint arXiv:2006.12569 (2020).
  * Zhao and Wang (2013) Y. Zhao and Z. Wang, Physical review letters 110, 240404 (2013).
  * Gardiner  _et al._ (2004) C. Gardiner, P. Zoller, and P. Zoller,  _Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics_ (Springer Science & Business Media, 2004).
  * Markos and Soukoulis (2008) P. Markos and C. M. Soukoulis,  _Wave propagation: from electrons to photonic crystals and left-handed materials_ (Princeton University Press, 2008).
  * Matlack  _et al._ (2018) K. H. Matlack, M. Serra-Garcia, A. Palermo, S. D. Huber, and C. Daraio, Nature materials 17, 323 (2018).
  * Köhl  _et al._ (2005) M. Köhl, H. Moritz, T. Stöferle, K. Günter, and T. Esslinger, Physical review letters 94, 080403 (2005).
  * Lidorikis  _et al._ (1998) E. Lidorikis, M. Sigalas, E. N. Economou, and C. Soukoulis, Physical review letters 81, 1405 (1998).
  * Bayindir  _et al._ (2000) M. Bayindir, B. Temelkuran, and E. Ozbay, Physical Review Letters 84, 2140 (2000).
  * Cerjan  _et al._ (2020) A. Cerjan, M. Jürgensen, W. A. Benalcazar, S. Mukherjee, and M. C. Rechtsman, Physical Review Letters 125, 213901 (2020).
  * Xiao  _et al._ (2015) M. Xiao, W.-J. Chen, W.-Y. He, and C. T. Chan, Nature Physics 11, 920 (2015).
  * Xue  _et al._ (2019) H. Xue, Y. Yang, G. Liu, F. Gao, Y. Chong, and B. Zhang, Physical review letters 122, 244301 (2019).
  * Sánchez-Burillo  _et al._ (2020) E. Sánchez-Burillo, D. Porras, and A. González-Tudela, Physical Review A 102, 013709 (2020).
  * Katznelson (2004) Y. Katznelson,  _An introduction to harmonic analysis_ (Cambridge University Press, 2004).
  * Brouder  _et al._ (2007) C. Brouder, G. Panati, M. Calandra, C. Mourougane, and N. Marzari, Physical review letters 98, 046402 (2007).
  * Alvarez (2017) O. Alvarez, arXiv preprint arXiv:1701.04446 (2017).
  * Gray (2012) A. Gray,  _Tubes_ , Vol. 221 (Birkhäuser, 2012).
  * Khudaverdian  _et al._ (2007) H. M. Khudaverdian _et al._ , Journal of Geometry and Symmetry in Physics 10, 29 (2007).
  * Choe and Koiso (2015) J. Choe and M. Koiso, Pacific Journal of Mathematics 280, 1 (2015).
  * (56) See supplemental material for details of the text derivations and more detailed discussion of the results presented in the main text .
  * (57) In this work an evanescent wave is understood in the sense of non-propagating: technically, we consider the wave decays superpolynomially .
  * Le Kien and Rauschenbeutel (2017) F. Le Kien and A. Rauschenbeutel, Physical Review A 95, 023838 (2017).
  * Pichler  _et al._ (2015) H. Pichler, T. Ramos, A. J. Daley, and P. Zoller, Physical Review A 91, 042116 (2015).
  * Stamnes (2002) J. J. Stamnes, Optics express 10, 740 (2002).
  * Zakharevich (2017) V. Zakharevich, Journal of Mathematical Physics 58, 083506 (2017).
  * Camley and Maradudin (1983) R. Camley and A. Maradudin, Physical Review B 27, 1959 (1983).
  * Chigrin (2004) D. N. Chigrin, Physical Review E 70, 056611 (2004).
  * Zhang and Mølmer (2019) Y.-X. Zhang and K. Mølmer, Physical review letters 122, 203605 (2019).
  * Van Loo  _et al._ (2013b) A. F. Van Loo, A. Fedorov, K. Lalumiere, B. C. Sanders, A. Blais, and A. Wallraff, Science 342, 1494 (2013b).
  * Asenjo-Garcia  _et al._ (2017) A. Asenjo-Garcia, M. Moreno-Cardoner, A. Albrecht, H. Kimble, and D. E. Chang, Physical Review X 7, 031024 (2017).
  * Kornovan  _et al._ (2019) D. Kornovan, N. Corzo, J. Laurat, and A. Sheremet, Physical Review A 100, 063832 (2019).
  * Eldredge  _et al._ (2016) Z. Eldredge, P. Solano, D. Chang, and A. V. Gorshkov, Physical Review A 94, 053855 (2016).
  * Bykov (1975) V. P. Bykov, Soviet Journal of Quantum Electronics 4, 861 (1975).
  * John and Wang (1990) S. John and J. Wang, Physical review letters 64, 2418 (1990).
  * Yang  _et al._ (2019) Y. Yang, W. Gao, L. Xia, H. Cheng, H. Jia, Y. Xiang, and S. Zhang, Physical review letters 123, 033901 (2019).
  * Cortes and Jacob (2017) C. L. Cortes and Z. Jacob, Nature communications 8, 1 (2017).
  * Gosselin  _et al._ (2006) P. Gosselin, F. Ménas, A. Bérard, and H. Mohrbach, EPL (Europhysics Letters) 76, 651 (2006).
  * Fang  _et al._ (2012) K. Fang, Z. Yu, and S. Fan, Nature photonics 6, 782 (2012).
  * Page (1930) L. Page, Physical Review 36, 444 (1930).
  * Leonforte  _et al._ (2021) L. Leonforte, A. Carollo, and F. Ciccarello, Physical Review Letters 126, 063601 (2021).
  * McGloin and Dholakia (2005) D. McGloin and K. Dholakia, Contemporary Physics 46, 15 (2005).
  * Broky  _et al._ (2008) J. Broky, G. A. Siviloglou, A. Dogariu, and D. N. Christodoulides, Optics express 16, 12880 (2008).
  * He  _et al._ (2019) P. He, C.-H. Hsu, S. Shi, K. Cai, J. Wang, Q. Wang, G. Eda, H. Lin, V. M. Pereira, and H. Yang, Nature communications 10, 1 (2019).
  * Rackauckas and Nie (2017) C. Rackauckas and Q. Nie, Journal of Open Research Software 5 (2017).
  * (81) S. G. Johnson, “QuadGK.jl: Gauss–Kronrod integration in Julia,” https://github.com/JuliaMath/QuadGK.jl.
  * Bachrathy and Stépán (2012) D. Bachrathy and G. Stépán, Periodica Polytechnica Mechanical Engineering 56, 81 (2012).
  * Des Cloizeaux (1964) J. Des Cloizeaux, Physical Review 135, A698 (1964).
  * Kato (2013) T. Kato,  _Perturbation theory for linear operators_ , Vol. 132 (Springer Science & Business Media, 2013).
  * Johnson (2015) S. G. Johnson, arXiv preprint arXiv:1508.04376 (2015).
  * Verpoort (2008) S. Verpoort, (2008).
  * Perczel and Lukin (2020) J. Perczel and M. D. Lukin, Physical Review A 101, 033822 (2020).
  * Katsura  _et al._ (1971) S. Katsura, T. Morita, S. Inawashiro, T. Horiguchi, and Y. Abe, Journal of Mathematical Physics 12, 892 (1971).
  * Izumiya  _et al._ (2016) S. Izumiya, M. d. C. Romero Fuster, M. A. S. Ruas, and F. Tari,  _Differential geometry from a singularity theory viewpoint_ (World Scientific, 2016).


aainstitutetext: Niels Bohr International Academy, Niels Bohr Institute,
University of Copenhagen,
Blegdamsvej 17, 2100, Copenhagen, Denmark††institutetext:  bInstituto Física
Teórica, UAM/CSIC, Calle Nicola’s Cabrera 13-15, Cantoblanco E-28049 Madrid,
Spain. & Research Center for Cosmic Neutrinos, Institute for Cosmic Ray
Research, University of Tokyo, Kashiwa, Chiba 277-8582, Japan ccinstitutetext:
Theoretical Physics Department, Fermi National Accelerator Laboratory, P. O.
Box 500, Batavia, IL 60510, USA

# Compact Perturbative Expressions for Neutrino Oscillations in Matter: II

Peter B. Denton b Hisakazu Minakata c Stephen J. Parke<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper we rewrite the neutrino mixing angles and mass squared
differences in matter given in our original paper, Denton:2016wmg , in a
notation that is more conventional for the reader. Replacing the usual
neutrino mixing angles and mass squared differences in the expressions for the
vacuum oscillation probabilities with these matter mixing angles and mass
squared differences gives an excellent approximation to the oscillation
probabilities in matter. Comparisons for T2K, NOvA, T2HKK and DUNE are also
given for neutrinos and anti-neutrinos, disappearance and appearance channels,
normal ordering and inverted ordering.

###### Keywords:

Neutrino Physics, Neutrino Oscillations in Matter, CP violation

††preprint: FERMILAB-PUB-18-014-T

##  1 Introduction

In this paper we rewrite the neutrino mixing angles and mass squared
differences in matter given in our original paper, Denton:2016wmg , in a
notation that is more conventional for the reader. Replacing the usual
neutrino mixing angles and mass squared differences in the expressions for the
vacuum oscillation probabilities with these matter mixing angles and mass
squared differences gives an excellent approximation to the oscillation
probabilities in matter. Higher orders are also easily calculated and provide
several orders of magnitude improvement per order.

In Section 2, we give the approximation to the mixing angles and mass squared
difference in matter and discuss how to use these to calculated the
oscillation probabilities in matter both at 0th order and 1st order. We also
give expansions of the mixing angles and mass squared differences in matter in
powers of $(a/\Delta m^{2})$. In Section 3, we make a detailed comparison
between the exact and the approximate oscillation probabilities in matter for
the T2K, NOvA, T2HKK and DUNE experiments. Section 4 is the Summary.

##  2 Mixing Angles and Mass Differences in Matter

###  2.1 Zeroth Order

In this section, a simple and accurate way to evaluate oscillation
probabilities, recently shown in Denton:2016wmg , is given.111In this note
$\phi$, $\psi$ and $\Delta\lambda_{jk}$ of Denton:2016wmg , are replaced with
the more traditional notation $\widetilde{\theta}_{13}$ and
$\widetilde{\theta}_{12}$ and $\Delta\,\widetilde{m^{2}}_{jk}$ respectively.
Details as to the why’s and how’s of this method are contained in that paper.

The mixing angles in matter, which we denote by a $\widetilde{\theta}_{13}$
and $\widetilde{\theta}_{12}$ here, can also be calculated in the following
way, using $\Delta m^{2}_{ee}\equiv\cos^{2}\theta_{12}\Delta
m^{2}_{31}+\sin^{2}\theta_{12}\Delta m^{2}_{32}$, as follows222Vacuum values
to be used in calculating $\Delta m^{2}_{ee}$., see Parke:2018brr :

$\displaystyle\cos 2\widetilde{\theta}_{13}$ $\displaystyle=$
$\displaystyle\frac{(\cos 2\theta_{13}-a/\Delta m^{2}_{ee})}{\sqrt{(\cos
2\theta_{13}-a/\Delta m^{2}_{ee})^{2}+\sin^{2}2\theta_{13}~{}}},$ (2.1.1)

where $~{}~{}a\equiv 2\sqrt{2}G_{F}N_{e}E_{\nu}~{}~{}$ is the standard matter
potential, and

$\displaystyle\cos 2\widetilde{\theta}_{12}$ $\displaystyle=$
$\displaystyle\frac{(\cos 2\theta_{12}-a^{\,\prime}/\Delta
m^{2}_{21})}{\sqrt{(\cos 2\theta_{12}-a^{\,\prime}/\Delta
m^{2}_{21})^{2}~{}+~{}\sin^{2}2\theta_{12}\cos^{2}(\widetilde{\theta}_{13}-\theta_{13})~{}~{}}},$
(2.1.2)

where $~{}~{}a^{\,\prime}\equiv a\,\cos^{2}\widetilde{\theta}_{13}+\Delta
m^{2}_{ee}\sin^{2}(\widetilde{\theta}_{13}-\theta_{13})~{}~{}$ is the
$\theta_{13}$-modified matter potential for the 1-2 sector. In these two
flavor rotations, both $\widetilde{\theta}_{13}$ and $\widetilde{\theta}_{12}$
are in range $[0,\pi/2]$.

$\theta_{23}$ and $\delta$ are unchanged in matter for this approximation.

The neutrino mass squared differences in matter, i.e. the $\Delta m^{2}_{jk}$
in matter, which we denote by $\Delta\,\widetilde{m^{2}}_{jk}$, are given by

$\displaystyle\Delta\,\widetilde{m^{2}}_{21}$ $\displaystyle=$
$\displaystyle\Delta m^{2}_{21}\,\sqrt{(\cos 2\theta_{12}-a^{\,\prime}/\Delta
m^{2}_{21})^{2}~{}+~{}\sin^{2}2\theta_{12}\cos^{2}(\widetilde{\theta}_{13}-\theta_{13})~{}~{}},$
$\displaystyle\Delta\,\widetilde{m^{2}}_{31}$ $\displaystyle=$
$\displaystyle\Delta
m^{2}_{31}+(\,a-\frac{3}{2}a^{\,\prime}\,)+\frac{1}{2}\,\left(\,\Delta\widetilde{m^{2}}_{21}-\Delta
m^{2}_{21}~{}\right),$ (2.1.3) $\displaystyle\Delta\,\widetilde{m^{2}}_{32}$
$\displaystyle=$
$\displaystyle\Delta\,\widetilde{m^{2}}_{31}-\Delta\,\widetilde{m^{2}}_{21}$

Note that the same square root333If $a=0$, then
$\widetilde{\theta}_{13}=\theta_{13}$ and since $a^{\prime}=0$ then
$\widetilde{\theta}_{12}=\theta_{12}$ and both $\sqrt{\cdots}=1$, also
$\Delta\,\widetilde{m^{2}}_{jk}=\Delta m^{2}_{jk}$ for all $(j,k)$ as
required. The identity $s^{2}_{\theta}=(1-\cos 2\theta)/2$ is useful for
calculating both $s_{\theta}$ and $c_{\theta}$.  appears in both
$\Delta\,\widetilde{m^{2}}_{21}$ and $\sin^{2}\widetilde{\theta}_{12}$. To see
that the $\Delta\,\widetilde{m^{2}}_{31}$ and $\Delta\,\widetilde{m^{2}}_{32}$
have the right asymptotic forms, use the fact that
$(\Delta\,\widetilde{m^{2}}_{21}-\Delta m^{2}_{21})=|a^{\,\prime}|+{\cal
O}(\Delta m^{2}_{21})$, for $|a|\gg\Delta m^{2}_{21}$.

In Fig. 1 and Fig. 2 the values of a, a${}^{\,\prime}$,
$\sin^{2}\widetilde{\theta}_{13}$, $\sin^{2}\widetilde{\theta}_{12}$,
$\widetilde{m^{2}}_{j}$ and $\Delta\widetilde{m^{2}}_{jk}$ as a function of
the neutrino energy for a density of 3.0 g.cm-3.

To calculate the oscillation probabilities, to 0th order, use the above
$\Delta\,\widetilde{m^{2}}_{jk}$ instead of $\Delta m^{2}_{jk}$ and replace
the vacuum MNS matrix as follows

$\displaystyle U^{0}_{MNS}\equiv
U_{23}(\theta_{23})\,U_{13}(\theta_{13},\delta)\,U_{12}(\theta_{12})$
$\displaystyle\Rightarrow$ $\displaystyle U^{M}_{MNS}\equiv
U_{23}(\theta_{23})\,U_{13}(\
\widetilde{\theta}_{13},\delta)\,U_{12}(\widetilde{\theta}_{12}).$

That is, replace

$\displaystyle\Delta m^{2}_{jk}$ $\displaystyle\rightarrow$
$\displaystyle\Delta\,\widetilde{m^{2}}_{jk}$ $\displaystyle\theta_{13}$
$\displaystyle\rightarrow$ $\displaystyle\widetilde{\theta}_{13}$
$\displaystyle\theta_{12}$ $\displaystyle\rightarrow$
$\displaystyle\widetilde{\theta}_{12},$ (2.1.4)

$\theta_{23}$ and $\delta$ remain unchanged, it is that simple. We call this
the 0th order DMP approximation.

These expressions are valid for both NO, $\Delta m^{2}_{ee}>0$ and IO, $\Delta
m^{2}_{ee}<0$. For anti-neutrinos, just change the sign of $a$ and $\delta$.
Our expansion parameter is

$\displaystyle\left|\,\sin(\widetilde{\theta}_{13}-\theta_{13})~{}\sin\theta_{12}\cos\theta_{12}~{}\frac{\Delta
m^{2}_{21}}{\Delta m^{2}_{ee}}\,\right|\leq 0.015,$ (2.1.5)

which is small and vanishes in vacuum, so that our perturbation theory
reproduces the vacuum oscillation probabilities exactly.

If $P_{\nu_{\alpha}\rightarrow\nu_{\beta}}(\Delta m^{2}_{31},\Delta
m^{2}_{21},\theta_{13},\theta_{12},\theta_{23},\delta)$ is the oscillation
probability in vacuum then
$P_{\nu_{\alpha}\rightarrow\nu_{\beta}}(\Delta\,\widetilde{m^{2}}_{31},\Delta\,\widetilde{m^{2}}_{21},\widetilde{\theta}_{13},\widetilde{\theta}_{12},\theta_{23},\delta)$
is the oscillation probability in matter, i.e. use the same function but
replace the mass squared differences and mixing angles with the matter values
given in eq. 2.1.1 \- 2.1.3. The resulting oscillation probabilities are
identical to the zeroth order approximation given in Denton, Minakata and
Parke, Denton:2016wmg .

Figure 1: In the normal ordering (NO): Top left, the matter potentials, $a$
and $a^{\,\prime}$, top right, sine squared of mixing angles in matter,
$\sin^{2}\widetilde{\theta}_{jk}$, bottom left, the mass squared eigenvalues
in matter, $\widetilde{m^{2}}_{j}$, and bottom right, the mass squared
differences in matter, $\Delta\,\widetilde{m^{2}}_{jk}$. $E_{\nu}\geq
0~{}(E_{\nu}\leq 0)$ is for neutrinos (anti-neutrinos). $E_{\nu}=0$ is the
vacuum values for both neutrinos and anti-neutrinos.

Figure 2: In the inverted ordering (IO): Top left, the matter potentials, $a$
and $a^{\,\prime}$, top right, sine squared of mixing angles in matter,
$\sin^{2}\widetilde{\theta}_{jk}$, bottom left, the mass squared eigenvalues
in matter, $\widetilde{m^{2}}_{j}$, and bottom right, the mass squared
differences in matter, $\Delta\,\widetilde{m^{2}}_{jk}$. $E_{\nu}\geq
0~{}(E_{\nu}\leq 0)$ is for neutrinos (anti-neutrinos). $E_{\nu}=0$ is the
vacuum values for both neutrinos and anti-neutrinos.

###  2.2 Higher Orders

If the 0th order is not accurate enough, going to 1st order is simple and
gives another two orders of magnitude in accuracy. First the
$\Delta\,\widetilde{m^{2}}_{jk}$ remain unchanged but the mixing matrix is
modified by

$\displaystyle U^{M}_{MNS}$ $\displaystyle\Rightarrow$ $\displaystyle V\equiv
U^{M}_{MNS}(1+W_{1}),$ (2.2.1)

where the matrix $W_{1}$ is given by

$\displaystyle W_{1}$ $\displaystyle=$
$\displaystyle\sin(\widetilde{\theta}_{13}-\theta_{13})~{}s_{12}c_{12}~{}\Delta
m^{2}_{21}$ (2.2.5)
$\displaystyle\left(\begin{array}[]{ccc}0&0&-\widetilde{s}_{12}\,e^{-i\delta}/\Delta\,\widetilde{m^{2}}_{31}\\\
0&0&+\widetilde{c}_{12}\,e^{-i\delta}/\ \Delta\,\widetilde{m^{2}}_{32}\\\
+\widetilde{s}_{12}\,e^{+i\delta}/\Delta\,\widetilde{m^{2}}_{31}&-\widetilde{c}_{12}\,e^{+i\delta}/\Delta\,\widetilde{m^{2}}_{32}&0\end{array}\right).$

where $\widetilde{s}_{12}=\sin\widetilde{\theta}_{12}$ and
$\widetilde{c}_{12}=\cos\widetilde{\theta}_{12}$ etc. The
$\Delta\,\widetilde{m^{2}}_{jk}$ and the $V$-mixing matrix can be used to
calculate the oscillation probabilities and improves the accuracy by two
orders of magnitude. We call this the 1st order DMP approximation. The next
highest order, 2nd order, is also discussed in Denton:2016wmg .

###  2.3 Expansions in $a/\Delta m^{2}$

If $|a|\ll|\Delta m^{2}_{ee}|$, that is when $E_{\nu}\ll
11\left(\,\rho\,/\,3~{}{\rm g\,cm^{-3}}\right)$ GeV, we have,

$\displaystyle\sin^{2}\widetilde{\theta}_{13}$ $\displaystyle\approx$
$\displaystyle s^{2}_{13}\biggr{[}1+2c^{2}_{13}(a/\Delta
m^{2}_{ee})+3(c^{2}_{13}-s^{2}_{13})c^{2}_{13}(a/\Delta m^{2}_{ee})^{2}+{\cal
O}(a/\Delta m^{2}_{ee})^{3})\biggr{]}$
$\displaystyle\sin^{2}(\widetilde{\theta}_{13}-\theta_{13})$
$\displaystyle\approx$ $\displaystyle s^{2}_{13}c^{2}_{13}(a/\Delta
m^{2}_{ee})^{2}\biggr{[}1+2(c^{2}_{13}-s^{2}_{13})(a/\Delta m^{2}_{ee})+{\cal
O}(a/\Delta m^{2}_{ee})^{2}\biggr{]}$ (2.3.1) $\displaystyle a^{\,\prime}$
$\displaystyle\approx$ $\displaystyle
ac^{2}_{13}\biggr{[}1-s^{2}_{13}(a/\Delta
m^{2}_{ee})-s^{2}_{13}(c^{2}_{13}-s^{2}_{13})(a/\Delta m^{2}_{ee})^{2}+{\cal
O}(a/\Delta m^{2}_{ee})^{3}\biggr{]}$

up to ${\cal O}(a/\Delta m^{2}_{ee})^{2}$. The expansion for $a^{\,\prime}$
can be used to calculate $\Delta\,\widetilde{m^{2}}_{31}$ as follows,

$\displaystyle\Delta\,\widetilde{m^{2}}_{31}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\Delta
m^{2}_{31}+~{}(\,a-a^{\,\prime}\,)+\frac{1}{2}\,\left[\,\Delta\widetilde{m^{2}}_{21}-\Delta
m^{2}_{21}-a^{\,\prime}~{}\right],&~{}~{}a,a^{\,\prime}>0\\\\[5.69054pt]
\Delta
m^{2}_{31}+(\,a-2a^{\,\prime}\,)+\frac{1}{2}\,\left[\,\Delta\widetilde{m^{2}}_{21}-\Delta
m^{2}_{21}+a^{\,\prime}~{}\right],&~{}~{}a,a^{\,\prime}<0\end{array}\right.$
(2.3.4)

where the quantities in $[\cdots]$ is of ${\cal O}(\Delta m^{2}_{21})$ for all
values of $E_{\nu}$.

As can be seen from Fig. 1 and Fig. 2, both $\Delta\,\widetilde{m^{2}}_{21}$
and $\sin^{2}\widetilde{\theta}_{12}$ make rapid changes in +150 MeV region.
Well away from this region, when $|a|\gg|\Delta m^{2}_{21}|$, that is
$E_{\nu}\gg 150\left(\,\rho\,/\,3~{}{\rm g\,cm^{-3}}\right)$ MeV, we can write

$\displaystyle\Delta\,\widetilde{m^{2}}_{21}$ $\displaystyle\approx$
$\displaystyle|\,a^{\,\prime}-\Delta m^{2}_{21}\cos 2\theta_{12}\,|\,,$
(2.3.5)

for $|\,~{}a^{\,\prime}-\Delta m^{2}_{21}\cos 2\theta_{12}\,|\gg\Delta
m^{2}_{21}$. This can be used to obtain the asymptotic values for neutrino
mass squareds in matter, which agree with the values given in Denton:2016wmg .

##  3 Oscillation Probabilities

###  3.1 Comparisons

Neutrino parameters relevant for oscillations:

$\displaystyle\Delta m^{2}_{ee}=\pm~{}2.5\times 10^{-3}~{}{\rm eV}^{2},$
$\displaystyle\Delta m^{2}_{21}=+~{}7.5\times 10^{-5}~{}{\rm eV}^{2}$
$\displaystyle\sin^{2}\theta_{12}=0.31,$
$\displaystyle\sin^{2}\theta_{23}=0.43$
$\displaystyle\sin^{2}\theta_{13}=0.022,$
$\displaystyle\delta=-72^{\circ}=-2\pi/5$ (3.1.1)

Where $\Delta m^{2}_{ee}>0$ gives a normal ordering (NO) neutrino spectrum and
$\Delta m^{2}_{ee}<0$ for inverted ordering (IO). Note we have avoided the
special points: $\theta_{23}=\pi/4$ as well as $\delta=0,\pm\pi/2,\pi$, so as
not to overestimate the precision.

We consider four experimental setups: to be comprehensive, the energy windows
are chosen to be wider than that accessible for a particular experiment.

  * •

T2K & T2HK: with baseline, L = 295 km, neutrino energy 0.2 $<E_{\nu}/{\rm
GeV}<$ 3.0,
and density, $\rho=$ 2.3 g.cm-3. See Fig. 3, 4.

  * •

NOvA: with baseline, L = 810 km, neutrino energy 0.6 $<E_{\nu}/{\rm GeV}<$
4.0,
and density, $\rho=$ 3.0 g.cm-3. See Fig. 5, 6.

  * •

T2HKK: with baseline, L = 1050 km, neutrino energy 0.3 $<E_{\nu}/{\rm GeV}<$
5.0,
and density, $\rho=$ 3.0 g.cm-3. See Fig. 7, Fig. 8.

  * •

DUNE: with baseline, L = 1300 km, neutrino energy 0.5 $<E_{\nu}/{\rm GeV}<$
7.0,
and density, $\rho=$ 3.0 g.cm-3. See Fig. 9, 10.

In these figures we have considered the channels $\nu_{\mu}$ disappearance and
$\nu_{e}$ appearance for neutrinos and anti-neutrinos and for both NO and IO.
Each figure consists of three panels: the top panel is the exact oscillation
probabilities in matter from Zaglauer and Schwarzer, Zaglauer:1988gz , as well
as the exact vacuum oscillation probability. The middle (bottom) panel shows
the difference (fractional difference) between

  1. 1.

the exact, Zaglauer:1988gz , and vacuum oscillation probabilities (black),

  2. 2.

the exact and the 0th order DMP approximation, Denton:2016wmg  (red),

  3. 3.

the exact and the 1st order DMP approximation, Denton:2016wmg  (green).

In the fractional differences, the denominator is the average of the two
probabilities being compared. In Fig, 3 to Fig. 10, the “dips” in middle and
bottom panels appear when $\Delta P$ changes sign.

In Table 1, we give the maximum difference and fractional difference of the
0th order approximation to the exact probability.

| T2K/HK | NOvA | T2HKK | DUNE
---|---|---|---|---
max $\Delta P$ | $10^{-5}$ | $10^{-4}$ | $10^{-4}$ | $10^{-4}$
max $\Delta P/P$ | $10^{-3}$ | $10^{-3}$ | $10^{-2.5}$ | $10^{-2}$
Table 1: The maximum $\Delta P$ and $\Delta P/P$ at 0th order in the DMP
approximation. The largest fraction difference occurs at oscillation maximum
for $\nu_{\mu}$ disappearance channel, where the oscillation probability is a
few %. None of the experiments included here, T2K, NOvA, T2HKK and DUNE will
be within an order a magnitude of being sensitive to any of these differences.

##  4 Summary

In summary, the simple 0th order approximation of DMP, Denton:2016wmg , is
sufficiently accurate for all of the accelerator based neutrino oscillation
experiments operating or planned: T2K, NOvA, T2HKK and DUNE. First order,
which is also simple to use, improves the accuracy by a further two orders of
magnitude for these experiments.

Figure 3: T2K, for normal ordering (NO): Top Left figure is $\nu_{\mu}$
disappearance, Top Right figure is $\bar{\nu}_{\mu}$ disappearance, Bottom
Left figure is $\nu_{\mu}\rightarrow\nu_{e}$ appearance, and Bottom Right is
$\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ appearance. In each figure, the top
panel is exact oscillation probability in matter , $P^{ex}_{mat}$, from
Zaglauer:1988gz , and in vacuum, $P_{vac}$. The Middle panel is difference
between exact oscillation probabilities in matter and vacuum (black), and the
difference between exact and 0th (red) and exact and 1st (green)
approximations to the matter probabilities using the DMP scheme,
Denton:2016wmg . Bottom panel is similar to middle panel but plotting the
fractional differences, $\Delta P/\overline{P}$. The density use is 2.3
g.cm-3.

Figure 4: T2K, for inverted ordering (IO): Top Left figure is $\nu_{\mu}$
disappearance, Top Right figure is $\bar{\nu}_{\mu}$ disappearance, Bottom
Left figure is $\nu_{\mu}\rightarrow\nu_{e}$ appearance, and Bottom Right is
$\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ appearance. In each figure, the top
panel is exact oscillation probability in matter , $P^{ex}_{mat}$, from
Zaglauer:1988gz , and in vacuum, $P_{vac}$. The Middle panel is difference
between exact oscillation probabilities in matter and vacuum (black), and the
difference between exact and 0th (red) and exact and 1st (green)
approximations to the matter probabilities using the DMP scheme,
Denton:2016wmg . Bottom panel is similar to middle panel but plotting the
fractional differences, $\Delta P/\overline{P}$. The density use is 2.3
g.cm-3.

Figure 5: NOvA, for normal ordering (NO): Top Left figure is $\nu_{\mu}$
disappearance, Top Right figure is $\bar{\nu}_{\mu}$ disappearance, Bottom
Left figure is $\nu_{\mu}\rightarrow\nu_{e}$ appearance, and Bottom Right is
$\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ appearance. In each figure, the top
panel is exact oscillation probability in matter , $P^{ex}_{mat}$, from
Zaglauer:1988gz , and in vacuum, $P_{vac}$. The Middle panel is difference
between exact oscillation probabilities in matter and vacuum (black), and the
difference between exact and 0th (red) and exact and 1st (green)
approximations to the matter probabilities using the DMP scheme,
Denton:2016wmg . Bottom panel is similar to middle panel but plotting the
fractional differences, $\Delta P/\overline{P}$. The density use is 3.0
g.cm-3.

Figure 6: NOvA, for inverted ordering (IO): Top Left figure is $\nu_{\mu}$
disappearance, Top Right figure is $\bar{\nu}_{\mu}$ disappearance, Bottom
Left figure is $\nu_{\mu}\rightarrow\nu_{e}$ appearance, and Bottom Right is
$\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ appearance. In each figure, the top
panel is exact oscillation probability in matter , $P^{ex}_{mat}$, from
Zaglauer:1988gz , and in vacuum, $P_{vac}$. The Middle panel is difference
between exact oscillation probabilities in matter and vacuum (black), and the
difference between exact and 0th (red) and exact and 1st (green)
approximations to the matter probabilities using the DMP scheme,
Denton:2016wmg . Bottom panel is similar to middle panel but plotting the
fractional differences, $\Delta P/\overline{P}$. The density use is 3.0
g.cm-3.

Figure 7: T2HKK, for normal ordering (NO): Top Left figure is $\nu_{\mu}$
disappearance, Top Right figure is $\bar{\nu}_{\mu}$ disappearance, Bottom
Left figure is $\nu_{\mu}\rightarrow\nu_{e}$ appearance, and Bottom Right is
$\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ appearance. In each figure, the top
panel is exact oscillation probability in matter , $P^{ex}_{mat}$, from
Zaglauer:1988gz , and in vacuum, $P_{vac}$. The Middle panel is difference
between exact oscillation probabilities in matter and vacuum (black), and the
difference between exact and 0th (red) and exact and 1st (green)
approximations to the matter probabilities using the DMP scheme,
Denton:2016wmg . Bottom panel is similar to middle panel but plotting the
fractional differences, $\Delta P/\overline{P}$. The density use is 3.0
g.cm-3.

Figure 8: T2HKK, for inverted ordering (IO): Top Left figure is $\nu_{\mu}$
disappearance, Top Right figure is $\bar{\nu}_{\mu}$ disappearance, Bottom
Left figure is $\nu_{\mu}\rightarrow\nu_{e}$ appearance, and Bottom Right is
$\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ appearance. In each figure, the top
panel is exact oscillation probability in matter , $P^{ex}_{mat}$, from
Zaglauer:1988gz , and in vacuum, $P_{vac}$. The Middle panel is difference
between exact oscillation probabilities in matter and vacuum (black), and the
difference between exact and 0th (red) and exact and 1st (green)
approximations to the matter probabilities using the DMP scheme,
Denton:2016wmg . Bottom panel is similar to middle panel but plotting the
fractional differences, $\Delta P/\overline{P}$. The density use is 3.0
g.cm-3.

Figure 9: DUNE, for normal ordering (NO): Top Left figure is $\nu_{\mu}$
disappearance, Top Right figure is $\bar{\nu}_{\mu}$ disappearance, Bottom
Left figure is $\nu_{\mu}\rightarrow\nu_{e}$ appearance, and Bottom Right is
$\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ appearance. In each figure, the top
panel is exact oscillation probability in matter , $P^{ex}_{mat}$, from
Zaglauer:1988gz , and in vacuum, $P_{vac}$. The Middle panel is difference
between exact oscillation probabilities in matter and vacuum (black), and the
difference, $\Delta P$, between exact and 0th (red) and exact and 1st (green)
approximations to the matter probabilities using the DMP scheme,
Denton:2016wmg . Bottom panel is similar to middle panel but plotting the
fractional differences, $\Delta P/\overline{P}$. The density use is 3.0
g.cm-3.

Figure 10: DUNE, for inverted ordering (IO): Top Left figure is $\nu_{\mu}$
disappearance, Top Right figure is $\bar{\nu}_{\mu}$ disappearance, Bottom
Left figure is $\nu_{\mu}\rightarrow\nu_{e}$ appearance, and Bottom Right is
$\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ appearance. In each figure, the top
panel is exact oscillation probability in matter , $P^{ex}_{mat}$, from
Zaglauer:1988gz , and in vacuum, $P_{vac}$. The Middle panel is difference
between exact oscillation probabilities in matter and vacuum (black), and the
difference between exact and 0th (red) and exact and 1st (green)
approximations to the matter probabilities using the DMP scheme,
Denton:2016wmg . Bottom panel is similar to middle panel but plotting the
fractional differences, $\Delta P/\overline{P}$. The density use is 3.0
g.cm-3.

##  5 Acknowledgements

This manuscript has been authored by Fermi Research Alliance, LLC under
Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of
Science, Office of High Energy Physics.

This project has received funding/support from the European Union’s Horizon
2020 research and innovation programme under the Marie Sklodowska-Curie grant
agreement No 690575. This project has received funding/support from the
European Union’s Horizon 2020 research and innovation programme under the
Marie Sklodowska-Curie grant agreement No 674896.

HM is supported by Instituto Física Teórica, UAM/CSIC in Madrid, via
“Theoretical challenges of new high energy, astro and cosmo experimental data”
project, Ref: 201650E082.

PBD acknowledges support from the Villum Foundation (Project No. 13164) and
the Danish National Research Foundation (DNRF91 and Grant No. 1041811001).

## References

  * (1) P. B. Denton, H. Minakata and S. J. Parke, “Compact Perturbative Expressions For Neutrino Oscillations in Matter,” JHEP 1606, 051 (2016) doi:10.1007/JHEP06(2016)051 [arXiv:1604.08167 [hep-ph]].
  * (2) S. J. Parke, P. B. Denton and H. Minakata, “Analytic Neutrino Oscillation Probabilities in Matter: Revisited,” arXiv:1801.00752 [hep-ph].
S.J. Parke and M.D. Messier, “Cross Check of NOvA Oscillation Probabilities”
(1/12/2018), NOVA-doc-25833, FERMILAB-FN-1049,
http://doi.org/10.5281/zenodo.1146747

  * (3) H. W. Zaglauer and K. H. Schwarzer, “The Mixing Angles in Matter for Three Generations of Neutrinos and the MSW Mechanism,” Z. Phys. C 40 (1988) 273.