Subradiant Bell states in distant atomic arrays

TL;DR

This paper demonstrates that two distant atomic arrays can support long-lived, subradiant Bell states that enable high-fidelity quantum information transfer and entanglement, with potential experimental realizations using cold atoms.

Contribution

It introduces a method to create and utilize subradiant Bell states in distant atomic arrays for quantum information processing, highlighting their long lifetime and non-local entanglement capabilities.

Findings

Supports long-lived Bell superposition states with strong subradiance

Enables high-fidelity quantum state transfer between arrays

Analyzes robustness against experimental imperfections

Abstract

We study collective `free-space' radiation properties of two distant single-layer arrays of quantum emitters as two-level atoms. We show that this system can support a long-lived Bell superposition state of atomic excitations exhibiting strong subradiance, which corresponds to a non-local excitation of the two arrays. We describe the preparation of these states and their application in quantum information as resource of non-local entanglement, including deterministic quantum state transfer with high fidelity between the arrays representing quantum memories. We discuss experimental realizations using cold atoms in optical trap arrays with subwavelength spacing, and analyze the role of imperfections.

Figures

Equations

$${\mathcal{H}}_{\mathbf{j},{\mathbf{j}}^{\prime }}\equiv -i({\gamma }_e/2){\mathbf{p}}^{\ast }\cdot \hat{\mathbf{G}}({\mathbf{r}}_{\mathbf{j}}-{\mathbf{r}}_{{\mathbf{j}}^{\prime }})\cdot \mathbf{p}$$

$${\gamma }_n=\Gamma [1+p_n\cos (k_{z}L)]\frac{k_z^2+|{\mathbf{q}}_n{|}^2/2}{k_0{q}_z},$$

$$|{\psi }_{d/b}⟩=\frac{1}{\sqrt{2N}}\sum _{{\mathbf{j}}_{\perp }}\left[{\sigma }_{({\mathbf{j}}_{\perp },1)}^{+}\mp (-1{)}^m{\sigma }_{({\mathbf{j}}_{\perp },2)}^{+}\right]|\mathcal{G}⟩$$

$$H_{\text{eff}}=\frac{\Omega }{\sqrt{2}}\left[{\sigma }_b^{+}(S_1^{-}+S_2^{-})+{\sigma }_d^{+}(S_1^{-}-S_2^{-})\right]+\text{h.c.}$$

$$\mathcal{F}\approx e^{-\pi \sqrt{2{\gamma }_d/{\gamma }_b}},$$

$$H_a=\sum _{\mathbf{j}}{\omega }_0{\sigma }_{\mathbf{j}}^{+}{\sigma }_{\mathbf{j}}^{-}+H_{\text{drive}},$$

$$V_{af}=-d\sum _{\mathbf{j}}({\sigma }_{\mathbf{j}}^{+}{\mathbf{p}}^{\ast }+{\sigma }_{\mathbf{j}}^{-}\mathbf{p})\cdot \hat{\mathbf{E}}({\mathbf{r}}_{\mathbf{j}})+\text{h.c.},$$

$$\hat{\mathbf{E}}(\mathbf{r})=i\int d\mathbf{k}\sum _{\lambda }{ϵ}_{|\mathbf{k}|}{b}_{\lambda ,\mathbf{k}}{e}^{i\mathbf{k}\cdot \mathbf{r}}{\mathbf{e}}_{\lambda ,\mathbf{k}},$$

$$\displaystyle\frac{d\rho}{dt}=-i\Big{[}H_{\text{drive}}+\text{Re}(H_{\text{dip}}),\rho\Big{]}-2$$

$${\mathcal{H}}_{\mathbf{j},{\mathbf{j}}^{\prime }}=-i({\gamma }_e/2)G({\mathbf{r}}_{\mathbf{j}}-{\mathbf{r}}_{{\mathbf{j}}^{\prime }}),$$

$$\displaystyle{\hat{\boldsymbol{G}}}({\boldsymbol{r}})=\frac{3e^{ik_{0}r}}{2i(k_{0}r)^{3}}\Big{[}$$

$$\displaystyle+\left(-(k_{0}r)^{2}-3ik_{0}r+3\right)\frac{{\boldsymbol{r}}\otimes{\boldsymbol{r}}}{r^{2}}\Big{]},$$

$$\nabla \times \nabla \times \hat{\mathbf{G}}(\mathbf{r})-k_0^2\hat{\mathbf{G}}(\mathbf{r})=-\frac{6\pi i}{k_0}\delta (\mathbf{r}),$$

$$\mathcal{L}(\rho )=\sum _{\mathbf{j},{\mathbf{j}}^{\prime }}\text{Im}(\mathcal{H}{)}_{\mathbf{j},{\mathbf{j}}^{\prime }}\left({\sigma }_{{\mathbf{j}}^{\prime }}^{-}\rho {\sigma }_{\mathbf{j}}^{+}-\frac{1}{2}\left\{{\sigma }_{\mathbf{j}}^{+}{\sigma }_{{\mathbf{j}}^{\prime }}^{-},\rho \right\}\right).$$

$$\rho (t)=e^{-{\gamma }_{n}t}|{\psi }_n⟩⟨{\psi }_n|+{\rho }^{\prime }(t),$$

$$\rho (t)={\rho }_{N_{\text{exc}}}(t)+{\rho }^{\prime }(t),$$

$$\frac{d{\rho }_{N_{\text{exc}}}}{dt}=-i{H}_{\text{dip}}{\rho }_{N_{\text{exc}}}+i{\rho }_{N_{\text{exc}}}{H}_{\text{dip}}^{†}.$$

$$\frac{e^{i{k}_{0}r}}{r}=\frac{i}{2\pi }\int d\mathbf{q}\frac{e^{i\mathbf{q}\cdot {\mathbf{r}}_{\perp }}{e}^{i{q}_z|z|}}{q_z},$$

$$\hat{\mathbf{G}}(\mathbf{r})=\frac{3}{4\pi k_0^3}\int d\mathbf{q}\left[k_0^{2}1-\overline{\mathbf{Q}}\otimes \overline{\mathbf{Q}}\right]\frac{e^{i\mathbf{q}\cdot {\mathbf{r}}_{\perp }}{e}^{i{q}_z|z|}}{q_z},$$

$$(c_n{)}_{\mathbf{j}}=e^{i{\delta }_{\perp }{\mathbf{j}}_{\perp }\cdot {\mathbf{q}}_n}{e}^{i\pi \frac{p_n-1}{2}(j_z-1)}/\sqrt{2N},$$

$$\sum _{{\mathbf{j}}_{\perp }}{e}^{i{\delta }_{\perp }{\mathbf{j}}_{\perp }\cdot \mathbf{q}}={\left(\frac{2\pi }{{\delta }_{\perp }}\right)}^2\sum _{\mathbf{g}}{\delta }^{(2)}(\mathbf{q}-\mathbf{g}),$$

$$\sum _{{\mathbf{j}}^{\prime }}{\mathcal{H}}_{\mathbf{j},{\mathbf{j}}^{\prime }}(c_n{)}_{{\mathbf{j}}^{\prime }}$$

$$=-\frac{3i{\gamma }_e}{8\pi k_0^3}\sum _{{\mathbf{j}}_{\perp }^{\prime },j_z^{\prime }}\int d\mathbf{q}{e}^{i{\delta }_{\perp }{\mathbf{j}}_{\perp }\cdot \mathbf{q}}\left[k_0^2-|\mathbf{q}\cdot \mathbf{p}{|}^2\right]\frac{e^{i{\delta }_{\perp }{\mathbf{j}}_{\perp }^{\prime }\cdot ({\mathbf{q}}_n-\mathbf{q})}{e}^{i{q}_z|z_{\mathbf{j}}-z_{{\mathbf{j}}^{\prime }}|}}{q_z\sqrt{2N}}{e}^{i\pi \frac{p_n-1}{2}(j_z^{\prime }-1)}$$

$$=-(c_n{)}_{\mathbf{j}}\frac{3(2\pi {)}^{2}i{\gamma }_e}{8\pi k_0(k_0{\delta }_{\perp }{)}^2}\sum _{j_z^{\prime }}\sum _{\mathbf{g}}\left[k_0^2-|({\mathbf{q}}_n-\mathbf{g})\cdot \mathbf{p}{|}^2\right]\frac{e^{i\sqrt{k_0^2-|{\mathbf{q}}_n-\mathbf{g}{|}^2}|z_{\mathbf{j}}-z_{{\mathbf{j}}^{\prime }}|}}{\sqrt{k_0^2-|{\mathbf{q}}_n-\mathbf{g}{|}^2}}{e}^{i\pi \frac{p_n-1}{2}(j_z^{\prime }-j_z)}$$

$$=-(c_n{)}_{\mathbf{j}}i(\Gamma /2)\sum _{\mathbf{g}}\frac{\left[k_0^2-|({\mathbf{q}}_n-\mathbf{g})\cdot \mathbf{p}{|}^2\right]}{k_0\sqrt{k_0^2-|{\mathbf{q}}_n-\mathbf{g}{|}^2}}\left(1+p_n{e}^{i\sqrt{k_0^2-|{\mathbf{q}}_n-\mathbf{g}{|}^2}L}\right),$$

$${\gamma }_n=\Gamma \overline{\sum _{\mathbf{g}}}$$

$$\left[1+p_n\cos \left(\sqrt{k_0^2-|{\mathbf{q}}_n-\mathbf{g}{|}^2}L\right)\right]$$

$$\overline{{\Delta }_n}=p_n(\Gamma /2)\sum _{\mathbf{g}}$$

$$\sin \left(\sqrt{k_0^2-|{\mathbf{q}}_n-\mathbf{g}{|}^2}L\right),$$

$${\text{TEM}}_{j,k}(\mathbf{r})=$$

$$e^{-(x^2+y^2)/w(z{)}^2}{e}^{i\left(k_0(x^2+y^2)/[2R(z)]-{\psi }_{j,k}(z)\right)},$$

$$w(z)=w_0\sqrt{1+(z/z_R{)}^2}$$

$$R(z)=z\left[1+(z_R/z{)}^2\right]$$

$${\psi }_{j,k}(z)=(j+k+1){\text{tan}}^{-1}(z/zR)$$

Full text

Introduction. — Recent advances in preparing regular arrays of atoms with optical traps Lester et al. (2015); Xia et al. (2015); Endres et al. (2016); Barredo et al. (2016) offer new opportunities to engineer strong collective coupling between atoms and light, with applications in quantum information science. In particular, a single layer of atoms loaded into a regular 2D array with sub-wavelength spacing has been proposed as an atomic mirror with high reflectivity García de Abajo (2007); Jenkins and Ruostekoski (2013); Bettles et al. (2016); Jenkins et al. (2017); Shahmoon et al. (2017), as quantum memory with efficient storage and retrieval Manzoni et al. (2018), and to implement topological quantum optics Perczel et al. (2017); Bettles et al. (2017); in addition, emission of single photons from bilayer atomic arrays can be engineered to be highly directional in free-space Grankin et al. (2018). Moreover, single-layered atomic arrays have been shown to support subradiant collective excitations Facchinetti et al. (2016); Plankensteiner et al. (2017); Asenjo-Garcia et al. (2017), which consist of excited superposition states of atoms decaying much slower than a single isolated excited atom, due to interference in spontaneous emission Dicke (1954); Scully (2015); Guerin et al. (2016); Sutherland and Robicheaux (2016); Jen et al. (2016); Solano et al. (2017); Moreno-Cardoner et al. .

Here we show that the composite quantum system consisting of two distant single-layered arrays of atoms [cf. Figs. 1(a-c)] can support an atomic Bell superposition state exhibiting strong subradiance. Remarkably, this non-radiating ‘dark’ state is a non-local entangled state, i.e. a superposition state of a collective excitation living in the first or second array, where the two arrays can be separated by a distance $L$ much larger than the transverse size $L_{\perp}$ of each individual array. This phenomenon relies on two ingredients. First, spontaneous emission from a collective atomic excitation in a single layer can be directional, with a proper phasing of the atomic dipoles, corresponding to light emission in both directions perpendicular to the atomic array, as in Fig. 1(a) García de Abajo (2007). Second, radiation from two distant atomic arrays can – provided the separation length $L$ is commensurate with half the optical wavelength [upper panel in Fig. 1(c)] – lead to destructive interference of light emitted to the left and to the right of the two arrays, corresponding to a subradiant state, i.e. this ‘dark’ state will show strongly suppressed radiative loss to the outside world. In contrast, the lower panel in Fig. 1(c) displays a ‘bright’ (i.e., radiating) state due to constructive interference.

Below we will show that these non-local subradiant atomic superposition states can be prepared naturally in setups involving two – or more – atomic arrays, and provide a source of entanglement shared between the two atomic arrays, with applications for quantum networking Northup and Blatt (2014). In particular, quantum information can be exchanged between the arrays representing ‘local’ quantum memories, in a coherent and deterministic process, with dark states acting as mediators.

Quantum optical model. — Our setup consists of two 2D arrays of $N=N_{\perp}\times N_{\perp}$ atomic emitters with lattice spacing $\delta_{\perp}$ and size $L_{\perp}\equiv N_{\perp}\delta_{\perp}$, separated by a distance $L$ along $z$. Each atom has a ground and an excited state, $\ket{g}_{\boldsymbol{j}}$ and $\ket{e}_{\boldsymbol{j}}$, and is coupled to free-space modes of the radiation field via a dipole transition with frequency $\omega_{0}=ck_{0}=2\pi c/\lambda_{0}$. Here the multi-index ${\boldsymbol{j}}=({\boldsymbol{j}}_{\perp},j_{z})$, where $j_{z}=1,2$ labels the arrays, while ${\boldsymbol{j}}_{\perp}=(j_{x},j_{y})$ label the atoms within each array, with $1\leq j_{x},j_{y}\leq N_{\perp}$. Atomic positions are denoted by ${\boldsymbol{r}}_{\boldsymbol{j}}=(x_{\boldsymbol{j}},y_{\boldsymbol{j}},z_{\boldsymbol{j}})$. We start by studying the dynamics of a single excitation with wave function $\ket{\psi(t)}=\sum_{{\boldsymbol{j}}}c_{\boldsymbol{j}}(t)\sigma_{\boldsymbol{j}}^{+}\ket{\mathcal{G}}\ket{0}+\int d{\boldsymbol{k}}\sum_{\lambda}\psi_{\lambda}({\boldsymbol{k}},t)\ket{\mathcal{G}}\ket{{\boldsymbol{k}},\lambda}$. Here $\sigma_{\boldsymbol{j}}^{+}=\ket{e}_{\boldsymbol{j}}\!\bra{g}$, $\ket{\mathcal{G}}=\otimes_{\boldsymbol{j}}\ket{g}_{\boldsymbol{j}}$, $\ket{0}$ is the photonic vacuum state and $\ket{{\boldsymbol{k}},\lambda}$ the state with a single photon with wave vector ${\boldsymbol{k}}$ and polarization $\lambda$. We extend our results below to states with multiple excitations.

The atomic dynamics, due to successive photon emissions and reabsorptions, is obtained by integrating out the dynamics of the radiation modes $\psi_{\lambda}({\boldsymbol{k}},t)$ in a Born-Markov approximation. Assuming the field initially in the vacuum state $\psi_{\lambda}({\boldsymbol{k}},0)=0$, this yields $\dot{c}_{\boldsymbol{j}}=-i\sum_{{\boldsymbol{j}}^{\prime}}\mathcal{H}_{{\boldsymbol{j}},{\boldsymbol{j}}^{\prime}}c_{{\boldsymbol{j}}^{\prime}}$, where in a frame rotating with $\omega_{0}$ Lehmberg (1970); James (1993); Novotny and Hecht (2006),

[TABLE]

is a non-hermitian effective Hamiltonian, whose hermitian part describes coherent exchanges of atomic excitations, while the non-hermitian part corresponds to dissipation accounting for radiation of photons. Here $\gamma_{e}$ is the spontaneous decay rate of each atom, and the dyadic Green’s tensor ${\hat{\boldsymbol{G}}}({{\boldsymbol{r}}})$, representing the electric field at position ${\boldsymbol{r}}$ generated by a dipole located at the origin, is the solution of $\boldsymbol{\nabla}\times\boldsymbol{\nabla}\times{\hat{\boldsymbol{G}}}({{\boldsymbol{r}}})-k_{0}^{2}{\hat{\boldsymbol{G}}}({{\boldsymbol{r}}})+(6\pi i/k_{0})\delta(\boldsymbol{r})=0$ with \hat{\boldsymbol{G}}\big{(}{\boldsymbol{0}}\big{)}\equiv{\boldsymbol{\mathbb{1}}} accounting for independent single-atom decay (see details in SM ). The atomic transition polarization ${\boldsymbol{p}}$ is taken circular, with $z$ as quantization axis.

Dark and bright eigenstates. — The dynamics of atomic excitations, including their radiative properties, can be understood by studying the spectrum of $\mathcal{H}$. Denoting its eigenvalues as $\epsilon_{n}=\Delta_{n}-i\gamma_{n}/2$ (with $n=1,\ldots,2N$), $\Delta_{n}$ is interpreted as the self-energy of the collective atomic excitation given by the corresponding eigenstate $c_{n}$, while $\gamma_{n}$ is its spontaneous emission rate. In particular, an eigenstate is subradiant (or ‘dark’) if spontaneous emission occurs with a rate suppressed below the single-atom decay rate $\gamma_{e}$. In view of the mirror symmetry of the system, all eigenstates have a definite parity, i.e., they can be written as $(c_{n})_{({\boldsymbol{j}}_{\perp},1)}=p_{n}(c_{n})_{({\boldsymbol{j}}_{\perp},2)}\equiv(v_{n})_{{\boldsymbol{j}}_{\perp}}/\sqrt{2}$, with parity $p_{n}=\pm 1$.

In Fig. 1(d) we plot the decay rates [for the setup of Fig. 1(c)], with parameters chosen as explained below. One of the eigenstates is remarkably subradiant, with a decay rate of $\gamma_{d}\sim 10^{-3}\gamma_{e}$. We also represent the mean absolute value of the transverse quasi-momentum $\overline{q}$, which is obtained from the discrete Fourier transform of the corresponding eigenvectors $(\tilde{v}_{n})_{{\boldsymbol{q}}}=\sum_{{\boldsymbol{j}}_{\perp}}(v_{n})_{{\boldsymbol{j}}_{\perp}}e^{i\delta_{\perp}{\boldsymbol{j}}_{\perp}\cdot{\boldsymbol{q}}}/\sqrt{N}$ as $\overline{q}=\sum_{\boldsymbol{q}}|(\tilde{v}_{n})_{{\boldsymbol{q}}}|^{2}|{\boldsymbol{q}}|$, with discrete quasi-momentum $\boldsymbol{q}=(q_{x},q_{y})$ where $q_{x,y}=-\pi/\delta_{\perp}+2\pi n_{x,y}/L_{\perp}$ ($n_{x,y}=0,1,...,N_{\perp}-1$). Two states have a distinctly low quasi-momentum $\overline{q}\ll k_{0}$: the dark state, as well as a ‘bright’ state, which radiates photons with a rate $\gamma_{b}$ comparable to $\gamma_{e}$. We contrast our dark states with the $\overline{q}>k_{0}$ subradiant states in single layer setups, studied e.g. in Refs. Asenjo-Garcia et al. (2017); Sutherland and Robicheaux (2016). In Fig. 1(e) we show the probability amplitude of the eigenvectors $|(v_{n})_{{\boldsymbol{j}}_{\perp}}|$ for the two states with lowest decay rates SM .

This pair of dark and bright states can be understood by considering first the situation where the arrays are infinite ($N_{\perp}\to\infty$), and the eigenstates are plane waves $(v_{n})_{{\boldsymbol{j}}_{\perp}}=e^{i\delta_{\perp}{\boldsymbol{j}}_{\perp}\cdot{\boldsymbol{q}}_{n}}/\sqrt{N}$ with continuous quasi-momentum ${\boldsymbol{q}}_{n}$. We now make two assumptions: First, the lattice spacing satisfies $\delta_{\perp}<\lambda_{0}$. Under this condition, we obtain, provided $|{\boldsymbol{q}}_{n}|\leq 2\pi/\delta_{\perp}-k_{0}$ SM ,

[TABLE]

with $\Gamma=3\pi\gamma_{e}/(k_{0}\delta_{\perp})^{2}$, $k_{z}=\sqrt{k_{0}^{2}-|{\boldsymbol{q}}_{n}|^{2}}$. Considering in particular the symmetric ($p_{n}=1$) and antisymmetric ($p_{n}=-1$) eigenstates with ${\boldsymbol{q}}_{n}={\boldsymbol{0}}$, we obtain a pair of states with decay rates $\gamma_{s/a}=\Gamma[1\pm\cos(k_{0}L)]$. Similarly, their self-energies are $\Delta_{s/a}=\pm(\Gamma/2)\sin(k_{0}L)+\Delta_{d}$, as depicted in Fig. 2(b), where $\Delta_{d}$ is a collective Lamb shift evaluated numerically. Our second assumption is that $k_{0}L=m\pi$ with integer $m$, so that either $\gamma_{s}$ or $\gamma_{a}$ vanishes due to interference in the emission of the two arrays, while the other reduces to $\gamma_{b}=2\Gamma$. The corresponding Bell states

[TABLE]

are thus respectively ‘dark’ and ‘bright’.

For finite-sized arrays, the eigenstates $(v_{n})_{{\boldsymbol{j}}_{\perp}}$ are confined, which has two consequences yielding a finite decay rate $\gamma_{d}$ for the dark state. First, photon emission in transverse directions is not perfectly cancelled. Second, photons emitted along $z$ have a finite spread of transverse momentum, and thus diffract when propagating between the two arrays, thereby hindering the interference of emission. This can be mitigated by curving the arrays according to the phase profile of a Gaussian mode $\mathcal{E}({\boldsymbol{r}})$ propagating along $z$ [as shown in Fig. 1(c)], in analogy to the mirrors of an optical cavity. As represented in Fig. 1(e), the spatial distribution of the dark state (as well as the bright state) is then $(v_{d})_{{\boldsymbol{j}}_{\perp}}\propto\mathcal{E}({\boldsymbol{r}}_{({\boldsymbol{j}}_{\perp},1)})$ par . Alternatively, one can add optical elements between the arrays, such as lenses or fibers.

The spatial profile of the electric field, generated by (virtual) photon exchanges between the atomic dipoles in the dark state, reads ${\boldsymbol{\psi}}({\boldsymbol{r}})\sim\sum_{\boldsymbol{j}}c_{\boldsymbol{j}}\hat{\boldsymbol{G}}({\boldsymbol{r}}-{\boldsymbol{r}}_{\boldsymbol{j}})\cdot{\boldsymbol{p}}$, and forms a standing wave [see Fig. 1(c)]. We emphasize that – although the system resembles a cavity with each array acting as a mirror – we are interested here in the quantum state of the atoms. More precisely, the ratio of atomic to photonic excitations in the dark state is given by $\Gamma L/(2c)$ with speed of light $c$ SM , which is assumed negligible when integrating the field dynamics above, amounting to neglecting retardation effects in the atomic dynamics. This is in analogy to atomic cavities built from strings of atoms coupled to a 1D waveguide Chang et al. (2012); Guimond et al. (2016).

We now discuss how the geometric parameters ($N_{\perp},L,\delta_{\perp}$) affect the spectral properties of the system. In Fig. 2(a) we show the scaling of $\gamma_{d}/\gamma_{b}$ as the relevant figure of merit, with the waist of $\mathcal{E}({\boldsymbol{r}})$ minimizing this ratio. Low ratios can be achieved for $L\lesssim L_{\perp}^{2}/\lambda_{0}$, a condition set by the diffraction limit, i.e. the spot size of the Gaussian mode must be smaller than the surface of the arrays. Remarkably, this condition allows to achieve strong subradiance even when the characteristic size of each array $L_{\perp}$ is much smaller than their separation $L$, i.e., the subradiant state is ‘non-local’. As an example, for $N_{\perp}=20$ and $\delta_{\perp}=0.8\lambda_{0}$ (i.e., $L_{\perp}=16\lambda_{0}$), we obtain $\gamma_{d}/\gamma_{b}\sim 10^{-2}$ for $L\sim 130\lambda_{0}$.

In Fig. 2(b) we observe that the interference mechanism is quite sensitive to the separation between arrays, as small deviations of $L$ compared to $\lambda_{0}$ will greatly increase the decay rate $\gamma_{d}$ [see Eq. (2)]. In Fig. 2(c) we show the effect of the lattice spacing on the saturation value of Fig. 2(a) for small $L$. The ratio of dark to bright state decay rates is minimal for $\delta_{\perp}=\lambda_{0}/2$, for which the emission in transverse directions is best cancelled, and scales with the atom number as $\gamma_{d}/\gamma_{b}\sim 1/N_{\perp}^{4}$. The collective shift $\Delta_{d}$ on the other hand is typically of the order of $\gamma_{e}$ [c.f. Fig. 2(d)]. It can be positive or negative depending on $\delta_{\perp}$, and vanishes around $\delta_{\perp}=0.2\lambda_{0}$ and $\delta_{\perp}=0.8\lambda_{0}$ (see also Ref. Shahmoon et al. (2017)).

Dark state preparation and quantum state transfer. — In order to prepare the atoms in the dark state, we consider the setup represented in Fig. 3(a), where the atomic level structure now includes a third state $\ket{s}$. We assume that the system is initially in a superposition state of the first array $S_{1}^{+}\ket{\mathcal{G}}$, with $S_{1}^{+}=\sum_{{\boldsymbol{j}}_{\perp}}(v_{d})_{{\boldsymbol{j}}_{\perp}}\ket{s}_{({\boldsymbol{j}}_{\perp},1)}\!\bra{g}$. This could be realized for instance using laser-dressed Rydberg-Rydberg interactions Grankin et al. (2018); Petrosyan and Mølmer (2018), or single photon pulses SM . Moreover, we assume a coherent field drives the $\ket{s}\to\ket{e}$ transition in the first array with Rabi frequency $\Omega$, resonantly with the collective shift $\Delta_{d}$. The atoms are thus driven from state $S_{1}^{+}\ket{\mathcal{G}}$ to a superposition of dark and bright states $\sum_{{\boldsymbol{j}}_{\perp}}(v_{d})_{{\boldsymbol{j}}_{\perp}}\sigma_{({\boldsymbol{j}}_{\perp},1)}^{+}\ket{\mathcal{G}}=(1/\sqrt{2})(\sigma_{b}^{+}+\sigma_{d}^{+})\ket{\mathcal{G}}$, where the operators $\sigma^{+}_{d/b}\equiv\sum_{{\boldsymbol{j}}_{\perp}}(v_{d})_{{\boldsymbol{j}}_{\perp}}(\sigma^{+}_{({\boldsymbol{j}}_{\perp},1)}\mp[-1]^{m}\sigma^{+}_{({\boldsymbol{j}}_{\perp},2)})/\sqrt{2}$ create a dark/bright atomic excitation. The decay rate of bright excitations can be orders of magnitude larger than for dark excitations, such that their contribution to the dynamics is vastly different. If $\gamma_{b}\gg\Omega$, the bright mode can be adiabatically eliminated, and contributes an effective loss with rate $\Omega^{2}/\gamma_{b}$, which can vanish in the spirit of a quantum Zeno effect. On the other hand, if $\Omega\gg\gamma_{d}$, the dynamics will yield oscillations between the initial state and the non-local dark state.

This mechanism can be exploited for quantum state transfer between the two arrays. Here, an initial qubit superposition state in the first array $\ket{\psi_{i}}=c_{g}\ket{\mathcal{G}}+c_{s}S_{1}^{+}\ket{\mathcal{G}}$ (with $|c_{g}|^{2}+|c_{s}|^{2}=1$) is transferred deterministically to the second array. That is, we realize the process $\ket{\psi_{i}}\to\ket{\psi_{f}}=c_{g}\ket{\mathcal{G}}+c_{s}S_{2}^{+}\ket{\mathcal{G}}$, where $S_{2}^{+}=\sum_{{\boldsymbol{j}}_{\perp}}(v_{d})_{{\boldsymbol{j}}_{\perp}}\ket{s}_{({\boldsymbol{j}}_{\perp},2)}\!\bra{g}$, with high fidelity $\mathcal{F}\approx 1$ foo . By driving atoms in both arrays with Rabi frequency $\Omega$, the state $S_{2}^{+}\ket{\mathcal{G}}$ is coupled to the opposite superposition $\sum_{{\boldsymbol{j}}_{\perp}}(v_{d})_{{\boldsymbol{j}}_{\perp}}\sigma_{({\boldsymbol{j}}_{\perp},2)}^{+}\ket{\mathcal{G}}=(1/\sqrt{2})(\sigma_{b}^{+}-\sigma_{d}^{+})\ket{\mathcal{G}}$, and we can write an effective model, where the system is described by four excitation modes: two ‘local’ modes, with creation operators $S_{1}^{+}$ and $S_{2}^{+}$, which represent quantum memories in $\ket{\psi_{i}}$ and $\ket{\psi_{f}}$; and two ‘non-local’ bright and dark modes, with creation operators $\sigma_{b}^{+}$ and $\sigma_{d}^{+}$, connecting the two memories. The dynamics can then be described by a Lindblad master equation for the density matrix of the atoms $\rho$, as $\dot{\rho}=-i[H_{\text{eff}},\rho]+\gamma_{d}\mathcal{D}[\sigma_{d}^{-}]\rho+\gamma_{b}\mathcal{D}[\sigma_{b}^{-}]\rho$, where $\mathcal{D}[a]\rho\equiv a\rho a^{\dagger}-(1/2)(a^{\dagger}a\rho+\rho a^{\dagger}a)$, and with an effective Hamiltonian

[TABLE]

The evolution of the system is shown in Fig. 3(b), demonstrating transfer at time $t=\pi/\Omega$ fn . We emphasize that our protocol does not require tailoring the temporal shape of exchanged photons, in contrast to deterministic quantum state transfer protocols with ‘flying’ photonic qubits Cirac et al. (1997); Grankin et al. (2018). Fig. 3(c) represents in red the optimal achievable fidelity for given $\gamma_{d,b}$, which reads

[TABLE]

showing the requirement $\gamma_{b}\gg\gamma_{d}$. The blue dots represent simulations for atomic arrays with the parameters of Fig. 2(a), with the optimal drive given by $\Omega=\sqrt{\gamma_{d}\gamma_{b}/8}$ SM . As noted above, our treatment neglects effects of retardation in atomic dynamics; Eq. (5) remains, however, valid even for large delay times, although at the cost of a slowdown of the dynamics SM .

Probing the dark state. — The existence of the dark state can be detected in the reflection of an external laser (see details in SM ). We consider here a weak probing field with frequency $\omega_{0}+\Delta_{d}$, propagating along $z$ in the Gaussian mode $\mathcal{E}({\boldsymbol{r}})$, and driving atoms prepared in the ground state $\ket{\mathcal{G}}$. Assuming the transition frequency of the atoms in each array is additionally detuned, by $\Delta$ for atoms in the first array and either $\Delta$ or $-\Delta$ for the second array, the dark and bright states are then revealed in the width of the resonance peak of the reflectivity $R(\Delta)$. We obtain $R=(\gamma_{b}-\gamma_{d})^{2}/(\gamma_{b}^{2}+4\Delta^{2})$ for symmetric detuning, and $R=(\gamma_{b}-\gamma_{d})^{2}/(\gamma_{b}+4\Delta^{2}/\gamma_{d})^{2}$ for opposite detuning, which both have a peak at $\Delta=0$ SM ; the widths of these peaks are given by $\gamma_{b}$ and $\sim\!\sqrt{\gamma_{d}\gamma_{b}}$, respectively, allowing for a direct probing of the dark state lifetime.

Experimental considerations. — The level structure can be implemented in neutral atoms using for instance stretched states of 87Rb for $\ket{g}=\ket{5S_{1/2},F=2,m_{F}=2}$ and $\ket{e}=\ket{5P_{3/2},F=3,m_{F}=3}$, along with a strong magnetic field to eliminate other hyperfine states from the dynamics. The level $\ket{s}$ needs to be coherently coupled to the excited state, while avoiding spontaneous decay from $\ket{e}$ to $\ket{s}$. This could be realized for example using a Rydberg state $\ket{s}=\ket{nS_{1/2},m=1/2}$, with higher energy Manzoni et al. (2018), or another ground state $\ket{s}=\ket{5S_{1/2},F=1,m_{F}=1}$, coupled to $\ket{e}$ via a two-photon transition Porras and Cirac (2008). Alternatively, one can use for the optical transition atoms with a $J=0\to J=1$ transition, e.g. 88Sr; while this introduces three excited states with orthogonal dipole matrix elements, our results for dark and bright state decay rates remain qualitatively similar SM .

The atomic trap is characterized by a finite temperature and Lamb-Dicke parameter $\eta$ Ludlow et al. (2015). The resulting spread of the atomic wavefunction yields a renormalization of the decay rates as $\gamma_{d/b}\to\gamma_{d/b}[1-\eta^{2}(2n_{\text{th}}+1)]+\gamma_{e}\eta^{2}(2n_{\text{th}}+1)$ SM , where $n_{\text{th}}$ is the thermal occupation number of trap states, and we assumed $\eta\sqrt{2n_{\text{th}}+1}\ll 1$ and $\gamma_{e}\eta\sqrt{2n_{\text{th}}+1}\ll\omega_{\nu}$, with $\omega_{\nu}$ the atomic motional frequency. We thus need $\eta^{2}(2n_{\text{th}}+1)\lesssim\gamma_{d}/\gamma_{e}$. The effect of missing atoms is similar SM ; for a defect probability $p$, we find $\gamma_{d/b}\to\gamma_{d/b}(1-p)+\gamma_{e}p+\mathcal{O}(p^{2})$, i.e. we require $p\lesssim\gamma_{d}/\gamma_{e}$.

Multiple excitations. — For states with multiple excitations, the dynamics can be studied again by analyzing the spectral properties of the non-hermitian effective Hamiltonian, which now takes the form $H_{\text{dip}}=\sum_{{\boldsymbol{j}},{\boldsymbol{j}}^{\prime}}\mathcal{H}_{{\boldsymbol{j}},{\boldsymbol{j}}^{\prime}}\sigma_{\boldsymbol{j}}^{+}\sigma_{{\boldsymbol{j}}^{\prime}}^{-}$ SM . Since each atom cannot be excited more than once, the doubly-excited state $(\sigma_{d}^{+})^{2}\ket{\mathcal{G}}$ cannot be an exact eigenstate of $H_{\text{dip}}$. An analytical expression for the resulting decay rates can, however, be obtained by treating the non-linearity as perturbation, where each excitation effectively acts as a defect for the other, with the ‘defect’ probability $p$ identified as the inverse participation ratio $p=\sum_{{\boldsymbol{j}}_{\perp}}|(v_{d})_{{\boldsymbol{j}}_{\perp}}|^{4}$ (see Ref. SM ). In Fig. 2(a) we show in red, for the eigenstate closest to $(\sigma_{d}^{+})^{2}\ket{\mathcal{G}}$, the ratio of the decay rate per excitation $\gamma_{d}^{(2)}$ and $\gamma_{b}$, which is well captured by this analytical approximation (dashed red curves).

For large $N_{\perp}$, we thus expect $\gamma_{d}^{(2)}\sim\gamma_{e}/N_{\perp}^{2}$, since $(v_{d})_{{\boldsymbol{j}}_{\perp}}\sim 1/N_{\perp}$. Two regimes can then be explored. First, for $\gamma_{d},\gamma_{d}^{(2)}\ll\gamma_{b}$ the system becomes effectively almost linear, and in particular the protocol for quantum state transfer above remains valid, with the replacement $\gamma_{d}\to\gamma_{d}^{(2)}$. This can be used to transfer states with more than one excitation, e.g. quantum error correcting states such as cat or binomial states Michael et al. (2016), allowing in principle to reach fidelities beyond Eq. (5). Second, if $\gamma_{d}\ll\gamma_{d}^{(2)},\gamma_{b}$, excitations of radiating two-excitation states can be adiabatically eliminated, exploiting again the quantum Zeno effect. This mechanism can be used to effectively block the transfer from the memories to the dark state, and thereby can operate as a controlled-phase gate Dzsotjan et al. (2010). Moreover, by the same principle, weakly driving the optical transition of atoms in one of the arrays generates Rabi oscillations between $\ket{\mathcal{G}}$ and $\sigma_{d}^{+}\ket{\mathcal{G}}$ as a two-level system, which can also be used to prepare the system in the dark state, e.g. for entanglement generation between memories, or as single-photon source.

Conclusion. — We have shown that distant single-layered arrays of two-level atoms can support subradiant (long-lived) states as collective excitations in the form of Bell superpositions. Our setup constitutes a building block for a modular quantum architecture, where quantum information, stored and processed in atomic arrays, is exchanged via dark modes. Moreover, the separation between arrays can be drastically increased by adding lenses or optical fibers to mediate photons between the arrays, although at the cost of adding decoherence channels. While we discussed here implementations with atoms in optical lattices, our results remain valid for other types of emitters, including for instance in solid-state platforms such as color centers in diamond Doherty et al. (2013), quantum dots Lodahl et al. (2015), or monolayers of transition metal dichalcogenides Zhou et al. .

I Quantum optical model

Here we provide details on our model and the definitions in Eq. (1). We first consider the full system comprising the atoms in the first array (labeled with $j_{z}=1$), in the second array ($j_{z}=2$), and including the electromagnetic field. Each atom has a ground state $\ket{g}$ and an excited state $\ket{e}$. The dynamics is governed by the Hamiltonian $H_{\text{tot}}=H_{a}+H_{f}+V_{af},$ where $H_{a}$ acts on the atoms, and reads ($\hbar=1$)

[TABLE]

with $\boldsymbol{j}=(\boldsymbol{j}_{\perp},j_{z})$, $\boldsymbol{j}_{\perp}=(j_{x},j_{y})$ and $1\leq j_{x},j_{y}\leq N_{\perp}$. Here $\omega_{0}$ is the atomic transition frequency, $\sigma_{\boldsymbol{j}}^{+}\equiv\ket{e}_{\boldsymbol{j}}\bra{g}$, and $H_{\text{drive}}$ is an additional term accounting for possible additional laser drivings. The electromagnetic field Hamiltonian reads $H_{f}=\int d{\boldsymbol{k}}\sum_{\lambda}\omega_{{\boldsymbol{k}}}b^{\dagger}_{\lambda,{\boldsymbol{k}}}b_{\lambda,{\boldsymbol{k}}}$, where $\omega_{{\boldsymbol{k}}}=c|{\boldsymbol{k}}|$ with $c$ the speed of light, $b_{\lambda,{\boldsymbol{k}}}$ is the annihilation operator for photons with helicity $\lambda=\pm 1$ satisfying $[b_{\lambda,{\boldsymbol{k}}},b^{\dagger}_{\lambda^{\prime},{\boldsymbol{k}}^{\prime}}]=\delta_{\lambda,\lambda^{\prime}}\delta({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime})$. Finally, the interaction Hamiltonian reads

[TABLE]

where $d$ is the atomic dipole, ${\boldsymbol{p}}$ the atomic transition polarization (we assume circular polarization), and the electric field operator expresses as

[TABLE]

with ${\boldsymbol{e}}_{\lambda,{\boldsymbol{k}}}$ the polarization unit vector, $\epsilon_{k}=\sqrt{ck/(2[2\pi]^{3}\varepsilon_{0})}$, and $\varepsilon_{0}$ the vacuum permittivity.

Assuming the electromagnetic field is initially in the vacuum state, the field dynamics can be integrated, to obtain a Lindblad master equation for the reduced system of the atoms, within a Born-Markov approximation. We obtain Lehmberg (1970); Novotny and Hecht (2006)

[TABLE]

The resulting non-hermitian dipole-dipole interaction Hamiltonian reads $H_{\text{dip}}=\sum_{{\boldsymbol{j}},{\boldsymbol{j}}^{\prime}}\mathcal{H}_{{\boldsymbol{j}},{\boldsymbol{j}}^{\prime}}\sigma_{\boldsymbol{j}}^{+}\sigma_{{\boldsymbol{j}}^{\prime}}^{-}$, where

[TABLE]

with $\gamma_{e}=k_{0}^{3}d^{2}/(3\pi\varepsilon_{0})$ the single-atom decay rate, $G({\boldsymbol{r}})={\boldsymbol{p}}^{*}\cdot\hat{\boldsymbol{G}}({\boldsymbol{r}})\cdot{\boldsymbol{p}}$, and the dyadic Green’s tensor taking the explicit form

[TABLE]

which represents the field at position ${\boldsymbol{r}}$ emitted by a dipole at the origin, with

[TABLE]

and we define $\hat{\boldsymbol{G}}({\boldsymbol{0}})={\boldsymbol{\mathbb{1}}}$. The last term in Eq. (6) reads

[TABLE]

In writing the master equation we moved to a rotating frame with the atomic transition frequency $\omega_{0}$, and we made the following assumptions. (i) Rotating wave approximation: counter-rotating terms (such as $\sigma_{\boldsymbol{j}}^{+}\sigma_{{\boldsymbol{j}}^{\prime}}^{+}$), which do not preserve the number of atomic excitations, are neglected. (ii) Markov approximation: retardation effects due to finite light velocity are also neglected.

We can notice from Eq. (6) that the collective emission properties of the arrays are determined by the spectrum of $H_{\text{dip}}$. Since $H_{\text{dip}}$ conserves the total number of atomic excitations $\sum_{\boldsymbol{j}}\sigma_{\boldsymbol{j}}^{+}\sigma_{\boldsymbol{j}}^{-}$, we can evaluate its eigenstates in each excitation subspace separately. In particular, for single excitation eigenstates this amounts to diagonalizing $\mathcal{H}$. Let us assume the atomic arrays are initially prepared in one of these eigenstates $\ket{\psi_{n}}$, with complex eigenvalue $\epsilon_{n}=\Delta_{n}-i\gamma_{n}/2$, containing $N_{\text{exc}}$ atomic excitations, i.e. $\rho(0)=\ket{\psi_{n}}\bra{\psi_{n}}$. Assuming here $H_{\text{drive}}=0$, the dynamics of Eq. (6) will yield

[TABLE]

where $\rho^{\prime}(t)$ contains strictly less than $N_{\text{exc}}$ atomic excitations, and $\gamma_{n}$ is thus interpreted as the decay rate of the eigenstate $\ket{\psi_{n}}$. More generally, starting from an initial mixture on the subspace with $N_{\text{exc}}$ excitations, we can write

[TABLE]

where $\rho_{N_{\text{exc}}}(t)$ is a density matrix with $N_{\text{exc}}$ excitations satisfying

[TABLE]

II Spectral analysis of $\mathcal{H}$

Here we discuss the spectrum of $\mathcal{H}$ (i.e., the spectrum of $H_{\text{dip}}$ in the single-excitation subspace), in the cases of infinite and finite arrays. We then explain how the system can be probed with a laser to measure the bright and dark state decay rates.

II.1 Infinite planar arrays

We first derive analytical expressions for the spectrum of $\mathcal{H}$ in the case of infinite planar arrays ($N_{\perp}\to\infty$). We use the identity Novotny and Hecht (2006)

[TABLE]

where ${\boldsymbol{r}}=({\boldsymbol{r}}_{\perp},z)$, and $q_{z}=\sqrt{k_{0}^{2}-|{\boldsymbol{q}}|^{2}}$. This allows us to rewrite the Green’s tensor as

[TABLE]

where ${\boldsymbol{\overline{Q}}}=({\boldsymbol{q}},q_{z}\text{sgn}(z))$.

Due to translational invariance and parity symmetry, the eigenstates of $\mathcal{H}$ are plane waves

[TABLE]

with ${\boldsymbol{q}}_{n}$ in the first Brillouin zone and $p_{n}=\pm 1$ the eigenstate parity. Next we use the relation

[TABLE]

where the sum on the right-hand side runs over vectors ${\boldsymbol{g}}$ of the reciprocal lattice, i.e., $g_{x,y}=(2\pi/\delta_{\perp})m_{x,y}$ with integer $m_{x,y}$. We thus obtain from Eq. (7)

[TABLE]

where $\Gamma=3\pi\gamma_{e}/(k_{0}\delta_{\perp})^{2}$.

From Eq. (11), a finite number of diffraction orders (i.e., vectors ${\boldsymbol{g}}$) contribute to the decay rate $\gamma_{n}$, as for $|{\boldsymbol{q}}_{n}-{\boldsymbol{g}}|\geq k_{0}$ the eigenvalue becomes purely real. We thus obtain

[TABLE]

where $\overline{\sum_{\boldsymbol{g}}}$ is restricted to vectors ${\boldsymbol{g}}$ of the reciprocal lattice satisfying $|{\boldsymbol{q}}_{n}-{\boldsymbol{g}}|<k_{0}$. In particular, for ${\boldsymbol{q}}_{n}={\boldsymbol{0}}$, a single order ($m_{x}=m_{y}=0$) contributes, provided $\delta_{\perp}<\lambda_{0}$, in which case we obtain Eq. (2). Moreover, this becomes valid for all ${\boldsymbol{q}}_{n}$ if $\delta_{\perp}<\lambda_{0}/2$.

The self-energies $\Delta_{n}$ can be similarly evaluated, however with a bit of caution. Indeed the real part in Eq. (11) diverges as all ${\boldsymbol{g}}$ with $|{\boldsymbol{q}}_{n}-{\boldsymbol{g}}|\geq k_{0}$ now contribute. We distinguish between two contributions, i.e. write $\sum_{\boldsymbol{g}}=\overline{\sum_{\boldsymbol{g}}}+\overline{\overline{\sum_{\boldsymbol{g}}}}$, where the second sum accounts for these vectors. Similarly, we write $\Delta_{n}=\overline{\Delta_{n}}+\overline{\overline{\Delta_{n}}}$ with

[TABLE]

while $\overline{\overline{\Delta_{n}}}$ is obtained numerically from the eigenvalues of $\mathcal{H}$. This last term is independent of $L$ as the exponential term in Eq. (11) vanishes and is thus identical to the self-energy of a single 2D array Shahmoon et al. (2017).

II.2 Finite-sized (curved) arrays

We now consider the case of finite atomic arrays. As discussed in the main text, in order to mitigate the spreading of wavepackets for photons propagating between the arrays, we assume the atoms in each array are located along the phase profile of a Gaussian mode $\mathcal{E}({\boldsymbol{r}})$. In the following we provide details on this Gaussian mode, derive analytical expressions for the spectrum of $\mathcal{H}$, and provide in the end a numerical study of the eigenstate distribution.

II.2.1 Definition of the Hermite-Gaussian modes and array curvature

Here we summarize the properties and notations of the Hermite-Gauss modes, which are solutions of the paraxial equation for light $[\partial_{z}-({i}/{2k_{0}})\nabla_{\perp}^{2}]\text{TEM}_{j,k}({\boldsymbol{r}})=0$ for modes propagating along $z$, with $\nabla_{\perp}^{2}\equiv\partial_{x}^{2}+\partial_{y}^{2}$. These modes represent a natural basis for the field generated by the atomic arrays. Assuming the focal point is here located at ${\boldsymbol{r}}={\boldsymbol{0}}$, these modes are defined by the waist $w_{0}$ as Novotny and Hecht (2006)

[TABLE]

with $j,k=(0,1,...)$, $H_{j}$ is the Hermite polynomial of order $j$,

[TABLE]

the mode width,

[TABLE]

the radius of curvature,

[TABLE]

the Gouy phase and $z_{R}=\pi w_{0}^{2}/\lambda_{0}$the Rayleigh length, and are normalized as

[TABLE]

In particular, the phase profile of the Gaussian mode, which determines the curvature of the atomic arrays, is taken as $\mathcal{E}({\boldsymbol{r}})\sim\text{TEM}_{0,0}({\boldsymbol{r}})e^{ik_{0}z}$. Specifically, for a given separation distance $L$ between arrays and mode waist $w_{0}$, the longitudinal position $z_{\boldsymbol{j}}$ of atom $\boldsymbol{j}$ satisfies

[TABLE]

where the $+$ and $-$ signs correspond respectively to atoms in the second array ($j_{z}=2$) and in the first array ($j_{z}=1$), such that the phase of $\mathcal{E}({\boldsymbol{r}}_{\boldsymbol{j}})$ only depends on $j_{z}$. From Eq. (15), the curvature radius for the arrays is maximal when $L\approx 2z_{R}$, yielding a displacement along $z$ of $\sim L_{\perp}^{2}/(2L)$ at the corners of the arrays, which can take values of the order of the wavelength $\lambda_{0}$. As we show in the next section, the condition of Eq. (16) allows us to construct non-local eigenstates of $\mathcal{H}$ with Gaussian distribution.

II.2.2 Analytical expressions

We now derive expressions for the eigenvalues of $\mathcal{H}$. Due to the finite array size, plane waves are not longer eigenstates of $\mathcal{H}$, however parity remains a symmetry of the system. We thus write the eigenstates as $(c_{n})_{({\boldsymbol{j}}_{\perp},1)}=p_{n}(c_{n})_{({\boldsymbol{j}}_{\perp},2)}\equiv(v_{n})_{{\boldsymbol{j}}_{\perp}}/\sqrt{2}$ with $p_{n}=\pm 1$. The matrix $\mathcal{H}$ can then be decomposed into 2 matrices of size $N$, namely

[TABLE]

accounting for the dipole-dipole interaction within each array, and

[TABLE]

accounting for the effective interaction between different arrays, with $v_{n}$ being an eigenstate of $\mathcal{H}_{0}+p_{n}\mathcal{H}_{1}$, with the same eigenvalue $\epsilon_{n}$.

The Green’s tensor in Eq. (7) can be formally decomposed as

[TABLE]

where $G_{\text{par}}({\boldsymbol{r}})$ is the Green’s function for paraxial modes, reading

[TABLE]

with

[TABLE]

for $z^{\prime}>z$.

We define the Fourier transform

[TABLE]

and the mean absolute quasi-momentum

[TABLE]

where $q_{x},q_{y}=-\pi/\delta_{\perp}+2\pi n_{x,y}/L_{\perp}$, with $n_{x,y}=0,1,...,N_{\perp}-1$. We consider in particular eigenstates with low quasi-momentum $\overline{q}$. As we saw in Sec. II.1, provided $\delta_{\perp}<\lambda_{0}$ only a single diffraction order contributes to the spontaneous emission, meaning that photons are emitted mostly in the direction normal to the arrays, and as such can be treated within a paraxial approximation. This motivates us to look for eigenvectors distributed as $(v_{(j,k)})_{{\boldsymbol{j}}_{\perp}}\sim\text{TEM}_{j,k}({\boldsymbol{r}}_{({\boldsymbol{j}}_{\perp},1)})e^{ik_{0}z_{({\boldsymbol{j}}_{\perp},1)}}$. From Eqs. (6), (19) and (21), we then have

[TABLE]

where $\Gamma=3\pi\gamma_{e}/(k_{0}\delta_{\perp})^{2}$, and we approximated the sum as an integral $\sum_{{\boldsymbol{j}}^{\prime}_{\perp}}\approx\int d{\boldsymbol{r}}^{\prime}_{\perp}/\delta_{\perp}^{2}$. The term $\gamma^{\prime}_{(j,k)}$ is a phenomenological decay by photon emission into non-paraxial modes, which we add as a perturbation accounting for $G^{\prime}({\boldsymbol{r}})$ in Eq. (19). In analogy to Sec. II.1, the self-energy $\overline{\overline{\Delta}}_{(j,k)}$ on the other hand diverges due to the divergence in Eq. (20) when $z\to 0$, and must be evaluated numerically. $v_{(j,k)}$ is thus approximately eigenstate of $\mathcal{H}_{0}$. Similarly, we get

[TABLE]

where the effect of $G^{\prime}({\boldsymbol{r}})$ is here neglected in a paraxial approximation for the photons exchanged between different arrays.

Using Eq. (16), we have

[TABLE]

such that $v_{(0,0)}$ is eigenstate of $\mathcal{H}_{0}\pm\mathcal{H}_{1}$ with eigenvalue

[TABLE]

This is the expression of the eigenvalues for the symmetric and anti-symmetric eigenstates of the main text, where we identify $\gamma_{d}\equiv\gamma^{\prime}_{(0,0)}$ and $\Delta_{d}\equiv\overline{\overline{\Delta}}_{0,0}$. For $k_{0}L=m\pi$ with integer $m$ one of these states is thus ‘dark’ (with minimal decay $\gamma_{d}$), while the other state is ‘bright’ (as it decays with rate $\gamma_{b}=2\Gamma+\gamma_{d}$). For $(j,k)\neq(0,0)$ on the other hand, we get similar expressions by considering that the Gouy phase $\psi_{j,k}$ in Eq. (13) is approximately constant for atoms within the same array. We then have

[TABLE]

where

[TABLE]

such that $v_{(j,k)}$ is eigenstate of $\mathcal{H}_{0}\pm\mathcal{H}_{1}$ with eigenvalue

[TABLE]

II.2.3 Numerical study

Here we provide details on the eigenstates $v_{n}$ and eigenvalues $\epsilon_{n}$ of $\mathcal{H}$ for finite arrays. The decay rates $\gamma_{d}$ and $\gamma_{b}$ are obtained by diagonalizing $\mathcal{H}$ and identifying the dark and bright states as the eigenstates with lowest quasi-momentum $\overline{q}$. We minimize the ratio $\gamma_{d}/\gamma_{b}$ by varying $w_{0}$, which sets the longitudinal atomic according to Eq. (16).

In Fig. 4(a) we show the optimal mode width $w$ [from Eq. (14)] for the parameters of Fig. 2(a) for a single excitation. At large $L\gtrsim L_{\perp}^{2}/\lambda_{0}$, we have $w=\sqrt{L\lambda_{0}/\pi}$, which is the minimal width achievable for fixed $L$ within the diffraction limit, where $z_{R}=L/2$. In this regime imperfections (i.e., finite $\gamma_{d}$) are mainly due to the array size being too small to fit a Gaussian mode connecting the arrays. At small $L\ll L_{\perp}^{2}/\lambda_{0}$, the width saturates to around $w\sim L_{\perp}/4$. This is a trade-off between having the Gaussian mode $\mathcal{E}({\boldsymbol{r}})$ fit the arrays, and increasing the number of participating atoms in order to minimize the emission to non-paraxial modes, both effects leading to a finite rate $\gamma_{d}$. In Fig. 4(b) we represent the overlap of the eigenvectors $v_{n}$ with the Gaussian mode $\mathcal{E}({\boldsymbol{r}})$. For each eigenvectors, this overlap is computed as $O_{n}=\left|\sum_{{\boldsymbol{j}}_{\perp}}\mathcal{E}({\boldsymbol{r}}_{({\boldsymbol{j}}_{\perp},1)})(v_{n})^{*}_{{\boldsymbol{j}}_{\perp}}\right|^{2}/\sum_{{\boldsymbol{j}}_{\perp}}|\mathcal{E}({\boldsymbol{r}}_{({\boldsymbol{j}}_{\perp},1)})|^{2}$. The overlaps for the dark and bright modes are represented as the upper points, and are close to $1$. Conversely, the sum of the overlaps of all other eigenstates is represented as the lower points, and takes very small values, vanishing for $L\lesssim L_{\perp}^{2}/\lambda_{0}$.

The requirement for curving the atomic arrays is studied in Fig. 4(c), where we compare the ratio of $\gamma_{d}/\gamma_{b}$ between flat arrays (for which $z_{\boldsymbol{j}}=\pm L/2$) and curved arrays (satisfying Eq. (16)). One sees that the curvature can improve this ratio by several orders of magnitude for $L\lesssim L_{\perp}^{2}/\lambda_{0}$.

In Fig. 5 we show the decay rates of all eigenstates $\gamma_{n}$ as well as their average quasi-momenta $\overline{q}$ for $N_{\perp}=10$. On the left, we represent the same situation as in Fig. 1. We note first that the dark state, labeled $1$, as well as the bright state above, have the distribution of a $\text{TEM}_{0,0}$ mode. The second most subradiant state, labeled $2$, corresponds to a $\text{TEM}_{1,1}$ mode. We note that its parity is opposite to that of $v_{1}$, which is due to the fact that here $L$ is large enough that $z_{R}=L/2$, and $\phi_{1,1}=\pi$ in Eq. (23). On the other hand, eigenvectors distributed according to $\text{TEM}_{1,0}$ and $\text{TEM}_{0,1}$ modes cannot be subradiant as they are out of phase, with $\phi_{1,0}=\phi_{0,1}=\pi/2$.

In the second column, we consider the situation with $L=2\lambda_{0}$, which corresponds to the opposite extreme regime where $z_{R}\gg L$, and $\phi_{j,k}\approx 0$, such that now the subradiant states have the same parity, and eigenvectors with $\text{TEM}_{1,0}$ and $\text{TEM}_{0,1}$ distributions can also be subradiant. The decay rate significantly increases with $\overline{q}$, which can be understood as a gradual breaking of the paraxial approximation, and can be seen from Eq. (12) as the interference between arrays becomes imperfect.

Finally, in the third column we show that subradiance can also appear in eigenvectors with large $\overline{q}$ when $\delta_{\perp}<\lambda_{0}/\sqrt{2}$. There, subradiance is due to the fact that these guided modes have their momentum larger than $k_{0}$, and as such reside outside of the light cone, and are studied e.g. in Refs. Sutherland and Robicheaux (2016); Asenjo-Garcia et al. (2017). These modes are localized in each array, and are thus degenerate in energy, in contrast to the non-local modes with low $\overline{q}$.

II.3 Dark and bright state probing

Finally, here we show how to probe the dark and bright state lifetimes in the reflectivity of a laser. The situation is represented in Fig. 6(a). We consider a weak laser with polarization ${\boldsymbol{p}}$ propagating along $z$ in the Gaussian mode $\mathcal{E}({\boldsymbol{r}})$ at frequency $\omega_{0}+\Delta_{d}$, and driving the system with the atoms in their ground state $\ket{\mathcal{G}}=\otimes_{\boldsymbol{j}}\ket{g}_{\boldsymbol{j}}$. Assuming the field is weak enough, each photon will be scattered by the system independently, and the reflectivity can be evaluated for single-photon pulses. We can thus write the state of the system for a single excitation as

[TABLE]

with here $c_{\boldsymbol{j}}(0)=0$. We assume moreover that atoms in each array $j_{z}$ additionally detuned by $\Delta_{j_{z}}$. We consider two situations, with either $\Delta_{j_{z}}=\Delta$ (i.e., with the same detuning for atoms in both arrays), or $\Delta_{z}=2\Delta(3/2-j_{z})$ (with opposite detuning between the two arrays).

As in Sec. I, the field dynamics can be integrated, yielding for the atoms

[TABLE]

where we moved to a frame rotating with $\omega_{0}+\Delta_{d}$, and defined the input field

[TABLE]

The resulting field on the other hand reads, neglecting retardation effects,

[TABLE]

For long pulses, i.e., varying over timescales much larger than the atomic response time $1/\gamma_{d}$, we can set $\dot{c}_{\boldsymbol{j}}\approx 0$ in Eq. (25), and get from Eq. (26)

[TABLE]

where $(\mathbb{\Delta}_{z})_{\boldsymbol{j},\boldsymbol{j}^{\prime}}=\Delta_{j_{z}}\delta_{\boldsymbol{j},\boldsymbol{j}^{\prime}}$.

The reflectivity $R$ is then obtained by taking $\psi^{\text{in}}({\boldsymbol{r}},t)=\mathcal{E}({\boldsymbol{r}})$ as the overlap between $\psi({\boldsymbol{r}})$ and a target mode $\psi^{\text{tar}}({\boldsymbol{r}})$

[TABLE]

where the target mode is the Gaussian mode propagating to the left, i.e., $\psi^{\text{tar}}({\boldsymbol{r}})=(\mathcal{E}({\boldsymbol{r}}))^{*}$. Within a paraxial approximation for the Gaussian mode we replace the Green’s function by its paraxial counterpart in Eq. (19), and apply Eq. (21) to obtain

[TABLE]

which can be evaluated numerically, and provides the green and red curves in Fig. 6(b). We stress that the expression of Eq. (28) is valid only within the paraxial approximation. This breaks down for configurations with $w_{0}\lesssim\lambda_{0}$, which can occur for small arrays with $L_{\perp}\lesssim 4\lambda_{0}$.

As we saw from Fig. 4(b), only the dark and bright states have significant overlap with the Gaussian distribution $\mathcal{E}({\boldsymbol{r}}_{\boldsymbol{j}})$, such that we can restrict the vector space of the matrix $\mathcal{H}+\mathbb{\Delta}_{z}-\Delta_{d}\mathbb{1}$ to these two states, and invert it on this subspace. Approximating $\sum_{\boldsymbol{j}}\approx\sum_{j_{z}}\int d{\boldsymbol{r}}_{\perp}/\delta_{\perp}^{2}$ and using Eq. (16), we obtain $R(\Delta)=(\gamma_{b}-\gamma_{d})^{2}/(\gamma_{b}^{2}+4\Delta^{2})$ for the case of symmetric detuning between the two arrays, and $R(\Delta)=(\gamma_{b}-\gamma_{d})^{2}/(\gamma_{b}+4\Delta^{2}/\gamma_{d})^{2}$ for the case of opposite detuning, with $R=1$ if $\Delta=0$. These expressions are represented in the dashed black curves of Fig. 6(b), which show excellent agreement with the numerical results.

III Photonic link between quantum memories

Here we write an effective model for the atomic dynamics, retaining four modes as expressed in Eq. (4). We then derive the expression for the fidelity of quantum state transfer in Eq. (5), and explain how to write and read from the local quantum memory states. We finally discuss how our results extend to non-markovian regimes, where retardation effects due to the finite speed of photons exchanged between arrays is no longer negligible.

III.1 Effective four mode model

We now consider each atom has a $\Lambda$ level structure, as represented in Fig. 3(a). We assume at most a single atom is in state $\ket{e}$ or $\ket{s}$ at a time, and a laser drives the $\ket{e}\to\ket{s}$ transition resonantly with the cooperative shift $\Delta_{d}$, with homogeneous Rabi frequency $\Omega$. In a rotating frame, the system is thus described by the master equation in Eq. (6), with

[TABLE]

with $s_{\boldsymbol{j}}^{+}=\ket{s}_{\boldsymbol{j}}\!\bra{g}$. We start from an initial state of the form

[TABLE]

where $S_{1}^{+}=\sum_{{\boldsymbol{j}}_{\perp}}(v_{d})_{{\boldsymbol{j}}_{\perp}}s_{(\boldsymbol{j}_{\perp},1)}^{+}$ creates an excitation of $\ket{s}$ in the first array, with $(v_{d})_{\boldsymbol{j}_{\perp}}\sim\mathcal{E}({\boldsymbol{r}}_{(\boldsymbol{j}_{\perp},1)})$ the probability amplitude corresponding to the dark (and bright) state. Similarly we define $S_{2}^{+}=\sum_{{\boldsymbol{j}}_{\perp}}(v_{d})_{{\boldsymbol{j}}_{\perp}}s_{(\boldsymbol{j}_{\perp},2)}^{+}$, which creates an excitation of $\ket{s}$ in the second array. The dynamics of Eq. (6) will then excite only the eigenstates of $\mathcal{H}$ with Gaussian distribution, i.e., the dark and bright states, as shown in Fig. 4(b). We thus obtain an effective model with only four modes: the two local modes, with creation operators $S_{1}^{+}$ and $S_{2}^{+}$, as well as the dark and bright modes, created by the operators

[TABLE]

The atomic dynamics thus follows

[TABLE]

with

[TABLE]

and $\mathcal{D}[a]\rho=a\rho a^{\dagger}-(1/2)(a^{\dagger}a\rho+\rho a^{\dagger}a)$. We note that Fig. 3 provides a numerical verification of this four modes model for the simulation of quantum state transfer, which we describe below.

III.2 Quantum state transfer

We now provide an analytical derivation of the fidelity for quantum state transfer. Assuming the system is initially prepared in the pure state $\ket{\psi_{i}}$ of Eq. (29), the atomic density matrix can be expressed as

[TABLE]

where

[TABLE]

with $c_{1}(0)=c_{s}$, $c_{2}(0)=c_{b}(0)=c_{d}(0)=P_{g}(0)=0$. We wish to transfer the quantum state $\ket{\psi_{i}}$ to the second array, i.e., have the system evolve to

[TABLE]

We define the fidelity of the state transfer as ${\cal F}\equiv\max_{t}\left|c_{2}\left(t\right)\right|^{2}$ for $c_{s}=1$. Eq. (30) then yields

[TABLE]

The general solution of Eq. (31) for $c_{2}\left(t\right)$ can be written in the form

[TABLE]

where, in the regime $\gamma_{d}\ll\Omega\ll\gamma_{b}$,

[TABLE]

with $g\equiv\sqrt{16\Omega^{2}-\gamma_{d}^{2}}$ the frequency of the oscillations. The corresponding amplitudes can be readily obtained using the initial conditions, and read

[TABLE]

Using Eqs. (32), (33) and (34), the first maximum of $c_{2}\left(t\right)$ is approximately at half the period of Rabi oscillations, i.e. $t_{\text{max}}\approx\left(\pi-\arctan\left[\gamma_{d}/g\right]\right)/g$. The fidelity $\mathcal{F}$, given by $c_{2}\left(t_{\text{max}}\right)$, depends on the drive $\Omega$. Expanding $c_{2}(t_{\text{max}})$ up to the first order in $\Omega/\gamma_{b}$ and $\gamma_{d}/\Omega$, we get

[TABLE]

The optimal drive then reads $\Omega_{\text{opt}}=\sqrt{\gamma_{b}\gamma_{d}/8}$ and the corresponding optimal fidelity of the state transfer is

[TABLE]

which is Eq. (5).

III.3 Write and read of quantum memory using single photon pulses

We now discuss how one can write and read from the quantum memories in the arrays. In particular, assuming the atoms are in state $\ket{\psi_{i}}$ as in Eq. (29) while the photonic field is in the vacuum state $\ket{0}$, we show that the atoms can be brought to their ground state $\ket{\mathcal{G}}$ while emitting a photonic qubit $c_{g}\ket{0}+c_{s}\ket{1}$, where $\ket{1}$ denotes a state with a single photon leaving the system in a well defined spatio-temporal mode, propagating in a given direction. The time-reversed process allows one to absorb a photonic qubit, thereby preparing the atoms in state $\ket{\psi_{i}}$.

We make the following two additional assumptions. First, the phase acquired by a photon propagating between the arrays $k_{0}L$ can be modified, e.g. by slightly changing the distance $L$ over a range of $\sim\lambda_{0}/2$. Second, the laser drive $\Omega$ can be turned off for the atoms in the second array. For convenience, we consider the system prepared in state $S_{1}^{+}\ket{\mathcal{G}}$, with the laser driving only the first array. The dynamics, and in particular the spatio-temporal shape of the emitted photon, can be obtained following the steps in Sec. II.3. We write here the state as

[TABLE]

with $\tilde{c}_{\boldsymbol{j}}(0)=(v_{d})_{\boldsymbol{j}_{\perp}}\delta_{j_{z},1}$ and $c_{\boldsymbol{j}}(0)=\psi_{\lambda}({\boldsymbol{k}},0)=0$. Integrating the field dynamics, we get

[TABLE]

where we moved to a frame rotating with $\omega_{0}+\Delta_{d}$. Assuming $k_{0}L\neq m\pi$ with integer $m$, neither the symmetric or anti-symmetric Gaussian states are dark, and their decay rate to paraxial modes is given by $\Gamma(1\pm\cos[k_{0}L])$ [see Eq. (22)]. Thus, provided $\Omega\ll|\Gamma(1\pm\cos[k_{0}L])|$ (ideally by setting $\cos[k_{0}L]=0$), the population of state $\ket{e}_{\boldsymbol{j}}$ can be adiabatically eliminated, i.e., we set $\dot{c}_{\boldsymbol{j}}\approx 0$ in Eq. (35), yielding

[TABLE]

Next we note that for the initial condition above, from Eq. (36) we have $\tilde{c}_{\boldsymbol{j}}(t)=\tilde{c}(t)(v_{d})_{\boldsymbol{j}_{\perp}}\delta_{j_{z},1}$. Restricting $\mathcal{H}-\Delta_{d}\mathbb{1}$ to the space spanned by the symmetric and anti-symmetric Gaussian states, we get

[TABLE]

i.e., the memory will spontaneous emit a photon with rate

[TABLE]

We remark that $\tilde{\gamma}$ is independent of $k_{0}L$ if $\gamma_{d}\approx 0$, and reduces to $\tilde{\gamma}=2\Omega^{2}/\Gamma$. The spatio-temporal shape of the outgoing photon is obtained from Eq. (26), as

[TABLE]

We note that the temporal distribution can be tailored by varying $\Omega$ in time.

The flux of photons emitted in the Gaussian mode $\mathcal{E}({\boldsymbol{r}})$ (i.e., propagating to the right) is obtained as the overlap

[TABLE]

where in the second line we replace the Green’s function by its paraxial counterpart as in Eq. (19) and used Eq. (21), and in the third line restricted $(\mathcal{H}-\Delta_{d}\mathbb{1})$ to the subspace of symmetric and anti-symmetric Gaussian states, replaced $\sum_{\boldsymbol{j}_{\perp}}\approx\int d{\boldsymbol{r}_{\perp}}/\delta_{\perp}^{2}$ and used the property of Eq. (16). We note that $P_{\rightarrow}$ vanishes for $\gamma_{d}\approx 0$. Similarly, the flux of photons emitted in the Gaussian mode $\mathcal{E}^{*}({\boldsymbol{r}})$ (i.e., propagating to the left), reads

[TABLE]

For $\gamma_{d}\approx 0$ this reduces to $P_{\leftarrow}(t)=\tilde{\gamma}|\tilde{c}(t)|^{2}$, showing that the photon is emitted in the left-propagating Gaussian mode. This allows to perform a transfer from the quantum memory state $\ket{\psi_{i}}$ to a propagating photonic qubit, as well as the time-reversed process.

III.4 Beyond the Markov approximation

Finally, we discuss the effects of time-delays in the atomic dynamics, arising from the finite propagation time of photons between arrays, which were neglected in the previous sections. These effects become relevant only when the arrays are separated by $L\gtrsim c/\Gamma\sim 10$m for $\Gamma$ in the MHz range, which can be realized by mediating photons exchanged between arrays with optical lenses or fibers. We first show how this affects the decay rates of dark and bright states, and then study its effect on the state transfer fidelity.

III.4.1 Effect on dark and bright states

Let us write the state of the system

[TABLE]

and integrate the field dynamics without neglecting retardation in photon propagation between different arrays, yielding Grankin et al. (2018), with the definitions of Eqs. (17) and (18)

[TABLE]

with $\tau=L/c$. Here we added an attenuation coefficient $\kappa$ accounting for additional decay channels induced by possible optical elements mediating the exchanged photons, but neglected any effect of dispersion. In particular, for the dark and bright state amplitudes we get

[TABLE]

Defining the Laplace transform variables $\tilde{c}(s)=\mathcal{L}[c(t)](s)$, we arrive at

[TABLE]

which cannot be analytically inverted directly. An analytical approximation can however be obtained by expanding $e^{-s\tau}=1-s\tau+\mathcal{O}(s\tau)^{2}$ to lowest order in $\Gamma\tau$, yielding

[TABLE]

From Eq. (37) increasing the retardation $\Gamma\tau$ decreases the decay of the dark state due to atomic losses with rate $\gamma_{d}$. This is a consequence of the dark state being now a superposition of field and atomic excitations, with only the atomic part decaying if $\kappa L=0$. If we include a finite attenuation $\kappa L$, the photonic component also induces losses.

The photonic field is expressed as

[TABLE]

providing for the photonic flux $E(z,t)=c\int d{\boldsymbol{r}_{\perp}}|\psi({\boldsymbol{r}},t)|^{2}$, within the paraxial approximation for the Green’s tensor,

[TABLE]

where $z_{1}=-L/2$ and $z_{2}=L/2$ denote the position of the first and second arrays along $z$. In particular, in the dark state the system can reach a quasi-equilibrium if the exponent in the decay of the dark state amplitude in Eq. (37) is much smaller than $1/\tau$, and we get (with $z_{R}\gg L$ for simplicity) $E(z,t)\approx{\Gamma}\sin(k_{0}z)^{2}|c_{d}(t)|^{2}$. While this quantity shows that the arrays continuously exchange photons at a rate $\sim\Gamma$ even when $\Gamma\tau\to 0$, the total number of photons between the arrays at any time, given by

[TABLE]

vanishes in that limit.

III.4.2 Effect on state transfer fidelity

The effect of retardation on the state transfer fidelity can be studied in a similar way. Following the same procedure, with the notations of Sec. III.2, we have

[TABLE]

with initial conditions $c_{1}(0)=1$, $c_{2}(0)=c_{d}(0)=c_{b}(0)=0$. The solution for the Laplace transform of $c_{2}$ then reads

[TABLE]

which is Laplace-inverted numerically.

The effect of retardation is represented for $\kappa=0$ in Fig. 7 with $\Gamma\tau\in[10^{-3},10^{3}]$ and $\gamma_{d}/\Gamma\in[10^{-5},10^{-1}]$. In all these figures we see that the effect is to rescale the parameters of the system. In Figs. 7(a,b) we see in particular that the fidelity for quantum state transfer is almost constant for all values of the retardation $\Gamma\tau$. This can be understood as while increasing $\Gamma\tau$ increases the state transfer time $t_{\text{max}}$, the decay rate of the dark state in Eq. (37) decreases, resulting in a constant overall loss probability. In Fig. 7(c,d) we see that the required optimal $\Omega$ decreases, such that the transfer time $t_{\text{max}}$ increases linearly with the delay at large $\Gamma\tau$. For finite attenuation $\kappa L$ this reduces the transfer fidelity.

IV Effect of finite Lamb-Dicke parameter and temperature

Here we discuss the effects of phononic degrees of freedom for atoms trapped in optical lattices with finite Lamb-Dicke parameter $\eta$, and thermal phonon distribution with mean number $n_{\text{th}}$. We derive a correction to $\mathcal{H}$, and in particular show that the spread of the atomic wavefunction leads to an additional individual decay of each atom of $\gamma_{e}\eta^{2}\left(2n_{\text{th}}+1\right)$.

IV.1 Model

Assuming that each atom is trapped with a harmonic potential with frequency $\omega_{v}$, the non-hermitian Hamiltonian $H_{\text{dip}}$, including now the coupling to motional degrees of freedom, reads

[TABLE]

where we performed a Taylor expansion for the Green’s tensor, with $\eta=k_{0}/\sqrt{2m\omega_{v}}$ the Lamb-Dicke parameter, $a_{\boldsymbol{j},\alpha}$ the annihilation operator of the motional excitation of atom ${\boldsymbol{j}}$ along axis $\alpha\in(x,y,z)$, $\hat{\boldsymbol{r}}_{\boldsymbol{j}}$ the quantized coordinates of atom $\boldsymbol{j}$ relative to its trap center position ${\boldsymbol{r}}_{\boldsymbol{j}}$, and where ${\boldsymbol{\nabla}}$ acts on ${G}({\boldsymbol{r}})$. We assumed the unperturbed density matrix of the system factorizes as $\rho^{\left(0\right)}=\rho_{\text{at}}\otimes\rho_{\text{th}}$ where $\rho_{\text{th}}$ stands for a thermal distribution of phononic modes with mean number $n_{\text{th}}$, while $\rho_{\text{at}}$ accounts for the internal atomic degrees of freedom.

IV.2 Elimination of phonon modes

We now perform an adiabatic elimination of the motional degrees of freedom in Eq. (38). We write $H_{\text{dip}}=H_{0}+V$ as a sum of a free Hamiltonian $H_{0}$ and an interaction term $V$. Assuming $\eta\sqrt{2n_{\text{th}}+1}\ll 1$, we can truncate the expansion in Eq. (38) to second order. We thus have

[TABLE]

Moving to an interaction picture with respect to $H_{0}$, we get an effective Hamiltonian $\tilde{H}$ to second order perturbation in $V$, assuming further that $\gamma\eta\sqrt{2n_{\text{th}}+1}\ll\omega_{v}$. We then obtain

[TABLE]

where we denote $\text{Tr}_{\text{ph}}$ for the trace over phononic degrees of freedom.

The first term in Eq. (39) is evaluated as

[TABLE]

Using the vector field identity ${\boldsymbol{\nabla}}\times(\boldsymbol{\nabla}\times\boldsymbol{V})=\boldsymbol{\nabla}\left(\boldsymbol{\nabla}\cdot\boldsymbol{V}\right)-(\boldsymbol{\nabla}\cdot\boldsymbol{\nabla})\boldsymbol{V}$ and Eq. (8), we further get

[TABLE]

The last term in this equation stands for the variation of the longitudinal part of the Green’s tensor, and introduces a renormalization of the near-field interaction. Dropping this near-field term, which decays as $|\boldsymbol{r}_{\boldsymbol{j}}-\boldsymbol{r}_{\boldsymbol{j}^{\prime}}|^{3}$, we finally have

[TABLE]

Using the identity $\int_{-\infty}^{t}ds\text{Tr}_{\text{ph}}\left[\hat{\alpha}_{\boldsymbol{j}}\left(t\right)\hat{\alpha}_{\boldsymbol{j}^{\prime}}^{\prime}\left(s\right)\rho_{\text{th}}\right]=\delta_{\alpha,\alpha^{\prime}}\delta_{\boldsymbol{j},\boldsymbol{j}^{\prime}}/\left(i\omega_{v}\right)$, with $\hat{\alpha}_{\boldsymbol{j}}$ the component of $\hat{\boldsymbol{r}}_{\boldsymbol{j}}$ along axis $\alpha$, we get for the second term in Eq. (39)

[TABLE]

This term can be evaluated in the limit of infinite arrays ($N_{\perp}\to\infty$), by splitting the sum as

[TABLE]

Using the representation of Eq. (9) for the Green’s tensor and the property of Eq. (10), we then find that the dark and bright states $(1/\sqrt{2N})\sum_{\boldsymbol{j}_{\perp}}(\sigma^{+}_{(\boldsymbol{j}_{\perp},1)}\pm\sigma^{+}_{(\boldsymbol{j}_{\perp},2)})\ket{\mathcal{G}}$ are eigenstates of the Hamiltonian in Eq. (41) with purely real eigenvalues, provided $\delta_{\perp}<\lambda_{0}$ and $k_{0}L=m\pi$. Therefore, for these states the term in Eq. (41) contributes only a frequency renormalization.

Thus, only the term of Eq. (40) contributes to the radiation, such that the decay rates are corrected as

[TABLE]

This shows an additional individual decay for each atom, and an additional rescaling of the interatomic interaction. In order have a given ratio for $\gamma_{d}/\gamma_{b}$, we must thus satisfy the condition $\eta^{2}(2n_{\text{th}}+1)\lesssim\gamma_{d}/\gamma_{b}$.

V Effect of missing atoms

Here we discuss how the presence of defects in the atomic arrays affect the system. The effect of holes can be accounted for by writing for the dipole-dipole interaction Hamiltonian of Eq. (7)

[TABLE]

where $\mathcal{H}^{\text{ideal}}$ is the matrix without defects, and $\mathcal{H}^{\text{holes}}_{\boldsymbol{j},\boldsymbol{j}^{\prime}}=\mathcal{H}^{\text{ideal}}_{\boldsymbol{j},\boldsymbol{j}^{\prime}}$ if atom $\boldsymbol{j}$ or $\boldsymbol{j}^{\prime}$ is missing and $\mathcal{H}^{\text{holes}}_{\boldsymbol{j},\boldsymbol{j}^{\prime}}=0$ otherwise.

We first consider the situation where a single atom, say atom $\boldsymbol{i}$, is missing from the arrays. The dipole-dipole interaction Hamiltonian can be expressed from Eq. (42) as

[TABLE]

For the state $\ket{\psi_{n}}\propto\sum_{\boldsymbol{j}\neq\boldsymbol{i}}(c_{n})_{\boldsymbol{j}}\sigma_{\boldsymbol{j}}^{+}\ket{\mathcal{G}}$, we then obtain

[TABLE]

Let us now consider a situation where each atom ${\boldsymbol{i}}$ has a probability $0\leq p\leq 1$ of being missing, and denote the associated random variable as $s_{\boldsymbol{i}}$. We then have

[TABLE]

Denoting $\overline{\phantom{aa}}$ for the statistical average, and using here $\overline{s_{\boldsymbol{i}}}=p$ and $\overline{s_{\boldsymbol{i}}s_{\boldsymbol{i}^{\prime}}}=p^{2}+\delta_{\boldsymbol{i},\boldsymbol{i}^{\prime}}(p-p^{2})$, we obtain for $\ket{\psi_{n}}\propto\sum_{\boldsymbol{j}}(c_{n})_{\boldsymbol{j}}(1-s_{\boldsymbol{j}})\sigma_{\boldsymbol{j}}^{+}\ket{\mathcal{G}}$

[TABLE]

In Fig. 8 we show the agreement between this expression and numerical simulations, where the dark and bright states decay rates are averaged over $100$ realizations of $s_{\boldsymbol{i}}$. In order to achieve a given ratio for $\gamma_{d}/\gamma_{b}$, we must thus have $p\lesssim\gamma_{d}/\gamma_{b}$.

VI Multiple excitations

For states with more than one atom in $\ket{e}$ or $\ket{s}$ the fact that each atom cannot support more than a single excitation generates an atomic non-linearity, which induces an additional decay rate. The state $\ket{\psi_{d}^{(2)}}\propto(\sigma_{d}^{+})^{2}\ket{\mathcal{G}}$ for instance is not an eigenstate of $H_{\text{dip}}$, however we can treat this non-linearity in first order perturbation theory. We get

[TABLE]

The decay rate per excitation, as represented in Fig. 2(a), is then obtained as

[TABLE]

From Eq. (43), we can interpret the result of Eq. (44) as one of the two excitations decays with rate $\gamma_{d}$, and acts as a defect for the other excitation with probability $p=\sum_{\boldsymbol{j}_{\perp}}|(v_{d})_{\boldsymbol{j}_{\perp}}|^{4}$ identified as twice the inverse participation ratio of the dark state.

VII Implementation with four-level atoms

In all the calculations above and in the main text we treated the atoms as two-level systems with a circular transition. Similar results can however also be obtained using instead atoms with a single ground state $\ket{g}_{\boldsymbol{j}}$ and three excited states $\ket{e_{i}}_{\boldsymbol{j}}$ ($i=x,y,z$), where $i$ denotes the dipole orientation axis. The non-hermitian dipole-dipole interaction Hamiltonian from Eq. (7) generalizes to

[TABLE]

where $\sigma_{\boldsymbol{j},i}^{-}=\ket{g}_{\boldsymbol{j}}\!\bra{e_{i}}$, which now mixes states with different polarizations. Diagonalizing this Hamiltonian, we obtain a degenerate pair of dark and bright states, polarized in the $x-y$ plane, with decay rates represented in Fig. 9. Notably, these decay rates remain close to the values obtained for two-level atoms.