Coherent and radiative couplings through 2D structured environments

TL;DR

This paper investigates how two or more quantum units interact through 2D structured environments, analyzing coherence, radiation, directionality, and entanglement decay in finite and infinite lattices, with implications for quantum technologies.

Contribution

It provides a comprehensive comparison of coherent and radiative interactions in 2D lattices, including finite-size effects, directionality disruption, and subradiant state configurations, advancing understanding of quantum emitter dynamics.

Findings

Directionality linked to isofrequency manifolds can be disrupted.

Finite lattices exhibit different mode scaling and van Hove singularities.

Subradiant states are more prominent in finite lattices, with zero induced coherence in dark states.

Abstract

We study coherent and radiative interactions induced among two or more quantum units, by coupling them to two-dimensional lattices acting as structured environments. This model can be representative of atoms trapped near photonic crystal slabs, trapped ions in Coulomb crystals or to surface acoustic waves on piezoelectric materials, cold atoms on state-dependent optical lattices, or even circuit QED architectures, to name a few. We compare coherent and radiative contributions for the isotropic and directional regimes of emission into the lattice, for infinite and finite lattices, highlighting their differences and existing pitfalls, e.g. related to long-time or large-lattice limits. We relate the phenomenon of directionality of emission with linear-shaped isofrequency manifolds in the dispersion relation, showing a simple way to disrupt it. For finite lattices, we study further details as the scaling of resonant number of lattice modes for the isotropic and directional regimes, and relate this behavior with known van Hove singularities in the infinite lattice limit. Further we export the understanding of emission dynamics with the decay of entanglement for two quantum, atomic or bosonic, units coupled to the 2D lattice. We analyze in some detail completely subradiant configurations of more than two atoms, which can occur in the finite lattice scenario, in contrast with the infinite lattice case. Finally we demonstrate that induced coherent interactions for dark states are zero for the finite lattice.

Equations

$${\omega }_{\vec{k}}={\omega }_0-2J(\cos k_x+\cos k_y)$$

$${\omega }_{\vec{k}}=\sqrt{{\omega }_0^2+8J({\sin }^2\frac{k_x}{2}+{\sin }^2\frac{k_y}{2})}$$

$$H_{\text{int}}=\lambda (a_{1}A({\vec{r}}_1{)}^{†}+a_{2}A({\vec{r}}_2{)}^{†}+h.c.)$$

$${\dot{\rho }}_S=-i[{\tilde{H}}_S(t),{\rho }_S]+\sum _{j,l=1}^2{\Gamma }_{jl}(\vec{r},t)(a_j{\rho }_S{a}_l^{†}-\frac{1}{2}\{a_l^{†}{a}_j,{\rho }_S\})$$

$${\lambda }_{LS}(\vec{r},t)=\frac{-{\lambda }^2}{(2\pi {)}^2}{\int }_{-\pi }^{\pi }{d}^2\vec{k}\frac{1-\cos [t(\Omega -{\omega }_{\vec{k}})]}{\Omega -{\omega }_{\vec{k}}}\cos (\vec{k}\vec{r})$$

$${\Gamma }_c(\vec{r},t)=\frac{2{\lambda }^2}{(2\pi {)}^2}{\int }_{-\pi }^{\pi }{d}^2\vec{k}\frac{\sin [t(\Omega -{\omega }_{\vec{k}})]}{\Omega -{\omega }_{\vec{k}}}\cos (\vec{k}\vec{r})$$

$${\lambda }_{LS}(\vec{r},\infty )=\frac{-{\lambda }^2}{(2\pi {)}^2}PV{\int }_{-\pi }^{\pi }{d}^2\vec{k}\frac{\cos (\vec{k}\vec{r})}{\Omega -{\omega }_{\vec{k}}}$$

$${\Gamma }_c(\vec{r},\infty )=\frac{2{\lambda }^2}{(2\pi {)}^2}{\int }_{-\pi }^{\pi }{d}^2\vec{k}\delta (\Omega -{\omega }_{\vec{k}})\cos (\vec{k}\vec{r})$$

$${\dot{\rho }}_S=\sum _{j=\pm }-i[H_j,{\rho }_S(t)]+{\Gamma }_j(\vec{r},t)(a_j{\rho }_S{a}_j^{†}-\frac{1}{2}\{a_j^{†}{a}_j,{\rho }_S\})$$

$$f_{\vec{n},\vec{k}}=\frac{2}{N+1}\sin \left(k_x{n}_x\right)\sin \left(k_y{n}_y\right)$$

$${\Gamma }_0=2{\lambda }^2\sum _{k_x,k_y=1}^N\frac{\sin [t(\Omega -{\omega }_{\vec{k}})]}{\Omega -{\omega }_{\vec{k}}}|f_{\vec{n},\vec{k}}{|}^2$$

$${\Gamma }_c=2{\lambda }^2\sum _{k_x,k_y=1}^N\frac{\sin [t(\Omega -{\omega }_{\vec{k}})]}{\Omega -{\omega }_{\vec{k}}}{f}_{\vec{n},\vec{k}}{f}_{{\vec{n}}^{\prime },\vec{k}}$$

$${\Gamma }_0(t)=\frac{2{\lambda }^{2}t}{(2\pi {)}^2}\sum _{\vec{k}}\mathrm{sinc}(\mathrm{t}(\mathrm{\Omega }-{\omega }_{\vec{\mathrm{k}}})).$$

Full text

Engineering the coherent interaction between many quantum units is a key ingredient to simulate quantum phases of many-body systems, whereas radiative interactions is instrumental to producing dark states and generating nonclassical states of light. One of the most advanced platforms for quantum simulators are e.g. 2D Coulomb crystals of trapped ions thompson ; Bollinger1 ; Bollinger2 , with vibrational quanta of the ions’ motion being used as a channel to effectively induce spin-spin interactions with tunable distance dependence. Excellent control of interactions between atoms and light in 1D waveguides has been recently shown painter2014 , leading to observation of many atoms superradiance painter2015 . This has spurred theoretical proposals to use the engineered properties of light to mediate strong long-range atom-atom interactions both in 1D chang2015 and 2D photonic crystal (PC) lattices tudela2015 ; tudela2016 ; galve2018 . Surface acoustic waves (SAW) on piezoelectric materials giedke2015 have been recently proposed as mediators between e.g. quantum dots, trapped ions, nitrogen-vacancy centers, or superconducting qubits. Control of mediators, be it phonons in the ion and SAW cases, or photons in the PC case, is thus of key importance to shape interactions to build future quantum simulators, as recently proposed for many-body physics chang2015 .

Structured environments, in addition to their usefulness to shape coherent interactions, are a good playground to understand effects of collective dissipation (radiative couplings) of distant units depending on the properties of the substrate. Understanding decoherence and dissipation in realistic scenarios is of utter importance for the successful realization of quantum technologies. This is for example the case when miniaturizing ion traps wine1 below the micrometer level, where details of decoherence sources, such as possible correlation lengths among adatoms’ dipoles dipoles ; Galve2017 could become accessible wine2 . Architectures of two or more units dissipating collectively can lead to super/sub-radiance Dicke ; Haroche and allow certain degrees of freedom to become noiseless noiseless ; pazroncaglia , which is of interest to realizing noise-free operations. Radiative coupling between bosonic probes enabling collective dissipation have been reported in cubic or triangular lattices, either isotropically decaying, or purely directional and long-range Galve2016 . This translates directly into sub/super-radiant dynamics, as compared to independent decay rates (no radiative interaction) between quantum units. Long range radiative couplings not only expand the toolbox of quantum optics, but have important practical implications for quantum networking qinternet1 ; qinternet2 or quantum memories qmem1 ; qmem2 ; qmem3 . Furthermore, understanding decoherence and dissipation in a broader class of platforms is of utter importance for the successful realization of quantum technologies, where environment microscopic details play a key role, like for instance for miniaturized ion traps wine1 ; dipoles ; Galve2017 .

Periodic structures can be described as tight-binding models for photons or phonons, disregarding the specific details of the materials used, and so allowing to study fundamental features. We focus on 2D structures because they allow richer features than 1D lattices, and can be implemented in e.g. PC lattices tudela2015 , 3D arrays of (evanescently-coupled) light waveguides whose propagation direction $z$ represents the flow of simulated time Sciarrino1 ; Sciarrino2 , 2D Coulomb crystals of trapped ions thompson ; Bollinger1 ; Bollinger2 or superconducting circuits cQED0 ; cQED1 ; cQED2 ; cQED3 to name a few.

In this work we characterize coherent and radiative couplings induced by a 2D periodic (squared lattice) structure, both infinite and finite, highlighting their similarities and differences. The paper is organized as follows: in section I we introduce the model and its description in terms of a time-dependent master equation for the infinite lattice case, where we show the behavior of coherent and radiative contributions for isotropic and directional regimes. We analyze finite-time coefficients, explaining their meaning, and also the relation of long-time decay rates with known van Hove singularities Tudela1 ; Tudela2 ; vanHove53 . We analyze the phenomenon of directional emission, showing how it can be broken, and highlight some pitfalls related to the ordering of large-lattice and long-time limits. An analysis on the validity of the infinite lattice limit is given. The decay dynamics of entanglement between atomic or bosonic units due to the 2D structure is then studied, showing that it maps with the behavior expected from the cross-talk prediction. In section II we introduce the finite lattice, and analyze how the divergence of decay rates due to van Hove singularities are related to the number of resonant modes in the lattice as a function of time. We study how to construct dark states of two, three and four atoms from the form of the dissipation matrix, and show that the form of the Lamb-shift matrix leaves dark states invariant. We finish by comparing exact dynamics of two-atom entanglement in the finite-lattice case, with the evolution under the master equation with time-dependent coefficients.

I Infinite lattice

We start by analyzing the interaction of a two-unit, bosonic or atomic, quantum system coupled to a structured substrate such as that of Fig. 1 (there a PC lattice is represented). Our analysis holds for periodic, bosonic, linear lattices (modeled by quadratic Hamiltonians). For a square-lattice symmetry dispersion is either

[TABLE]

for coupled cavities or tight-binding models, or

[TABLE]

for harmonic (spring-like) couplings. It can be seen that both dispersions have almost the same shape if they are scaled properly to fit in the range of values. The tight-binding dispersion ranges from $\omega_{0}-4J$ to $\omega_{0}+4J$, whereas the harmonic one ranges from $\omega_{0}$ to $\sqrt{\omega_{0}^{2}+16J}$, so when we use the latter we set $J=3\omega_{0}^{2}/16$ so that $\omega_{\vec{k}}\in[1,2]\omega_{0}$ (as in Fig. 1). Of course, we can offset all energies so that the frequency band lies e.g. $\omega_{\vec{k}}\in[-1,1]\omega_{0}$ for both dispersions, being only important the detuning between emitters and the resonant pseudomomentum manifold in the lattice. The cases we study in this paper behave equally for both dispersion, because the implied resonant manifolds are the same, so we will interchangeably use both, making the choice explicit only when we use the tight-binding one.

The two identical and independent emitters can be either bosons, $H_{S}=\Omega(a_{1}^{\dagger}a_{1}+a_{2}^{\dagger}a_{2})$, or two-level atomic systems (TLS), $H_{S}=\Omega(\sigma_{1}^{+}\sigma_{1}^{-}+\sigma_{2}^{+}\sigma_{2}^{-})$. We consider local coupling to the structure

[TABLE]

or $H_{\text{int}}=\lambda(\sigma_{1}^{-}A(\vec{r}_{1})^{\dagger}+\sigma_{2}^{-}A(\vec{r}_{2})^{\dagger}+h.c.)$ for TLS.The annihilation operators $A(\vec{r})$ for the light field (atoms interacting with PC) or surface displacement field (SAW case), are evaluated at the positions of the atoms (bosons) $\vec{r}_{1}$ and $\vec{r}_{2})$. Interactions extended in space over a range $\sigma$ (instead of local) are known to introduce frequency cut-off functions which limit interaction strength for lattice’s frequencies above $\sim 1/\sigma$ Galve2016 . The interaction between system and lattice $\lambda$ is assumed to be weak, allowing for a Born-Markov treatment breuer . The dynamics, for environment at T=0, is given by

[TABLE]

with $\tilde{H}_{S}(t)=H_{S}+H_{LS}(t)$ the Lamb-Shift corrected Hamiltonian, $H_{LS}(t)=\Omega_{LS}(t)(a_{1}^{\dagger}a_{1}+a_{2}^{\dagger}a_{2})+\lambda_{LS}(\vec{r},t)(a_{1}a_{2}^{\dagger}+h.c)$ and $\vec{r}=|\vec{r}_{1}-\vec{r}_{2}|$. We define the cross-damping $\Gamma_{c}\hat{=}\Gamma_{12}=\Gamma_{21}$ and the self-damping $\Gamma_{0}\hat{=}\Gamma_{11}=\Gamma_{22}$ rates. Using $A(\vec{r})=(2\pi)^{-1}\int_{-\pi}^{\pi}d^{2}\vec{k}\ \exp(-i\vec{k}\cdot\vec{r})A(\vec{k})$ (i.e. infinite-lattice), the spatial dependent coefficients are

[TABLE]

[TABLE]

and the onsite coefficients are just $\Omega_{LS}(t)=\lambda_{LS}(0,t)$, $\Gamma_{0}(t)=\Gamma_{c}(0,t)$. They are the sum of counter-propagating plane waves, weighted by their resonance with $\Omega$. The emission of only one unit is proportional to the same integral (note that due to the $\vec{k}\to-\vec{k}$ symmetry we can replace cos by an exponential), and thus the spatial features of $\Gamma_{c}$ give a precise idea of the radiation pattern of one emitter. We retain, and analyze later, the time dependence of the coefficients in the master equation, which account for the build-up of a communication channel between quantum units Benedetti . This accounts to having performed Born and first-Markov approximations breuer . Usually in the literature the second-Markov approximation (long-time $t\to\infty$ limit) is used:

[TABLE]

[TABLE]

The dynamics for identical uncoupled probes can be diagonalized at all times, in the common and relative coordinates $a_{\pm}=(a_{1}\pm a_{2})/\sqrt{2}$, resulting in two completely independent dynamics:

[TABLE]

with $H_{\pm}=(\Omega+\Omega_{LS}\pm\lambda_{LS})a_{\pm}^{\dagger}a_{\pm}$ and $\Gamma_{\pm}=\Gamma_{0}\pm\Gamma_{c}$. The latter expression indicates that those modes can become noiseless when $\Gamma_{c}=\pm\Gamma_{0}$, which is known to lead to preservation of entanglement at long times pazroncaglia ; Galve2016 . When $\Gamma_{c}=+\Gamma_{0}$ we have $\Gamma_{-}=0$ and thus the mode $a_{-}$ is not dissipating, what is typically called a common bath situation (CB). With $\Gamma_{c}=0$ both modes dissipate, the separate baths scenario (SB). Otherwise, if $\Gamma_{c}=-\Gamma_{0}$ so we have $\Gamma_{+}=0$ and thus it is the center of mass motion which is noiseless Galve2017 ; zoller2008 . We termed this situation anti-common bath (aCB) Galve2017 .

The equivalent expressions for the two TLS case are obtained by substituting $a_{j}\to\sigma_{j}^{-}$ and $a_{j}^{\dagger}\to\sigma_{j}^{+}$. The dynamics is then diagonal in the operator basis $\sigma_{\pm}^{-}=(\sigma_{1}^{-}\pm\sigma_{2}^{-})/\sqrt{2}$, which has as fixed points $|00\rangle$ and $|\Psi_{-}\rangle=$ for CB ($\Gamma_{c}/\Gamma_{0}=1$) with definition $|\Psi_{\pm}\rangle=(|01\rangle\pm|10\rangle)/\sqrt{2}$. The state $|00\rangle$ is trivial in the sense that the lattice is a $T=0$ reservoir. However, if the system starts with $|\Psi_{-}\rangle$, entanglement is conserved asymptotically. For the aCB case ($\Gamma_{c}/\Gamma_{0}=-1$) it is $|\Psi_{+}\rangle$ which is conserved, and none is conserved asymptotically when SB ($|\Gamma_{c}|<\Gamma_{0}$).

I.1 Finite time coefficients

The coefficients $\Omega_{LS}$ and $\Gamma_{0}$, equal to their single-unit counterparts, can be seen in figure 1. Both settle to constant values rather fast, at $\omega_{0}t\sim 10$, justifying the broadly used long-time limit master equations. The single-unit decay rate $\Gamma_{0}$ also settles quickly to a constant, except for the specific splitting $\Omega=\sqrt{5/2}\omega_{0}$ (from now on we will work with the harmonic dispersion, unless otherwise estated), where it has its highest value (for $J=3/16\omega_{0}^{2}$). The decay rate at this particular $\Omega$ is divergent, see Fig. 1d right, in the long time limit $\Gamma_{0}(t\to\infty)\propto\int_{-\pi}^{\pi}d^{2}\vec{k}\ \delta(\Omega-\omega_{\vec{k}})=\int_{\vec{k}_{\Omega}}d^{2}\vec{k}\ |\vec{\nabla}\omega(\vec{k}_{\Omega})|^{-1}$. It is the sum of all inverse group velocities (density of states, DoS, see fig. 1c) of the manifold of excitations $\vec{k}_{\Omega}$ resonant with $\Omega$. For $\Omega=\sqrt{5/2}\omega_{0}$ (red line) it picks up points (black regions in plot) where the density of states diverges, known as van Hove singularities vanHove53 . Due to the periodic structure of reciprocal space these singularities are known to appear generically Morse1925 and lead to logarithmic divergences in 2D. We also note that other points (not relevant for our work) have zero DoS, at the extrema of the band; here bound states form and perturbative approaches are no longer valid wang90 ; quang94 ; intJ ; boundCirac ; boundCicca .

The most important consequence of retaining time-dependent coefficients in the master equation, instead of just using their long-time limits (the usual Markovian master equation), comes from the cross-talk. From figures 2, 3 and 4 we see that it takes a non-negligible time for the coherent and radiative interactions to become nonzero for distant units. During that time, the two quantum units see independent environments because $\Gamma_{c}=0$, and only when it becomes nonzero they start ‘seeing’ a correlated environment. Take for example two units with $\Omega=1.01\omega_{0}$, as in figure 2, at a distance $\sim 30$: before $\omega_{0}t=1000$ they have $\Gamma_{c}=0$ and thus decay independently, whereas from that time on they have $\Gamma_{c}\sim\Gamma_{0}/3$, leading to $\Gamma_{-}/\Gamma_{+}\sim 1/2$. This means that an initial state $|\psi_{+}\rangle$ will decay twice faster than $|\psi_{-}\rangle$, but only from that time onwards.

If the system frequency is increased reaching the middle of the band, fig. 3 ($\Omega=\sqrt{5/2}\omega_{0}$), it was predicted in 2Dphotonic ; Galve2016 that it leads to ‘diagonal-only’ propagation of excitations, and thus to decay of cross-damping in all directions except for the diagonals. This behavior has been checked against the exact (atom-lattice case) dynamics (see galve2018 ; Tudela1 ; Tudela2 ).

We plot next the Lamb shift LS (coherent) coupling $\lambda_{LS}(\vec{r})$ and the cross-damping coefficient $\Gamma_{c}(\vec{r})$, without the irrelevant factor $\lambda^{2}/[(2\pi)^{2}]$ for the 3 cases highlighted in Fig.1c, corresponding to increasing frequency of the units interacting with infinite planar squared lattices. In the first case ($\Omega\simeq\omega_{0}$) the system is resonant with the low-momentum manifold of excitations in the lattice, corresponding to an approximately isotropic dispersion relation $\omega_{\vec{k}}\simeq\sqrt{\omega_{0}^{2}+2g|\vec{k}|^{2}}$, as can be seen in the spatial shape of the propagation in Fig. 2. The case of high*-k* ($\Omega=1.999\omega_{0}$) in Fig. 4 resonates with the ‘corners’ of the dispersion relation (Fig.1c) and, although apparently anisotropic, can be seen to lead also to isotropic decay of the master equation, with a superimposed alternating oscillation. Writing for one of the four corners $\vec{k}=(\pi-\delta_{x},\pi-\delta_{y})$, we have $\cos(\vec{k}\vec{r})=\cos[\pi(x+y)]\cos(x\delta_{x}+y\delta_{y})$. Integration of variables $\delta_{x,y}$ is again isotropic and thus leads to a $J_{0}(|\delta||\vec{r}|)$ decay, times the oscillating factor $\cos[\pi(x+y)]$. Finally, the ‘diagonal’ case ($\Omega=\sqrt{5/2}\omega_{0}$) is a peculiar case where the excitations conspire to produce a purely diagonal propagation galve2018 as displayed in Fig. 3.

At long times distant points are both coherently and radiatively coupled with different strengths (orange curves respectively in Fig. 2, 3,4). It is to be noted that the spatial dependence of coherent and incoherent contributions (for each case) are rather similar although not equal, see e.g. the beating of long time Lamb-shift as compared to cross-damping (top fig. 3), or their different decay behaviors (top fig. 4), since both consist in integrating $\cos(\vec{k}\vec{r})$ but with the different weights $[1-\cos(\Omega-\omega_{\vec{k}})]/(\Omega-\omega_{\vec{k}})$ and $\text{sinc}(\Omega-\omega_{\vec{k}})$. Both weights are most important when $\Omega\simeq\omega_{\vec{k}}$ but have different functional behavior, leading to these differences.

I.2 Radiative directional coupling

The ‘diagonal’ ($\Omega=\sqrt{5/2}\omega_{0}$ or $\Omega=\omega_{0}$ for the tight-binding dispersion relation) was highlighted in Galve2016 as displaying purely directional radiative coupling (also radiation by one emitter). This happens for cubic and triangular lattices Galve2016 , and also their higher-dimensional analogues, as well as for graphene graphene , being then present in common geometric configurations. It has been recently explained as a consequence of linear-shape iso-frequency manifolds in the dispersion relation galve2018 : a rotation into that momentum coordinate selects a specific value of that momentum, which produces a plane wave (or sine) in that direction, and is thus non-decaying. In the orthogonal coordinate we still sum many momenta, which gives generically an incoherent sum of waves and thus a decaying function along the corresponding spatial coordinate. A simple mechanism breaking the linear shape for a square-symmetric lattice was proposed in galve2018 , consisting in adding a hopping term between diagonal neighbors in the lattice. This results in e.g. a modified tight-binding dispersion relation $\omega_{\vec{k}}=\omega_{0}-2J(\cos k_{x}+\cos k_{y})-4\tilde{J}\cos k_{x}\cos k_{y}$, whose last term gives curvature to the isofrequency manifold. This is easily seen in Fig. 5.

Because of the van Hove singularities in this manifold, the case of directional emission allows inspection of some possible pitfalls related to the ordering of the large-lattice and long-time limits. If we first obtain the analytic shape of $\Gamma_{c}(\vec{r})$ for periodic boundary conditions by taking a finite lattice with $t\to\infty$ (which exactly selects $N/2$ modes on each of the red straight lines in fig. 1c) and taking then the large-lattice limit $N\to\infty$, which yields $\Gamma_{c}(\vec{r})/\Gamma_{0}=(x\sin\pi x-y\sin\pi y)/[\pi(x^{2}-y^{2})]$, i.e. $[{\rm sinc}(\pi|\vec{r}|)+\cos(\pi|\vec{r}|)]/2$ at the diagonals and $\Gamma_{c}(\vec{r})/\Gamma_{0}=\text{sinc}(\pi r_{y})$ at $r_{x}=0$. This function has a maximum value $1$ at $|\vec{r}|=0$ and after the sinc decays, it has $1/2$ at the diagonals and [math] elsewhere. It is a peculiar fact that interchanging the limits gives a different result. By setting $t<N$ and doing $N\to\infty$ in fact gives $\Gamma_{c}(\vec{r})/\Gamma_{0}=1$. At intermediate $N$ we see a shape similar to directional radiation but with a region near $\vec{r}=0$ with almost maximum cross-talk (see figure  6), i.e. higher than the $1/2$ value obtained when $t\gg N$. In addition, its absolute value is logarithmic-divergent, as is well-known from van-Hove singularities vanHove53 . The very different behavior can be explained as follows: in the first case, $t\gg N$, only $N/2$ modes are resonant, because the spacing in momentum space goes as $\Delta k\sim 1/N$ and $t$ is so large that the $sinc$ function selects a momentum manifold always smaller than $\Delta k$. In that way, we only sum the contribution of a finite number of modes. In spite of that, see later discussion on the scaling of the decay rate, the cross-talk diverges as $N\ t$. In the second case, it is even worse: for $t\ll N$ the amount of resonant modes selected by the $sinc$ is even bigger, and does not have only contributions from the momentum manifold yielding directionality; it has more contributions, and this is reflected in the modified cross-talk shape in Fig. 6.

I.3 Timescales and modes choice

A comment is in order about timescales. The possibility to use exponential waves as ansatz for guided modes in a lattice totally disregards that lattices are always finite, an assumption that can be perfectly valid within a range of parameters. Take for example the tight-binding dispersion and note that the highest group velocity $d\omega_{\vec{k}}/d\vec{k}$ (see Fig. 1c) is $2\sqrt{2}J$. Once any signal reaches the boundary of the lattice, bouncings and revivals of the system dynamics are expected to occur. This happens for times $J\ t\sim N$. The decay dynamics of the system is however governed by the system-lattice coupling $\lambda$, yielding a decay rate $\Gamma_{0}\sim\lambda^{2}$, and thus significant decay is expected to happen at times $\sim 1/\lambda^{2}$. Thus the exponential wave ansatz is valid whenever the characteristic decay time is at least smaller than the revival of signals at the boundaries. This implies $1/\lambda^{2}<N/J$, and thus requires $\lambda^{2}>J/N$. Taking $J$ as energy scale, and a lattice of $N\sim 100$ sites, already implies $\lambda\sim 0.1$, within the strong coupling regime, as in Tudela1 . For this reason it is advisable to use finite-lattice mode functions, as we do in later sections.

I.4 Entanglement decay

We analyze next how an initial entangled state of the two bosonic/TLS probes evolves as a result of the interaction with an infinite lattice. We consider a two-mode squeezed vacuum 2mode1 ; 2mode2 ; 2mode3 for the oscillator case, and one of the two maximally entangled Bell states $|\Psi_{\pm}\rangle=(|01\rangle\pm|10\rangle)/\sqrt{2}$ for the two TLS case. To quantify entanglement for two TLS, we use the concurrence $C$ wootters and for the bosonic case the logarithmic negativity vidalwerner . We start with the dynamics when the atom-lattice interaction is so weak that the master equation coefficients can be assumed to reach their constant asymptotic values, corresponding physically to the regime where waves between atoms have enough time to become stationary (i.e. first and second Markov approximations). We show results for a couple of TLS for different times $\omega_{0}t$ (to be compared with the system-bath interaction time $1/\lambda^{2}=10^{4}/\omega_{0}$) and distances $\vec{r}$ in figure 7, where we plot the surviving concurrence for initial states $|\Psi_{\pm}\rangle$. They follow the qualitative behavior of the cross-damping’s spatial dependence, including the alternating behavior of $|\Psi_{\pm}\rangle$ as the sign of $\Gamma_{c}$ changes.

There is noiseless behavior (total subradiance) only if we are able to make $\vec{r}\ll 1/|\vec{k}_{\Omega}|$ (low k), for $|\Psi_{-}\rangle$, or with $\Omega$ at the high-k edge of the band placing the TLS at odd ($|\Psi_{-}\rangle$) or even ($|\Psi_{+}\rangle$) distances from each other. It must be noted though, that when the units are near the band edges they form atom-photon (or ion-phonon) bound states which cannot be treated under Born-Markov conditions wang90 ; quang94 ; intJ . The interesting ‘diagonal’ case features $\Gamma_{c}(\vec{r})/\Gamma_{0}=1/2$ as maximum values, so it never produces noiseless dynamics.

The case of a bosonic two-mode state is, in comparison, insensitive to the distinction between CB and aCB, as can be seen in fig. 8 (in some sense each degree of freedom acts as each Bell state, but now both are present in the two-particle dynamics). This leads to similar although faster decay. In order to see a behavior similar to atomic Bell states, we need to initialize the state as $|r\rangle_{+}\otimes|0\rangle_{-}$ (or $|0\rangle_{+}\otimes|r\rangle_{-}$), the equivalent of $|\psi_{+}\rangle$ (or $|\psi_{-}\rangle$, respectively), with the mode corresponding to operator $a_{+}$ squeezed and other in vacuum (or vice versa). In comparison, the two-mode squeezed state is $|r\rangle_{+}\otimes|-r\rangle_{-}$. In this asymmetrically squeezed case, the dynamics is qualitatively similar to atomic Bell states, as can be seen in Fig. 9.

II Finite size lattices

Finite size PC lattices typically need finite element calculations to compute the guided modes bookPHOT , and are leaky at their boundaries. However we can use guided modes coupled to the atoms: these modes have evanescent tails at the outside of the crystal. Also, we could enforce Dirichlet boundary conditions for all allowed modes (i.e. also guided resonance modes which are within the light cone and thus can interact with the outer modes, allowing extraction of light JoanoGUIDED ) with Bragg mirrors at the boundaries Bragg . Assuming that in practice this can be achieved with enough precision, traveling waves are then replaced by standing waves

[TABLE]

with $k_{x,y}=\pi l_{x,y}/(N_{x,y}+1)$, $n_{x,y}\in[1,N_{x,y}]$ and the same for $l_{x,y}$. $N_{x,y}$ is given by the amount of dielectric function changes in each direction, typically associated with holes made in a dielectric slab, or pillars above it bookPHOT . For the finite size lattice we will take the tight-binding dispersion relation, which is as before $\omega_{0}-2J(\cos k_{x}+\cos k_{y})$. A unitary transformation of the Hamiltonian removes the irrelevant offset $\omega_{0}$, which from now on we take as zero. Now the damping coefficients become

[TABLE]

[TABLE]

and accordingly the LS coherent contributions.

II.1 Number of resonant modes

Once a discrete spectrum is assumed, consistent with the finite size of the lattice, one can understand the divergence in the decay rate as a scaling of the number of resonant modes as follows. Let us take for simplicity the case of periodic boundaries, so $f_{\vec{n},\vec{k}}=e^{-i\vec{k}\cdot\vec{n}}/2\pi$ and $|f_{\vec{n},\vec{k}}|^{2}=(2\pi)^{-2}$. The decay rate in this case is

[TABLE]

For long times, this sum approximately counts how many modes in the crystal are resonant with $\Omega$, because for them the sinc function is 1, and 0 elsewhere. In figure 10 we can see that for frequencies $\Omega\neq 0$, at the beginning all modes are resonant, but as time grows, the number of resonant modes quickly reaches a scaling $t^{-1}$, making the decay rate tend to a constant already for $J\cdot t\sim 10$, consistent with figure 1d. This is a common situation, justifying the use of long-time master equations (Markov 2 approximation breuer with constant coefficients. However, for $\Omega=0$ we see that the number of resonant modes scales as $~{}t^{-0.8}$, and thus the decay rate scales as $t^{0.2}$. At long times it saturates at a value $4\times N/2$ and does not fall more, because the resonant manifold defined by the red square in fig. 1b ($N/2$ modes in each of the four lines) perfectly fits the $sinc$ function in eq. (7), whatever $t$. From this time on, the decay rate scales as $N\cdot t$, diverging for $N\to\infty$ as expected from a van Hove singularity.

II.2 Multi-atom dark states

Dark states (perfect subradiance) can arise among distant atoms in 1D but not in 2D isotropic environments. In the following we show that the ‘diagonal’ case with its linear radiation pattern can be used to construct subradiant states, in analogy with the 1D case where $\Gamma_{c}^{(1D)}/\Gamma_{0}=\cos(k_{\Omega}x)$ painter2015 ; Galve2016 , i.e. it makes $|\Psi_{\pm}\rangle$ subradiant for specific distances between atoms due to coherent cancellation of waves. It is easy to see geometrically in an open 2D configuration that, because both atoms radiate with a cross pattern, only one of the lines (the one which joins their positions) can cancel. As noted recently Tudela1 ; Tudela2 , no subradiant 2D configuration of two atoms exists in infinite planar settings. Geometrically again, we see that if we can join the remaining two lines, total cancellation can occur galve2018 . Indeed this is the case for reflecting boundaries, as seen in Fig. 11a. Only four positions for the two atoms will yield total cancellation (equivalently, $\Gamma_{c}/\Gamma_{0}=1$, red points): placing two atoms on any of those points with a state $|\Psi_{-}\rangle$ yields total subradiance, up to retardation effects. We stress that these dark states can be created at arbitrary distances, depending on the off-set from the center of the PC slab.

The reflection of (photonic/phononic) waves can be used to build peculiar radiation patterns of multi-unit architectures. In Fig. 11b,c,d we show the cross-damping $\Gamma_{c}/\Gamma_{0}$ of one atom with respect to other atom at position $(r_{x},r_{y})$: one unit (green) shares a partial subradiant condition ($\Gamma_{c}/\Gamma_{0}=1/2$) with another unit (red), the red unit shares partial subradiance with the blue unit ($\Gamma_{c}/\Gamma_{0}=1/2$), but the green and blue units radiate independently. We have plotted the cross-talk of each of the atoms with the rest of the possible atomic positions in the crystal, to make it more understandable. In this particular case the decay rate matrix would be $\Gamma_{i,j}=\{1,1/2,0\},\{1/2,1,1/2\}\{0,1/2,1\}\}$, with the slowest decay channel corresponding to state $|\psi\rangle=(|100\rangle-\sqrt{2}|010\rangle+|001\rangle)/\sqrt{2}$ (with atom notation $|$green, red, blue$\rangle$) radiating  5.8 and  3.4 slower than the other two orthogonal states of three atoms.

Another interesting case is to put four atoms, e.g. each in the points highlighted in fig 11c (rhomboid configuration). This leads to $\Gamma_{i,j}=1$, $\forall i,j$ and thus to one superradiant state with decay $4\Gamma_{0}$, $|\psi_{\text{super}}\rangle=(|1000\rangle+|0100\rangle+|0010\rangle+|0001\rangle)/2$ or simply $|\psi_{\text{super}}\rangle=\frac{1}{2}(1,1,1,1)$ (expressed in the basis $\{|1000\rangle,|0100\rangle,|0010\rangle,|0001\rangle\}$, with atom notation $|$up,bottom,left,right$\rangle$). Also to completely subradiant states of two (e.g. $(1,0,0,-1)/\sqrt{2}$), three (e.g. $(-1,0,2,1)/\sqrt{3}$) and four atoms (e.g. $(1,-3,1,1)/2\sqrt{3}$; see also galve2018 ). Notice that one can keep several of the atoms in the ground state and they will not affect the subradiance of the remaining excited atoms. In contrast, in infinite lattices the decay rates we have the same superradiant state with lower decay rate $2\Gamma_{0}$, a totally subradiant four atom state Tudela1 ; Tudela2 $|\psi_{1}\rangle=\frac{1}{2}(-1,-1,1,1)$, and the states of two or three atoms now radiate with $\Gamma_{0}$. These examples show the potential of this platform to explore the exotic scenarios arising when all sites are filled with atoms in different architectures.

II.3 Lamb-shift coupling

One could worry about the stability of dark states because of the Lamb-shift coupling, which in principle could modify the phases of dark states turning them into bright. For the most interesting ‘diagonal’ case, and in contrast with the infinite lattice scenario, the values of $\lambda_{LS}$ are far smaller than the radiative coupling $\Gamma_{c}$. In spite of both having the same spatial shape (see Fig. 12, now taking the tight-binding dispersion), we observe (not shown) that the ratio $\lambda_{LS}/\Gamma_{c}$ grows linearly with a rate $\sim 10^{-6}$, i.e. it takes $J~{}t\sim 10^{5}$ to reach $\lambda_{LS}/\Gamma_{c}\sim 0.1$. In addition, let us place four atoms in the four special points of maximum $\Gamma_{c}$ in a subradiant state. Because the matrix of Lamb-shift couplings $\lambda_{LS}(\vec{r}_{i},\vec{r}_{j})$ (with $\lambda_{LS}(\vec{r}_{i},\vec{r}_{i})=\Omega_{LS,i}$) has the same signs and structure as the cross-talk matrix $\Gamma_{c}(\vec{r}_{i},\vec{r}_{j})$, it is easy to show that it does not perturb dark states galve2018 . This is because dark states are eigenstates of the cross-talk matrix with eigenvalue 0, and thus also of the Lamb-shift matrix.

II.4 Entanglement decay: perturbative vs. exact

We finish by comparing the validity of our master equation with time-dependent coefficients and the exact dynamics of the model in the one-excitation sector galve2018 ; Tudela1 . Take two atoms in the configuration of Fig. 12, at positions $\vec{r}_{1}=(25,20)$, $\vec{r}_{2}=(25,30)$, with initial Bell states $|\psi_{\pm}\rangle$. Since at those points $\Gamma_{c}/\Gamma_{0}=+1$, we should, and do, observe no decay of the concurrence for $|\psi_{-}\rangle$ and fast decay for $|\psi_{+}\rangle$. While the shapes are not exactly the same, the master equation faithfully captures the time it takes for the cross-talk to build-up. During this time, both Bell states decay because $\Gamma_{c}$ has not yet reached its final value. From that time on, only one of the Bell states decays further. It is also interesting to see that, once the concurrence of $|\psi_{+}\rangle$ has fallen to zero, it can never go up. This is a clear effect of the Born approximation which takes the state of the lattice to be always the T=0 thermal state (vacuum). Because of this, no excitation can ever go back to the atoms. We also observe that our master equation underestimates decay, as compared to the exact dynamics which gives oscillatory behavior (revivals). Finally, for system-lattice strong coupling $\lambda=0.1J$ we start observing more pronounced differences in the dynamics. This is a very interesting subject for future studies.

III Conclusions

We have studied the interaction of two probes induced by a 2D structured periodic environment. The physics of this problem is generic enough that it can be applied to e.g. 2D Coulomb crystal of trapped ions thompson ; Bollinger1 ; Bollinger2 , PCs with nearby trapped atoms painter2014 ; painter2015 ; chang2015 ; tudela2015 , 3D printed photonic circuits Sciarrino1 ; Sciarrino2 , trapped-ions near piezoelectric substrates giedke2015 or even circuit QED architectures cQED0 ; cQED1 ; cQED2 ; cQED3 , only requiring translational invariance of the lattice and weak system-lattice coupling. Coherent (Lamb Shift) and incoherent (radiative) couplings are studied both for finite and infinite lattices, analyzing many details and possible pitfalls. We show the usefulness of the master equation description, in terms of radiative couplings between quantum units, by showing that it directly translates into the behavior of entanglement decay of two units. Further, for more than two units, we show how to use the radiative coupling matrix to build dark states.

IV Acknowledgements

We acknowledge fruitful discussions with A. González-Tudela and J.I. Cirac. This work has been supported by the EU through the H2020 Project QuProCS (Grant Agreement 641277), by MINECO/AEI/FEDER through projects NoMaQ FIS2014-60343-P, QuStruct FIS2015-66860-P and EPheQuCS FIS2016-78010-P. FG acknowledges funding from ‘Vicerectorat d’Investigació i Postgrau’ of the UIB.