Extended range of dipole-dipole interactions in periodically structured photonic media

TL;DR

This paper demonstrates that structured photonic media can significantly extend the range of dipole-dipole interactions between quantum systems, enabling long-range quantum effects through the manipulation of isosurfaces.

Contribution

It introduces a general relation linking interaction range to isosurface properties in structured photonic media, revealing a new method for mediating long-range quantum interactions.

Findings

Extended interaction range is robust and independent of dipole orientation.

Diminishing momentum isosurfaces enable long-range quantum interactions.

A general relation between isosurface properties and interaction range is established.

Abstract

The interaction between quantum two-level systems is typically short-range in free space and most photonic environments. Here we show that diminishing momentum isosurfaces with equal frequencies can create a significantly extended range of interaction between distant quantum systems. The extended range is robust and does not rely on a specific location or orientation of the transition dipoles. A general relation between the interaction range and properties of the isosurface is described for structured photonic media. It provides a new way to mediate long-range quantum behavior.

Equations

$$\begin{split}\mathcal{H}=&\sum_{i=\mathrm{1,2}}{\omega}_{0}{{\hat{\sigma}}^{\dagger}_{i}\hat{\sigma}}_{i}+\sum_{\bm{\mathrm{k}}}{{\omega}_{\bm{\mathrm{k}}}{\hat{a}}^{\dagger}_{\bm{\mathrm{k}}}{\hat{a}}_{\bm{\mathrm{k}}}}\\ +&i\sum_{i\mathrm{=1,2}}\sum_{\mathbf{k}}\Big{[}ig_{\mathbf{k}}(\mathbf{r}_{i})\left(\hat{\sigma}_{i}^{{\dagger}}+\hat{\sigma}_{i}\right)\hat{a}_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{r}_{i}}+\mathrm{H.c.}\Big{]},\end{split}$$

$${\mathrm{\Gamma }}_{\mathrm{Re}}={\iint }_{S_{{\omega }_0\left(\mathrm{k}\right)}}{{\rho }_{\mathbf{k}}e}^{i\mathrm{k}\cdot \mathrm{R}}d{S}_{\mathrm{k}}.$$

$$\mathcal{H}={\mathcal{H}}_{\mathrm{ph}}+{\mathcal{H}}_{\mathrm{tls}}+{\mathcal{H}}_{\mathrm{int}}.$$

$$\begin{split}\mathcal{H}_{\mathrm{ph}}=&\sum_{\mathbf{k},\eta}\omega_{\mathbf{k}}\hat{a}_{\mathbf{k},\eta}^{{\dagger}}\hat{a}_{\mathbf{k},\eta}\\ \mathcal{H}_{\mathrm{tls}}=&\sum_{i=1,2}\omega_{0}\hat{\sigma}_{i}^{{\dagger}}\hat{\sigma}_{i}\\ \mathcal{H}_{\mathrm{int}}=&\sum_{i=1,2}\sum_{\mathbf{k},\eta}\Big{[}ig_{\mathbf{k},\eta}(\mathbf{r}_{i})\left(\hat{\sigma}_{i}^{{\dagger}}+\hat{\sigma}_{i}\right)\hat{a}_{\mathbf{k},\eta}e^{i\mathbf{k}\cdot\mathbf{r}_{i}}+\mathrm{H.c.}\Big{]},\end{split}$$

$${\mathcal{M}}_{\mathrm{FI}}=⟨\mathrm{F}|{\mathcal{H}}_{\mathrm{int}}|\mathrm{I}⟩+\sum _{\alpha }\frac{⟨\mathrm{F}|{\mathcal{H}}_{\mathrm{int}}|{\mathrm{R}}_{\alpha }⟩⟨{\mathrm{R}}_{\alpha }|{\mathcal{H}}_{\mathrm{int}}|\mathrm{I}⟩}{{\mathcal{E}}_{\mathrm{I}}-{\mathcal{E}}_{{\mathrm{R}}_{\alpha }}}+\cdots .$$

$$\begin{split}\mathcal{M}_{\mathrm{FI}}=\sum_{{\mathbf{k}},\eta}\Big{(}&g_{\mathbf{k},\eta}\left(\mathbf{r}_{1}\right)^{\ast}g_{\mathbf{k},\eta}\left(\mathbf{r}_{2}\right)\frac{e^{i\mathbf{k}\cdot\mathbf{R}}}{\omega_{\mathbf{k}}-\omega_{0}}\\ +&g_{\mathbf{k},\eta}\left(\mathbf{r}_{1}\right)g_{\mathbf{k},\eta}\left(\mathbf{r}_{2}\right)^{\ast}\frac{e^{-i\mathbf{k}\cdot\mathbf{R}}}{\omega_{\mathbf{k}}+\omega_{0}}\Big{)},\end{split}$$

$$\begin{split}\Gamma=&\ \Gamma_{\mathrm{Re}}+i\Gamma_{\mathrm{Im}}\\ \\ =&\sum_{\eta}\iiint_{\mathcal{V}_{\mathbf{k}}}d^{3}\mathbf{k}\frac{v_{g}(\mathbf{k})}{\pi}\Bigg{[}\pi\delta\left(\omega_{\mathbf{k}}-\omega_{0}\right)\rho_{\mathbf{k},\eta}e^{i\bm{\mathrm{k}}\cdot{{\mathbf{R}}}}\\ &\ \ \ \ \ \ \ \ \ +i\mathbb{P}\left(\rho_{\mathbf{k},\eta}\frac{e^{i\bm{\mathrm{k}}\cdot{{\mathbf{R}}}}}{\omega_{\mathbf{k}}-\omega_{0}}+\rho_{\mathbf{k},\eta}^{\ast}\frac{e^{-i\mathbf{k}\cdot\mathbf{R}}}{\omega_{\mathbf{k}}+\omega_{0}}\right)\Bigg{]},\end{split}$$

$${\rho }_{\mathbf{k},\eta }=\frac{\omega }{16{\pi }^2{\epsilon }_0{v}_g(\mathbf{k})}{\left({\text{muup}}_1\cdot {\text{epsilonup}}_{\mathbf{k},\eta }\right)}^{\ast }\left({\text{muup}}_2\cdot {\text{epsilonup}}_{\mathbf{k},\eta }\right).$$

$${\Gamma }_{\mathrm{Re}}\left({\omega }_0\right)=\sum _{\eta }{\iint }_{{\mathcal{S}}_{{\omega }_0(\mathbf{k})}}{\rho }_{\mathbf{k},\eta }\left({\omega }_0\right)e^{i\mathbf{k}\cdot \mathbf{R}}d{\mathcal{S}}_{\mathbf{k}},$$

$$\begin{array}{rl}{\Gamma }_{\mathrm{Im}}\left({\omega }_0\right)=& \frac{1}{\pi }\mathbb{P}{\int }_0^{\infty }d\omega \sum _{\eta }{\iint }_{{\mathcal{S}}_{\omega (\mathbf{k})}}d{S}_{\mathbf{k}}\\ & \times \left(\frac{{\rho }_{\mathbf{k},\eta }}{\omega -{\omega }_0}{e}^{i\mathrm{k}\cdot \mathbf{R}}+\frac{{\rho }_{\mathbf{k},\eta }^{\ast }}{\omega +{\omega }_0}{e}^{-i\mathrm{k}\cdot \mathbf{R}}\right)\\ & \\ =& \frac{1}{\pi }\mathbb{P}{\int }_0^{\infty }d\omega \left(\frac{{\Gamma }_{\mathrm{Re}}(\omega )}{\omega -{\omega }_0}+\frac{{\Gamma }_{\mathrm{Re}}^{\ast }(\omega )}{\omega +{\omega }_0}\right).\end{array}$$

$${\omega }_{\mathbf{k}}=c|\mathbf{k}|=ck,$$

$$\begin{array}{rl}{\Gamma }_{\mathrm{Re}}=& \frac{{\mu }_1{\mu }_2{k}^3}{16{\pi }^2{\epsilon }_0}\sum _{\eta }{\int }_0^{2\pi }d\theta {\int }_0^{\pi }d\phi \sin \phi \\ \times & {\left({\text{muup}}_1\cdot {\text{epsilonup}}_{\mathbf{k},\eta }\right)}^{\ast }\left({\text{muup}}_2\cdot {\text{epsilonup}}_{\mathbf{k},\eta }\right)e^{ikR{\xi }_{\mathbf{k}}}.\end{array}$$

$$\sum _{\eta }{ϵ}_{\mathbf{k},i}^{(\eta )}{ϵ}_{\mathbf{k},j}^{(\eta )}={\delta }_{12}-{\hat{k}}_1{\hat{k}}_2$$

$$\begin{split}\Gamma_{\mathrm{Re}}=&\frac{\mu_{1}\mu_{2}k^{3}}{16\pi^{2}\varepsilon_{0}}\int_{0}^{2\pi}d\theta\int_{0}^{\pi}d\varphi\left(\delta_{12}-{\hat{k}}_{1}{\hat{k}}_{2}\right)\sin{\varphi}e^{ikR\xi_{\mathbf{k}}}\\ =&\frac{\mu_{1}\mu_{2}k}{16\pi^{2}\varepsilon_{0}}\left(-\nabla^{2}\delta_{12}+\nabla_{1}\nabla_{2}\right)\int_{0}^{2\pi}d\theta\int_{0}^{\pi}d\varphi\sin{\varphi}e^{ikR\xi_{\mathbf{k}}}\\ =&\frac{\mu_{1}\mu_{2}k}{4\pi\varepsilon_{0}}\left(-\nabla^{2}\delta_{12}+\nabla_{1}\nabla_{2}\right)\frac{\sin{kR}}{R}\\ =&\frac{\mu_{1}\mu_{2}k^{3}}{4\pi\varepsilon_{0}}\Bigg{(}\left(\delta_{12}-\hat{{R}}_{i}\hat{{R}}_{j}\right)\frac{\sin{kR}}{kR}\\ &\ \ \ \ \ \ \ \ \ \ +\left(\delta_{12}-3\hat{{R}}_{i}\hat{{R}}_{j}\right)\left(\frac{\cos{kR}}{(kR)^{2}}+\frac{\sin{kR}}{(kR)^{3}}\right)\Bigg{)},\end{split}$$

$$\begin{array}{rl}{\Gamma }_{\mathrm{Im}}=& \frac{{\mu }_1{\mu }_2}{16{\pi }^3{\epsilon }_0}\left(-{\nabla }^2{\delta }_{12}+{\nabla }_1{\nabla }_2\right){\int }_0^{2\pi }d\theta {\int }_0^{\pi }d\phi \sin \phi \\ & \times \mathbb{P}{\int }_0^{\infty }dkk\left(\frac{e^{ikR{\xi }_{\mathbf{k}}}}{k-k_0}+\frac{e^{-ikR{\xi }_{\mathbf{k}}}}{k+k_0}\right).\end{array}$$

$$\begin{split}\Gamma_{\mathrm{Im}}=&\frac{\mu_{1}\mu_{2}k_{0}}{4\pi\varepsilon_{0}}\left(-\nabla^{2}\delta_{12}+\nabla_{1}\nabla_{2}\right)\frac{\cos{k_{0}R}}{R}\\ =&\frac{\mu_{1}\mu_{2}k_{0}^{3}}{4\pi\varepsilon_{0}}\Bigg{(}\left(\delta_{12}-\hat{{R}}_{i}\hat{{R}}_{j}\right)-\frac{\cos{k_{0}R}}{k_{0}R}\\ &\ \ \ \ \ +\left(\delta_{12}-3\hat{{R}}_{i}\hat{{R}}_{j}\right)\left(\frac{\sin{k_{0}R}}{(k_{0}R)^{2}}+\frac{\cos{k_{0}R}}{(k_{0}R)^{3}}\right)\Bigg{)}.\end{split}$$

$${\mathcal{H}}_{\mathrm{wp}}\left(\mathbf{k}\right)=\sum _{i=x,y,z}{v}_i{q}_i{\sigma }_i,$$

$${\omega }_{\mathbf{q}}={\omega }_{\mathrm{wp}}\pm v|\mathbf{q}|.$$

$$\begin{array}{rl}{\Gamma }_{\mathrm{Re}}(\omega )=& \frac{{\mu }_1{\mu }_2\omega }{16{\pi }^2{\epsilon }_{0}v}{e}^{i{\mathbf{k}}_{\mathrm{c}}\cdot \mathbf{R}}\\ & \times {\iint }_{{\mathcal{S}}_{\omega \left(\mathbf{q}\right)}}d{\mathcal{S}}_{\mathbf{q}}{\left({\widehat{\text{muup}}}_1\cdot {\text{epsilonup}}_{\mathbf{k},\eta }\right)}^{\ast }\left({\widehat{\text{muup}}}_2\cdot {\text{epsilonup}}_{\mathbf{k},\eta }\right)e^{i\mathbf{q}\cdot \mathbf{R}}.\end{array}$$

$${\rho }_{\mathbf{k},\eta }=\frac{{\mu }_1{\mu }_2\omega }{16{\pi }^2{\epsilon }_{0}v}{\left({\widehat{\text{muup}}}_1\cdot {\text{epsilonup}}_{\mathbf{k},\eta }\right)}^{\ast }\left({\widehat{\text{muup}}}_2\cdot {\text{epsilonup}}_{\mathbf{k},\eta }\right)\simeq \rho .$$

$$\begin{array}{rl}{\Gamma }_{\mathrm{Re}}\simeq & \rho e^{i{\mathbf{k}}_{\mathrm{c}}\cdot \mathbf{R}}{\iint }_{{\mathcal{S}}_{{\omega }_0\left(\mathbf{q}\right)}}d{\mathcal{S}}_{\mathbf{q}}{e}^{i\mathbf{q}\cdot \mathbf{R}}\\ =& 4\pi q^2\rho e^{i{\mathbf{k}}_{\mathrm{c}}\cdot \mathbf{R}}\frac{\sin qR}{qR}.\end{array}$$

$$\begin{split}\Gamma_{\mathrm{Re}}=&B^{\left(3D\right)}\frac{k^{2}}{c}\Bigg{(}\left({\delta}_{12}-{\hat{R}}_{1}{\hat{R}}_{2}\right)\ \frac{{\mathrm{\ sin}kR\ }}{kR}+\\ &\left(\delta_{12}-3{\hat{R}}_{1}{\hat{R}}_{2}\right)\ \left(\frac{{\mathrm{\ cos}kR\ }}{{\left(kR\right)}^{2}}+\frac{{\mathrm{\ sin}kR\ }}{{\left(kR\right)}^{3}}\right)\ \Bigg{)},\end{split}$$

$${\Gamma }_{\mathrm{Re}}\left(R>\lambda \right)\cong B^{\left(3D\right)}\left({\delta }_{12}-{\hat{R}}_1{\hat{R}}_2\right)\frac{k^2}{c}\frac{\sin kR}{kR}.$$

$$\begin{array}{rl}{\mathrm{\Gamma }}_{\mathrm{Re}}& \cong \frac{{\mu }_1{\mu }_2\omega }{16{\pi }^2{\epsilon }_0}\frac{1}{v}{e}^{i{\mathrm{k}}_c\cdot \mathrm{R}}{\iint }_{S_{\omega \left(\mathrm{q}\right)}}d{S}_{\mathrm{q}}{e}^{i\mathrm{q}\cdot \mathrm{R}}\\ & =B^{\left(3D\right)}\frac{q^2}{v}\frac{\sin qR}{qR}{e}^{i{\mathrm{k}}_c\cdot \mathrm{R}}.\end{array}$$

$${\omega }_{\mathbf{k}}=\beta {\left|\mathrm{k}-{\mathrm{k}}_{\mathrm{c}}\right|}^2+{\omega }_c=\beta q^2+{\omega }_c,$$

$$\begin{split}\Gamma_{\mathrm{Re}}=&\frac{\omega}{16{\pi}^{2}{\varepsilon}_{0}}\frac{q}{2\beta}\mathop{\int\!\!\!\!\int}\nolimits_{S_{\omega\left(\bm{\mathrm{q}}\right)}}{dS_{\bm{\mathrm{q}}}{\left({\bm{\mathrm{\muup}}}_{1}\cdot{\bm{\mathrm{\epsilonup}}}_{\bm{\mathrm{k}}}\right){\left({\bm{\mathrm{\muup}}}_{2}\cdot{\bm{\mathrm{\epsilonup}}}_{\bm{\mathrm{k}}}\right)}^{*\ }e}^{i\bm{\mathrm{q}}\cdot\bm{\mathrm{R}}}}\\ =&B^{\left(3D\right)}\frac{q}{2\beta}\Bigg{(}\left({\delta}_{12}-{\hat{R}}_{1}{\hat{R}}_{2}\right)\ \frac{{\mathrm{\ sin}qR}}{qR}+\\ &\left(\delta_{12}-3{\hat{R}}_{1}{\hat{R}}_{2}\right)\ \left(\frac{{\mathrm{cos}qR\ }}{{\left(qR\right)}^{2}}+\frac{{\mathrm{sin}qR\ }}{{\left(qR\right)}^{3}}\right)\ \Bigg{)}.\end{split}$$

$${\Gamma }_{\mathrm{Re}}\left(R>\lambda \right)\cong B^{\left(3D\right)}\left({\delta }_{12}-{\hat{R}}_1{\hat{R}}_2\right)\frac{q}{2\beta }\frac{\sin qR}{qR}.$$

$${\Gamma }_{\mathrm{Re}}\left(\omega \right)=\frac{\omega }{8\pi {\epsilon }_0\hslash }{\int }_{{\ell }_{\omega \left(\mathrm{k}\right)}}d{\ell }_{\mathrm{k}}\left({\text{muup}}_1\cdot {\text{epsilonup}}_{\mathrm{k}}\right){\left({\text{muup}}_2\cdot {\text{epsilonup}}_{\mathrm{k}}\right)}^{\ast }\frac{e^{ik\cdot R}}{v_g(\mathrm{k})},$$

$$\begin{array}{rl}{\Gamma }_{\mathrm{Re}}\left(\omega \right)& =\frac{{\mu }_1{\mu }_2\omega k}{4{\epsilon }_{0}c}{\int }_0^{2\text{piup}}d\theta e^{ikR\cos \theta }\\ & =B^{\left(2D\right)}\frac{k}{c}{J}_0\left(kR\right),\end{array}$$

$${\omega }_{\mathrm{k}}=v\left|\mathrm{k}-{\mathbf{k}}_{\mathrm{c}}\right|+{\omega }_c=vq+{\omega }_c.$$

Full text

The resonant dipole-dipole interaction between two quantum two-level systems (TLS) is typically short-range. There has been strong interest in realizing long-range interactions to exploit collective physics such as superradiance Scully and Svidzinsky (2009); Solano et al. (2017), collective frequency shift Meir et al. (2014), F$\ddot{\mathrm{o}}$rster resonance energy transfer Clegg (1995); Garcia-Vidal and Feist (2017), and quantum entanglement Van Loo et al. (2013); Burkard and Imamoglu (2006); Petrosyan and Fleischhauer (2008); Mingaleev et al. (2000); Gonzalez-Tudela et al. (2011); Shahmoon and Kurizki (2013); Hood et al. (2016). The ability to modulate the distance dependence of these processes could have potential applications in quantum information processing Imamog et al. (1999); Petrosyan and Fleischhauer (2008) and energy conversion Maxwell et al. (2013). Two components contribute to the interaction: the evanescent near fields and the propagating far fields (Fig. 1a&b). To enable long-range interaction from the evanescent fields, one could use evanescent fields with a long tail, such as defect modes in the photonic bandgap Kurizki (1990); Douglas et al. (2015); Notararigo et al. (2018). However, it is less obvious how to engineer propagating far fields to enable long-range interaction. It is the goal of this letter to provide a new perspective to understand the general physical mechanism that is responsible for long-range interaction induced by propagating far fields, and identify photonic structures that are capable of extending the interaction range.

In free space, the range of far-field interaction is limited to the wavelength scale. When the wavelength is long, such as in index-near-zero materials Fleury and Alu (2013); Mahmoud et al. (2017); Liberal and Engheta (2018); Gundogdu et al. (2015); Serebryannikov et al. (2019), the interaction range can increase proportionally. However, there are a few intriguing examples where the interaction range extends beyond the effective wavelength. These include low-dimensional spaces, such as photonic crystal waveguides and fibers Sato et al. (2012); Le Kien and Rauschenbeutel (2017); Solano et al. (2017); Lecamp et al. (2007); Hughes (2007); Yao and Hughes (2009); Minkov and Savona (2013); Hung et al. (2013); Vasco et al. (2014), or hyperbolic materials in selected directions Biehs et al. (2016); Cortes and Jacob (2017). These interesting but isolated examples heavily rely on very specific configurations. Thus, it is difficult to generalize the theoretical treatments to identify the underlying physics, which unfortunately remains elusive. In this letter, we show the deep connection between the interaction range and the size and shape of the isofrequency surface in momentum space. It can be generalized to a broad range of physical systems and can reveal new systems capable of realizing long-range interactions.

We begin by examining the interaction between two TLSs over a long distance. The TLSs are embedded in a photonic environment that can be described by a dispersion relation $\omega=\omega\left(\mathbf{k}\right)$. For example, in free space, $\omega=c|\mathbf{k}|=ck$, where $c$ is the speed of light. Other dispersion relations can be seen in metamaterials, photonic crystals or waveguides. In general, the Hamiltonian of the TLSs and the photonic modes is given by Bay et al. (1997)

[TABLE]

where $\omega_{0}$ is the resonant transition frequency of TLSs. $\hat{\sigma}^{\dagger}_{i}(\hat{\sigma}_{i})$ is the raising (lowering) operator of $i$th TLS. $\omega_{\mathbf{k}}$ and $\hat{a}^{\dagger}_{\mathbf{k}}\left({\hat{a}}_{\mathbf{k}}\right)$ are the frequency and creation (annihilation) operator of photons, respectively. $g_{\mathbf{k}}\left(\mathbf{r}_{i}\right)=\sqrt{\omega_{k}/2{\varepsilon}_{0}V}{\bm{\mathrm{\muup}}}_{i}\cdot{\bm{\mathrm{\epsilonup}}}_{\mathbf{k}}$ is the coupling between the $i$th TLS and the photonic mode $\mathbf{k}$, where ${\bm{\mathrm{\muup}}}_{i}$ is the transition dipole moment of the $i$th TLS and ${\bm{\mathrm{\epsilonup}}}_{\bm{\mathrm{k}}}$ is the polarization direction of the photonic mode $\mathbf{k}$ . One can derive the radiative interaction $\Gamma=\Gamma_{\mathrm{Re}}+i\Gamma_{\mathrm{Im}}$ between two TLSs based on the above Hamiltonian. The real and imaginary parts describe the cooperative decay rate and cooperative energy shift, respectively. The focus of this letter will be the cooperative decay rate. Similar conclusions can be drawn for the cooperative energy shift.

We first provide a graphic illustration of why the interaction between TLSs is short-range in free space. Unlike most theoretical treatments used in the literature Cortes and Jacob (2017), we do not use the Green’s function method to describe the radiative environment. Instead, we try to keep all radiative modes in their explicit forms in order to gain a more intuitive picture. As shown in Section I of Supplementary Material (SM), the real part of the radiative interaction between TLSs can be expressed in the following form:

[TABLE]

The integral is performed on an isosurface in momentum space, i.e. all wavevectors $\mathbf{k}$ that satisfy $\omega\left(\bm{\mathrm{k}}\right)=\omega_{0}$. The integrand includes two terms. The first term is simply a polarization factor $\rho_{\mathbf{k}}=\frac{\omega_{0}}{16\varepsilon_{0}\pi^{2}v_{g}\left(\mathbf{k}\right)}\left({\bm{\mathrm{\muup}}}_{1}\cdot{\bm{\mathrm{\epsilonup}}}_{\bm{\mathrm{k}}}\right)^{\ast}{\left({\bm{\mathrm{\muup}}}_{2}\cdot{\bm{\mathrm{\epsilonup}}}_{\mathbf{k}}\right)}$, which describes the relative orientation of the transition dipole $\bm{\mathrm{\muup}}$ and the polarization of the electric field $\bm{\mathrm{\epsilonup}}$. Here $v_{g}(\mathbf{k})$ is the group velocity of mode $\mathbf{k}$. For degenerate polarization states, the integration should also include all polarizations. Since the polarization factor ${\rho}_{\mathbf{k}}$ is independent of the inter-TLS distance, it does not affect the interaction range. It is the second term, $e^{i\mathbf{k}\cdot\mathbf{R}}$, that plays the critical role in the physics of the interaction range. Here $\mathbf{R}=\mathbf{r}_{1}-\mathbf{r}_{2}$ is the distance vector between the two TLSs. The integrand $\rho_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{R}}$ is a fast oscillating function, which generally results in cancellation of the integration when the inter-TLS distance $R$ is large. Therefore, the interaction is always short-range. We can see this effect in Fig. 2a. Here we consider two TLSs in free space. The spherical isosurface has a radius of $k=|\mathbf{k}|={\omega}_{0}/c$. The real part of $\rho_{k}e^{i\mathbf{k}\cdot\mathbf{R}}$ is plotted on the isosurface. When $R=10\lambda$, there are rapid oscillations as $\mathbf{k}$ varies on the isosurface. The resulting value of the integral is small, and therefore the interaction is weak at this long distance. When the inter-TLS distance is small, for example $R=0.3\lambda$, the oscillation is slow (Figure 2c), leading to a sizeable value of the integral and thus a strong interaction. The interaction decays as the distance $R$ grows (Fig. 2d).

The graphic illustration also indicates that the interaction range is inversely proportional to the size of the isosurface in momentum space. A large inter-TLS distance $R$ on a large isosurface leads to a fast oscillating integrand on the isosurface that results in a small value of the integral. One way to counteract this effect is to substantially reduce the isosurface size. Small isosurfaces can save the integral from cancellation even for a fast-oscillating function. Figure 2c shows the real part of the integrand $\rho_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{R}}$ with a long inter-TLS distance $R=10\lambda$ on an isosurface that has a radius that is $0.03$ times that of the free-space isosurface. While the oscillation is still fast, the small isosurface cannot accommodate many oscillations, yielding a sizable value of the integral. Figure 2e shows that this strong interaction is sustained over a long distance if the isosurface is small. Specifically, for an isosurface with a radius of $q$, the real part of interaction $\Gamma_{\mathrm{Re}}$ scales as ${\sin{\left(qR\right)}}/{qR}$. As the isosurface radius approaches zero $q\to 0$, the range becomes infinite. Here, we use a polarization factor $\rho_{\mathbf{k}}$ based on plane waves, which, although a simplification, is sufficient for estimating the scaling.

The size of isosurface is fixed in free space. But there are many structured photonic environments that offer smaller isosurfaces. Here, we use Weyl photonic crystals as an example to demonstrate the inverse relationship between the interaction range and the isosurface. Weyl photonic crystals Lu et al. (2013); Yang et al. (2018) exhibit a conic dispersion relation in three-dimensional space, similar to Dirac dispersion relations in two-dimensional space. The isosurface gradually reduces to a point around the apex of the conic dispersion, i.e. the Weyl point. Observation of this small isosurface suggests that we could expect long-range interactions around isolated Weyl points. Specifically, we consider a double gyroid structure described by $g\left(\bm{\mathrm{r}}\right)\mathrm{=}{\mathrm{sin}(2\pi x/a)\ }{\mathrm{cos}(2\pi y/a)\ }+{\mathrm{sin}(2\pi y/a)\ }{\mathrm{cos}(2\pi z/a)\ }+{\mathrm{sin}(2\pi z/a)\ }{\mathrm{cos}(2\pi x/a)\ }$, where $a$ is the lattice constant. A material with a dielectric constant $\varepsilon_{r}=13$ fills the regions defined by $|g\left(\bm{\mathrm{r}}\right)|>1.1$. Four air spheres are placed in the dielectric material as defects to break parity symmetry yielding two pairs of Weyl points at identical frequencies Wang et al. (2016). The unit structure is shown Fig. 3(a). The dispersion relation on the momentum plane of $k_{z}=0$ is shown in Fig. 3b with two pairs of Weyl points at the frequency ${\omega}_{\mathrm{wp}}=0.55096[2\pi c/a]$. The isosurface becomes infinitesimally small at the Weyl point.

Using these isosurfaces, we numerically calculate the interaction between two TLSs placed inside the Weyl crystal. The photonic modes are simulated using the MPB software package Johnson and Joannopoulos (2001). The details of the calculation are shown in SM. Figure 3c shows the interaction as a function of the inter-TLS distance for three different transition frequencies, which are also labeled by white lines in Fig. 3b. The isosurfaces have four lobes because there are four Weyl points, as shown in Fig. 3c (i-iii). As the TLS transition frequency approaches the Weyl point, the isosurface size decreases, causing the interaction extends to extend to a significantly greater range. When the transition frequency is $0.00024[2\pi c/a]$ away from the Weyl point (panel iii in Fig. 3c), the interaction shows a negligible decay even at $180$ wavelengths (Fig. 3c bottom).

The decaying and oscillating patterns in these curves are attributed to a few different origins. At the largest scale, the envelop scales as $\sin{\left(\bar{q}R\right)}/\bar{q}R$, where we use $\bar{q}$ to roughly characterize the size of the isosurface (we will discuss the impact of the shape of isosurface later). The medium-range oscillation is due to the interplay of four Weyl points at the same frequency. The fastest oscillation is due to the modulation of the nonuniform field within a unit cell of the photonic crystal. The long-range interaction observed here is robust in that it does not rely on the orientation of the dipole direction or the spatial placement of TLSs (See more discussion in SM).

We can quantitatively characterize the interaction range by numerically fitting the envelope. These envelops are shown by the dashed line in Fig. 3c. We further define a range $\ell_{D}$ as the distance when the envelop drops to half of its maximum value. We calculate this range for TLSs at different transition frequencies near the Weyl points, corresponding to different isosurface sizes. The results are shown in Fig. 3d. A clear linear relationship is demonstrated between $\ell_{D}$ and the inverse of the isosurface size $1/\bar{q}$. Because the isosurfaces are not spherical, we use $\bar{q}=\sqrt{S_{\omega}/4\pi}$ to define the isosurface size, where $S_{\omega}$ is the surface area of isosurfaces.

Thus far, we have shown that the size of the isosurface plays a critical role in the interaction range. Next, we will discuss the role of the shape of the isosurface. A spherical isosurface leads to an isotropic interaction range. On the other hand, a non-spherical isosurface generally creates an anisotropic interaction range: the interaction range depends on the direction of the inter-TLS distance vector $\hat{\mathbf{R}}$. There is a general reciprocal relationship between the interaction range and the size of the isosurface when projected along $\hat{\mathbf{R}}$.

Let us take the example of an ellipsoidal isosurface in an anisotropic media. The interaction range is longer when the two TLSs are placed along the direction of the short axis of the ellipsoid $\hat{\bm{s}}$, than when they are along the long axis $\hat{\bm{l}}$. We can easily see this effect by observing the oscillation pattern of $\rho_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{R}}$ on an ellipsoidal isosurface as shown in Fig. 4a. When $\hat{\mathbf{R}}$ is parallel to the long axis $\hat{\bm{l}}$, we have many oscillations and strong cancellation of the integration. On the other hand, when $\hat{\mathbf{R}}$ is parallel to the short axis $\hat{\bm{s}}$, we have fewer oscillations and weaker cancellation.

To demonstrate this effect in Weyl photonic crystals, we plot the isosurface at frequency $\omega={\omega}_{\mathrm{wp}}+0.00404[2\pi c/a]$, where the isosurface has a flat edge-softened rectangular geometry (Fig. 4b). We plot the real part of the integrand in Eq. (2) on the isosurface for three different ${\mathbf{R}}$. Here the magnitude of $\mathbf{R}$ is fixed, but its direction $\hat{\mathbf{R}}$ varies from the short axis  $\hat{\bm{s}}$ to the long axis $\hat{\bm{l}}$. The cancellation effect is weaker when $\mathbf{R}$ is aligned with the short axis and stronger along the long axis. We also calculate the interaction as a function of the distance for the three directions shown in Fig. 4b. The range is conspicuously longer for TLSs placed along the short axis of the isosurface than that for the long axis as shown in Fig. 4c. In the case shown in Fig. 4, the frequency is greatly detuned from the Weyl point, and thus, the interaction range is not as long as those shown in Fig. 3.

The extended range of the dipole-dipole interaction extends beyond quantum systems. In the microwave regime, where Weyl photonic crystals have been experimentally realized on a printed circuit board Yang et al. (2018), the resonant dipole-dipole interaction range can also be extended. The range will also be limited by the propagation length of the waves inside such systems due to finite absorption by metallic materials.

We also emphasize that the relation between the interaction range and the isosurface is not unique to Weyl photonic crystals. It is generally applicable to periodically structured media. For example, in two-dimensional space, the scaling of the interaction range follows $J_{0}(kR)$, where $J_{0}$ is the Bessel function of the first kind. For a two-dimensional photonic crystal, a spherical isosurface with a radius of $q$ creates a different scaling law that follows $J_{0}(qR)$. More examples are discussed in Sec. II of SM.

We have discussed that the interaction range. Another important aspect is the strength of the interaction. We chose the linear dispersion near Weyl points because it makes it easy to separate the effect of the isosurface from other effects such as group velocity and density of states. However, the shrinking isosurface combined with a finite group velocity also decreases the interaction strength. At the Weyl point, the interaction strength is zero. The linear dispersion near a Weyl point results in a trade-off between the interaction range and strength. Such a trade-off can be alleviated in two-dimensional crystals and with a high order dispersion relation. We discussed the scaling of the interaction strength in Sec. II in SM.

Visual inspection of the isosurface provides a convenient tool to understand a broad class of long-range interaction phenomena. We now comment on the connection between our approach and the existing literature. The behavior of index-near-zero materials Fleury and Alu (2013) was explained by a long effective wavelength. Alternatively, it can also be conveniently explained by our method: the index-near-zero material also has an ultra-small isosurface. In addition to these examples, we can envision that Dirac points in two-dimensional photonic crystals also provide small ‘isosurfaces’ (isofrequency contours) for long-range interaction. Ref. Douglas et al. (2015) shows that inside the photonic bandgap, long tails of evanescent fields can induce long-range interaction. Here we can also see that outside the photonic bandgap but near the band edge, the propagating far fields have small isosurfaces, offering a different mechanism for long-range interaction. A hyperbolic material, where long-range interactions were allowed along specific directions, was treated using the Green’s function method Cortes and Jacob (2017). Using our graphic interpretation allows one to intuitively see that only special directions allow long-range interactions (see the visualization in SM).

To conclude, we show the deep connection between the interaction range and the isosurface in momentum space. Both the size and shape of the isosurface affect the interaction range. The method introduced here provides an intuitive understanding of underlying physics that is somewhat buried in traditional treatments, and we were able to use our method to help understand several photonic systems from the existing literature. It also provides a general recipe to search for new photonic systems that support long-range interactions.

This work was supported by the National Science Foundation (NSF) through the University of Wisconsin Materials Research Science and Engineering Center DMR-1720415. L.Y. and Z.Y. were also supported by the Defense Advanced Research Projects Agency (DARPA) (DETECT program). L. Y. also acknowledges the financial support from NSF EFRI Award-1641109.

Supplementary Material

I Theory of long-range interaction

I.1 General resonant interaction theory between two quantum two-level systems

The Hamiltonian of quantum two-level systems (TLSs) in an arbitrary photonic environment is given by

[TABLE]

They are explicitly written as Bay et al. (1997) ($\hbar=1$)

[TABLE]

where $\omega_{0}$ and $\hat{\sigma}_{i}^{\dagger}(\hat{\sigma}_{i})$ are the transition frequency and raising (lowering) operator of $i$th two-level system (TLS). $\omega_{\mathbf{k}}$ and $\hat{a}_{\mathbf{k},\eta}^{{\dagger}}(\hat{a}_{\mathbf{k},\eta})$ are the frequency and creation (annihilation) operator of photon. $g_{\mathbf{k},\eta}(\mathbf{r}_{i})=\sqrt{\omega_{\mathbf{k}}/2\varepsilon_{0}V}\bm{\mathrm{\muup}}_{i}\cdot\bm{\mathrm{\epsilonup}}_{\mathbf{k},\eta}$ is the coupling between $i$th TLS and photonic mode $\mathbf{k}$. $\bm{\mathrm{\muup}}_{i}$ is the transition dipole moment of $i$th quantum TLS and $\bm{\mathrm{\epsilonup}}_{\mathbf{k},\eta}$ is the polarization of photonic mode $\mathbf{k}$ with polarization index $\eta$.

The transition probability from initial to final states is given by the Fermi’s Golden rule $2\pi/\hbar|\mathcal{M}_{\mathrm{FI}}|^{2}\delta\left(\mathcal{E}_{\mathrm{F}}-\mathcal{E}_{\mathrm{I}}\right)$, where the transition matrix element $\mathcal{M}_{\mathrm{FI}}$ can describe the resonant dipole-dipole interaction between two TLSs. For the weak light-matter interaction, it can be written as the second-order form:

[TABLE]

Here, $|\mathrm{I}\rangle=|e_{1},g_{2};0\rangle$ and $|\mathrm{F}\rangle=|g_{1},e_{2};0\rangle$ denote initial and final states, where‘$e$’ and ‘$g$’ in the Dirac bracket notions represent excited and ground states, respectively, and the number ‘[math]’ or ‘$1$’ is the photon number in the photonic environment. The intermediate state $|\mathrm{R}_{\alpha}\rangle$ has two options: $|g_{1},g_{2};1_{\mathbf{k},\eta}\rangle$ with energy $\mathcal{E}_{\mathrm{R}_{1}}=\mathcal{E}_{g}^{(1)}+\mathcal{E}_{g}^{(2)}+\hbar\omega_{\mathbf{k}}$ and $|e_{1},e_{2};1_{\mathbf{k},\eta}\rangle$ with energy $\mathcal{E}_{\mathrm{R}_{2}}=\mathcal{E}_{e}^{(1)}+\mathcal{E}_{e}^{(2)}+\hbar\omega_{\mathbf{k}}$. The energy of the initial state is $\mathcal{E}_{\mathrm{I}}=\mathcal{E}_{g}^{(1)}+\mathcal{E}_{g}^{(2)}+\hbar\omega_{\mathbf{k}}$. Since two identical TLSs are considered, we have $\mathcal{E}_{e}^{(1,2)}-\mathcal{E}_{g}^{(1,2)}=\hbar\omega_{0}$. Then, Eq. (3) can be explicitly given by

[TABLE]

where $\mathbf{R}=\mathbf{r}_{2}-\mathbf{r}_{1}$ and $\omega_{\mathbf{k}}=\omega\left(\mathbf{k}\right)$. The summation over $\mathbf{k}$ can be written as an integral as $\sum_{\bm{\mathrm{k}}}=[{V}/{{\left(2\pi\right)}^{3}}]\iiint_{\mathcal{V}_{\mathbf{k}}}d\mathbf{k}^{3}$. The transition matrix element can also be expressed as Kurizki (1990)$\mathcal{M}_{\mathrm{FI}}=i\Gamma^{\ast}$, where $\Gamma$ is the radiative interaction. Utilizing the relation $\int_{0}^{\infty}[f(x)/\left(x-x_{0}\right)]dx=\int_{0}^{\infty}\left[\mathbb{P}\left(f(x)/\left(x-x_{0}\right)\right)+i\pi\delta(x-x_{0})f(x)\right]dx$, we can write the radiative interaction as

[TABLE]

where $\mathbb{P}$ denotes the Cauchy principal value and the polarization factor is

[TABLE]

The real part of $\Gamma$ is the cooperative decay rate and its explicit expression is

[TABLE]

where $\mathcal{S}_{\omega_{0}(\mathbf{k})}$ is the isosurface of $\omega=\omega_{0}$ in momentum space and $v_{g}\left(\mathbf{k}\right)=|\nabla_{\mathbf{k}}\omega_{\mathbf{k}}|$ is the group velocity of mode $\mathbf{k}$. The cooperative energy shift is

[TABLE]

I.2 Interaction in free space vacuum

In the free space, the dispersion relation is given by

[TABLE]

where $c$ is the speed of light. Because the isosurface is isotropic, we have $\iint_{\mathcal{S}_{\omega(\mathbf{k})}}d\mathcal{S}_{\mathbf{k}}=k^{2}\iint d\Omega_{\mathbf{k}}=k^{2}\int^{\mathrm{2}\pi}_{0}d\theta\int^{\pi}_{0}d\varphi\mathrm{sin}\varphi$ and group velocity $v_{g}\left(\mathbf{k}\right)=c$. Assuming $\bm{\mathrm{k}}\mathrm{\cdot}{\bm{\mathrm{R}}}\mathrm{=}{\xi}_{\bm{\mathrm{k}}}kR=\mathrm{cos}\varphi kR$, the cooperative decay rate in free space is written as ($c=1$)

[TABLE]

The polarization sum rule is given by

[TABLE]

with $\epsilon^{(\eta)}_{\mathbf{k},i}=\hat{\bm{\mathrm{\muup}}}_{i}\cdot\bm{\mathrm{\epsilonup}}_{\mathbf{k},\eta}$, $\delta_{12}=\hat{\bm{\mathrm{\muup}}}_{1}\cdot\hat{\bm{\mathrm{\muup}}}_{2}$, and ${\hat{k}}_{i}=\hat{\mathbf{k}}\cdot\hat{\bm{\mathrm{\muup}}}_{i}$. Then, we have

[TABLE]

where $\hat{{R}}_{1,2}=\hat{\bm{\mathrm{\muup}}}_{1,2}\cdot\hat{\mathbf{R}}$. Also, with the relation in Eq. (13), the cooperative energy shift is given by

[TABLE]

After calculating the Cauchy principal integral and integral over isosurface andrews2004virtual, we have

[TABLE]

I.3 Interaction near Weyl points

At first, we only consider a single Weyl point, as shown in Fig. 2c of the main text. The Hamiltonian for the continuum around the Weyl point is given by Lu et al. (2013)

[TABLE]

where $\sigma_{x,y,z}$ are Pauli matrices and $\mathbf{q}=\left(q_{x},\ q_{y},\ q_{z}\right)=\mathbf{k}-\mathbf{k}_{c}$ is the distance to the Weyl point in momentum space. The Weyl point is at $\mathbf{k}_{c}$ when $\mathbf{q}=0$. $v_{x,y,z}$ are the $x,y,z$ components of the group velocity. For simplicity, we assume the isosurface is isotropic, i.e. $v_{x}=v_{y}=v_{z}=v$. The dispersion relation near the Weyl point is given by

[TABLE]

Then, the cooperative decay rate is given by

[TABLE]

Here, we use a polarization based on plane waves. Although this approximation is simplified, it is sufficient for calculating the scaling. At frequencies near the Weyl point, $q\ll k_{\mathrm{c}}$ and the polarization factor term is a constant

[TABLE]

Consequently, the cooperative decay rate is

[TABLE]

The imaginary part of the radiative interaction (the cooperative energy shift) is an integral over frequencies from zero to infinity, as shown in Eq.  (10). For Weyl photonic crystals, the dispersion relation is quite different from Eq. (18) at frequencies far from the Weyl point Lu et al. (2013). Thus, we do not show an analytic estimation of the cooperative energy shift here, but numerical details will be discussed in Sec. III.

II The scaling of the interaction strength in 3D and 2D photonic environments

The interaction range increases as the isosurface decreases. However, at the same time, the density of states also decreases, particularly when the group velocity does not scale to zero. The consequence is that the interaction strength reduces while the range extends unless the group velocity scales in a way to cancel the effect. Here below, we discuss how the strength scales in different photonic environments.

II.1 Three-dimensional photonic media

In the 3D case, the derivation starts from the definition of the real part of ${\Gamma}$ as shown in Eq. (9). Here, we will show the real part of radiative interaction in different 3D photonic media.

II.1.1 3D free space

For the free space case, the dispersion relation is ${\omega}_{\bm{\mathrm{k}}}=c\left|\bm{\mathrm{k}}\right|=ck$, where $c$ is the speed of light. The radiative interaction strength is given by (see detailed in Sec. I B)

[TABLE]

where $B^{\left(3D\right)}={{\mu}_{1}{\mu}_{2}\omega}/{4\pi{\varepsilon}_{0}}$. In the long-distance regime, we have

[TABLE]

II.1.2 3D Weyl photonic environment

The linear dispersion relation near a Weyl point is given by Eq. (18). Assuming $\left(\hat{\bm{\mathrm{\muup}}}_{1}\cdot\bm{\mathrm{\epsilonup}}_{\mathbf{k},\eta}\right)^{\ast}\left(\bm{\hat{\mathrm{\muup}}}_{2}\cdot\bm{\mathrm{\epsilonup}}_{\mathbf{k},\eta}\right)\approx 1$, we have the radiative interaction (also see details in Sec. IC)

[TABLE]

As the radius of the small isosurface $q$ shrinks to zero, the interaction strength accordingly diminishes to zero. However, if the group velocity $v$ and the radius of isosurface $q$ are finite and small, the interaction can still be much larger and longer than that in vacuum (see comparison in Fig. S1(d)).

II.1.3 3D quadratic photonic environment

For a quadratic dispersion relation such as

[TABLE]

the group velocity is $v_{g}=2\beta q$. Supposing ${\bm{\mathrm{k}}}_{\bm{\mathrm{c}}}=0$, we have

[TABLE]

In the long-distance regime, the radiative interaction is written as

[TABLE]

The interaction strength reduces to zero when q = 0, which is similar to Weyl point. However, when the group velocity $v$ and the radius of isosurface $q$ are finite and small, the radiative interaction can be significantly larger than that in vacuum.

II.2 Two-dimensional photonic media

In 2D case, the real part of the radiative interaction is written as

[TABLE]

where $\ell_{\omega\left(\mathbf{k}\right)}$ is the isofrequency contour.

II.2.1 2D free space

We assume the dipole direction is normal to the 2D plane, the radiative decay rate is written as

[TABLE]

where $J_{0}\left(x\right)$ is the zero-order of the Bessel functions of the first kind and $B^{\left(2D\right)}={{\mu}_{1}{\mu}_{2}\omega}/{4{\varepsilon}_{0}}$ is a frequency-related parameter.

II.2.2 2D linear dispersion near Dirac point

The dispersion relation near a Dirac point is written as

[TABLE]

Then the group velocity is written as $v_{g}=v$. The radiative interaction is given by

[TABLE]

The 2D Dirac points provide stronger and longer-range interaction than Weyl points because of the reduced dimensionality.

II.2.3 2D quadratic dispersion near band edges.

If the 2D dispersion relation is quadratic such as

[TABLE]

and the group velocity is $v_{g}=2\beta q$. Then,

[TABLE]

In this case, the strength do not be limited by small $q$.

III The Weyl Photonic Crystal

III.1 Numerical method and details

The index of polarization $\eta$ is replaced by the index of band $n$. If the transition frequency of TLSs $\omega_{0}$ is in the $n_{0}$th band, the cooperative decay rate is given by

[TABLE]

The cooperative energy shift is written as

[TABLE]

Numerically, we use the MPB software package Johnson and Joannopoulos (2001) to calculate eigen-modes of the Weyl photonic crystal in Fig. 3 of the main text. We set the resolution in the unit cell as $30\times 30\times 30$. Then, the frequency of the Weyl points is $\omega_{\mathrm{wp}}=0.55096[2\pi c/a]$, which falls in between $4$th and $5$th bands.

III.2 Dipole orientations and spatial placements of TLSs

The details of the gyroid photonic crystal are shown in Fig. S1. To demonstrate the effect of TLS dipole orientations, we fix the direction of the distance vector as $\hat{\mathbf{R}}=\left(-1,\ 1,\ 1\right)/\sqrt{3}$. Fig. S2 shows the real part of radiative interaction $\Gamma_{\mathrm{Re}}$ as a function of distance $R$ with arbitrary dipole orientations. Although the the dipole orientations affect the oscillation patterns of $\Gamma_{\mathrm{Re}}$ curves, the envelopes of all curves show negligible decay even at 30 wavelengths. In addition, the amplitudes of $\Gamma_{\mathrm{Re}}$ with different dipole orientations stay on a same order of magnitudes.

To show the effect of distance vector, the dipole orientations are fixed as $\hat{\bm{\mathrm{\muup}}}_{1,2}=\left(0,\ 1,\ 0\right)$. Fig. S3 shows the real part of radiative interaction $\Gamma_{\mathrm{Re}}$ as a function of distance $R$ with different distance vectors. Similarly to the effect of dipole orientations, the variation of $\hat{\mathbf{R}}$ only influences the oscillation patterns of $\Gamma_{\mathrm{Re}}$ curves. However, the envelopes of all curves exhibit negligible decay at 30 wavelengths. Similary, the variation of the first TLS location $\mathbf{r}_{1}$ in a unit cell does not affect the interaction range.

IV Interactions in other photonic environments

Our theory is also applicable for understanding the dipole-dipole interactions in other photonic environments, such as near the bandedge of photonic crystals, index-near-zero materials, hyperbolic materials, etc. Here, we describe the interaction in hyperbolic media as an especially interesting case. Its isosurface is infinitely large, which normally leads to a very short interaction range. This is indeed the case. However, because of its unique shape, for the direction $\mathbf{R}$ that is normal to the isosurface, the interaction range can be very long. This can also be clearly seen in the illustration of the integrand on the surface as we will show now. The dispersion relation of hyperbolic materials can be described by

[TABLE]

Here, we choose the second type hyperbolic material with $\varepsilon_{z}=-\varepsilon_{x}=1$. The isosurface and the real part of integrand $\rho_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{R}}$ are shown in Fig. S4 (a) and (b) with $\theta=45^{\circ}$ and $80^{\circ}$, respectively. Here, $\theta$ is the included angle between $\hat{\mathbf{R}}$ and $x$-axis. When $\hat{\mathbf{R}}$ is normal to the isosurface, i.e. $\theta=45^{\circ}$,there is a thick red (positive) strip in the oscillating pattern and it results in a large integral value. If $\hat{\mathbf{R}}$ is aligned with a different direction, say $\theta=80^{\circ}$, the fast oscillation of $e^{i\mathbf{k}\cdot\mathbf{R}}$ results in cancellation of the integral and thus a weak interaction strength. In Fig. S4(c), we plot the cooperative decay rate at $R=10\lambda$ as a function of $\theta$. The left- and right-hand red dots correspond to the case in Fig. S4 (a) and (b), respectively. Our results agree with the exact numerical result in Ref. Cortes and Jacob (2017).

Furthermore, we also plot the radiative interaction $\Gamma_{\mathrm{Re}}$ as a function of distance for the hyperbolic material compared to the case in free space, as shown in Fig. S4 (d). It greatly agrees with the result in Ref. Cortes and Jacob (2017).