Optical waveguiding by atomic entanglement in multilevel atom arrays

TL;DR

This paper demonstrates that atomic entanglement in multilevel atom arrays can enable optical waveguiding, revealing collective quantum effects that persist beyond simple atomic models and advancing understanding of realistic quantum emitter systems.

Contribution

It shows that waveguiding arises from atomic entanglement in multilevel arrays, extending previous models to more realistic atomic structures with hyperfine levels.

Findings

Waveguiding can result from atomic entanglement in multilevel atom arrays.

Necessary conditions for waveguiding involve specific entanglement and atomic structure.

The work advances understanding of collective effects in realistic atomic systems.

Abstract

The optical properties of sub-wavelength arrays of atoms or other quantum emitters have attracted significant interest recently. For example, the strong constructive or destructive interference of emitted light enables arrays to function as nearly perfect mirrors, support topological edge states, and allow for exponentially better quantum memories. In these proposals, the assumed atomic structure was simple, consisting of a unique electronic ground state. Within linear optics, the system is then equivalent to a periodic array of classical dielectric particles, whose periodicity supports the emergence of guided modes. However, it has not been known whether such phenomena persist in the presence of hyperfine structure, as exhibited by most quantum emitters. Here, we show that waveguiding can arise from rich atomic entanglement as a quantum many-body effect, and elucidate the necessary conditions. Our work represents a significant step forward in understanding collective effects in arrays of atoms with realistic electronic structure.

Equations

$${\hat{\Sigma }}_{iq}^{†}=\sum _{m_g=-F_g}^{F_g}{C}_{m_g,q}{\hat{\sigma }}_{F_e{m}_g-q,F_g{m}_g}^i,$$

$$C_{m_g,q}=(-1{)}^{F_g-m_g}\left(\begin{array}{ccc}{F}_g& 1& F_e\\ -m_g& q& m_g-q\end{array}\right),$$

$$\mathcal{H}=\hslash \sum _{i,j=1}^N\sum _{q,q^{\prime }=-1}^1{J}_{ijq{q}^{\prime }}{\hat{\Sigma }}_{iq}^{†}{\hat{\Sigma }}_{j{q}^{\prime }},$$

$$\mathcal{L}[\rho ]=\sum _{i,j=1}^N\sum _{q,q^{\prime }=-1}^1\frac{{\Gamma }_{ijq{q}^{\prime }}}{2}(2{\hat{\Sigma }}_{j{q}^{\prime }}\rho {\hat{\Sigma }}_{iq}^{†}-{\hat{\Sigma }}_{iq}^{†}{\hat{\Sigma }}_{j{q}^{\prime }}\rho$$

$$-\rho {\hat{\Sigma }}_{iq}^{†}{\hat{\Sigma }}_{j{q}^{\prime }}),$$

$$J_{ijq{q}^{\prime }}=-\frac{{\mu }_0{\omega }_0^2}{\hslash }|\wp {|}^2{\hat{\text{e}}}_q\cdot \text{Re}\mathbf{G}({\mathbf{r}}_i,{\mathbf{r}}_j,{\omega }_0)\cdot {\hat{\text{e}}}_{q^{\prime }}^{\ast },$$

$${\Gamma }_{ijq{q}^{\prime }}=\frac{2{\mu }_0{\omega }_0^2}{\hslash }|\wp {|}^2{\hat{\text{e}}}_q\cdot \text{Im}\mathbf{G}({\mathbf{r}}_i,{\mathbf{r}}_j,{\omega }_0)\cdot {\hat{\text{e}}}_{q^{\prime }}^{\ast },$$

$$\mathbf{G}(\mathbf{r},{\omega }_0)=\frac{e^{\mathrm{i}{k}_{0}r}}{4\pi k_0^2{r}^3}[(k_0^2{r}^2+\mathrm{i}{k}_{0}r-1)1+$$

$$+(-k_0^2{r}^2-3\mathrm{i}{k}_{0}r+3)\frac{\mathbf{r}\otimes \mathbf{r}}{r^2}],$$

$${\mathcal{H}}_{\text{eff}}=\hslash \sum _{i,j=1}^N\sum _{q,q^{\prime }=-1}^1\left(J_{ijq{q}^{\prime }}-\mathrm{i}\frac{{\Gamma }_{ijq{q}^{\prime }}}{2}\right){\hat{\Sigma }}_{iq}^{†}{\hat{\Sigma }}_{j{q}^{\prime }}.$$

$${\hat{\mathbf{E}}}^{+}(\mathbf{r})={\mu }_0{\omega }_0^2\sum _{j=1}^N\sum _{q=-1}^1\mathbf{G}(\mathbf{r},{\mathbf{r}}_j,{\omega }_0)\cdot {\hat{\text{e}}}_q^{\ast }\wp {\hat{\Sigma }}_{jq}.$$

$${\mathcal{H}}_{\text{eff}}=\hslash \sum _{i,j=1}^N\sum _{q=-1}^1\left(J_{ijq}-\mathrm{i}\frac{{\Gamma }_{ijq}}{2}\right){\hat{\Sigma }}_{iq}^{†}{\hat{\Sigma }}_{jq}.$$

$${\tilde{\mathcal{H}}}_{\text{q}}=\hslash \sum _k\left(J_{k,q}-\mathrm{i}\frac{{\Gamma }_{k,q}}{2}\right){\hat{S}}_{k,q}^{†}{\hat{S}}_{k,q},$$

$$J_{k,q}=-\frac{{\mu }_0{\omega }_0^2}{\hslash }|\wp {|}^2{\hat{\text{e}}}_q\cdot \text{Re}\tilde{\mathbf{G}}(k,{\omega }_0)\cdot {\hat{\text{e}}}_q^{\ast },$$

$${\Gamma }_{k,q}=\frac{2{\mu }_0{\omega }_0^2}{\hslash }|\wp {|}^2{\hat{\text{e}}}_q\cdot \text{Im}\tilde{\mathbf{G}}(k,{\omega }_0)\cdot {\hat{\text{e}}}_q^{\ast }.$$

$$|\psi ⟩=\sum _{i,j=1}^N{c}_{ij}{\hat{\sigma }}_{20}^i{\hat{\sigma }}_{10}^j{|0⟩}^{\otimes N}.$$

$${\mathcal{C}}_{mn}(i,j)=⟨{\hat{\sigma }}_{mm}^i{\hat{\sigma }}_{nn}^j⟩-⟨{\hat{\sigma }}_{mm}^i⟩⟨{\hat{\sigma }}_{nn}^j⟩$$

$$H=\hslash \sum _{q=-1}^1\sum _k\left(J_{k,q}-\mathrm{i}\frac{{\Gamma }_{k,q}}{2}\right){\hat{S}}_{k,q}^{†}{\hat{S}}_{k,q},$$

$$|{\psi }_{\tilde{k}}⟩=\mathcal{N}\sum _{j=1}^N{e}^{\mathrm{i}\tilde{k}{z}_j}(|2_j⟩|{\mathcal{D}}_{3/2}^j⟩+\beta |3_j⟩|{\mathcal{D}}_{1/2}^j⟩$$

$$+\beta |4_j⟩|{\mathcal{D}}_{-1/2}^j⟩+|5_j⟩|{\mathcal{D}}_{-3/2}^j⟩),$$

$$|{\mathcal{D}}_{\alpha }^j⟩=\sum _{m\ne j}{\mathcal{P}}_m\left\{{|0⟩}^{\otimes n_0}\otimes {|1⟩}^{\otimes n_1}\right\}$$

$$H|{\psi }_{\tilde{k}}⟩=\hslash \sum _{q=-1}^1\left(J_{\tilde{k},q}-\mathrm{i}\frac{{\Gamma }_{\tilde{k},q}}{2}\right){\hat{S}}_{\tilde{k},q}^{†}{\hat{S}}_{\tilde{k},q}|{\psi }_{\tilde{k}}⟩={\lambda }_{\tilde{k}}|{\psi }_{\tilde{k}}⟩,$$

$$J_{ij,q=\pm 1}^{\prime }(\text{r})=-\frac{360\Delta }{7}{\left(\frac{d}{\pi r}\right)}^4.$$

Full text

Realizing efficient atom-light interactions is a major goal in quantum optics. Due to the intrinsically weak coupling between photons and atoms in free space, atomic ensembles have risen as one of the workhorses of the field, as the interaction probability with photons is enhanced due to the large number of atoms in the cloud Hammerer et al. (2010). Atomic ensembles have broad potential applications, which include, among others, photon storage and retrieval Hau et al. (1999); Fleischhauer et al. (2005); Hammerer et al. (2010), few-photon non-linear optics Pritchard et al. (2013); Peyronel et al. (2012); Gorniaczyk et al. (2014); Tiarks et al. (2014), and metrology Wasilewski et al. (2010); Leroux et al. (2010); Sewell et al. (2012). The fidelity of an atomic ensemble in carrying any of these applications is fundamentally limited by the so-called optical depth, which is a product of the interaction probability between a single atom and a photon in a given optical mode and the total number of atoms. While the important role of optical depth is ubiquitously stated in literature Gorshkov et al. (2007a, b); Chang et al. (2012); Paulisch et al. (2016), the underlying arguments in fact rely on one crucial assumption: that the atoms do not interact with each other and thus that photon emission happens at a rate given by that of single atoms. It is clear, however, that this approximation breaks down when atoms are close to each other, as photon emission is a wave phenomenon and interference and multiple scattering effects will be relevant at short distances.

In dense and ordered atomic arrays Lester et al. (2015); Endres et al. (2016); Barredo et al. (2016); Labuhn et al. (2016); Barredo et al. (2018); Bernien et al. (2017); Bakr et al. (2009, 2010); Weitenberg et al. (2011), strong constructive or destructive interference of light emitted by excited atoms allows one to exceed the fidelities predicted by these simple optical depth arguments in applications Asenjo-Garcia et al. (2017a). For example, it has been theoretically shown that interference can impact communication and metrology applications: it enables both an exponential improvement in the fidelity of a quantum memory Asenjo-Garcia et al. (2017a); Manzoni et al. (2018) and an improvement of the signal-to-noise ratio in optical lattice clocks Ostermann et al. (2013); Henriet et al. (2019). More generally, interference in arrays can give rise to exotic phenomena Jenkins and J.Ruostekoski (2012); Plankensteiner et al. (2015); Bettles et al. (2016); Facchinetti et al. (2016); Sutherland and Robicheaux (2016); Zoubi and Ritsch (2010, 2011); Grankin et al. (2018); Guimond et al. (2019), which have no counterpart in disordered atomic gases. These include perfect reflection of light Bettles et al. (2015); Shahmoon et al. (2017) or the existence of guided topological edge states of light in two-dimensional arrays Perczel et al. (2017); Bettles et al. (2017).

In these previous theoretical works, the atoms were assumed to have a unique electronic ground state. For two level atoms, and within the single-excitation manifold, multiple scattering enables a process where an excited atom $i$ can interact and exchange its excitation with another atom $j$ in its ground state [shown in Fig. 1(a)]. The resulting dynamical equations are exactly equivalent to $N$ classical polarizable dipoles interacting via their radiated fields. In particular, it is well known that ordered arrays of dielectric particles can support lossless guided modes Joannopoulos et al. (2008); Blaustein et al. (2007); Halir et al. (2015). Within the context of infinite atomic arrays, waveguiding manifests itself in the form of perfectly “subradiant” single-excitation states with zero decay rate Sutherland and Robicheaux (2016); Asenjo-Garcia et al. (2017a), a key idea underlying the previously proposed phenomena. In reality, though, most atoms display a rich hyperfine structure – which arises from the coupling between the total electron and nuclear angular momenta – and have more than one ground state. Given the growing body of theoretical and experimental literature about atomic arrays, it is critical to understand the underlying physics of collective optical phenomena for atoms with non-trivial internal structure.

The complexity introduced by hyperfine structure is illustrated in Fig. 1, where one can see that light-mediated dipole-dipole interactions generally do not allow the atomic dynamics to be confined to a two-level subspace. In particular, even if atoms are initialized in such a subspace, emitted photons can drive other atoms out of the two level manifold, as photons do not have a uniform polarization in space. Once an atom is excited out of this subspace, the possibility to decay into unoccupied ground states cannot be suppressed by interference. Thus, even for a single excitation, the mechanism of subradiance, if it exists, could involve some many-body phenomenon. Indeed, the condition for subradiance to exist has already been investigated in the “Dicke” limit Hebenstreit et al. (2017), where all atoms are located at a single point and thus effectively interact with a single, common electromagnetic mode. Interestingly, it was found that subradiance required a specific entanglement structure within the ground state manifold.

In this manuscript, we tackle the problem of collective effects in extended arrays of atoms with hyperfine structure. In particular, using a generalized “spin model” describing dipole-dipole interactions in the presence of hyperfine structure, we identify and analyze different classes of subradiant single-excitation states in a 1D atomic array. We find that the classical waveguiding effect still underlies the vast majority of subradiant states. Here, over large spatial regions, atoms in the array essentially live within a two-level subspace, interrupted by local “defect” states or domain walls that divide “phase separated” regions. However, we also describe a new, truly many-body mechanism, where waveguiding is enabled through rich, long-range entanglement within the ground-state manifold, and we elucidate the necessary conditions for its existence. These results are an important step forward in understanding collective effects in atoms with realistic electronic structure.

I Spin model and minimal toy atom

Here, we introduce a spin model to describe the photon-mediated quantum interactions between atoms. We will consider atoms whose electronic structure consists of a ground- and excited-state manifold, with total hyperfine angular momenta quantum numbers $F_{g}$ and $F_{e}$, respectively. A complete basis can be obtained by labeling states according to the projection of angular momentum along the $z$-axis, $\mathinner{|{F_{g/e}\,m_{g/e}}\rangle}$, where $m_{g/e}\in[-F_{g/e},F_{g/e}]$ are Zeeman sublevels. We thus can describe the state of an atom by its quantum numbers, i.e., $\mathinner{|{F_{e}\,m_{e}}\rangle}$ if it is excited or $\mathinner{|{F_{g}\,m_{g}}\rangle}$ if it is in the ground state manifold. The ground and excited states couple to light via well-defined selection rules, such that $m_{e}=m_{g}+q$, with $q=\{0,\pm 1\}$ denoting the units of angular momentum that can be transferred by a photon. We can define an atomic raising operator that depends on $q$ as

[TABLE]

where $\hat{\sigma}^{i}_{F_{e}m_{g}-q,F_{g}m_{g}}=\mathinner{|{F_{e}\,m_{g}-q}\rangle}_{i}\mathinner{\langle{F_{g}\,m_{g}}|}_{i}$ is the atomic coherence operator between the ground and excited states of atom $i$. This operator conveniently groups all the possible transitions which transfer $q$ units of angular momentum along $z$. These are weighted by the Clebsch-Gordan coefficients

[TABLE]

written here in terms of a Wigner 3j-symbol, which reflect the different strengths that the possible transitions couple to light of a given polarization.

Spin model for multilevel atoms– Intuitively, the interaction with light allows for processes of photon-mediated emission and re-absorption, where the excitation of one atom decays and another is excited. We describe such dynamics by means of a spin model Gruner and Welsch (1996); Buhmann (2012); Asenjo-Garcia et al. (2017b); Gross and Haroche (1982); Evers et al. (2006); Kiffner et al. (2007); Munro et al. (2018); HOR where we integrate out the photons and find an all-atomic density matrix that only depends on the internal degrees of freedom of the atoms. Specifically, the evolution of the atomic density matrix $\hat{\rho}_{\rm A}$ obeys $\dot{\hat{\rho}}_{\rm A}=-({\rm i}/\hbar)\,[\mathcal{H},\hat{\rho}_{\rm A}]+\mathcal{L}[\hat{\rho}_{\rm A}]$, where $\mathcal{H}$ is the Hamiltonian, and $\mathcal{L}[\hat{\rho}_{\rm A}]$ is the Lindblad operator. For atoms with hyperfine structure and resonance frequency $\omega_{0}$, these operators read

[TABLE]

where we have defined the polarization-dependent spin-exchange and decay rates as

[TABLE]

with $\bm{\wp}=\mathinner{\langle{F_{g}||e\hat{{\bf r}}||F_{e}}\rangle}$ being the reduced matrix element associated with the transition. In the above equations, $\hat{\textbf{e}}_{\pm 1}=\mp\frac{1}{\sqrt{2}}(\hat{x}\pm{\rm i}\hat{y})$, and $\hat{\textbf{e}}_{0}=\hat{z}$ are spherical basis vectors. Physically, the spin-exchange and decay rates are proportional to the classical field amplitude projected along polarization $q$ at position ${\bf r}_{i}$, due to a classical oscillating dipole of polarization $q^{\prime}$ at ${\bf r}_{j}$. Naturally, they are given in terms of the free-space electromagnetic Green’s tensor, ${\bf G}({\bf r}_{i},{\bf r}_{j},\omega_{0})\equiv{\bf G}({\bf r}_{ij},\omega_{0})$, with ${\bf r}_{ij}={\bf r}_{i}-{\bf r}_{j}$, which reads

[TABLE]

where $r\equiv|{\bf r}|$ and $k_{0}=2\pi/\lambda_{0}=\omega_{0}/c$ is the wave number corresponding to the atomic transition energy. Note that the Green’s function ${\bf G}_{\alpha\beta}$ is a tensor quantity ($\{\alpha,\beta\}=\{x,y,z\}$), as both the electromagnetic field and the atomic transition have specific polarizations. These equations are derived within the Markovian approximation, which is highly justified for our system under consideration Chang et al. (2012); Shi et al. (2015); Guimond et al. (2016).

The dynamics under the master equation can analogously be described in the quantum jump formalism Gardiner and Zoller (2004), where the last two terms in the parentheses of Eq. (3b) are combined with $\mathcal{H}$ to form a non-Hermitian Hamiltonian that characterizes the deterministic evolution. Within this formalism, the effective non-Hermitian Hamiltonian is readily given by

[TABLE]

While this equation is quite general, we now highlight several points, which will be key to understanding collective decay:

• •

First, for a single, isolated atom, we obtain the total decay rate from excited sublevel $m_{e}$ as $\Gamma_{0}=\sum_{q}\Gamma_{iiqq}C_{m_{e}+q,q}^{2}$, which is equal for all $m_{e}$. The branching fraction into a specific ground state $m_{g}=m_{e}+q$ is simply given by $C_{m_{e}+q,q}^{2}$.

• •

Second, the form of $\mathcal{H}_{\text{eff}}$ illustrates that, in general, lowering one atom angular momentum projection by $q$ does not necessarily imply that another atom’s projection is raised by $q$ ($q^{\prime}\neq q$), as $q$ refers to the dipole matrix element and not the global polarization of the emitted photon [see Fig. 1(c)]. Thus, a system initialized in a two-level subspace does not generally remain in that space under coherent dynamics. In particular, even if the initially involved excited state decays to a unique ground state, a different excited state can become populated later in time [see atoms $i$ and $k$ in Fig. 1(c)]. We note that one way to suppress such undesired excitation pathways (e.g., the $\sigma^{+}$ transition of atom $k$) is by applying large Zeeman shifts. While this might indeed constitute a practical way to maintain the desired effects of subradiance to some extent given hyperfine structure, our focus here will be on identifying new mechanisms for subradiance, which fundamentally persist in the presence of multiple pathways.

• •

Finally, we consider the situation where an excited atom has an allowed decay channel into a ground state that is unoccupied by any atoms [illustrated by atom $k$ in Fig. 1(c) following absorption of a $\sigma^{+}$ photon]. In that case, there is no interference; the excited state decays into that ground-state level with the full strength of a single, isolated atom. Presumably, any mechanism for strong subradiance must then prevent this process from occurring.

It will also be useful to calculate the emitted quantum field associated with any given atomic state. It can be shown that the positive frequency component reads

[TABLE]

From the previous equation, we can find the negative-frequency component $\hat{{\bf E}}^{-}({\bf r})$ by taking the Hermitian conjugate of $\hat{{\bf E}}^{+}({\bf r})$. The field also in principle contains a vacuum noise component Asenjo-Garcia et al. (2017a), which will not affect our quantities of interest and is thus not explicitly written.

One-dimensional chain– From now on, we will focus on the problem of a 1D chain. We will see that this system both encodes the classical waveguiding effect previously identified for two-level atoms, and allows for new highly correlated subradiant states due to multi-level structure. Moreover, 1D allows for numerics on relatively large systems (not dominated by boundaries), and we can take advantage of additional conserved quantities that allow us to simplify the problem. For the sake of simplicity, we choose the direction of the chain to align to the quantization axis of the atoms (i.e., along $z$). This guarantees the conservation of polarization of the emitted photons, as ${\bf G}_{\alpha\beta}=0$ along the axis if $\alpha\neq\beta$ (i.e., ${\bf G}_{\alpha\beta}$ is a diagonal tensor in the polarization indices). Therefore, the only non-zero interactions are those preserving polarization (i.e., $q=q^{\prime}$) and the effective Hamiltonian simplifies to

[TABLE]

where $J_{ijq}\equiv J_{ijqq}$ and $\Gamma_{ijq}\equiv\Gamma_{iiqq}$.

Until now, our spin model captures the dynamics of atoms with any kind of hyperfine structure. Hereafter, we focus on what we consider to be the minimal toy model that captures all the relevant physics: 6-level atoms with $F_{g}=1/2$ and $F_{e}=3/2$ [see Fig. 1(b)]. This specific hyperfine structure displays closed transitions, where an excited state with maximum angular momentum projection $m_{e}=\pm 3/2$ can only decay to only one ground state, also with maximum angular momentum projection $m_{g}=\pm 1/2$ in the ground state manifold by emitting a circularly-polarized photon ($q=\mp 1$). It also exhibits excited states with $m_{e}=\pm 1/2$ where decay into two different ground states $m_{g}=\pm 1/2$ is allowed, which involve emission of photons of different polarizations ($q=\pm 1$ or $q=0$). To simplify notation, in what follows we will label ground states with $m_{g}=\{-1/2,1/2\}$ as $\{\mathinner{|{0}\rangle},\mathinner{|{1}\rangle}\}$, respectively, and excited states with $m_{e}=\{-3/2,-1/2,1/2,3/2\}$ as $\{\mathinner{|{2}\rangle},\mathinner{|{3}\rangle},\mathinner{|{4}\rangle},\mathinner{|{5}\rangle}\}$, respectively [see Fig. 1(b)].

The conservation of angular momentum projection along the direction of the chain allows us to diagonalize the Hamiltonian by blocks with well-defined $F_{z}=\sum_{i=1}^{N}m_{i}$, where the sum over Zeeman sublevels includes both ground and excited states (see Fig. 2). Note that subspaces with equal magnitude of angular-momentum projection but opposite sign (i.e., $\pm F_{z}$) display the same energy and decay spectra (i.e. the physics is identical if we simultaneously flip $m_{g}$ to $-m_{g}$ and $m_{e}$ to $-m_{e}$). In what follows, we study single-excitation subradiant states of various angular momentum manifolds, starting by maximum $|F_{z}|$ (as we will show, the most “classical” manifold) and ending by $F_{z}=0$ (the most complex, and where subradiant states enabled by rich entanglement live).

Motivated by previous work on 1D and 2D arrays with simple atomic structure Asenjo-Garcia et al. (2017a), here we seek to elucidate the properties of eigenstates $\mathinner{|{\psi_{\xi}}\rangle}$ of $\mathcal{H}_{\text{eff}}$ in the single-excitation manifold. In particular, such states will have complex eigenvalues $\omega_{\xi}=J_{\xi}-i\Gamma_{\xi}/2$ characterizing the energy shifts and decay rates. Of specific interest will be the identification of states where the decay rate approaches zero as $N\rightarrow\infty$, which implies that the states decouple from radiation fields and correspond to guided modes; and to elucidate the properties of the eigenstates that enable the waveguiding phenomenon. Importantly, while for simple two-level atoms, the Hilbert space of the single-excitation manifold increases as $\sim N$ and encodes classical linear optics, here, the manifold size increases exponentially, raising the possibility for waveguiding through entanglement.

II Two-level subradiance: $|F_{z}|=F_{z}^{\text{max}}$

In each subspace $F_{z}$, it is first helpful to consider the full set of possible basis states. For example, within the single-excitation manifold, the maximum allowed value of $F_{z}$ is given by $F_{z}^{\text{max}}=\pm|F_{e}+(N-1)F_{g}|=\pm(N/2+1)$, which is achieved when one atom is excited in the state $\mathinner{|{2}\rangle}$, with $m_{e}=-3/2$, while all other atoms are in the ground state $\mathinner{|{0}\rangle}$, with $m_{g}=-1/2$. These states are connected through a $\sigma^{-}$ transition (this is, $q=1$). Since $\mathcal{H}_{\text{eff}}$ is block-diagonalizable, $|F_{z}|=F_{z}^{\text{max}}$ thus corresponds to an effective two-level atom subspace.

Numerically, we diagonalize $\mathcal{H}_{\text{eff}}$ within this subspace as a function of atom number $N$ and for selected different lattice constants $d$. For each $N$ and $d$, we then find the eigenstate with the minimum decay rate $\Gamma$, which we plot in Fig. 3(a). In agreement with our previous results Asenjo-Garcia et al. (2017a), we find subradiant states if the interatomic distance $d$ is such that $d<\lambda_{0}/2$, $\lambda_{0}$ being the wavelength of the resonant transition. The decay rate of the most subradiant eigenstate scales as $\Gamma/\Gamma_{0}\sim 1/N^{3}$, as also found previously for two level atoms. We calculate the field intensity radiated by the most subradiant eigenstate $\mathinner{|{\psi}\rangle}$ by numerically finding the expectation value of the intensity operator $I({\bf r})=\mathinner{\langle{\psi|\hat{{\bf E}}^{-}({\bf r})\cdot\hat{{\bf E}}^{+}({\bf r})|\psi}\rangle}$, according to Eq. (7). The intensity pattern in Fig. 3(b) demonstrates strong emitted intensity as the ends of the chain and absence of intensity in the middle, as would be expected for a finite-size waveguide.

As previously remarked, the guiding effect is entirely classical. Because the ground state is unique and there is only one excitation in the system (i.e., $N$ degrees of freedom), one can alternatively deduce this state by considering $N$ coupled harmonic oscillators. In that case, the stationary states under full master equation evolution are coherent states $\mathinner{|{\psi_{\xi}^{ho}}\rangle}=\mathinner{|{\alpha_{1}e^{-i\omega_{\xi}t}}\rangle}\mathinner{|{\alpha_{2}e^{-i\omega_{\xi}t}}\rangle}...\mathinner{|{\alpha_{N}e^{-i\omega_{\xi}t}}\rangle}$. Here, the coherent state amplitude $\alpha_{i}\propto c_{i}$ of harmonic oscillator $i$ is proportional to the spin wavefunction amplitude of the corresponding single-excitation spin eigenstate $\mathinner{|{\psi_{\xi}}\rangle}=\sum_{i}c_{i}\sigma_{20}^{i}\mathinner{|{0}\rangle}^{\otimes N}$ found earlier. These coherent states amplitudes evolve in time as $e^{-i\omega_{\xi}t}$ with a frequency $\omega_{\xi}$ corresponding to the spin state eigenvalue. These harmonic oscillator states physically describe classical resonant dipole arrays, interacting with each others’ radiated fields. In particular, it is well-known that an infinite array can support guided modes. Starting from such states, we can apply an operator $\mathcal{P}$ that projects these states into the single-excitation manifold, which then exactly reproduces the single-excitation spin eigenstates $\mathinner{|{\psi_{\xi}}\rangle}$. This implies that the resulting physics can be understood purely classically – even if the spin eigenstates themselves are technically entangled – as the entanglement simply arises from the projection at the end. As we will see later on, there are entangled states that do not support any analogy with those of classical dipoles.

The waveguiding concept can be additionally confirmed by considering an infinite chain. In this case, we can diagonalize the Hamiltonian using spin waves with well-defined momentum $k$ along $z$, and readily find $\mathcal{H}_{\text{eff}}=\sum_{q=-1}^{1}\tilde{\mathcal{H}}_{\text{q}}$, where

[TABLE]

with $\hat{S}^{\dagger}_{k,q}=N^{-1/2}\sum_{j}e^{{\rm i}kdj}\hat{\Sigma}^{\dagger}_{jq}$ and

[TABLE]

In the above equations, $\tilde{{\bf G}}(k)=\sum_{j}e^{-{\rm i}kdj}{\bf G}(z_{j})$ is the discrete Fourier transform of the free-space Green’s tensor. Eigenstates in the single-excitation manifold can be generated by applying a spin raising operator $S^{\dagger}_{k,1}$ to the product ground state $\mathinner{|{0}\rangle}^{\otimes N}$. For $k>\omega_{0}/c$ – such that the spin wavevector exceeds the wavevector of free-space radiation – one finds that $\Gamma_{k,q}=0$, indicating the decoupling of the spin wave from radiation fields and thus the guided nature of these excitations. We note that the ability to restrict dynamics to a two-level subspace is unique to a 1D chain, as in higher dimensions the Hamiltonian of Eq. (6) is not restricted to $q=q^{\prime}$. This motivates the deeper investigation of subradiance in subspaces of other total $F_{z}$, to find a true many-body mechanism.

III Defect states: $|F_{z}|=|F_{z}^{\text{max}}|-1$

Reducing the angular momentum projection by one unit increases the complexity of the problem. This is manifest in the larger size of the Hilbert subspace, of dimension $dim(\mathcal{H})=N^{2}$. The basis states are of two types: (i) “defect” states where one atom is excited in $\mathinner{|{2}\rangle}$ ($m_{e}=-3/2$), $N-2$ atoms are in $\mathinner{|{0}\rangle}$ ($m_{g}=-1/2$) (reminiscent of the two-level subspace of the previous sub-section), and one “defect” atom is in the ground state $\mathinner{|{1}\rangle}$ ($m_{g}=1/2$); (ii) states where all atoms except one are in the ground state $\mathinner{|{0}\rangle}$, and the excited state involves level $\mathinner{|{3}\rangle}$ ($m_{e}=-1/2$), a state without maximal angular momentum projection. From the considerations of Sec. II, it is already clear that any highly subradiant states must exhibit the following features. First, because the basis states containing excited state $\mathinner{|{3}\rangle}$ have an allowed transition (via a $q=1$ photon) to state $\mathinner{|{1}\rangle}$, which is unoccupied, such states will decay into $\mathinner{|{1}\rangle}$ at the full strength $\Gamma_{0}/3$ of a single, isolated atom. Thus, these basis states must constitute a vanishingly small weight of a subradiant state. Second, for basis states composed of a single atom $j$ in state $\mathinner{|{1}\rangle}$, one can conclude that coherent interactions with an atom $i$ in state $\mathinner{|{2}\rangle}$ allow for an exchange $\mathinner{|{2_{i}\,1_{j}}\rangle}\rightarrow\mathinner{|{0_{i}\,3_{j}}\rangle}$, while we argued already that state $\mathinner{|{3}\rangle}$ is undesirable. Thus, a subradiant state must find some mechanism to “hide” the atom in state $\mathinner{|{1}\rangle}$ from the dynamics.

This intuition agrees with our numerical findings. For the most subradiant eigenstates, the population of state $\mathinner{|{3}\rangle}$ decreases as $\sum_{j}\mathinner{\langle{\hat{\sigma}_{33}^{j}}\rangle}\sim 1/N^{3}$, i.e., in the thermodynamic limit the most subradiant eigenstate is mostly formed by “defect” states and is of the form

[TABLE]

Figure 4(a) shows the spatial profile of the population in levels $\mathinner{|{1}\rangle}$ and $\mathinner{|{2}\rangle}$ of the most subradiant eigenstate of a chain of $N=30$ atoms. The population in the excited state $\mathinner{|{2}\rangle}$ is spatially smooth, as it is distributed among all atoms. On the other hand, the population in the defect level $\mathinner{|{1}\rangle}$ is localized in a superposition between both edges of the chain (in particular, the defect atom is the third one, counting from both edges). The population of level $\mathinner{|{3}\rangle}$ is spatially correlated with that of $\mathinner{|{1}\rangle}$, as $\mathinner{\langle{\sigma_{33}^{j}}\rangle}\propto\mathinner{\langle{E^{-}({\bf r}_{j})E^{+}({\bf r}_{j})}\rangle}\mathinner{\langle{\sigma_{11}^{j}}\rangle}$ (see SI). The population in this level is also localized around the third atom, but negligible in the scale of the figure. As shown in Fig. 4(b), the decay rate of the most subradiant state with atom number scales as $\Gamma/\Gamma_{0}\sim 1/N^{3}$. Figure 4(c) shows the scaling of the field intensity at the position of the central atom of the chain versus that at the defect atom. The latter scales as $I\sim 1/N^{3}$. This confirms our previous intuition, that the defect atom is efficiently “hidden” from the rest of the chain. In particular, the low intensity seen by the defect atom ensures that it cannot be efficiently excited to level $\mathinner{|{3}\rangle}$, and that the large decay rate of this excited state does not destroy subradiance. The precise spatial location of the defect atom depends on microscopic details, and changes depending on the length of the chain (it appears just at the edge for short chains and in atoms further away from the boundaries for longer chains), but the precise location does not qualitatively alter the physics (see SI).

Furthermore, when $F_{z}=|F_{z}^{max}|-n_{d}$ is further reduced, (with $n_{d}\ll|F_{z}^{max}|$), the most subradiant states are of similar character, with long chains of “two-level” atoms interrupted by $n_{d}$ defects, whose potentially detrimental effects are reduced as they are largely decoupled (see SI). Importantly, the character of such states, dominated by the waveguiding of effective two-level atoms, do not yet point to any non-trivial mechanism for subradiance.

IV Phase-separated and symmetric states: $F_{z}=0$

For an even number of atoms, there is a subspace with zero angular momentum projection. This manifold is the largest one, with dimension that scales exponentially with atom number. In this manifold, we have approximately the same number of atoms in each of the degenerate ground states, allowing for many different combinations. We are able to find two different kinds of subradiant states: phase separated states and symmetric states. In the following, we study these two families in detail.

IV.1 Phase-separated states

Decreasing angular momentum increases the number of defects, which tend to localize at the chain edges. In $F_{z}=0$, approximately half of the chain is made of defects, and phase separation and domain walls emerge, which we discuss in this section. Phase-separated states occur when the number of defects is so large that half of the chain is in a “defect state” with respect to the other half. The “active” phase involves one atom excited in $\mathinner{|{2}\rangle}$ (or $\mathinner{|{5}\rangle}$), and $N/2-2$ in $\mathinner{|{0}\rangle}$ ($\mathinner{|{1}\rangle}$), which participate in two-level subradiance: the atoms exchange photons with $q=1(-1)$ and get excited and decay via closed transitions. The “inactive” or defect phase consists of $N/2+1$ atoms in $\mathinner{|{1}\rangle}$ ($\mathinner{|{0}\rangle}$). It should be noted that the eigenstate is a spatial superposition, where the active part can equally occupy the left and right side of the chain, with a domain wall in between the phases, as shown in Fig. 5(a).

To observe signatures of phase separation, we calculate connected correlation functions between different Zeeman sublevels $m,n$. These are defined as

[TABLE]

where the expectation value is calculated over the most subradiant eigenstate, and $i,j$ are atomic indices. Figure 5(b) shows the correlations $\mathcal{C}_{00}(i,j)$ between $\mathinner{|{0}\rangle}$-states in different locations for the most subradiant state of a chain of $N=14$ atoms. We observe $\mathinner{|{0}\rangle}$-states clustering around each other, and “repelling” atoms in $\mathinner{|{1}\rangle}$-states, displaying anti-correlations. The orange line in Fig. 5(c) shows the scaling of the decay rate with respect to atom number. Due to the exponentially large size of the Hilbert space, we can only exactly diagonalize up to $N=14$ atoms. However, based on our intuition of phase separation, we expect the decay rate to closely coincide with that obtained for the most subradiant eigenstate of a chain of $N/2-1$ atoms in the $|F^{\text{max}}_{z}|$ subspace. This calculation is plotted in Fig. 5(c) (red curve) and seems to agree with the scaling of the phase-separated states.

As in our previous analysis of other manifolds, thus far these results do not point to some fundamentally new mechanism for waveguiding: these states can be understood as the chain separating into two domains where “two-level” subradiance occurs independently from the other.

IV.2 Symmetric states

We also numerically find evidence of another type of subradiant state, with significant population in states $\mathinner{|{3}\rangle}$ and $\mathinner{|{4}\rangle}$. These states also show no particular length scales or features in pairwise correlations, suggesting a qualitatively different mechanism for subradiance. To better grasp the underlying physics, we first introduce a “toy model” Hamiltonian, which is inspired by Eq. (9). In particular, while the Hamiltonian can exactly be written in terms of pure spin wave operators for an infinite system, we consider a hypothetical finite system whose Hamiltonian also takes the same form,

[TABLE]

where $\hat{S}^{\dagger}_{k,q}=N^{-1/2}\sum_{j}e^{{\rm i}kdj}\hat{\Sigma}^{\dagger}_{jq}$ and $J_{k,q}$ and $\Gamma_{k,q}$ are given in Eq. (10). In particular, as these quantities correspond to the results for an infinite system, we have that $\Gamma_{k,q}=0$ when $|k|>k_{0}$. The wavevector $k$ is now taken to be a discrete index, with $k=2\pi n/N$ (with $n\in[0,N-1]$) to ensure periodic boundary conditions.

A perfectly subradiant eigenstate $\mathinner{|{\psi_{\tilde{k}}}\rangle}$, of well-defined quasi-momentum $\tilde{k}$, with zero decay rate fulfills $\text{Im}\{H\mathinner{|{\psi_{\tilde{k}}}\rangle}\}=0$. It should be noted that despite the apparent simplicity of Eq. (13), it is in general challenging to diagonalize, as the involved operators have complicated commutation relations. To proceed, we first find a state that fulfills the less demanding condition of being a zero-eigenvalue eigenstate of all spin operators $\hat{S}_{k,q}$ if $k\neq\tilde{k}$. Such a state is found to be

[TABLE]

where $\beta=(C_{-1/2,1}/C_{1/2,1})$, $\mathcal{N}$ is a normalization constant, and

[TABLE]

is a Dicke state where all the atoms except $j$ are in ground states $\mathinner{|{0}\rangle}$ or $\mathinner{|{1}\rangle}$ such that the total angular momentum projection of the Dicke state is $\alpha=\{\pm 1/2,\pm 3/2\}$. In the above expression, n_{0(1)}=(N-1\mathrel{\raise 2.0pt\hbox{\tiny(}\!\raise 1.0pt\hbox{+}\hskip-7.7778pt\raise-2.4pt\hbox{-}\!\raise 2.0pt\hbox{\tiny)}}2\alpha)/2 is the number of atoms in state $\mathinner{|{0}\rangle}(\mathinner{|{1}\rangle})$, and $\sum_{m}\mathcal{P}_{m}$ denotes the sum over all distinct permutations of the ground states. It should be noted that Dicke states are known to exhibit significant multipartite entanglement Toth (2007); Duan (2011). More importantly, the entanglement of this state cannot be obtained simply by applying a projection operator on a coherent state, as was the case for two-level subradiance. Remarkably, the limit where $H$ is composed of only one $k$-vector corresponds to Ref. Hebenstreit et al. (2017), where they show the existence of states that are zeros of a single jump operator. Our states, however, must satisfy many relationships involving all spin wave operators $\hat{S}_{k\neq\tilde{k},q}$.

For the state in Eq. (14) to be a lossless eigenstate of $H$ it needs to fulfill

[TABLE]

with $\lambda_{\tilde{k}}\in\mathbb{R}$. In the equation above, we have utilized the property that $S_{k\neq\tilde{k},q}\mathinner{|{\psi_{\tilde{k}}}\rangle}=0$, to eliminate all but a single spin-wave operator from $H\mathinner{|{\psi_{\tilde{k}}}\rangle}$. For $|\tilde{k}|>k_{0}$ (such that $\Gamma_{\tilde{k},q}=0$), we find that the eigenstate equation (16) is satisfied provided that the dispersion relations for both polarizations $|q|={0,1}$ coincide at $\tilde{k}$, i.e. $J_{\tilde{k},q}\equiv J_{\tilde{k}}$, with a corresponding eigenvalue $\lambda_{\tilde{k}}=C_{-1/2,1}^{2}\hbar J_{\tilde{k}}$. We subsequently show that an intersection of the dispersion relations at $|\tilde{k}|>k_{0}$ gives rise to waveguiding not only in our toy model, but in the original physical system.

In free space the dispersion relations $J_{k,\pm 1}$ and $J_{k,0}$ are in general different Asenjo-Garcia et al. (2017a), and only intersect for a given $\tilde{k}>k_{0}$ for distances $d\lesssim 0.17\lambda_{0}$. Numerically, however, it is difficult to confirm directly that an intersection in such a case leads to a decay rate approaching zero as $N\rightarrow\infty$. We attribute this to the fact that for the limited $N$ that we can simulate, the small lattice constant $d\ll\lambda_{0}$ still results in significant finite-size effects due to the short length of the chain, and because the intersection of the dispersion relations occurs close to radiative wavevectors $|\tilde{k}|\sim k_{0}$ (this is, for a short system, significant components of the many-body wavefunction still have wavevectors that couple to radiation). Thus, to better confirm our hypotheses, we add an additional short-range term to the $q=\pm 1$ interaction rates, $J_{ij,q=\pm 1}\rightarrow J_{ijq=\pm 1}+J^{\prime}_{ijq=\pm 1}$, given by

[TABLE]

Here, $\Delta=J_{k=\pi/d,q=0}-J_{k=\pi/d,q=1}$ is the difference between the original free-space dispersion relations of Eq. (10), evaluated at the edge of the Brillouin zone $k=\pi/d$. It can readily be shown that this extra term guarantees that both the dispersion relations and the slopes of the $q=0,\pm 1$ polarizations coincide at the Brillouin zone edges, $|k|=\pi/d$. In the following, we will not focus on the details of how such dispersion engineering is implemented, although we note that it can potentially be realized by introducing some dielectric structure to change the Green’s function itself, or by dressing of atomic levels.

Figure 6(a) shows the spatial profile of the most subradiant eigenstate of a finite chain of $N=16$ atoms obtained with dispersion engineering. The population in levels $\mathinner{|{2}\rangle}$ and $\mathinner{|{5}\rangle}$ is identical, and so are the populations in $\mathinner{|{3}\rangle}$ and $\mathinner{|{4}\rangle}$. In Fig. 6(b), we plot the decay rate of the most subradiant state vs. atom number. In particular, we observe a scaling of the decay rate with atom number of $\Gamma/\Gamma_{0}\sim 1/N^{3}$, identical to that observed for “classical” waveguiding. Figure 6(c) shows the fidelity between the most subradiant eigenstate of the finite chain and the toy-model state $\mathinner{|{\psi_{\tilde{k}}}\rangle}$ given by Eq. (14). There is extremely good agreement at $\tilde{k}=\pi/d$. In particular the infidelity at $\tilde{k}=\pi/d$ scales as $\epsilon=1-|\mathinner{\langle{\psi|\psi_{\tilde{k}}}\rangle}|^{2}\sim 1/N^{2}$, which indicates that our toy model accurately captures the physics of the actual system. Finally, Fig. 6(d) shows the field intensity emitted by the most subradiant eigenstate of the chain with dispersion engineering. The intensity pattern is consistent with waveguiding along the chain, where most of the field is radiated through the edges. We note that these states are robust against classical fluctuations in the atomic positions, as shown in the SI.

We emphasize the importance of realizing an intersection of the dispersion relations $J_{k,\pm 1}$ and $J_{k,0}$, in order to yield this new class of highly entangled, waveguiding states. In particular, in the gray points and curves of Figs. 6(b) and (c), respectively, we plot the corresponding results where only the free-space Green’s function is used, without the additional term of Eq. (17). It can be seen that although a moderate reduction of decay rate $\Gamma/\Gamma_{0}\sim 0.1$ can be achieved, this rate apparently does not decrease with increasing $N$. Furthermore, the overlap fidelity between the numerically obtained eigenstate and the ansatz state does not approach 1, as the differing coherent interactions associated with $q=0,\pm 1$ mix in additional contributions to the eigenstate. We note that, in contrast with the dispersion-engineered situation, these are not the most subradiant eigenstates in the $F_{z}=0$ manifold (as they are given by the phase-separated states). Thus, we have found these states by exactly diagonalizing the Hamiltonian and filtering out states such that $\sum_{j}\mathinner{\langle{\sigma^{j}_{33}}\rangle}>\sum_{j}\mathinner{\langle{\sigma^{j}_{22}}\rangle}$, and then selecting the most subradiant states that satisfy that condition. In this situation, we need to diagonalize the full Hamiltonian matrix, instead of simply finding the eigenvalue with the smallest decay rate by using iterative sparse-matrix diagonalization algorithms, which limits the size of the chain that we can successfully diagonalize to $N=12$.

V Discussion

To summarize, we have shown that arrays of atoms with hyperfine structure can support highly subradiant, waveguiding states, and we have elucidated the conditions for their existence. In contrast to the “classical” effect that occurs in atoms with a single ground state, here, the waveguiding is fundamentally enabled by rich, many-body correlations within the ground-state manifold.

Having shown the existence of true waveguiding states, one interesting question going forward is how novel phenomena or applications, previously identified for simple atoms, can be encoded into these highly entangled states. For example, it would be interesting to investigate if such states support more powerful quantum memory protocols, or if generalizations of such states exist in higher dimensions, such as to support topological edge states. In the case of higher dimensions, while numerics will likely be highly challenging, a promising approach could be the generalization of toy models such as Eq. (13). It would also be interesting to investigate multi-excitation subradiant states, e.g., to see if they exhibit the same “fermionic” correlations as multi-excitation states in atoms with simple (two-level) structure Asenjo-Garcia et al. (2017a).

An important associated question is how the necessary many-body entanglement in the ground-state manifold can be generated in the first place. We speculate that the same correlated dissipation processes encoded in the dipole-dipole interactions of Eq. (3) might be used to generate the necessary entanglement, such as through the steady state obtained under constant driving (i.e., correlated optical pumping).

More generally, the insight developed in this work could give rise to broader opportunities. For example, under constant driving, one could investigate whether correlated dissipation can give rise to useful many-body correlations within the ground-state manifold, such as for quantum-enhanced metrology Leroux et al. (2010); Wasilewski et al. (2010); Bohnet et al. (2014); Hosten et al. (2016); Masson and Parkins (2019). It would also be interesting to more systematically understand the forms and range of entanglement that can arise, as the Green’s function is varied (e.g. through a dielectric structure) and/or the atomic level structure is altered. Finally, as an optical phenomenon, it is intriguing that we have identified a mechanism for waveguiding that explicitly relies on entanglement. It would be interesting to more broadly search for novel mechanisms of waveguiding and other optical effects BPP , which can only emerge through quantum correlations.

Acknowledgments.– A. A.-G. work was supported as part of Programmable Quantum Materials, an Energy Frontier Research Center funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), under award DE-SC0019443. D.E.C. acknowledges support from Fundación Ramon Areces, Fundacio Privada Cellex, Spanish MINECO Severo Ochoa Program SEV-2015-0522, Spanish Plan Nacional Grant ALIQS (funded by the MCIU, AEI, and FEDER), CERCA Programme/Generalitat de Catalunya, ERC Starting Grant FOQAL, and AGAUR Grant 2017 SGR 1334. H.J.K. acknowledges funding from the Office of Naval Research (ONR) Grant N00014-16-1-2399, the ONR Multidisciplinary University Research Initiative (MURI) Quantum Opto-Mechanics with Atoms and Nanostructured Diamond Grant N00014-15-1-2761, and the Air Force Office of Scientific Research MURI Photonic Quantum Matter Grant FA9550-16-1-0323.

Supplementary Information

Here we present further details on the physics of defect states and we show that subradiant states are robust against (classical) fluctuations in atomic positions.

Appendix A Defect states

In Fig. 7(a) we plot the decay rate of the most subradiant state in the manifold $|F_{z}|=|F_{z}^{\text{max}}|-1$ as a function of atom number $N$, and classified by the position of the defect atom ($n=1,2,3,4$). It can be seen that although the position of the defect atom $n$ changes for the most subradiant state overall (minimized over all possible positions $n$), all of these defect states exhibit qualitatively similar behavior. Considering the most subradiant defect state, we plot in Fig. 7(b) the total population of state $\mathinner{|{3}\rangle}$, which is seen to decrease with atom number as $\sum_{j}\mathinner{\langle{\hat{\sigma}_{33}^{j}}\rangle}\sim 1/N^{3}$. In manifolds with lower angular momentum, defect atoms pile up at the edges of the chain. This is illustrated in Fig. 7(c), where we plot the individual atomic state populations for the most subradiant state of $N=20$ atoms, in the $|F_{z}|=|F_{z}^{\text{max}}|-2$ manifold. Here, the two defect atoms, characterized by a large population in state $\mathinner{|{1}\rangle}$, are positioned as the two outermost atoms at either of the chain ends.

Appendix B Effect of disorder

In Fig. 8 we plot the decay rate of the most subradiant state in the (a) $|F_{z}|=|F_{z}^{\text{max}}|$ and (b) $|F_{z}|=0$ manifolds as a function of atom number $N$ for different degrees of spatial disorder. Specifically, here we consider classical disorder, where each atom is randomly displaced around its equilibrium position by a normal distribution of standard deviation $\sigma$. We plot the decay rate of the most subradiant eigenstate after performing an average over disorder realizations. The chain lies along $z$ and the disorder is introduced only along this dimension (i.e. there is no disorder in the positions along $x,y$), as disorder in the perpendicular directions would break the conservation of angular momentum projection, which would make the calculations unmanageable. While increasing disorder leads to an average increase in the decay rate, it is clear from the figure that there are still deeply subradiant states even in the presence of significant fluctuations ($\sigma/d=0.1$). The number of disorder realizations is 200. For (b), the numerical difficulty of the problem prevents us from going beyond $N=10-12$ for significant disorder. In particular, sparse diagonalization methods seem to fail, which requires one to resort to exact diagonalization instead. Nevertheless, there does not seem to be any drastic difference in how it is affected by disorder, as compared to classical subradiance.

We now briefly discuss the case of quantum disorder. Specifically, we consider the that each atom is tightly trapped in the Lamb-Dicke regime, and in the motional ground state of its trap. When an atom is excited, there is the possibility of an inelastic process, where a photon is emitted and a phonon is created at the same time. As shown in Ref. Guimond et al. (2019), this leads to an additional decay rate of each atom of $\Gamma^{\prime}\sim\Gamma_{0}\eta^{2}$, where $\eta$ is the Lamb-Dicke parameter. This decay process is independent (does not depend on correlations), as the final atomic motional state (with one phonon) is distinguishable.