Quantum Emitters in Two-dimensional Structured Reservoirs in the Non-Perturbative Regime

TL;DR

This paper explores the complex quantum dynamics of quantum emitters coupled to a two-dimensional structured reservoir, revealing non-perturbative effects like directional emission and perfect subradiance, which are not captured by traditional models.

Contribution

It introduces a non-perturbative analysis of quantum emitters in 2D reservoirs, uncovering novel phenomena such as directional emission and subradiance.

Findings

Exponential relaxation with non-Golden rule rate

Directional emission into the reservoir

Perfect subradiance for four emitters

Abstract

We show that the coupling of quantum emitters to a two-dimensional reservoir with a simple band structure gives rise to exotic quantum dynamics with no analogue in other scenarios and which can not be captured by standard perturbative treatments. In particular, for a single quantum emitter with its transition frequency in the middle of the band we predict an exponential relaxation at a rate different from that predicted by the Fermi's Golden rule, followed by overdamped oscillations and slow relaxation decay dynamics. This is accompanied by directional emission into the reservoir. This directionality leads to a modification of the emission rate for few emitters and even perfect subradiance, i.e., suppression of spontaneous emission, for four quantum emitters.

Equations

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

$$G_e(z)=\frac{1}{z-\Delta -{\Sigma }_e(z)},G_{\mathbf{k}}(z)=g\frac{G_e(z)}{z-\omega (\mathbf{k})},$$

$${\Sigma }_e(E+i{0}^{+})\approx \frac{g^2}{4J}\left[\mathrm{sgn}(E)-\frac{2i}{\pi }\log \left(\frac{|E|}{16J}\right)\right].$$

$${\bar{\Gamma }}_e\approx \frac{g^2}{\pi J}\log \left(\frac{32\pi J^2}{g^2}\right)$$

$$\underset{t\to \infty }{\lim }{C}_{\mathbf{k}}(t)=\frac{g{e}^{-i\omega (\mathbf{k})t}}{\omega (\mathbf{k})-\Delta -{\Sigma }_e(\omega (\mathbf{k}))},$$

$${\Sigma }_{4,-}(z)=\frac{4g^2}{{\pi }^2}{\iint }_0^{\pi }d\mathbf{q}\frac{{\sin }^2(2q_{x}n){\sin }^2(2q_{y}n)}{z+4J\cos (q_x)\cos (q_y)},$$

$$C_{4,-}(\infty)=\frac{1}{1-\partial_{z}\Sigma_{4,-}(z)}\Big{|}_{z=0}=\frac{1}{1+\frac{g^{2}n^{2}}{J^{2}}}\,.$$

Full text

The interaction of quantum emitters (QEs) with propagating bosonic particles, e.g., photons, lies at the core of Quantum Optics Cohen-Tannoudji et al. (1992). This interaction leads, for example, to collective interactions between QEs Lehmberg (1970a, b) which can be harnessed for both quantum information and simulation applications. New avenues in the integration of QEs with nanophotonic structures Vetsch et al. (2010); Huck et al. (2011); Hausmann et al. (2012); Laucht et al. (2012); Thompson et al. (2013); Goban et al. (2014); Lodahl et al. (2015); Goban et al. (2015); Sipahigil et al. (2016); Bermúdez-Ureña et al. (2015); Béguin et al. (2014) provide us with systems in which the QEs interact with low dimensional bosonic modes, with complicated energy dispersions in the case of engineered dielectrics Hausmann et al. (2012); Laucht et al. (2012); Thompson et al. (2013); Goban et al. (2014); Lodahl et al. (2015); Goban et al. (2015); Sipahigil et al. (2016). Despite originally the main motivation of such implementations was to exploit the small sizes to enhance light-matter interactions, it was soon realized that intriguing phenomena arise because of the reduced dimensionality. One particular aspect is the possibility of realizing chiral emission Mitsch et al. (2014); Söllner et al. (2015); Lodahl et al. (2016), which can display very uncommon features Ramos et al. (2014); Pichler et al. (2017). Another one is the possibility of exploiting the phenomena of sub and superradiance Goban et al. (2015); Corzo et al. (2016); Sørensen et al. (2016); Sipahigil et al. (2016); Solano et al. (2017), e.g., to enhance the coupling to the emitter Chang et al. (2012), to generate QE entanglement González-Tudela and Porras (2013); Ramos et al. (2014), to produce non-classical light González-Tudela et al. (2015, 2016), or even to perform quantum computation Paulisch et al. (2016).

The dynamics of QEs in 1D reservoirs is relatively simple, specially when their transition frequency, $\omega_{e}$, lies within a band. Perturbative treatments predict that a single QE initially excited decays at a rate, $\Gamma$, given by the Fermi’s Golden Rule (FGR), i.e. proportional to the density of states of the bath at $\omega_{e}$. The emission mostly occurs in the bath modes that are resonant with that frequency. Typically, there are two such modes of associated momentum $\pm k_{e}$, leading to a symmetric left/right emission. When two (or more) QEs are present, the existence of only two such modes leads to super/subradiant states, where the emission is enhanced or suppressed by interference. In higher dimensions for structureless baths, a single QE will also decay at a rate given by the FGR. However, the emission takes place in different directions as there are many resonant modes in the bath. For two QEs, the interference in the emission cannot occur in all those modes at the same time 111Unless their separation is smaller than the associated wavelength Dicke (1954), something we are not considering here and thus, the phenomena of sub and superradiance are generically absent.

In this manuscript we use non-perturbative methods to analyze the dynamics of QEs interacting with a simple two-dimensional (2D) structured reservoir. By structured we mean with a periodic structure giving rise to a dispersion relation containing frequency bands. In particular, we contrast our results with the predictions of perturbative treatments based on a markovian master equation approach. First, we show how a single excited QE with its energy tuned in the middle of the band shows an exponential decay at a rate different from that predicted by FGR. Moreover, at longer times, this exponential relaxation is followed by an oscillation and subexponential dynamics. These QE dynamics are followed by a directional emission into the bath in two orthogonal quasi-1D directions, as also predicted for classical sources of light Mekis et al. (1999) and sound Langley (1996). As a consequence, when several QEs are coupled to the bath, this directional emission induces anisotropic collective dissipation. For two QEs we observe a modification of the spontaneous emission rate as a consequence of this directional emission when the QE lie on a line at $45$ degrees. A related behavior has been predicted in Galve et al. (2017) using perturbative master equations. We find that a total suppression of the emission is not possible for two QEs, and explain this fact in terms of a partial interference effect. In contrast, we show how to design a perfect subradiant state with four QEs which survives even in the non-perturbative regime, since in that case the interference can be fully destructive,

We assume a 2D bath with a square-like symmetry described by $N\times N$ bosonic modes with energy $\omega_{a}$ and with nearest neighbour coupling $J$. The Hamiltonian is given by (using $\hbar=1$) $H_{B}=-J\sum_{\langle\mathbf{n},\mathbf{m}\rangle}\left(a_{\mathbf{m}}^{\dagger}a_{\mathbf{n}}+\mathrm{h.c.}\right)$, where $\mathbf{n}=(n_{x},n_{y})$ is a vector indicating the bosonic mode position within the lattice. We have used a rotating frame at a frequency $\omega_{a}$, such that the zero energy corresponds to the center of the band structure (see below). Despite the simplicity of the model, we expect it to describe more complex materials such as photonic crystals, in the same way that the 1D tight-binding model does for structured waveguide QED. The bath Hamiltonian can be diagonalized in $\mathbf{k}$ space by introducing periodic boundary conditions and the operators $a_{\mathbf{k}}=\frac{1}{N}\sum_{\mathbf{n}}e^{-i\mathbf{k}\mathbf{n}}a_{\mathbf{n}}$, where $\mathbf{k}=(k_{x},k_{y})$ and $k_{x},k_{y}=\frac{2\pi}{N}(-\frac{N}{2},\dots,\frac{N}{2}-1)$, such that $H_{B}=\sum_{\mathbf{k}}\omega(\mathbf{k})a_{\mathbf{k}}^{\dagger}a_{\mathbf{k}}$ with $\omega(\mathbf{k})=-2J\left[\cos(k_{x})+\cos(k_{y})\right]$. In Fig. 1, we plot the energy dispersion of the band $\omega(\mathbf{k})$ together with the associated density of states in the limit $N\rightarrow\infty$, $D(E)=\frac{1}{(2\pi)^{2}}\iint d\mathbf{k}\delta[E-\omega(\mathbf{k})]$. The band extends from $[-4J,4J]$. At the band edges, the $D(E)$ is nearly constant as predicted for isotropic dispersions González-Tudela et al. (2015). At the middle of the band it displays a divergence associated to the saddle point appearing at the energy dispersion $\omega(\mathbf{k})$, as it also happens for real materials Mekis et al. (1999). As we show below, this has important consequences in the dynamics beyond purely enhancing the emission. We also consider one (or several) QEs described as two-level systems $\{\left|g\right\rangle_{j},\left|e\right\rangle_{j}\}$, with transition frequency, $\omega_{e}$, whose Hamiltonian reads: $H_{S}=\Delta\sum_{j}\sigma^{j}_{ee}$. We use the notation $\sigma^{j}_{\alpha\beta}=\left|\alpha\right\rangle_{j}\left\langle\beta\right|$ for the spin operators and $\Delta=\omega_{e}-\omega_{a}$ represents the detuning with respect to the middle of the band. Finally, we assume a local coupling of the QEs to the bath modes, described by $H_{\mathrm{int}}=g\sum_{j}\left(a_{\mathbf{n}_{j}}\sigma_{eg}^{j}+\mathrm{h.c.}\right)$. We assumed to be in a parameter regime where the QEs are coupled to a single band of the bosonic reservoir.

Along this manuscript we consider the QE(s) to be initially excited in certain QE state $\left|\Phi_{0}\right\rangle_{S}$, whereas the bath starts initially empty, i.e., $\left|\mathrm{vac}\right\rangle_{B}=\left|0\right\rangle^{\otimes N^{2}}$. Then, we let the system free to evolve under the total Hamiltonian $H=H_{S}+H_{B}+H_{\mathrm{int}}$, and study the QEs relaxation. The situation when the QE energies lie outside the band, $|\Delta\pm 4J|\gg g$, where the dynamics are dominated by the presence of a bound state, has been explored extensively in other works (see, e.g., Bykov (1975); John and Wang (1990); Kurizki (1990)). We focus here instead on the situation when the QE energies lie within the band, $\Delta\in[-4J,4J]$, with emphasis on what happens in the middle of the band. In our illustrations, we use relatively large values of $g$ in order to emphasize the features in the figures; however, the conclusions can be extended to smaller $g$’s that is the common situation in the optical regime. Nevertheless, we make a complete discussion about the whole range of parameters and give more details about the calculation in the accompanying paper González-Tudela and Cirac (2017).

We first consider that a single QE is coupled to the bath, such that we drop the index $j$ in the QE operators, and assume $\left|\Phi_{0}\right\rangle_{S}=\left|e\right\rangle$. As the initial state contains a single excitation, the state at any time, $t$, can be written as follows:

[TABLE]

In Fig. 2 we show the results of the numerical integration of the dynamics of the QEs for several detunings for $\Delta/J=-3,-2,-1,0$ and $g=0.1J$, which we complement with the state of the bath excitations, i.e., $|C_{\mathbf{n}}|$, at $tJ=100$ for $\Delta=-3J,0$. From Fig. 2(a), it seems that the decay of the QE is basically exponential, with an enhanced decay rate as its energy is tuned closer to the center of the band. Naively, this is what one expects from perturbative approaches, which predict $|C_{e}(t)|^{2}\approx e^{-\Gamma_{e}(\Delta)t}$, with a decay rate given by FGR, $\Gamma_{e}(\Delta)=2\pi g^{2}D(\Delta)$, where $D(\Delta)$ is the density of modes evaluated at the QE energy $\Delta$. For $\Delta/J=-3,-2,-1$ such prediction works well; however, for $\Delta=0$ an straightforward application of FGR predicts an infinite decay rate that we do not observe in Fig. 2(a). Moreover, when plotting the population dynamics in logarithmic scale (inset) we observe that the relaxation is actually non-monotonic showing both an oscillation and a subexponential decay for long times.

To gain analytical insight into this exotic relaxation, we apply standard techniques Cohen-Tannoudji et al. (1992) to rewrite the probability amplitudes in terms of their Fourier transform $C_{e,\mathbf{k}}(t)=\frac{i}{2\pi}\int_{-\infty}^{\infty}dEG_{e,\mathbf{k}}(E+i0^{+})e^{-iEt}$, where:

[TABLE]

and $\Sigma_{e}(z)=\frac{g^{2}}{N^{2}}\sum_{\mathbf{k}}\frac{1}{z-\omega(\mathbf{k})}$ is the so-called QE self-energy that captures the effect of the coupling to the bath on the QE dynamics. From now on we focus on what happens for energies, $E$, around the middle of the band ($|E|\ll J$) where the self-energy can be expanded as Morita (1971); Economou (1983); González-Tudela and Cirac (2017)

[TABLE]

In this expression we observe that around $E=0$ the real part, that we denote as $\delta\omega_{e}(E)$, has a discontinuous jump, whereas the imaginary one, denoted as $\Gamma_{e}(E)$, has a logarithmic divergence. For the latter, we use the same notation as the FGR decay rate because they are connected. More concretely, the standard perturbative approaches, such as the Markov approximation Gardiner and Zoller (2000), assume the self-energy to smoothly vary around $\Delta$ and replace $\Sigma_{e}(E+i0^{+})\approx\Sigma_{e}(\Delta+i0^{+})$, recovering the exponential relaxation of $|C_{e}(t)|^{2}$ with decay rate, $\Gamma_{e}(\Delta)$, as given by the FGR.

The divergence appearing in the middle of the band, however, forces us to go beyond the perturbative treatment to unravel the results. In the standard quantum optical scenario Cohen-Tannoudji et al. (1992) one calculates the exact Fourier transform of $C_{\alpha}(t)$ by closing the contour in the lower half complex plane $\mathrm{Im}(E)<0$, taking detours at the band edges of $\omega(\mathbf{k})$ because of the presence of branch cuts in the self-energy. The poles in the real axis describe bound states, and give rise to fractional decay John and Quang (1994). The complex ones lead to an exponential decay. Finally, branch cuts associated to the band-edges give rise to power-law decays, which are, however, typically hidden by the fractional decay induced by bound states John and Quang (1994). In the case studied here, an additional branch cut appears in the (negative) imaginary axis, associated to the divergence of the $D(E)$ in the middle of the band. This forces us to take an extra detour in the integration contour which has two visible consequences in the inset of Fig. 2(a): i) the slow relaxation dynamics at long times scaling as $O[\left(t\log(16Jt)^{2}\right)^{-2}]$; ii) there exists a regime of $\Delta\in[-\frac{g^{2}}{2J},\frac{g^{2}}{2J}]$, in which two unstable poles appear in the analytical continuation of $G_{\alpha}(z)$. For $\Delta=0$, their imaginary part coincides and is given by González-Tudela and Cirac (2017)

[TABLE]

obtained in the limit $\frac{32\pi J^{2}}{g^{2}}\gg 1$. Their real part is given by $\approx\pm\frac{g^{2}}{2J}$. This explains both the finite time scale observed in Fig. 2(a), and the oscillation observed in the inset, which has a frequency proportional to the difference of the real part of the unstable poles ($\sim\frac{g^{2}}{J}$). Both behaviours emerge from the failure of perturbative treatments due to the divergent density of states. The oscillations can also be intuitively understood as part of the $\mathbf{k}$-modes propagate very slowly allowing them to be reabsorbed by the QE. Remarkably, for short times one observes an exponential relaxation with a timescale given by Eq. 4 rather than the one expected from perturbative arguments that would be $O(g^{2}/J)$.

We also plot the emission into the bath in Figs. 2(b-c). The bath population can also be obtained exactly in the long-time limit, as the probability amplitude $C_{\mathbf{k}}(t)$ is still dominated by the pole contribution of $G_{\mathbf{k}}(E+i0^{+})$ of Eq. 2 at $E=\omega(\mathbf{k})$, which yields:

[TABLE]

This expression tells us that the $\mathbf{k}$-modes around $\omega(\mathbf{k})\approx\Delta$ for $g\ll J$ are the ones dominating the emission. This explains why when $\Delta$ is away from zero, the $\mathbf{k}$ modes are isotropically populated, as $\omega(\mathbf{k})\approx f(|\mathbf{k}|^{2})$. For $\Delta\approx 0$, we have that the modes that dominate fulfill $k_{x}\pm k_{y}=\pi$, as sketched in Fig. 1(a). This explains why the emission is anisotropic in this case: the only resonant modes fulfill $k_{x}\pm k_{y}\approx\pm\pi$, and thus propagate along the diagonals.

Now, we include several QEs in the discussion and take $\Delta=0$. We explore the interplay between the anisotropic emission and the relative position, $\mathbf{n}_{12}$ of the QEs, to check up to which point super/subradiance phenomena survive within a non-perturbative picture. We first study the scenario with two QEs prepared in a (anti)symmetric superposition $\left|\Phi_{0}\right\rangle_{S}=\left|\Phi_{2,\pm}\right\rangle=\frac{1}{\sqrt{2}}\left(\sigma_{eg}^{1}\pm\sigma_{eg}^{2}\right)\left|g\right\rangle^{\otimes 2}$. Interestingly, when $\omega(\mathbf{k})=\omega(-\mathbf{k})$ the dynamics of the symmetric/antisymmetric state separate because they couple to orthogonal bath modes. Consequently, their relaxation dynamics is calculated analogously to the one of a single QE, but replacing: $\Sigma_{e}\rightarrow\Sigma_{\pm}=\Sigma_{e}\pm\Sigma_{12}$, where $\Sigma_{12}=\frac{g^{2}}{N^{2}}\sum_{\mathbf{k}}\frac{e^{i\mathbf{k}\cdot\mathbf{n}_{12}}}{z-\omega(\mathbf{k})}$ is the collective QE interaction induced by the environment, which can be obtained recursively Morita (1971); Economou (1983); González-Tudela and Cirac (2017). Along this manuscript, we focus on a situation when the two QEs lie along a diagonal, $\mathbf{n}_{12}=(n,n)$, as this is where most of the emission occurs [see Fig. 2(c)], and may lead to modifications of collective decay Galve et al. (2017). We find that that perturbative Markov approaches, i.e., $\Sigma_{12}(i0^{+};\mathbf{n})/\Sigma_{e}(i0^{+})=(-1)^{n}$, which point to the possibility of perfect super/subradiance. However, from our single QE study we know at $\Delta=0$ perturbative approaches may fail.

To go beyond perturbative treatments, we first numerically integrate the dynamics and show that if $n$ is even (odd), the states $\left|\Phi_{\pm(\mp)}\right\rangle$ have enhanced (suppressed) decay rates, as expected from the propagation phases of the bath modes satisfying $k_{x}\pm k_{y}=\pi$, which are the ones dominating the dynamics when $\Delta=0$. We show an example of such dynamics in Figs. 3(a-b) for $\mathbf{n}_{12}=(6,6)$, where we observe collective effects leading to enhancement/suppression of spontaneous emission. Notice, however, that neither the suppression nor the enhancement is perfect as it occurs in 1D systems. The reason behind that can be intuitively understood from the bath population that we plot in Figs. 3(b-c). For the superradiant state, i.e., $\left|\Phi_{2,+}\right\rangle$, the emission occurs in 3 quasi-1D modes, one collective along the diagonal where the QEs are placed and two independent ones along the orthogonal directions. On the contrary, the subradiant state $\left|\Phi_{2,-}\right\rangle$ suppresses the emission in the diagonal where they are placed, while emitting into (two) quasi-1D modes along the orthogonal diagonals. These independent decay channels forbid finding a perfectly subradiant state with only two QEs. This enhancement/suppression can also be explained in terms of constructive/destructive interference of the bath modes emitted by the two QEs at $\omega(\mathbf{k})=0$, which are mainly given $k_{x}+k_{y}=\pm\pi$, such that their phases: $e^{i(k_{x}+k_{y})n}\pm e^{i(k_{x}+k_{y})n}=e^{\pm i\pi n}\pm e^{\pm i\pi n}$, add up constructively/destructively depending on the relative phase between QEs. Apart from that, we observe: i) retardation effects, as they also occur in 1D systems; ii) other non-Markovian effects introduced by the divergence of the density of states, which lead to a logarithmic correction of the super/subradiant decay with the distance González-Tudela and Cirac (2017).

The intuition obtained with two QEs allows us to build perfect subradiant states with four QEs. Let us consider four QEs at positions $(0,2n)$, $(2n,0)$, $(2n,4n)$ and $(4n,2n)$ and assume they are in a state $\left|\Phi_{4,\pm}\right\rangle=\frac{1}{2}\left(\sigma_{eg}^{1}\pm\sigma_{eg}^{2}+\sigma_{eg}^{3}\pm\sigma_{eg}^{4}\right)\left|g\right\rangle^{\otimes 4}$, The decay of $\left|\Phi_{4,+}\right\rangle$ is collectively enhanced with an emission in $8$ orthogonal directions as shown in Fig. 3(d). Starting out in $\left|\Phi_{4,-}\right\rangle$, the emission is completely suppressed (up retardation effects) because perfect destructive interference occurs in the eight emission directions, trapping the light between the four QEs [see Fig. 3(d)]. The analysis to show that the $\left|\Phi_{4,-}\right\rangle$ is perfectly subradiant relies on the fact that in this position configuration, the QE modes defined by $\left(\sigma_{eg}^{1}\pm\sigma_{eg}^{2}+\sigma_{eg}^{3}\pm\sigma_{eg}^{4}\right)$ couple to orthogonal bath modes, that allow us to consider their dynamics as those of a single QE, but with a modified self-energy. In particular, for the subradiant case the associated self-energy can be shown to be González-Tudela and Cirac (2017):

[TABLE]

where we have used the rotated axis coordinates $k_{x,y}=q_{x}\pm q_{y}$. If $\Sigma_{4,-}(z)$ vanishes at $z=0$, this implies that a real bound state emerges within the band. This can be shown to be the case because the integrand at $z=0$ is separable in $q_{x,y}$, and each integrand satisfies: i) $I(q_{x,y})=-I(\frac{\pi}{2}-q_{x,y})$ and ii) its divergence appearing because of the zeros in the denominator for the $\mathbf{q}$ modes satisfying $q_{x,y}=\pi/2$ is canceled by the one in the numerator. The last thing to show is that, indeed, its associated residue, connected to the steady-state population, is not zero. We can explicitly calculate it obtaining:

[TABLE]

which is therefore very close to $1$ as long as retardation effects are small $\frac{g^{2}n^{2}}{J^{2}}\ll 1$. In any case, we obtain that independent of the distances there always remains some excitation within the QEs.

Finally, let us briefly comment how the Hamiltonian $H$ can be obtained on a platform beyond engineered dielectrics Hausmann et al. (2012); Laucht et al. (2012); Thompson et al. (2013); Goban et al. (2014); Lodahl et al. (2015); Goban et al. (2015); Sipahigil et al. (2016), namely, cold atoms in state dependent optical lattices de Vega et al. (2008); Navarrete-Benlloch et al. (2011). The simplest scenario consists of a bosonic atom with two metastable states, $a$/$b$, trapped in a shallow/deep optical potential. In that situation, only the $a$-atoms tunnel to nearest neighbour at a rate, $J$, of the order of $10$ KHz in typical experiments Bloch et al. (2008). Thus, the $a/b$ atoms play the role of the propagating bath modes/QEs respectively. Their coupling is obtained through either two-photon Raman transition or, in the case of Alkaline-Earth atoms, through direct optical coupling and can be $g\sim O(J)$, which allows one to tune $g$ and $\Delta$, from [math] to several MHz, and thus investigate all the parameter regimes. Atomic motion effects can be suppressed by cooling the atom initially to the ground state and operating in the Lamb-Dicke regime. The scattering losses induced by the trapping potential, finite lifetime of the optically metastable state in Alkaline-Earth atoms, or other imperfections lead to decoherence in both the QE/bath modes which can be as small as $\Gamma_{\mathrm{loss}}\sim$ Hz Daley et al. (2008); Schreiber et al. (2015), such that $g/\Gamma_{\mathrm{loss}}\sim 10^{4}$. The impact of losses on the observation of the phenomenology has been considered in Ref. González-Tudela and Cirac (2017). The possibility of addressing single sites Weitenberg et al. (2011); Endres et al. (2013) and distinguishing the internal states, makes this setup ideal to observe our predictions.

To sum up, we have explored the non-perturbative QEs dynamics emerging from their coupling to a two-dimensional bosonic bath with a square-like geometry when their energies lie in the middle of the band. For a single QE, we predict an exponential relaxation at short times, but with a timescale that escapes Fermi’s Golden Rule description, which is followed by reversible and slow relaxation dynamics at longer times. Moreover, such phenomena are accompanied by strongly directional emission into the bath along two quasi-1D orthogonal directions, which lead to super/subradiant states when many QEs are coupled to the bath. For two QEs, the perturbative predictions are corrected due to the divergent density of states in the middle of the band, which forbids perfect super/subradiance. We characterize mathematically the phenomena and give an intuitive explanation in terms of interference. This understanding allows us to build perfect subradiant states with four QEs where the emission is trapped within them in spite of the 2D character of the bath.

Both the directionality and the divergent density of states responsible for the phenomena we describe, are associated to the saddle points in $\omega(\mathbf{k})$, and not its interplay with polarization as in 1D chiral systems Lodahl et al. (2016). Since those points are ubiquitous in 2D reservoirs, we conjecture that our findings will be relevant in more general situations beyond the simple model employed here Mekis et al. (1999). Apart from the platforms mentioned, the predictions can also be tested in circuit QED Astafiev et al. (2010); Hoi et al. (2013); Van Loo et al. (2013); Liu and Houck (2017).