Purely Long-Range Coherent Interactions in Two-Dimensional Structured Baths

TL;DR

This paper investigates how quantum emitters interact coherently over long distances within a two-dimensional hexagonal bosonic bath, revealing unique relaxation and interaction behaviors near the Dirac point.

Contribution

It demonstrates the emergence of long-range coherent interactions and logarithmic relaxation laws for quantum emitters in a 2D structured bath with Dirac points.

Findings

Single emitter relaxation follows a logarithmic law in time.

Emitters exhibit inverse-distance decay of coherent interactions.

Long-range interactions occur without exponential attenuation.

Abstract

In this work we study the quantum dynamics emerging when quantum emitters exchange excitations with a two-dimensional bosonic bath with hexagonal symmetry. We show that a single quantum emitter spectrally tuned to the middle of the band relaxes following a logarithmic law in time due to the existence of a singular point with vanishing density of states, i.e., the Dirac point. Moreover, when several emitters are coupled to the bath at that frequency, long-range coherent interactions between them appear which decay inversely proportional to their distance without exponential attenuation. We analyze both the finite and infinite system situation using both perturbative and non-perturbative methods.

Figures

Equations

$$H_B=\sum _{\mathbf{k}}\left(f(\mathbf{k}){\hat{a}}_{\mathbf{k}}^{†}{\hat{b}}_{\mathbf{k}}+h.c.\right)=\sum _{\mathbf{k}}\omega (\mathbf{k})({\hat{u}}_{\mathbf{k}}^{†}{\hat{u}}_{\mathbf{k}}-{\hat{l}}_{\mathbf{k}}^{†}{\hat{l}}_{\mathbf{k}}),$$

$$H_{\mathrm{int}}=g\sum _{j,\alpha }\left(d_{{\mathbf{n}}_j,\alpha }{\sigma }_{eg}^{j_{,}\alpha }+h.c.\right)$$

$$|\Phi (t)⟩=\left[C_e(t){\sigma }_{eg}+\sum _{{\mathbf{n}}_{,}\alpha =a,b}{C}_{{\mathbf{n}}_{,}\alpha }(t)d_{{\mathbf{n}}_{,}\alpha }^{†}\right]|g⟩\otimes {|\mathrm{vac}⟩}_B$$

$$C_e(t)=-\frac{1}{2\pi i}{\int }_{-\infty }^{\infty }dE{G}_e(E+i{0}^{+})e^{-iEt},$$

$$G_e(z)=\frac{1}{z-\Delta -{\Sigma }_e(z)},$$

$${\Sigma }_e(z)=\frac{g^2}{2N^2}\sum _{\mathbf{k}}\left(\frac{1}{z-|f(\mathbf{k})|}+\frac{1}{z+|f(\mathbf{k})|}\right)=\frac{g^2}{N^2}\sum _{\mathbf{k}}\frac{z}{z^2-|f(\mathbf{k}){|}^2}.$$

$${\Sigma }_e(z)$$

$$C(z)$$

$$k(z)$$

$$K^{\mathrm{I}}(m)$$

$$K(m)={\int }_0^{\pi /2}d\theta \frac{1}{\sqrt{1-m{\sin }^2(\theta )}}.$$

$$\Sigma (E+i{0}^{+})\approx \frac{g^2}{\pi \sqrt{3}{J}^2}\left[E\log \left(\frac{E^2}{9J^2}\right)-i\pi |E|\right],$$

$$C_e(t)=R_0+\sum _{jBS}{R}_{jBS}{e}^{-i{z}_{jBS}}+\sum _{jUP}{R}_{jUP}{e}^{-i{z}_{jUP}}+\sum _{\alpha }{C}_{\alpha BC}(t),$$

$$R_{\mathrm{BS}}=\frac{1}{1-\partial_{z}\Sigma_{e}(z)}\Big{|}_{z=z_{\mathrm{BS}}}\,.$$

$$R_0={\left[\frac{1}{1-{\partial }_z{\Sigma }_e(z)}\right]}_{z=i{0}^{+}}=\frac{1}{1+\frac{g^2}{J^2}g(N)},$$

$$\underset{t\to \infty }{\lim }{C}_{\mathrm{MBC}}(t)\approx -\frac{\pi \sqrt{3}{J}^2}{2g^2}{\int }_0^{\infty }dy{e}^{-yt}\frac{1}{y\log {\left(\frac{y}{3J}\right)}^2}.$$

$$\underset{t\to \infty }{\lim }{\int }_0^{\infty }dy{e}^{-yt}{y}^{\alpha -1}(-\log (y){)}^{\beta }=$$

$$\displaystyle=\frac{1}{t^{\alpha}}\sum_{k=0}^{\infty}(-1)^{k}\binom{\beta}{k}\Gamma^{(k)}(\alpha)\big{(}\log(t))^{\beta-k}\,.$$

$$\frac{d}{dt}\left[{\int }_0^{\infty }dy{e}^{-yt}{y}^{-1}(\log (y){)}^{-2}\right]=-{\int }_0^{\infty }dy{e}^{-yt}(\log (y){)}^{-2},$$

$$\underset{t->\infty }{\lim }{C}_{\mathrm{MBC}}(t)\approx -\frac{\pi \sqrt{3}{J}^2}{2g^2\log (3Jt)}.$$

$$H_{\mathrm{int}}$$

$$+\frac{g}{N\sqrt{2}}\sum _j\sum _{\mathbf{k}}\left({\hat{l}}_{\mathbf{k}}\left(e^{-i\mathbf{k}\cdot {\mathbf{n}}_j}{\sigma }_{eg}^{j,A}-e^{-i(\mathbf{k}\cdot {\mathbf{n}}_j+\varphi (\mathbf{k}))}{\sigma }_{eg}^{j,B}\right)+h.c.\right).$$

$$H_{\mathrm{int}}=\frac{g}{\sqrt{2}N}\sum _j\sum _{\mathbf{k}}\left[e^{-i\mathbf{k}\cdot {\mathbf{n}}_j}{\sigma }_{eg}^j\left({\hat{u}}_{\mathbf{k}}+{\hat{l}}_{\mathbf{k}}\right)+h.c.\right],$$

$$H_{\mathrm{int}}$$

$$\displaystyle\left(\hat{u}_{\mathbf{k}}-\hat{l}_{\mathbf{k}}\right)e^{-i(\mathbf{k}\cdot\mathbf{n}_{2}+\phi(\mathbf{k}))}\sigma_{eg}^{2}+\mathrm{h.c.}\Big{]}\,.$$

$$H_{\mathrm{int}}^{\mathrm{AA}}$$

$${\hat{u}}_{\mathbf{k},\pm }$$

$$\displaystyle\left(e^{i\mathbf{k}\cdot\mathbf{n}_{1}}\pm e^{i\mathbf{k}\cdot\mathbf{n}_{2}}\right)\hat{u}_{-\mathbf{k}}\Big{]}\,,$$

$${\hat{l}}_{\mathbf{k},\pm }$$

$$\displaystyle\left(e^{i\mathbf{k}\cdot\mathbf{n}_{1}}\pm e^{i\mathbf{k}\cdot\mathbf{n}_{2}}\right)\hat{l}_{-\mathbf{k}}\Big{]}\,.$$

$$H_{\mathrm{int}}^{\mathrm{AB}}$$

$$\displaystyle\sqrt{1\mp\cos(\mathbf{k}\cdot\mathbf{n}_{12}+\phi(\mathbf{k}))}\hat{l}_{\mathbf{k},\pm}\sigma_{\pm}^{\dagger}+\mathrm{h.c.}\Big{]}\,,$$

$${\hat{u}}_{\mathbf{k},\pm }=\frac{1}{2\sqrt{1\pm \cos (\mathbf{k}\cdot {\mathbf{n}}_{12}+\varphi (\mathbf{k}))}}$$

$$\displaystyle\left(e^{i\mathbf{k}\cdot\mathbf{n}_{1}}\pm e^{i(\mathbf{k}\cdot\mathbf{n}_{2}+\phi(\mathbf{k}))}\right)\hat{u}_{-\mathbf{k}}\Big{]}\,,$$

Full text

I Introduction

Recently, the interest on the dynamics of quantum emitters (QEs) coupled to a structured baths has revived due to the prospects offered by new experimental platforms Vetsch et al. (2010); Huck et al. (2011); Hausmann et al. (2012); Laucht et al. (2012); Thompson et al. (2013); Goban et al. (2014); Béguin et al. (2014); Lodahl et al. (2015); Bermúdez-Ureña et al. (2015); Sipahigil et al. (2016); Corzo et al. (2016); Sørensen et al. (2016); Sipahigil et al. (2016); Solano et al. (2017); de Vega et al. (2008); Navarrete-Benlloch et al. (2011); Houck et al. (2012); Astafiev et al. (2010); Hoi et al. (2013); Van Loo et al. (2013); Liu and Houck (2017). Of particular interest are the way in which the dynamics of the QE is modified by the bath (see Purcell et al. (1946); Nakazato et al. (1996); Lambropoulos et al. (2000); Woldeyohannes and John (2003); Giraldi and Petruccione (2011); Breuer et al. (2016); de Vega and Alonso (2017) and references therein), the way the latter can mediate interactions or collective dissipation Le Kien et al. (2005); Dzsotjan et al. (2010); Gonzalez-Tudela et al. (2011); Chang et al. (2012); Shahmoon and Kurizki (2013); Douglas et al. (2015); González-Tudela et al. (2015a); Shahmoon et al. (2016) or how the QEs can induce entanglement in the bath González-Tudela et al. (2015b, 2016); Douglas et al. (2016).

For a single QE, unconventional relaxations dynamics have been predicted when the QE energies lie close to band-edges John and Quang (1994) or close to 2D Van-Hove singularities in the middle of the band González-Tudela and Cirac (2017a, b). The former is associated to the emergence of atom-photon bound states Bykov (1975); John and Wang (1990); Kurizki (1990), which leads to fractional decay, whereas the latter is also accompanied by a highly directional emission into the bath González-Tudela and Cirac (2017a, b); Galve et al. (2017); Mekis et al. (1999); Langley (1996), and very slow relaxation dynamics scaling with $\sim 1/(t^{2}\log(t)^{2})$.

For two (or more) QEs, the possibility of obtaining long-range coherent interactions between quantum emitters can be harnessed, e.g., to generate long-distance entanglement Shahmoon and Kurizki (2013), but also to simulate spin-models with long-range interactions. Non-conventional phenomena occur when the interactions decay with the distance, $r$, as $1/r^{a}$, with $a<D$ and $D$ being the spatial dimension on the bath. For instance, non-local transmission of correlations Hauke and Tagliacozzo (2013); Richerme et al. (2014); Maghrebi et al. (2016), violation of the area law Koffel et al. (2012), fast equilibration Kastner (2011) and quantum state transfer Vodola et al. (2014); Eldredge et al. (2016).

The coupling to the three-dimensional structureless photonic bath ($D=3$) induces coherent interactions decaying with the distance as a power law Lehmberg (1970a, b). However, they are unavoidably accompanied by dissipative terms of the same order which compete with them. A way of avoiding such dissipative terms consists in using structured reservoirs, such as photonic band-gap materials Joannopoulos et al. (1995), and exploit the emergence of atom-photon bound states Bykov (1975); John and Wang (1990); Kurizki (1990) to mediate the interactions Douglas et al. (2015, 2016). The price to pay in this case for the cancellation of the dissipative terms, however, is that the interactions are ultimately exponentially attenuated by the length of the bound states. Even though this length can be large, it is of fundamental and practical interest to know whether there is a bath that can combine the best of both scenarios, that is, having no dissipative terms as photonic band-gap materials, while keeping unattenuated interactions like 3D structureless baths. Inspired by pioneering works on graphene (see Pereira et al. (2006); Castro Neto et al. (2009) and references therein), in this manuscript we provide a positive answer to this question. In particular, we show how by coupling the QEs to a bosonic bath with a Dirac cone dispersion relation gives rise to such scenario.

In the context of graphene research Castro Neto et al. (2009), transport of energies around the Dirac cone in the presence of classical scatterers has been thoroughly analyzed. In that problem, a localized state emerging around the scatterers plays a crucial role. Its nature is very special given that its wavefunction decays very slowly with the distance, $r$, to the scatterer [$\sim 1/r$]. This function is (marginally) not square-integrable in 2D, and thus it is labeled as quasi-bound state (qBS) Pereira et al. (2006). As we will show here, such state also has crucial consequences when instead of scatterers we have QEs, and instead of a fermionic bath, we replace it by a bosonic one.

In this manuscript, we solve exactly the quantum dynamics of initially excited QEs with their transition frequency exactly at the Dirac point. We compare the results with those predicted by perturbative methods based on Born-Markov (or single-pole approximations) Lehmberg (1970a, b) that are normally used in quantum optics. For a single QE, Markov approximation predicts no decay because the density of states is zero, whereas the exact dynamics predicts a logarithmic relaxation to the ground state in the thermodynamic limit; that is, when the number of modes, $N$, in the bath diverges. This very slow relaxation quenches at a certain time for finite systems, leading to a fractional decay of excitations. For two QEs coupled to different sublattices, Born-Markov approximation predicts purely long-range coherent interactions with $a=1$ and without exponential attenuation in the thermodynamic limit, whereas the exact treatment predicts no interactions at all. We characterize the convergence to the thermodynamic limit showing that the exact dipole-dipole interactions decay with a $1/\log(N)$ dependence. Thus, for practical purposes these interactions can still be harnessed in large systems as we show numerically. We provide both analytical and numerical understanding of the phenomena and the interpolation between finite and infinite systems, which allows us to evaluate to which extent these interactions can be observed.

The manuscript is structured as follows: in Section II we describe the system under study by writing explicitly all the terms of the Hamiltonian and explain the bath properties. Then, in Section III, we characterize the dynamics of a single QE in both finite and infinite systems. In Section IV we study the physics emerging when two QEs are coupled to the bath, generalizing the conclusions to many QEs. In Section V we give a brief overview of several experimental limitations, and how they will affect to the observation of the phenomen. Finally, we conclude the paper in Section VI by summarizing the main results and pointing to future directions of work.

II System

For the bath, we use a two-band model formed by two different lattices, that we denote as A/B lattices. The annihilation operators of the A/B lattice modes are described in terms of bosonic operators, $a_{\mathbf{n}}$ and $b_{\mathbf{n}}$, respectively. We assume both lattices to be degenerate in energy $\omega_{a}=\omega_{b}$. The A/B lattices are connected through nearest-neighbor coupling $J$, such that the bath Hamiltonian reads as follows (using $\hbar=1$): $H_{B}=J\sum_{\mathbf{n}}\sum_{\mathbf{m}=0,\mathbf{a}_{1},\mathbf{a}_{2}}\left(a^{\dagger}_{\mathbf{n}+\mathbf{m}}b_{\mathbf{n}}+\mathrm{h.c.}\right)$. The positions within the A/B lattices are given by two integer numbers $n_{1,2}=1,\dots,N$, such that $\mathbf{n}=(n_{1},n_{2})\equiv n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2}$, where $\mathbf{a}_{1,2}$ are the primitive lattice vectors, e.g., for hexagonal lattices $\mathbf{a}_{1,2}=(\frac{3}{2},\pm\frac{\sqrt{3}}{2})$ [and $0=(0,0)$]. We have written $H_{B}$ in a frame rotating with $\omega_{a}$ that we use for the rest of the Hamiltonians of the manuscript. This bath is the bosonic analogue of graphene, which can be implemented with different systems such as cold atoms Polini et al. (2013) or circuit QED Houck et al. (2012). Despite the simplicity of the model, we also expect it to capture the most important features of more complex photonic realizations where Dirac cone dispersions emerge Bravo-Abad et al. (2012); Xie et al. (2014).

By imposing periodic boundary conditions, we define the operators $\hat{c}_{\mathbf{k}}=\frac{1}{N}\sum_{\mathbf{n}}c_{\mathbf{n}}e^{-i\mathbf{k}\cdot\mathbf{n}}$ for $c=a,b$ bath modes, where $\mathbf{k}=(k_{1},k_{2})\equiv k_{1}\mathbf{b}_{1}+k_{2}\mathbf{b}_{2}$, with $\mathbf{b}_{1,2}$ satisfying $\mathbf{a}_{i}\cdot\mathbf{b}_{j}=\delta_{ij}$ and $k_{1,2}=\frac{2\pi}{N}(-\frac{N}{2},\dots,\frac{N}{2}-1)$. Using those operators in momentum space, the bath Hamiltonian reads:

[TABLE]

where the eigenoperators $\hat{u}_{\mathbf{k}}(\hat{l}_{\mathbf{k}})=\frac{1}{\sqrt{2}}\left(\hat{a}_{\mathbf{k}}+(-)\hat{b}_{\mathbf{k}}e^{i\phi(\mathbf{k})}\right)$ represent the annihilation operators of upper(lower) band modes, respectively, and $f(\mathbf{k})=J(1+e^{ik_{1}}+e^{ik_{2}})=\omega(\mathbf{k})e^{i\phi(\mathbf{k})}$ 111Notice, we use a different convention than the one used in traditional graphene studies Castro Neto et al. (2009).. The dispersion relation extends from $[-3J,3J]$ and has two bands that touch at $\mathbf{K}_{\pm}=\frac{2\pi}{3}(\pm 1,\mp 1)$, the so-called Dirac points corresponding to $\omega(\mathbf{K}_{\pm})\equiv 0$ [see Fig 1(a)]. At these points, the energy dispersion can be linearized, i.e., $f(\mathbf{K}_{\pm}+\mathbf{q})\approx J\mathbf{h}_{\pm}\cdot\mathbf{q}$, with $\mathbf{h}_{\pm}=i\left(e^{\frac{\pm 2\pi i}{3}},e^{\frac{\mp 2\pi i}{3}}\right)$.

The density of states of this bath in the limit $N\rightarrow\infty$ is plotted in Fig 1(b). Apart from discontinuities and divergences also appearing in other structured baths, it possesses a singular point at the Dirac point, i.e., $E=0$. This is the region that we focus along this manuscript as, up our knowledge, has no analogue in other quantum optical scenarios, and it is the source of many interesting behaviour in other contexts Pereira et al. (2006); Castro Neto et al. (2009).

We are interested in the dynamics of $N_{e}$ QEs described by two-level systems, $\{\left|g\right\rangle_{j,\alpha},\left|e\right\rangle_{j,\alpha}\}$, whose free Hamiltonian reads: $H_{S}=\Delta\sum_{j,\alpha}\sigma_{ee}^{j,\alpha}$, where $\Delta=\omega_{e}-\omega_{a}$ is the detuning between the transition frequency of the QEs, $\omega_{e}$, and the reference energy of the bath modes, $\omega_{a}$. We also defined the spin operator of the $(j,\alpha)$-th QE as $\sigma_{\alpha\beta}^{j,\alpha}=\left|\alpha\right\rangle_{j,\alpha}\left\langle\beta\right|$. The indices $(j,\alpha)$ denote both the position, $\mathbf{n}_{j}$, and the sublattice $\alpha=A,B$, of the bath mode that the QE is coupled to. Thus, the interaction Hamiltonian generally reads:

[TABLE]

with $d_{j,A/B}=a_{\mathbf{j}}/b_{\mathbf{j}}$ and $g$ is the coupling strength of the interaction. We study the dynamics of the QEs when they are initialized on a given state $\left|\Phi_{0}\right\rangle_{S}$ containing a single excitation and the bath is in the vacuum, i.e., $\left|\mathrm{vac}\right\rangle_{B}=\left|0\right\rangle_{A}^{\otimes N^{2}}\otimes\left|0\right\rangle_{B}^{\otimes N^{2}}$. Then, we study the evolution of the combined system under the total Hamiltonian $H=H_{S}+H_{B}+H_{\mathrm{int}}$.

III Single quantum emitter dynamics

It is instructive to first study the single excited QE, i.e., $\left|\Phi_{0}\right\rangle_{S}=\left|e\right\rangle$, as the structure of the bath has already remarkable consequences on its dynamics. We assume the QE to be coupled to the mode $a_{\mathbf{n}_{0}}$ (the coupling to the B lattice leads to similar behaviour). As the total Hamiltonian $H$ conserves the number of excitations, the wavefunction at any time $t>0$ can be written:

[TABLE]

with initial condition $C_{e}(0)=1$ and $C_{\mathbf{n},\alpha}(0)=0$ . In Fig. 2 we show the result of a numerical simulation of the dynamics with $g=0.1J$ and for several illustrative $\Delta$’s, together with the bath population in the A/B lattices at $tJ=200$. For $\Delta$’s close to the band edges ($\Delta=2.5J$) the emission is mostly isotropic and equally distributed among the A and B sublattices, as shown in Figs. 2(b-c). The isotropic emission can be traced back to the isotropic character of $\omega(\mathbf{k})$ around the band edges. For $\Delta=J$ the emission is highly anisotropic Mekis et al. (1999); Langley (1996); Galve et al. (2017); González-Tudela and Cirac (2017a, b), predominantly in three directions in both the A/B lattices, as shown in Fig. 2(d-e). This is also accompanied by overdamped oscillations and power law decay dynamics [see inset of Fig. 2(a)]. These non-perturbative dynamics associated to divergences of $D(E)$ were explored in Ref. González-Tudela and Cirac (2017a, b). Finally, for regions close to the Dirac point, the dynamics get slower because of the reduction of the density of modes until it seems to get quenched at $\Delta=0$. This behaviour appears together with an asymmetry in the population of the A/B lattices, as shown in Fig. 2(f-g), where the B lattice shows a larger bath population around the impurity than the A sublattice that the QE is coupled to. This will have implications when more than one QE are coupled to the lattice, as we show below.

To gain more intuition about the dynamics at $\Delta=0$, we plot in Fig. 3 the excited state population, $|C_{e}(t)|^{2}$, for a larger timescale and different ratios $g/J$. We observe that the dynamics at this point have two unexpected behaviours: i) for long times the emission appears to be oscillating around a constant value which is smaller the larger $g/J$; ii) For short times [Fig. 3(b)], $|C_{e}(t)|^{2}$ shows a logarithmic relaxation until a time $t_{0}$ related to the finite size, $N$, of the numerical simulation.

To explain this unconventional behaviour, we use the resolvent operator technique Cohen-Tannoudji et al. (1992), in which the probability amplitude $C_{e}(t)$ is calculated as the Fourier transform:

[TABLE]

of the single QE Green function, $G_{e}(z)$. The latter is given by:

[TABLE]

where $\Sigma_{e}(E)$ is the so-called self-energy which contains the effect that the bath produces in the QE. In contrast to standard single band models, the self-energy contains the contribution of the two baths (upper/lower band), leading to:

[TABLE]

This function can be calculated analytically in the limit $N\rightarrow\infty$, in terms of elliptic integrals Horiguchi (1972), obtaining

[TABLE]

for $z\notin[-3,3]J$, and where $K^{\mathrm{I}}(m)$ has to be defined in a piecewise manner to guarantee analyticity in its definition domain:

[TABLE]

being $K(m)$ the complete elliptic integral of the first kind defined by:

[TABLE]

Evaluating $\Sigma_{e}(z)$ slightly above the real axis, $\Sigma_{e}(E+i0^{+})=\delta\omega_{e}(E)-i\frac{\Gamma_{e}(E)}{2}$, we obtain the functions $\delta\omega_{e}(E)$ and $\Gamma_{e}(E)$, that we plot in Fig. 1(b). There, we observe that the function $\Gamma_{e}(E)$ [$\delta\omega_{e}(E)$] has a discontinuity [divergence] at $E=\pm 3J$ appearing at the border of the band edges of the upper/lower band, $\pm\omega(\mathbf{k})$. Moreover, at $E=\pm J$ the function also displays a divergence [discontinuity] at $E=\pm J$. Close to the Dirac point, for $|E|\ll J$, the self-energy can be expanded as:

[TABLE]

such that both $\Gamma_{e}(E)$ and $\delta\omega_{e}(E)$ are continuous functions at $E=0$, but have discontinuous derivatives.

Perturbative approaches simply replace $\Sigma_{e}(E+i0^{+})\approx\Sigma_{e}(\Delta+i0^{+})$ in Eq. 4, and assumes that the dynamics of $C_{e}(t)$ are just given by the contribution of the pole at $z_{M}=\Delta+\Sigma_{e}(\Delta+i0^{+})=\Delta+\delta\omega_{M}-i\frac{\Gamma_{M}}{2}$. Thus, the population dynamics within that approximation are given by $|C_{e}(t)|^{2}\approx e^{-\Gamma_{M}t}$, with $\Gamma_{M}$ being the one predicted by Fermi’s Golden Rule. As $z_{M}=0$ at $\Delta=0$, this perturbative approach predicts no decay, i.e., $C_{e}(t)=1$, and fails therefore to describe both the logarithmic relaxation and fractional decay that we observe for finite systems in Fig. 3.

To obtain the exact dynamics of the probability amplitude $C_{e}(t)$ one can calculate the Fourier Transform of Eq. 4 by closing the contour of integration in the lower plane of the complex plane ($\mathrm{Im}[z]<0$). However, when the $\Sigma_{e}(z)$ have regions where it is non-analytic, as the one of Eq. 7, one must be careful when closing the contour to avoid them. One possible choice for the contour of integration consists of taking 5 detours of the integration contour at $E=\pm 3J,\pm J$ and [math] as schematically depicted in Fig. 4(a). In the regions II-V, we must analytically continue the function $\Sigma_{e}(z)$ into other Riemann sheets (see Appendix for details) to ensure the continuity of $\Sigma_{e}(z)$ along the integration contour. Thus, the total dynamics can be shown to be given by a sum of different contributions:

[TABLE]

where we have

• •

The contribution from real bound states (BS) of the single QE Green function $G_{e}(E)$ that appear outside the continuum. Their energies are given by the solution $z_{\mathrm{BS}}-\Delta-\Sigma_{e}(z_{\mathrm{BS}})=0$, and which corresponding Residue is obtained by using Residue Theorem yielding:

[TABLE]

• •

Apart from the real BS, we can also find complex or unstable poles (UPs) when $\mathrm{Re}[z]\in[-3J,3J]$ appearing in the analytic continuation of the Green Function to other Riemann sheets. Since our Green function is extended into different sheets it is possible to find more than one UP for the same value of $\Delta$. These UPs are obtained from the same pole equation than the real ones, but replacing $\Sigma_{e}$ by the corresponding analytic continuation of the region. The corresponding residue is also calculated in an analogue way to the BS case.

• •

We must also take into account the contributions from the detour because of the BCs. Notice, that in each side of the detour one must use the corresponding $\Sigma_{e}$ depending on the region where the integral appears.

• •

Finally, we have separated the contribution corresponding to the quasi-bound state (qBS) at $z=0$, with residue $R_{0}$, as it will play an important role in the discussion of finite systems. The residue $R_{0}$ can be calculated as follows:

[TABLE]

where $g(N)=\frac{J^{2}}{N^{2}}\sum_{\mathbf{k}}\frac{1}{|f(\mathbf{k})|^{2}}\sim\frac{2}{\pi\sqrt{3}}\log(N)$ for $N\gg 1$ as we show in the Appendix. This residue can be understood as the overlap of the initial wavefunction with the qBS emerging at $\Delta=0$ Pereira et al. (2006), $\left|\Phi_{\mathrm{qBS}}\right\rangle$, i.e., $R_{0}=\langle\Phi_{0}|\Phi_{\mathrm{qBS}}\rangle$ and is ultimately responsible of the excitation remaining in the QE as $C_{e}(\infty)=R_{0}$. For a finite system, the qBS is indeed square integrable which is why $R_{0}$ can be finite. However, the very slow decay ($1/r$) with the distance makes it not square integrable in the thermodynamic limit, which is why $R_{0}\rightarrow 0$ when $N\rightarrow\infty$.

• •

The rest of the contours can be shown to give no contribution.

To illustrate the weight in the thermodynamic limit of the different processes to the dynamics depending on $\Delta$, we plot the absolute value of each contribution at $t=0$ for fixed $g=0.6J$ and different $\Delta$’s in Fig. 4(b). For the sake of simplicity, we plot the sum of all the branch cut contributions using the same color, although we calculate them independently. We observe several interesting regimes:

• •

When $\Delta$ is tuned very far from the band, the main contributions is either from the LBS or UBS. This predicts that the QE will not decay as there are no modes energetically resonant with the QE. As expected from the divergence of $\delta\omega_{e}(E)$ at $E=\pm 3J$ Shi et al. (2016), both the UBS/LBS survive for all $\Delta$, including those which lie within the band.

• •

As the QE transition frequency gets closer to the band edges, $|\Delta\pm 3J|\ll g,J$, the contribution of the BS decreases and the BC contribution becomes more important, until the weight moves to the UPs.

• •

As we already showed for a square lattice González-Tudela and Cirac (2017a, b), the divergence of the imaginary part of $\Sigma_{e}(E)$, in this case at $E=\pm J$, leads to the coexistence of two UPs for a range of $\Delta$’s given by the jump of $\delta\omega_{e}(E)$ around that value. Moreover, the imaginary part of the pole saturates to a finite value, as shown in Fig. 4(c), compared to the infinite value predicted by Fermi’s Golden Rule. As we already explored extensively this phenomenon in Refs. González-Tudela and Cirac (2017a, b), we will not investigate it further in this manuscript.

• •

Finally, around the Dirac point we observe that as $\Delta$ gets closer to [math] the BC contribution increases, making a discontinuous jump at $\Delta=0$, where the contribution of the left/right UPs vanishes. The middle BC contribution can be calculated by using the analytical expansions of the elliptic functions around $E\approx 0$ which yield:

[TABLE]

The asymptotic expansion of Laplace transform of functions with logarithmic singularities was studied in Ref. Wong and Lin (1978), where it was shown that:

[TABLE]

where $\mathrm{Re}\alpha$ must be $>0$. Unfortunately, our $C_{\mathrm{MBC}}(t)$ lies at the border of validity of this formula as $\alpha=0$. However, one can realize that by taking the first derivative of the integral once that we obtain a formula that lies within the region of validity:

[TABLE]

which has the shape of Eq. • ‣ III, but with $\alpha=1$. Applying now the formula, the leading term of the derivative at long times is $\sim 1/(t\log(t)^{2})$. Now, integrating back this leading term one can find:

[TABLE]

such that $|C_{\mathrm{MBC}}|^{2}\propto 1/(\log(t))^{2}$. We also checked numerically that the convergence to the leading term is slow such that intermediate times one may find other scalings for the relaxation.

Summarizing, for a single QE we have shown that in the vicinity of the Dirac point, perturbative methods fail to capture the exact behaviour of both finite and infinite systems. In finite systems, it does not capture the existence of fractional decay of excitations ($R_{0}<1$), whereas in infinite ones it completely neglects the logarithmic relaxation which dominates the dynamics.

IV Two quantum emitters dynamics

Now we will study the consequences of such exotic behaviour when two QEs are coupled to the bath. In particular, we aim at discerning whether it is possible to obtain long-range dipole-dipole interactions by exploiting the existence of the qBS. To explore that situation, we assume that two QEs initially start in $\left|\Psi_{0}\right\rangle_{S}=\left|e\right\rangle_{1}\otimes\left|g\right\rangle_{2}$ and study the populations of the QEs, $|C_{1,2}(t)|^{2}$, to see whether there are coherent oscillations or not.

IV.1 Rewriting interaction Hamiltonian for two QEs

Before moving on to the solution of the dynamics $|C_{1,2}(t)|^{2}$ it is instructive for this Section, to rewrite $H_{\mathrm{int}}$ in terms of the eigenmodes of the $H_{B}$ operator, which reads:

[TABLE]

Let us particularize for the two situations of interest of this Section, namely,

• •

For two QEs at positions $\mathbf{n}_{1,2}$ coupled within the same sublattice, e.g., A, the interaction Hamiltonian reads:

[TABLE]

where we dropped the index $\alpha$ denoting the sublattice the QE is coupling to. Apart from the coupling of the QEs to two independent baths (upper/lower modes), the only difference between the coupling of the two QEs to each bath is the phase difference introduced by $e^{i\mathbf{k}\cdot\mathbf{n}_{12}}$, where we use the notation $\mathbf{n}_{12}=\mathbf{n}_{2}-\mathbf{n}_{1}$ for the vector connecting the two impurities. This is analogue to what happens with standard open quantum systems.

• •

Different from the previous situation, when two QEs at positions $\mathbf{n}_{1,2}$ are coupled to different sublattices, e.g., the $\mathbf{n}_{1/2}$ one to be coupled the A/B lattice respectively. In that case, the interaction Hamiltonian is given:

[TABLE]

which apart from the propagation phase ($e^{i\mathbf{k}\cdot\mathbf{n}_{12}}$) contains an extra $\mathbf{k}$-dependent phase, $e^{i\phi(\mathbf{k})}$, that has consequences on the dynamics.

It is also interesting to consider the bath modes that are coupled to the QE symmetric/antisymmetric operators $\sigma^{\dagger}_{\pm}=\frac{1}{\sqrt{2}}\left(\sigma_{eg}^{1}\pm\sigma_{eg}^{2}\right)$ and $\sigma_{\pm}=\left(\sigma_{\pm}^{\dagger}\right)^{\dagger}$. With that definition, it is possible to rewrite $H_{\mathrm{int}}$ in such a way that $\sigma^{\dagger}_{\pm}$ couple to two orthogonal bath modes for the upper/lower band operators. This allows us to solve the symmetric/antisymmetric dynamics independently. In particular, for the case of two QEs coupled to the AA lattice, $H_{\mathrm{int}}$ is rewritten as:

[TABLE]

where $u_{\mathbf{k},\pm},l_{\mathbf{k},\pm}$ are two orthogonal bath modes defined as:

[TABLE]

For the case of two QEs, $\mathbf{n}_{1}$ and $\mathbf{n}_{2}$, coupled to the A/B lattice respectively, the interaction Hamiltonian, $H_{\mathrm{int}}$, changes to:

[TABLE]

where $\hat{u}_{\mathbf{k},\pm},\hat{l}_{\mathbf{k},\pm}$ are two orthogonal bath modes defined as:

[TABLE]

The important consequence of this separation is that for the case of two QEs we will be able to work separately with the symmetric & antisymmetric subspaces. This means that the Green functions $G_{1,2}(z)$ associated to the probability amplitudes $C_{1,2}(t)$ can be obtained from Green Functions $G_{\pm}(z)$ associated to the symmetric/antisymmetric combination of excitations as follows: $G_{1,2}(z)=\left[G_{+}(z)\pm G_{-}(z)\right]/2$. Remarkably, $G_{\pm}(z)$ are given by the same expression than the single QE, $G_{e}(z)$, just replacing $\Sigma_{e}(z)\rightarrow\Sigma^{\beta}_{\pm}(z;\mathbf{n}_{12})=\Sigma_{e}(z)\pm\Sigma^{\beta}_{12}(z;\mathbf{n}_{12})$. The collective self-energy, $\Sigma^{\beta}_{12}(z;\mathbf{n}_{12})$, reads:

[TABLE]

where $\beta=AA,AB,BA,BB$ is an index that denotes to which sublattices the QEs are coupled to. Depending on $\beta$, we have that: $D_{AA}=D_{BB}=z$, whereas $D_{AB}=D_{BA}^{*}=f^{*}(\mathbf{k})$.

Perturbative approaches predict a dynamics in the symmetric and antisymmetric subspace $C_{\pm}(t)\approx e^{-i(J_{M,\pm}-i\frac{\Gamma_{M,\pm}}{2})t}$, where we define $J_{M,\pm}-i\frac{\Gamma_{M,\pm}}{2}=\Sigma_{\pm}(\Delta+i0^{+})$. Thus, in the original basis they turn into:

[TABLE]

Therefore, purely coherent oscillations of $|C_{1,2}(t)|^{2}$ are associated to the existence of two real bound states[$\Gamma_{M,\pm}=0$], $\left|\Psi_{\pm}\right\rangle$, of the symmetric/antisymmetric component with different real component, $J_{M,+}-J_{M,-}\neq 0$.

IV.2 Two QEs coupled to the AA/BB lattices

At this point, it is easy to show that $\Sigma_{\pm}^{\mathrm{AA/BB}}(0)=0$, such that $z=0$ satisfies the pole equation for both the symmetric and antisymmetric wavefunction. Within a perturbative treatment, this will imply that these states are subradiant and should not relax to the ground state. However, as it occurred for a single QE when doing the exact calculation one finds different result. In particular, it can be shown that the residue associated to $z=0$ reads:

[TABLE]

where:

[TABLE]

Reminding that the divergence with $N$ occurs because $f(\mathbf{K}_{\pm})=0$ at the $\mathbf{K}_{\pm}$ points, we see that there are two interesting regimes depending on the relative position between the QEs, $\mathbf{n}_{12}=(n_{1},n_{2})$.

• •

When $\mathbf{n}_{12}$ is such that $e^{i\mathbf{K}_{\pm}\cdot\mathbf{n}_{12}}\neq 1$, then the numerator is always finite at the $\mathbf{K}_{\pm}$ points, yielding $g_{\pm}(\mathbf{n}_{12},N)\propto\log(N)$. The dynamics in the symmetric/antisymmetric subspaces is then very similar to the single QE situation, that is, an initial logarithmic decay for short times, until it quenches for a finite system or completely decays for an infinite one.

• •

On the other hand when $e^{i\mathbf{K}_{\pm}\cdot\mathbf{n}_{12}}=1$, which occurs when $n_{1}-n_{2}=3m$, with $m\in\mathbb{Z}$, then $g_{-}(\mathbf{n}_{12},N)\rightarrow C(\mathbf{n}_{12})$ when $N\rightarrow\infty$. This occurs because the divergence of the denominator cancels with the zero of the numerator. This means that for such positions, the antisymmetric component is indeed a real pole of the Green Function in the continuum limit, that is, a perfect subradiant state for all parameter regimes. As it occurs with subradiant states in other reservoirs, the only effect of the distance is to decrease the residue of the pole due to retardation effects. We can estimate these retardation effects, for example, for the situation $n_{1}=n_{2}=n$ by assuming that most of the contribution is coming from the $\mathbf{k}$ modes around $\mathbf{K}_{\pm}$, where the divergence of the denominator occurs. Then, we can displace the integral to $\mathbf{k}=\mathbf{K}_{\pm}+\mathbf{q}$ to obtain:

[TABLE]

where the $\mathbf{p}$-coordinates are defined by $q_{1,2}=\frac{3}{2}\left(p_{1}\pm\frac{1}{\sqrt{3}}p_{2}\right)$ and they lead to an isotropic $\omega(\mathbf{p})$. In the last approximation of the integral we expanded for $n\gg 1$. The constant $D$ depends on the momentum cut-off $q_{c}$ and we numerically calculate it using the exact sum of Eq. 32 for $\mathbf{n}_{12}=(1,1)$ obtaining $D\approx 0.6$. Thus, the final residue of these subradiant states is given by:

[TABLE]

Notice, the very slow decay of the residue of such subradiant states with $n$, i.e., $R_{\mathrm{sb}}\propto 1/\log(n)$, which is considerably slower than the one of the subradiant states for the square lattice González-Tudela and Cirac (2017a, b) which decay with $1/n^{2}$ or even slower than the 1D situation $\propto 1/n$.

Even though these subradiant states are interesting on its own, their associated real part, which is the one responsible for coherent interactions, is also zero. As in this manuscript we are interested in obtaining coherent interactions, we will not discuss them further, and leave for future work a more detailed study about them.

IV.3 Two QEs coupled to different sublattices.

The situation with two QEs coupled to different sublattices is radically different from the single QE or when the two QEs are coupled to the same sublattice. To illustrate it, we plot in Fig. 5(a), the numerical evolution of $C_{1,2}(t)$ for two QEs coupled to the A/B lattice respectively, with $g=0.1J$, for bath sizes of $N=2^{6}$ (blue) and $2^{10}$ (red), and relative position $\mathbf{n}_{12}=\mathbf{n}^{A}_{1}-\mathbf{n}^{B}_{2}=(1,1)$. We observe (almost) complete oscillations with a very slight dependence of both the frequency and amplitude on the system size, $N$.

Let us first try to explain these numerical results through perturbative approaches. In this case, the collective self-energy within the Born-Markov approximation is finite and real in the limit $N\rightarrow\infty$, i.e., $\Sigma^{AB}_{\pm}(i0^{+};\mathbf{n}_{12})=J_{M,\pm}(\mathbf{n}_{12})=\pm J_{AB,M}(\mathbf{n}_{12})\in\mathbb{R}$, which predicts therefore coherent oscillations. We use Eq. 29 with $z=0$ to calculate numerically $J_{AB,M}(\mathbf{n}_{12})$ that we plot in solid green in Fig. 5(b). An asymptotic expression of $J_{AB,M}(\mathbf{n}_{12})$ for large distances $|\mathbf{n}_{12}|\gg 1$ reads Basko (2008); Bena (2009) (see Appendix):

[TABLE]

where $\mathbf{m}_{12}=(m_{1},m_{2})=\left(\frac{3}{2}\left(n_{1}+n_{2}\right),\frac{\sqrt{3}}{2}\left(n_{1}-n_{2}\right)\right)$ is a rescaled vector in terms of the original coordinates of $\mathbf{n}_{12}=(n_{1},n_{2})$. For the situation we are plotting in Fig. 5, i.e., $\mathbf{n}_{12}=(n,n)$, reads:

[TABLE]

Summing up, perturbative approaches predict the QEs should experience purely coherent interactions even in the limit $N\rightarrow\infty$. However, as we show in the single QE situation, the perturbative approach may introduce artifacts at points where the self-energy is not analytical. Thus, it is relevant to analyze to what extent the previously mentioned results survive after an exact analysis.

When we take into account the analytical structure of $\Sigma_{\pm}^{\beta}$, we find that the dynamics in the symmetric/antisymmetric subspace is given by:

[TABLE]

where $R_{\pm}$ is the residue associated to the pole of $G_{\pm}(z)$ at $z_{\pm}=J_{\pm}-i\frac{\Gamma_{\pm}}{2}$, whereas $C_{\pm,\mathrm{MBC}}(t)$ is the contribution of the detours associated to the branch cut in the middle of the band. The rest of the contributions can be shown to be small for $g/J\ll 1$. In Fig. 5(b), we plot the results obtained by solving the exact pole equation for $N=10^{2},10^{3}$ and $10^{4}$, showing that indeed the coherent dipole-dipole interactions, $z_{\pm}=\pm J_{AB}\in\mathbb{R}$, survive even for large systems. Interestingly, one can show that the exact dipole-dipole interactions now depend on the system size, $N$, renormalizing to:

[TABLE]

This expression has a very intuitive meaning: the qBS is able to mediate coherent interactions as long as the overlap with a QE excitation, $R_{0}(N)$, is finite. In the thermodynamic limit, these coherent interactions vanish as $R_{0}=0$. However, they decrease so slowly with the system size, $J_{AB}(\mathbf{n}_{12},N)\propto 1/\log(N)$, that in practical situations they can be harnessed even for very large systems, as shown in Fig. 5(a-b). Another physical picture to understand the decrease of the interactions with the system size is that the mode volume of the qBS, which is the one mediating the interactions, grows with the system size. As it occurs in cavity QED Ritsch et al. (2013), the strength of the light-mediated interactions is inversely proportional to the mode volume of the photon mediating the interactions. As the mode volume for the qBS is infinite in the thermodynamic limit, the mediated interactions are strictly zero in that limit.

The last remaining thing to prove is that the associated residue for the exact pole, i.e., the overlap with the initial wavefunction, is finite. We calculate it numerically using the expression:

[TABLE]

and show the results in Fig. 5(c). Moreover, we are also able to obtain an asymptotic expansion of $R_{\pm}$ in the limit $g/J\ll 1$, showing that $R_{\pm}\approx R_{0}$ (see Appendix), which tell us that the overlap can be very large as long as $\frac{g^{2}}{J^{2}}\log(N)\ll 1$. To sum up, the perturbative results for two QEs coupled to the AB lattices get corrected in two ways: i) the dipole-dipole interactions depend on the system size as $J_{AB}(\mathbf{n}_{12};N)\propto J_{AB,M}(\mathbf{n}_{12})/\log(N)$, and therefore vanish in the limit $N\rightarrow\infty$; ii) the residues of the associated poles also depend on the system size $R_{\pm}\propto 1/\log(N)$, also vanishing in the thermodynamic limit.

Finally, we predict that the interaction of QEs of baths with Dirac Cone dispersion relations leads to very exotic models if one considers many QEs: apart from the long-range character of the interactions, they occur only within different sublattices. The predicted many-body spin Hamiltonian can be written as:

[TABLE]

where $\mathbf{n}_{A},\mathbf{m}_{B}$ denotes that the coupling is only induced between the QEs coupled to different sublattices and $J_{AB}$ decreases very slowly with the system size, as $1/\log(N)$. This will certainly give rise to very rich many-body behavior without analogue in previously studied models.

V Experimental considerations

Along this paper we have discussed a very general model, which can be implemented in very different platforms, ranging from atoms coupled to the guided photons in photonic crystal structures Bravo-Abad et al. (2012); Xie et al. (2014) or even purely atomic implementations using cold atoms trapped in state-dependent optical lattices de Vega et al. (2008); Navarrete-Benlloch et al. (2011). In the latter, two atomic states must be trapped in deep/shallow potentials to play the role of the QE/photons respectively. This can be achieved, for example, with Alkaline-Earth atoms in which an optically excited state exists with very long lifetime (of the order of seconds Schreiber et al. (2015); Daley et al. (2008); Snigirev et al. (2017)) which allows one to design very different trapping potentials for ground/excited states. This is a very attractive platform because it combines: i) large flexibility in the design of the trapping potentials, and therefore, of the bath geometry; ii) low decoherence rates $\sim$ Hz Bloch et al. (2008); Schreiber et al. (2015); Daley et al. (2008); Snigirev et al. (2017) compared to the tunneling ($J$) and coupling ($g$) rates which can be of the order of $10$ KHz Bloch et al. (2008); iii) the possibility of single-site addressing and detection Sherson et al. (2010); Weitenberg et al. (2011).

Irrespective of the realization, the coupling to external reservoirs will induce extra decoherence in the bath/QE modes at rates $\kappa/\Gamma^{*}$ respectively. To illustrate the effect of the losses on the observed phenomena, we will consider a particular situation where the solution to the dynamics can be obtained analytically. As was shown in Ref. González-Tudela and Cirac (2017b) when $\kappa=\Gamma^{*}\equiv\Gamma_{\mathrm{loss}}$, the solution of the density matrix of the combined system of the QE and bath is given by:

[TABLE]

where $\left|\Psi(t)\right\rangle=e^{-iHt}\left|\Psi(0)\right\rangle$, i.e., the evolved state in the absence of losses. Thus, the main effect of the losses is to exponentially attenuate with time of the dynamics induced by the system bath coupling in $\left|\Psi(t)\right\rangle$. With that solution it is possible to estimate the conditions to observe the emergent phenomena described in this paper.

Observing the complete logarithmic relaxation will be challenging due to its slow decay. However, by going to situations where $g\sim O(J)$, which can be obtained in the cold atoms realization, one can make the dynamics faster. For example, in Fig. 3 we see how for $g=0.5J$, the QE has relaxed to approx. 20% from its initial value at times $tJ\approx 1000$, which will be within the reach of state-of-the-art parameters for the cold atom setup. Regarding the oscillations between two QEs at a given distance, $R$, they will be visible as long as $J_{AB}(\mathbf{n}_{12})/\Gamma_{\mathrm{loss}}\sim\frac{g^{2}}{JR}\gg 1$. In that case, the finite lifetime $\Gamma_{\mathrm{loss}}^{-1}$ sets a critical distance where these oscillations can still be observed.

VI Conclusions

To sum up, we have studied the quantum dynamics emerging when QEs are interacting though a reservoir with hexagonal lattice symmetry, focusing on the physics emerging because of the Dirac cone dispersion relation of the bath modes. For a single QE, we show that an exact treatment predicts the fractional decay of excitations and a logarithmic relaxation which can not be captured by perturbative treatments. Moreover, when several QEs are coupled to the reservoir purely long range coherent dipole interactions emerge between QEs coupled to different sublattices, which leads to very exotic many-body spin Hamiltonians. Further flexibility to engineer the QEs coupling can be obtained by introducing multi-level atoms and magnetic field gradients Douglas et al. (2015); González-Tudela et al. (2015a); Hung et al. (2016). Apart from the photonic realizations Bravo-Abad et al. (2012); Xie et al. (2014), we foresee that these type of 2D structured reservoirs can be obtained in a wide variety of systems ranging from cold atoms Polini et al. (2013) in state dependent optical lattices de Vega et al. (2008); Navarrete-Benlloch et al. (2011) or circuit QED Houck et al. (2012).

Acknowledgements.

Acknowledgements. The work of AGT and JIC was funded by the EU project SIQS and by the DFG within the Cluster of Excellence NIM. AGT also acknowledges support from Intra-European Marie-Curie Fellowship NanoQuIS (625955). We also acknowledge discussions with T. Shi, Y. Wu, S.-P. Yu, J. Muñiz, and H.J. Kimble.

Appendix 1 Single QE dynamics

As we show in the main text, the honeycomb like structure makes that finite size effects play a relevant role in the calculations. For that reason, we first show in Section 1.1 the results associated to the continuum limit, $N\rightarrow\infty$, where one can obtain an analytical expression of $\Sigma_{e}(z)$. Afterwards, in Section 1.2 we study how finite size effects correct the continuum predictions.

1.1 Continuum limit

The key point of this Section is that the analytical expression of $\Sigma_{e}(z)$ can be obtained by taking the continuum limit, that is,

[TABLE]

In this case, the $\mathbf{k}$ integrals can be calculated exactly Horiguchi (1972), first by demanding that $z\in\mathbb{R}$ and $z\notin[-3J,3J]$ and then, analytically continuing it to the rest of the complex plane. By doing so, one obtains the following expression:

[TABLE]

where $K^{\mathrm{I}}(m)$ has to be defined in a piecewise manner to guarantee analyticity in its definition domain:

[TABLE]

being $K(m)$ the complete elliptic integral of the first kind defined by:

[TABLE]

The regions with different signs of $\mathrm{Im}[z^{2}]$ are $\mathrm{Im}[k(z)]$ plotted in Figs.A1(a-b), together with the ones for the sign of $\mathrm{Re}[k(z)]$, which gives us the complete information about of the phase of the argument of the elliptic function [as $\mathrm{Im}[k^{2}(z)]=2\mathrm{Re}[k(z)]\mathrm{Im}[k(z)]$]. Evaluating $\Sigma_{e}(z)$ slightly above the real axis, $\Sigma_{e}(E+i0^{+})=\delta\omega_{e}(E)-i\frac{\Gamma_{e}(E)}{2}$, we obtain the functions $\delta\omega_{e}(E)$ and $\Gamma_{e}(E)$, that we plot in Fig. 1(b). There, we observe that the function $\Gamma_{e}(E)$ [$\delta\omega_{e}(E)$] has a discontinuity [divergence] at $E=\pm 3J$ appearing at the border of the band edges of the upper/lower band, $\pm\omega(\mathbf{k})$. Moreover, at $E=\pm J$ the function also displays a divergence [discontinuity] at $E=\pm J$. Finally, at $E=0$ both functions are continuous, but have discontinuous derivative, therefore, being also non-analytical at that point.

From these functions we can already obtain the Markov prediction, that tells us that a single QE excitation decays at a rate $\Gamma_{M}=\Gamma_{e}(\Delta)$ to the bath. Thus, it predicts that: i) outside of the bands, i.e., $\Delta\notin[-3J,3J]$, or exactly at the Dirac point, $\Delta=0$, the QE should not decay because $\Gamma_{M}=0$; ii) Inside of the band, $\Delta\in[-3J,0)\cup(0,3J]$, the QE must decay with a rate $\Gamma_{M}$, which is infinite at $\Delta=\pm J$. However, we have seen in the main text, that these predictions get corrected when calculating the exact dynamics, so let us explain how these corrections appear within this formalism and which is their origin.

A typical way to obtain the exact dynamics of the probability amplitude $C_{e}(t)$ consists in calculating the Fourier Transform by closing the contour of integration in the lower plane of the complex plane ($\mathrm{Im}[z]<0$). However, when the $\Sigma_{e}(z)$ have regions where it is non-analytic, as the one of Eq. A2, one must be careful when closing the contour to avoid them. One possible choice for the contour of integration is depicted in Fig. A1(d), where we take 5 detours of the integration contour at $E=\pm 3J,\pm J$ and [math]. Moreover, in the regions II-V, we analytically continue the function $\Sigma_{e}(z)$ into other Riemann sheets changing $K^{\mathrm{I}}(m)\rightarrow K^{\alpha}(m)$. The functions $K^{\alpha}(m)$ are obtained by imposing that $\Sigma_{e}^{\alpha}$ is continuum along the integration contour of Fig. A1. We exploit the fact that the elliptic functions are doubly periodic in which the different branches of $K(m)$ are given by combinations $pK(m)+qiK(1-m)$, with $p,q\in\mathbb{Z}$ Abramowitz et al. (1966). Using that, we obtain:

[TABLE]

To make it more clear, in Fig. A1(e) we make a diagram with the different combinations of $(p,q)$ in all the regions relevant for the integration contour defined in Fig. A1(d). The results of the complete integration are shown in the main text.

1.2 Finite system

Due to the very slow decay of the excitations in the infinite case, it is expected that finite size effects play a relevant role in the discussion. In fact, we already show in Fig. 3 of the main text how the logarithmic spontaneous emission quenches after a finite time, $N$, and oscillates around a constant value which depends on the ratio $g/J$. In Fig. A2(a), we evidence that this effect depends strongly on the system size by plotting the decay of a QE with a fixed $g=0.1J$ for different number $N$ of bath modes. There, we observe that the quenching time appears at a longer timescale the larger the system size. To further characterize this phenomenon, we go back to the finite system self-energy expression, which reads:

[TABLE]

It is clear that $z=0$ is a solution from the pole equation, because $\Sigma_{e}(0)=0$. Its residue, which can be interpreted as the overlap of the initial wavefunction with the quasi-bound state (qBS) emerging around the QE, is given by:

[TABLE]

where $g(N)$ is a function that only depends on the system size $N$. This function can be approximated by going to the continuum limit and changing to the $\mathbf{p}$-coordinates where the linear dispersion is isotropic in $\mathbf{p}$ [$q_{1,2}=\frac{3}{2}\left(p_{1}\pm\frac{1}{\sqrt{3}}p_{2}\right)$], which leads:

[TABLE]

where we include a minimum/maximum value for the region of integration to regularize the divergences of the integral. The minimum momentum is naturally provided by the discreteness of the lattice, i.e., $p_{\mathrm{min}}\propto\frac{1}{N}$. Thus, the $g(N)$ depends on the system size as follows:

[TABLE]

where $C$ is a constant that depends on the momentum cut-off $p_{c}$ and that we determine numerically by evaluating numerically the sum at a finite $N$, finding $C\approx 0.2$. Thus, the final residue of the qBS can be approximated by:

[TABLE]

which goes to [math] logarithmically with the systems size, $N$. The intuitive explanation for this is that for any finite system, the qBS is indeed square integrable as the rest of the eigenstates of the total Hamiltonian. Thus, the overlap with the initial wavefunction can be finite. However, for an infinite system the very slow decay of the qBS [$1/r$] makes it [marginally] not square integrable. Thus, it belongs to the continuous spectrum which has zero overlap with the initial wavefunction. In any practical situation, the relaxation dynamics will be a combination of a logarithmic decay for initial times, together with an absence of decay for long times due to the existence of the qBS.

Appendix 2 Two QE dynamics

As we reviewed in the main text, the two QE dynamics simplifies because the symmetric and antisymmetric subspaces decouple. This is convenient because their corresponding probability amplitudes can be calculated through a similar Fourier transform than the one of a single QE but with a modified self-energy: $\Sigma^{\beta}_{\pm}(z;\mathbf{n}_{12})=\Sigma_{e}(z)\pm\Sigma^{\beta}_{12}(z;n_{12})$, where the collective self-energy reads:

[TABLE]

where $D_{\beta}$ is a function which depends on whether the QE coupled to the $\beta=$AA/BB/AB sublattices. As we show in the main text $D_{\mathrm{AA}}=D_{\mathrm{BB}}=z$, whereas $D_{\mathrm{AB}}=D_{\mathrm{BA}}^{*}=f^{*}(\mathbf{k})$. In this case, it is difficult to obtain the analytical expression of the self-energy for every $z$ and $\mathbf{n}_{12}$, and the expressions not be very enlightening Horiguchi (1972). Thus, we decide to follow the intuition developed from the single QE analysis to obtain some general expression for the $\Delta=0$ situation for both finite and infinite systems.

2.1 Two QEs coupled to the AB sublattices

The more attractive situation for coherent interactions is the case where the two QEs coupled to different sublattices, where:

[TABLE]

To find an exact expression of $\Sigma_{12}^{AB}(z;\mathbf{n}_{12})$ is hard because of the complicated shape of the energy dispersion $|f(\mathbf{k})|$. However, as in the previous Section one can obtain a closed simple expression for $\mathrm{Re}(z)/J\ll 1$ by taking it as the sum around the $\mathbf{K}_{\pm}$ points and making the transformation to the $\mathbf{p}$-coordinates where the dispersion is isotropic. The contribution coming from the $\mathbf{K}_{\pm}$ point can be finally writen:

[TABLE]

where $\mathbf{h}_{p,\pm}=(-1,\pm i)$. Notice, we have redefined the $\bar{z}=2z/(3J)$ to make it adimensional and written a rescaled real space coordinates, $\mathbf{m}_{12}=(m_{1},m_{2})=\left(\frac{3}{2}\left(n_{1}+n_{2}\right),\frac{\sqrt{3}}{2}\left(n_{1}-n_{2}\right)\right)$. From here, a final (assymptotic) expression can be obtained:

[TABLE]

where $H_{1}^{1}(x)$ is the Hankel Function of the first kind. For small arguments, the Hankel function satisfies: $H_{1}(x)\approx-\frac{2i}{\pi x}+O(x)$, such that:

[TABLE]

For example, for $n_{1}=n_{2}=n$ as we are choosing for Fig. 5 of the main manuscript, where $\mathbf{m}_{12}=(3n,0)$ and $|\mathbf{m}_{12}|=3n$.

[TABLE]

For $n_{1}=-n_{2}=n$, then $\mathbf{m}_{12}=(0,\sqrt{3}n)$ and $|\mathbf{m}_{12}|=\sqrt{3}n$, we have:

[TABLE]

where $\sin\left(\frac{4\pi n}{3}\right)=-\frac{\sqrt{3}}{2},\frac{\sqrt{3}}{2},0,-\frac{\sqrt{3}}{2},\frac{\sqrt{3}}{2},0,\dots$.

Thus, within the Markov approximation the theory predicts coherent interactions decaying as a power law without any exponential attenuation, and with no imaginary component. However, we know from our previous experience with the single QE that one must be careful when extracting conclusions as non-Markovian corrections can be important around the Dirac point. Thus, in the next Section we deal with the exact problem and try to extract general conclusions.

2.1.1 Exact integration

We know that the pole equation for the symmetric/antisymmetric component for two QEs with $\Delta=0$ is given by:

[TABLE]

In the main text, we show how for finite systems the solution to this pole equation is strictly real, that is, it is a real BS with no decay and which can give rise to oscillations [see Fig. 5]. In order to further characterize it, we use the intuition developed from the previous Section which tell us that $\Sigma^{AB}_{\pm}(z;\mathbf{n}_{12})$ goes to a function independent of $z$ as $z\rightarrow 0$. Thus, we can try to solve the previous Equation self-consistently by expanding the right-hand side of Eq. A24 around $z=0$, arriving to:

[TABLE]

which yields:

[TABLE]

where $R_{0}$ is the residue for the (quasi)-bound state that appears for a single QE at $\Delta=0$. Thus, the Markovian coherent interactions get renormalized by a factor $R_{0}\propto 1/\log(N)$, which goes to [math] in the thermodynamic limit. However, the convergence to zero with the system size is so slow that for any practical situation one can indeed observe the coherent interactions as shown in Fig. 5 of the main text. The last remaining thing is to calculate the associated residue of the symmetric/antisymmetric pole, which is given by:

[TABLE]

This is what we show in Fig. 5 of the main text, showing how it can still be large for large systems. To estimate the conditions under which this residue can be large, we can expand $\partial_{z}\Sigma^{AB}_{\pm}(z)$ for small $z$ [as we know $z_{\pm}\sim O\left(\frac{g^{2}}{J\log(N)}\right)$]:

[TABLE]

It can be easily shown that $\partial_{z}\Sigma^{AB}_{\pm}(0)=-\frac{g^{2}}{J^{2}}g(N)$, that is the same function that controls the single QE residue. The only thing left to calculate is the contribution associated to the second derivative at $z=0$, which is given by:

[TABLE]

Following similar recipes as the ones we have used for the rest of the results, that is, calculating the contribution around $\mathbf{K}_{\pm}$, changing to the $\mathbf{p}$ coordinates where $f(\mathbf{p})$ is isotropic, we arrive to:

[TABLE]

As the minimum momentum is set by: $p_{\mathrm{min}}=\frac{2\pi}{N}$ and $z_{\pm}\sim O\left(\frac{g^{2}}{J\log(N)}\right)$, then we find $\partial^{2}_{z}\Sigma^{AB}_{\pm}(z)|_{z=0}z_{\pm}\sim O\left(\frac{g^{4}}{J^{4}}\right)$. Thus, we can neglect the second term in the expansion of Eq. A28 as long as $g^{2}/J^{2}\ll 1$. Finally, the residue $R_{\pm}$ will be large as long as $g^{2}\log(N)/J^{2}\ll 1$.