\section{\label{app:long-range}Derivation of the long-range behaviour of the Green's function}

\subsection{\label{app:long-range:nodegen}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 \cite{bloch_analytic} 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 
\begin{equation}
    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}')}{\Delta - \omega_{\nu}(\mathbf{k})}
\end{equation}
becomes non-smooth \cite{katoperturbation2013}, 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 \cite{sanchez_limits}. Recall the unitary transformation between eigenfunctions $\psi_{i}$ and the Bloch wave corresponding to an individual sublattice:

\begin{equation}
    |\psi_{\nu}(\mathbf{k})\rangle = \sum_{j}U_{\nu j}(\mathbf{k})|\mathbf{k}_{j}\rangle,
\end{equation}
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:

\begin{equation}
    \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}'}{(2\pi)^d}e^{i\mathbf{k}\cdot \mathbf{r}_{i}}U_{\nu j}(\mathbf{k})\langle \mathbf{k}'_{i}|\mathbf{k}_{j}\rangle = e^{i\mathbf{k}\cdot \mathbf{r}_{i}}U_{\nu i}(\mathbf{k}).
\end{equation}
    Assuming $\mathbf{r}$ and $\mathbf{r}'$ belong to lattices $i$ and $i'$ respectively, we rewrite the integrand
    
\begin{equation}
    I = \sum_{\nu}\int_{BZ}\frac{d^{d}\mathbf{k}}{(2\pi)^d}\frac{e^{i\mathbf{k}\cdot \bm{\rho}}U_{i' \nu}^{*}(\mathbf{k})U_{\nu i}(\mathbf{k})}{\Delta - \omega_{\nu}(\mathbf{k})},
\end{equation}
or more succinctly in matrix notation
\begin{equation}
    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' i}.
\end{equation}
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 

\begin{equation}
       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' i} =  \int_{BZ}\frac{d^{d}\mathbf{k}}{(2\pi)^d},
\end{equation}
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 \cite{katoperturbation2013}, 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.
 
 
 
 \subsection{\label{app:long-range:coherent-part} 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:

\begin{equation}
\label{appeq:greenfunction}
	{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'})}{\Delta - \omega(\mathbf{k}) + i0},
\end{equation}
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:

\begin{eqnarray}
\label{appeq:om}{\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'})}{\Delta - \omega(\mathbf{k})},
\end{eqnarray}
and proceed with an asymptotic large $\bm{\rho} = \mathbf{r} - \mathbf{r}'$ treatment, where the direction of $\bm{\rho}$ is kept fixed. Having shown analyticity of the Kernel around degeneracies in Appendix  \ref{app:long-range:nodegen}, we now address the  (near) resonant contributions from $S$. Under the assumption of locality in Appendix \ref{app:long-range:nodegen}, 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 \cite{2015BumpFourier} 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:

\begin{equation}
{\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})},
\end{equation}
where we increase the separation with $\mathbf{r}, \mathbf{r}'$ fixed modulo the lattice, so that $\Phi(\mathbf{k}) = \psi(\mathbf{k},\mathbf{r})\psi^{*}(\mathbf{k},\mathbf{r'})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 \textit{tubular neighbourhood} of $S$:

\begin{equation}
	\mathbf{k} = (t,\bm{\xi}) = \bm{\xi} + t\hat{\mathbf{n}}({\bm{\xi}}),
\end{equation}
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$ \cite{WeylVolume2017}, 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:

\begin{equation}
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})},
\end{equation}
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$:

\begin{equation}
    \phi_{\delta S}(\mathbf{k}) = \phi_{\epsilon}(t),
\end{equation}

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

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

We now change coordinates systems, giving
the volume element obtained in the derivation of\textit{Weyl's tube formula} \cite{WeylVolume2017,gray_tubes,tube_backup,tube_backup_two}:
\begin{equation}
	d^{d}\mathbf{k} = \det (\mathbf{I} + t \mathbf{K})d\bm{\xi}dt,
\end{equation}

where $\mathbf{K}$ is the \textit{second fundamental form} of the surface $S$. In the case of an implicit surface we have \cite{secondform}

\begin{equation}
	\mathbf{K} = \mathbf{K}(\bm{\xi}) = -\frac{\mathbf{H}_{T}(\bm{\xi})}{v(\bm{\xi})},
\end{equation}

for the Hessian 

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

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

\begin{equation}
	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}),
\end{equation}


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

\begin{equation}
	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}).
\end{equation}

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

\begin{equation}
	\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}),
\end{equation}
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 phase\cite{StatPhase2017}. 


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

\begin{equation}
	\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'})d^{d-1} \mathbf{k}}{v(\mathbf{k})(2\pi)^{d-1}} + O(\rho^{-\infty}),
\end{equation}

where we have relabelled $S \ni \bm{\xi}$ with $\mathbf{k} \in S$ and where we have the phase function 
\begin{equation}
\label{appeq:phasefunc}
	\Theta_{\rho }(x) = \int_{-\infty}^{\infty} \frac{dt \sin(tx)\phi_{\rho \epsilon}(t)}{\pi t}.
\end{equation}
\begin{figure}[tb]
    \centering
    \includegraphics[scale = 0.4]{app_fig_1.eps}
    \caption{The Fourier transform $\Theta_{\rho}(x)$ for different values of $\rho$. Due to the identical scaling of $\Theta_{\rho}$ with $\epsilon$ (in Equation \eqref{appeq:phasefunc}), 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)).$}
    \label{app:fig:bump}
\end{figure}
$\Theta_{\rho}$ acts as a sigmoid function, varying smoothly from -1 at $-\infty$ to 1 at $\infty$ (see Figure \ref{app:fig:bump}), and as the integral of the Fourier transform of the bump function \cite{2015BumpFourier}, has the asymptotic behaviour

\begin{equation}
\label{app:bumpasymp}
    \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)
\end{equation}
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:$
\begin{equation}
	\label{appeq:bestapprox}
	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'})d^{d-1} \mathbf{k}}{v(\mathbf{k})(2\pi)^{d-1}},
\end{equation}
where the amplitude of the main text is given by

\begin{equation}
    \Pi_{\rho}(
	\hat{\bm{\rho}}\cdot\hat{\mathbf{n}}) = \frac{(1+\Theta_{\rho}(
	\hat{\bm{\rho}}\cdot\hat{\mathbf{n}}))}{2}
\end{equation}

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 \cite{tunable_tudela,tudela_dirac,lukin_green_1,cirac_noperturb2017}.

    
\subsection{\label{app:long-range:ibp} Calculation of the exact Green's function}
The strongly singular manifold present in \eqref{eq:greenfunction} of the main text presents significant calculational difficultly, so that to obtain exact results either recurrence relations \cite{horiguchi} 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 \eqref{eq:gam}. In one dimension repeated integration by parts on periodic functions $f, \omega$ gives

\begin{equation}\label{appeq:ibp1}
        \mathcal{P}\int \frac{dk f(k)}{\Delta - \omega(k)} = -\mathcal{P}\int dk\log|\Delta - \omega(k)| \bigg(\frac{f(k)}{v(k)}\bigg)' = \mathcal{P}\int dk(\Delta - \omega(k))\big[\log|\Delta - \omega(k)| - 1\big]
        \Bigg(\frac{1}{v(k)}\bigg(\frac{
    f'(k)}{v(k)}\bigg)'\Bigg)',
\end{equation}
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 \eqref{appeq:ibp1}, 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

\begin{equation}\label{appeq:ibp2}
        \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),
\end{equation}
and applying integration by parts again for the final result. The bump function may be used here to regularize if necessary.

\section{\label{app:anisotropicemission}Characteristics of anisotropic emission}

\subsection{\label{app:anisotropicemission:ddos}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 \textit{Gauss map}, whose Jacobian is given as the Gaussian curvature. The change of coordinates is then
\begin{equation}
\label{appeq:area_element}
    d^{d-1}\mathbf{k} = \frac{d\Omega}{|K(\mathbf{k})|},
\end{equation}
for the intrinsic Gaussian curvature of $S$ at $\mathbf{k}.$ We may then decompose the formula for density of states:
\begin{equation}
\label{appeq:doscrosssection}
    D(\Delta) = \int_{S^{d-1}}\frac{d\Omega}{(2\pi)^{d}v(\mathbf{k})|K(\mathbf{k})|},
\end{equation}

so that the \textit{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):

\begin{equation}
\label{appeq:curvature_dos}
    D(\hat{\mathbf{v}},\Delta)  = \frac{\mathcal{A}}{(2\pi)^{d}v(\mathbf{k})|K(\mathbf{k})|}.
\end{equation}
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
\begin{equation}
	[\mathbf{m}^{-1}]_{ij} = \frac{\partial^2 \omega(\mathbf{k})}{\partial k_{i} \partial k_{j}}.
\end{equation}
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 \begin{equation}
    \mathbf{J}_{\hat{\mathbf{n}}}(\mathbf{k}) = \frac{\mathbf{H}(\mathbf{k})[\mathbf{I} - \hat{\mathbf{n}}\hat{\mathbf{n}}]}{v(\mathbf{k})},    
\end{equation}
which equals
\begin{equation}
    \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})},    
\end{equation}
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 
\begin{eqnarray}
    K(\mathbf{k}) = -\det[\mathbf{J}_{\hat{\mathbf{n}}}(\mathbf{k})],
\end{eqnarray}
leading to the state density  with respect to the cross section:
\begin{equation}
\label{appeq:crosssection}
    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|,
\end{equation}

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

\label{app:long-time}
We begin with the long time population of bath modes as $t \to \infty$\cite{cirac_noperturb2017,scatter_weyl_two}:
\begin{equation}
\label{appeq:bath_long}
    C_{\mathbf{k}} \to -\frac{ge^{-i\omega(\mathbf{k})t}}{\omega(\mathbf{k}) - \Delta' + i\gamma/2}\sim -\frac{ge^{-i\omega(\mathbf{k})t}}{\omega(\mathbf{k}) - \Delta' + i\gamma/2},
\end{equation}
where the frequency shift has been absorbed into $\Delta \to \Delta'$, and $\gamma = g^2 \Gamma$. The bath population is maximized when $\omega(\mathbf{k}) = \Delta'$. In particular, the $\mathbf{k}$-space cross section is obtained by \cite{scatter_weyl_two,scatter_weyl}
\begin{equation}
    \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')^2 + \frac{\gamma^2}{4}}
\end{equation}

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') \sim S(\Delta)$, obtained through the coordinate transform:

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

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




\subsection{\label{app:anisotropic-emission:stationary-phase}Stationary phase approximation on the resonant level set}


To evaluate the asymptotics 

\begin{equation}
	\label{app:bestapprox}
	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'})d^{d-1} \mathbf{k}}{v(\mathbf{k})(2\pi)^{d-1}},
\end{equation}

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 manifold\cite{StatPhase2017}:
\begin{equation}
    \label{appeq:manifoldstatphase}
    \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}),
\end{equation}
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 \eqref{appeq:bestapprox}, 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 \cite{height2016}, 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 \eqref{appeq:bestapprox} gives the final long-range approximation:

\begin{equation}
    G_{R} = i\mathcal{A}\frac{\psi(\mathbf{k}_{0},\mathbf{r})\psi^{*}(\mathbf{k}_{0},\mathbf{r'})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},
\end{equation}

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

\begin{equation}
	\label{app:stat_phase}
	\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}).
\end{equation}
    
    
\subsection{Asymptotics beyond the critical angle\label{app:anisotropic-emission: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 
    \begin{equation}
        \label{appeq:gam}{\Gamma}  = \mathcal{A} \int_{S}\frac{d\mathbf{k}}{(2\pi )}\frac{\psi(\mathbf{k},\mathbf{r})\psi^{*}(\mathbf{k},\mathbf{r'})}{{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})}.
    \end{equation}
    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
    
    \begin{eqnarray}
        &\hat{\bm{\epsilon}} \cdot \mathbf{v}(\mathbf{k}_{\infty}) = 0, \\
        &{m}^{T}(\mathbf{k}_{\infty}) \epsilon = 0, \\
        &\bm{\rho}_{\infty} \parallel \mathbf{v}(\mathbf{k}_{\infty}).
    \end{eqnarray}
    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
    
    \begin{eqnarray}
    (\hat{\bm{\rho}} + \hat{\bm{\epsilon}}) \cdot \mathbf{k}'(t) = 0. 
    \end{eqnarray}

    For small $\epsilon,$ we expect the stationary point to comprise a correspondingly small correcting complex part $t'$ in $t$ and thus $\mathbf{k}.$ We expand 
    
    \begin{equation}
        0 = \hat{\bm{\rho}} \cdot \mathbf{k}'(t) = \hat{\bm{\rho}} \cdot \mathbf{k}'(t_{\infty} + t') = (\hat{\bm{\rho}}_{\infty} + \hat{\bm{\epsilon}}) \cdot \left[\mathbf{k}'(t_{\infty}) + \mathbf{k}''(t_{\infty}) t' + \frac{t'^2}{2} \mathbf{k}'''(t_{\infty}) \right].
    \end{equation}
        To lowest order in ${\epsilon}$ and $t',$ we obtain 
        
    \begin{equation}
        0 = \hat{\bm{\epsilon}} \cdot \mathbf{k}'(t_{\infty}) + \frac{t'^2}{2}  \hat{\bm{\rho}}_{\infty} \cdot \mathbf{k}'''(t_{\infty}). 
    \end{equation}
    
    Under assumption of simple caustic, the third derivative is non-zero and the above equation is balanced for purely imaginary correction $t'$:  $i\sqrt{\epsilon} \sim t'$. The normalized argument of the exponent at the new stationary point has a real and imaginary part to leading order
    
  
    \begin{eqnarray}
    &\mathfrak{R}\big[ i\hat{\bm{\rho}} \cdot \mathbf{k}(t_{\infty} + t')\big] \sim -\epsilon^{3/2}, \\
    &\mathfrak{I}\big[i\hat{\bm{\rho}} \cdot \mathbf{k}(t_{\infty} + t')\big] \sim \hat{\bm{\rho}} \cdot \mathbf{k}(t_{\infty}).
    \end{eqnarray}

    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:
    
    \begin{eqnarray}
    \Gamma \sim \frac{\exp\big[-\epsilon^{3/2}\rho + i\bm{\rho} \cdot \mathbf{k}_{\infty}\big]}{\sqrt{\rho}}.
    \end{eqnarray}
    
    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:
    
    \begin{eqnarray}
    \Gamma \sim \frac{\exp\big[-\epsilon^{3/2}\rho + i\bm{\rho} \cdot \mathbf{k}_{\infty}\big]}{\rho^{(d-1)/2}}.
    \end{eqnarray}
    
    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.
    
  

\subsection{Supplementary numerics}
\subsubsection{Green's function elements in the honeycomb lattice}
\label{app:additional-numerics:honeycomb}

\begin{figure}
    \centering
    \includegraphics[scale=0.25]{app_fig_2.eps}
    \caption{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 \eqref{eq:medrange} together with the classical approximation obtained by \eqref{eq:longrange}. Relative error is shown in the inset.}
    \label{appfig:honey_compare}
\end{figure}

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 \ref{appfig:anisotropy-honey}(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:
\begin{equation}
    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})},
\end{equation}
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) \cite{tudela_dirac}. We take $\Delta = 2J,$ where the Born-Markov approximation for two spins is accurate\cite{tudela_dirac}. The results with $T=0$ can be seen in Fig \ref{appfig:honey_compare}(b), whereby the calculation of the exact expression \eqref{eq:om} is compared with the real part of \eqref{eq:medrange}. 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 coherent\cite{tudela_dirac}, or even perfectly coherent\cite{tunable_tudela} 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.

\begin{figure}[ht!]
    \centering
    \includegraphics[scale=0.40]{app_fig_3.eps}
    \caption{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}.$}
\label{appfig:anisotropy-honey}
\end{figure}


The investigation corresponding to Figure \ref{fig:anisotropy_square} of the main text is given in Figure \ref{appfig:anisotropy-honey}, 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 \ref{appfig:anisotropy-honey}(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.


\clearpage
\newpage

\end{widetext}
