Exponential improvement in photon storage fidelities using subradiance and "selective radiance" in atomic arrays

TL;DR

This paper demonstrates that ordered atomic arrays can significantly suppress spontaneous emission through subradiance and introduces the concept of selective radiance, leading to exponential improvements in photon storage fidelity.

Contribution

It provides a comprehensive theoretical framework for understanding subradiance and introduces selective radiance, enabling efficient photon storage with exponentially improved fidelity.

Findings

Subradiant states in atomic arrays are guided modes emitting only from array ends.

Selective radiance states can radiate into desired channels while suppressing undesired scattering.

Photon storage error scales exponentially better with atom number using these states.

Abstract

A central goal within quantum optics is to realize efficient interactions between photons and atoms. A fundamental limit in nearly all applications based on such systems arises from spontaneous emission, in which photons are absorbed by atoms and then re-scattered into undesired channels. In typical treatments of atomic ensembles, it is assumed that this re-scattering occurs independently, and at a rate given by a single isolated atom, which in turn gives rise to standard limits of fidelity in applications such as quantum memories or quantum gates. However, this assumption can be violated. In particular, spontaneous emission of a collective atomic excitation can be significantly suppressed through strong interference in emission. Thus far the physics underlying the phenomenon of subradiance and techniques to exploit it have not been well-understood. In this work, we provide a comprehensive treatment of this problem. First, we show that in ordered atomic arrays in free space, subradiant states acquire an interpretation in terms of optical modes that are guided by the array, which only emit due to scattering from the ends of the finite chain. We also elucidate the properties of subradiant states in the many-excitation limit. Finally, we introduce the new concept of selective radiance. Whereas subradiant states experience a reduced coupling to all optical modes, selectively radiant states are tailored to simultaneously radiate efficiently into a desired channel while scattering into undesired channels is suppressed, thus enabling an enhanced atom-light interface. We show that these states naturally appear in chains of atoms coupled to nanophotonic structures, and we analyze the performance of photon storage exploiting such states. We find that selectively radiant states allow for a photon storage error that scales exponentially better with number of atoms than previously known bounds.

Equations

$$\nabla \times \nabla \times \mathbf{G}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )-\frac{{\omega }^2}{c^2}ϵ(\mathbf{r},\omega )\mathbf{G}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )=\delta (\mathbf{r}-{\mathbf{r}}^{\prime })1,$$

$${\hat{\mathbf{E}}}^{+}(\mathbf{r})$$

$$\mathcal{H}=\hslash {\omega }_0\sum _{i=1}^N{\hat{\sigma }}_{ee}^i+\hslash \sum _{i,j=1}^N{J}^{ij}{\hat{\sigma }}_{eg}^i{\hat{\sigma }}_{ge}^j,$$

$$\mathcal{L}[{\hat{\rho }}_{\mathrm{A}}]=\sum _{i,j=1}^N\frac{{\Gamma }_{ij}}{2}\left(2{\hat{\sigma }}_{ge}^j{\hat{\rho }}_{\mathrm{A}}{\hat{\sigma }}_{eg}^i-{\hat{\sigma }}_{eg}^i{\hat{\sigma }}_{ge}^j{\hat{\rho }}_{\mathrm{A}}-{\hat{\rho }}_{\mathrm{A}}{\hat{\sigma }}_{eg}^i{\hat{\sigma }}_{ge}^j\right).$$

$$J^{ij}$$

$${\Gamma }^{ij}=\frac{2{\mu }_0{\omega }_0^2}{\hslash }{\wp }^{\ast }\cdot \text{Im}\mathbf{G}({\mathbf{r}}_i,{\mathbf{r}}_j,{\omega }_0)\cdot \wp ,$$

$${\mathcal{H}}_{\mathrm{eff}}=-{\mu }_0{\omega }_0^2\sum _{i,j=1}^N{\wp }^{\ast }\cdot \mathbf{G}({\mathbf{r}}_i,{\mathbf{r}}_j,{\omega }_0)\cdot \wp {\hat{\sigma }}_{eg}^i{\hat{\sigma }}_{ge}^j,$$

$${\mathbf{G}}_0(\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}],$$

$$\frac{J_{\text{k}}}{{\Gamma }_0}$$

$$\frac{{\Gamma }_{\text{k}}}{{\Gamma }_0}$$

$$\frac{J_{k_z}^{||}}{{\Gamma }_0}$$

$$-\mathrm{i}{k}_{0}d{\text{Li}}_2(e^{\mathrm{i}(k_0+k_z)d})-\mathrm{i}{k}_{0}d{\text{Li}}_2(e^{\mathrm{i}(k_0-k_z)d})],$$

$$\frac{J_{k_z}^{\perp }}{{\Gamma }_0}$$

$$-\mathrm{i}{k}_{0}d{\text{Li}}_2(e^{\mathrm{i}(k_0+k_z)d})-\mathrm{i}{k}_{0}d{\text{Li}}_2(e^{\mathrm{i}(k_0-k_z)d})$$

$$+k_0^2{d}^2\text{Ln}(1-e^{\mathrm{i}(k_0+k_z)d})+k_0^2{d}^2\text{Ln}(1-e^{\mathrm{i}(k_0-k_z)d})],$$

$$\frac{{\Gamma }_{k_z}^{||}}{{\Gamma }_0}$$

$$\frac{{\Gamma }_{k_z}^{\perp }}{{\Gamma }_0}$$

$$\frac{{\Gamma }_{\text{k}}^{||}}{{\Gamma }_0}$$

$$\frac{{\Gamma }_{\text{k}}^{\perp }}{{\Gamma }_0}$$

$$c_{\mathrm{ans},k_n}^j$$

$${\mathcal{H}}_{1D}=-\mathrm{i}\frac{\hslash {\Gamma }_{1D}}{2}\sum _{i,j=1}^N{e}^{\mathrm{i}{k}_{1D}|z_i-z_j|}{\hat{\sigma }}_{eg}^i{\hat{\sigma }}_{ge}^j,$$

$${\mathcal{H}}^{\prime }=-\frac{3\pi \hslash {\Gamma }_0}{k_0}\sum _{i,j=1}^N{G}_{\rho \rho }^{\prime }({\mathbf{r}}_i,{\mathbf{r}}_j,{\omega }_0){\hat{\sigma }}_{eg}^i{\hat{\sigma }}_{ge}^j$$

$${\mathcal{H}}_{\text{indep}}^{\prime }=\hslash (J^{\prime }-\mathrm{i}{\Gamma }^{\prime }/2)\sum _{j=1}^N{\hat{\sigma }}_{ee}^j.$$

$${\mathcal{H}}_{\text{drive}}=-\hslash \Delta \sum _{i=1}^N{\hat{\sigma }}_{ee}^i-\hslash \sum _{i=1}^N\left(\Omega e^{\mathrm{i}{k}_{1D}{z}_i}{\hat{\sigma }}_{eg}^i+\text{h.c.}\right).$$

$${\dot{c}}_e^j(t)=$$

$$-$$

$${\hat{\mathbf{E}}}^{+}(\mathbf{r})$$

$${\hat{E}}_{\text{1D,R}}^{+}(z)$$

$$+\mathrm{i}\frac{{\Gamma }_{1D}}{2}\sum _{j=1}^N{e}^{\mathrm{i}{k}_{1D}(z-z_j)}\Theta (z-z_j){\hat{\sigma }}_{ge}^j,$$

$${\hat{E}}_{\text{1D,L}}^{+}(z)$$

$$+\mathrm{i}\frac{{\Gamma }_{1D}}{2}\sum _{j=1}^N{e}^{\mathrm{i}{k}_{1D}(z_j-z)}\Theta (z_j-z){\hat{\sigma }}_{ge}^j,$$

$$T_{\text{indep}}\simeq \text{exp}\left[\frac{-D}{1+4(\Delta -J^{\prime }{)}^2/{\Gamma }^{\prime 2}}\right],$$

$${\dot{c}}_e^j(t)=$$

$$+\mathrm{i}\Omega e^{\mathrm{i}{k}_{1D}{z}_j}+\mathrm{i}\frac{3\pi {\Gamma }_0}{k_0}\sum _{i=1}^N{G}_{\rho \rho }^{\prime }({\mathbf{r}}_i,{\mathbf{r}}_j,{\omega }_{\text{p}})c_e^i(t).$$

$${\mathcal{H}}_{\mathrm{c}}=-\hslash \sum _{i=1}^N{\Delta }_{\mathrm{s}}{\hat{\sigma }}_{ss}^i-\hslash \sum _{i=1}^N{\Omega }_c(z_i)\left({\hat{\sigma }}_{es}^i+{\hat{\sigma }}_{se}^i\right),$$

$${\dot{c}}_e^i(t)$$

Full text

I Introduction

The ability to achieve controlled, deterministic interactions between photons and atomic media constitutes an important resource in applications ranging from quantum information processing to metrology. As single photons and atoms typically do not interact efficiently, a common approach has been to employ atomic ensembles, where the interaction probability with a given optical mode is enhanced via a large number of atoms Hammerer et al. (2010). Atomic ensembles have enabled a number of spectacular proof-of-principle demonstrations of quantum protocols, such as coherent photon storage and quantum memories for light Hau et al. (1999); Fleischhauer et al. (2005); Hammerer et al. (2010), entanglement generation between light and atomic spins Julsgaard et al. (2001), nonlinear interactions between photons at the level of individual quanta Pritchard et al. (2013); Peyronel et al. (2012); Gorniaczyk et al. (2014); Tiarks et al. (2014), and quantum-enhanced metrology Wasilewski et al. (2010); Leroux et al. (2010); Sewell et al. (2012). It has also been proposed that such systems could lead to exotic many-body physics, such as strongly correlated photon “gases” Noh and Angelakis (2016).

A fundamental limitation in nearly all such possibilities arises from spontaneous emission, wherein photons in a desired optical mode (e.g., a Gaussian input beam) that facilitate the process are absorbed by the atoms and then re-scattered into other inaccessible modes or channels. Within the context of quantum light-matter interfaces based upon atomic ensembles, it is typically assumed that spontaneous emission occurs independently, and at the same rate given by a single, isolated atom. In that case, the infidelity or error arising from spontaneous emission for a desired process typically decreases with the “optical depth” $D$ of the medium as $1/D$ or slower. The optical depth is given by $D\sim(\lambda_{0}^{2}/A_{\text{eff}})N$, where $N$ is the atom number, and $\lambda_{0}^{2}/A_{\text{eff}}$ represents the interaction probability between a single atom and a single photon in the preferred optical mode ($\lambda_{0}$ being the wavelength associated with the atomic transition and $A_{\text{eff}}$ the beam area). Intuitively, the $1/D$ (or $1/N$) scaling directly reflects the fact that a given atom is assumed to succeed or fail independently, and that the success is enhanced by the number of atoms involved.

Technically, however, the assumption of independent emission cannot strictly be correct. In particular, as scattering is a wave phenomenon, the emission into other directions may exhibit collective interference. In fact, the possibility that an atomic ensemble can experience a significantly enhanced radiation rate via interference (“superradiance”) was already pointed out in the seminal work of Dicke Dicke (1954), and has been thoroughly studied for decades Gross and Haroche (1982). The complementary phenomenon of subradiance, in which photon emission becomes highly suppressed, has also been theoretically studied Prasad and Glauber (2000); Svidzinsky et al. (2010); Jenkins and J.Ruostekoski (2012); Plankensteiner et al. (2015); Bettles et al. (2015, 2016a); Facchinetti et al. (2016); Sutherland and Robicheaux (2016); Zoubi and Ritsch (2010, 2011); Kornovan et al. (2016), and even observed in recent experiments DeVoe and Brewer (1996); Guerin et al. (2016); SBF . Clearly the possibility to enhance atom-light interfaces by suppressing unwanted emission is a tantalizing one, and has started to gain theoretical interest Plankensteiner et al. (2015). However, finding protocols where subradiance clearly improves the scaling of errors remains an elusive task, in part because techniques to efficiently address subradiant states remain poorly developed.

In this paper, we provide a comprehensive description of subradiance in the case where atoms form ordered arrays. We also present an explicit construction of a protocol exploiting suppressed emission into undesired channels, which enables an exponential improvement in infidelity as a function of atom number over previously known bounds. Our main results are summarized as follows:

• •

We first consider infinite 1D or 2D arrays of atoms, which consist of an electronic ground state $\mathinner{|{g}\rangle}$ and excited state $\mathinner{|{e}\rangle}$ that couple to light through a dipole transition. Examining the case of a single collective excitation, we find that a set of perfectly subradiant states with zero decay rate emerge, which can be interpreted as optical “guided modes.” Specifically, in exact analogy to guided modes of conventional optical fibers or photonic structures, the spin-wave excitations that constitute these subradiant states have associated wave vectors that are mismatched from free-space radiation fields, which consequently prevents the decay of energy from these states.

• •

In the case of a finite array, a set of single-excitation collective atomic modes can exhibit decay rates which are polynomially suppressed with atom number $N$. The finite decay rate can be understood as emerging from scattering of guided excitations into free space through the boundaries of the array.

• •

We go beyond the most frequently studied case of subradiance within a single-excitation manifold, and investigate the nature of multi-excitation subradiant modes. Specifically, we show that subradiance is largely destroyed when excitations spatially overlap, as the scattering of two excitations generates many wave vectors that couple to free-space radiation due to the “hard core” nature of spins. In 1D arrays, we find that a “fermionic” ansatz works well to describe multi-excitation subradiant modes, where these multi-excitation states are constructed from anti-symmetric combinations of single-excitation subradiant modes in order to enforce a spatial repulsion (i.e., “Pauli exclusion”) of excitations. These states preserve the same polynomial suppression of decay rate with atom number for any low density of excitations.

• •

Having elucidated the salient properties of subradiant states in free-space atomic arrays, we introduce the new concept of “selectively radiant” states. In particular, while subradiant states couple weakly to all electromagnetic modes, to realize an efficient atom-light interface it is instead desirable to construct states that are simultaneously superradiant to a preferred photonic mode and subradiant to all the others. We show that one natural way to achieve such a scenario is by coupling an atomic array to the guided modes of a nanophotonic structure, such as an optical nanofiber Vetsch et al. (2010); Goban et al. (2012); Grover et al. (2015); Sørensen et al. (2016); Corzo et al. (2016). As the wave vectors of the guided modes of the structure itself are mismatched from free-space radiation, it becomes possible for a set of atomic spin waves to efficiently couple to these guided modes, while retaining a suppressed coupling to free-space modes. We analyze the specific protocol of photon storage Lukin (2003); Fleischhauer et al. (2005) using an atomic array coupled to a nanofiber Sayrin et al. (2015); Gouraud et al. (2015), and find numerically a storage infidelity that is exponentially small in the atom number or optical depth, $\sim\exp(-D)$. This scaling represents an exponential improvement over the best previously established error bound of $\sim 1/D$ Gorshkov et al. (2007a, b), derived assuming that emission into undesired modes is independent.

This article is structured as follows. In Sec. II we begin by introducing a theoretical framework for atom-light interactions that does not invoke the typical assumption of independent atomic emission, and instead formally and exactly describes collective emission and interactions of atoms via photon fields. In Sec. III we apply this formalism to investigate single- and multi-excitation subradiant states in atomic arrays, with the main results having already been summarized above. In Sec. IV we present the idea of selectively radiant states, and analyze the efficiency of a quantum memory consisting of a chain of atoms close to a nanofiber. Finally, in Sec. V we discuss possible implementations and other photonic platforms for observing subradiant physics. An outlook is provided in Sec. VI.

II Spin model

Here we introduce a theoretical formalism to describe the fully quantum interaction between atoms and radiation fields, which is valid in the presence of any linear, isotropic, dielectric media. This rather general formalism will enable us to equally treat the case of atomic arrays in free space (Sec. III), or interacting via the guided modes of an optical fiber (Sec. IV). In particular, we present a model in which the field is integrated out and the dynamics of the atomic internal (“spin”) degrees of freedom follow a master equation that only depends on atomic operators. Moreover, once the time evolution of the atoms is solved for, one can recover the field at any point in space by means of a generalized input-output equation.

The first step to describe how atoms couple to radiation is to quantize the electromagnetic field. The traditional approach involves explicitly finding a normal mode decomposition of the fields, and associating bosonic annihilation and creation operators to each mode. This is well-suited to cases where a limited number of modes are assumed to be relevant (such as a high-Q cavity). In our case, though, as we want to exactly capture collective effects in spontaneous emission involving all modes, such an approach becomes unwieldy (as in free space) Gross and Haroche (1982); Bhat and Sipe (2006) or impossible, such as for complex dielectric structures. We require a more general technique that allows us to treat these situations. Such a framework was developed by Welsch and coworkers Gruner and Welsch (1996); Dung et al. (2002); Buhmann and Welsch (2007); Buhmann (2012), and is based on the classical electromagnetic Green’s function (or Green’s tensor).

The Green’s function ${\bf G}({\bf r},{\bf r}^{\prime},\omega)$ is the fundamental solution of the electromagnetic wave equation, and obeys Novotny and Hecht (2006):

[TABLE]

where $\epsilon({\bf r},\omega)$ is the position- and possibly frequency-dependent relative permittivity of the medium. The Green’s function physically describes the field at point ${\bf r}$ due to a normalized, oscillating dipole at ${\bf r}^{\prime}$. ${\bf G}_{\alpha\beta}$ is a tensor quantity ($\{\alpha,\beta\}=\{x,y,z\}$), as $\alpha$ and $\beta$ refer to the possible orientations of the field and dipole, respectively. Here, we will deal with cases where the Green’s function can be solved analytically, but for more complex structures it is also possible to obtain it numerically  Oskooi et al. (2010); Manga Rao and Hughes (2007); Hood et al. (2016). In the following, we introduce a prescription of how to write down an equation that relates the field and the atomic coherence operators, built upon the intuition provided by classical physics. For a more formal derivation of the field quantization, we refer the reader to Refs. Gruner and Welsch (1996); Dung et al. (2002); Buhmann and Welsch (2007); Buhmann (2012); Asenjo-Garcia et al. (2017).

In the frequency domain, the analogous classical problem that one would like to solve is to find the total field ${\bf E}({\bf r},\omega)$ at point ${\bf r}$, given a known input field ${\bf E}_{\rm p}({\bf r},\omega)$ and a collection of $N$ polarizable dipoles ${\bf p}_{j}(\omega)$ located at ${\bf r}_{j}$, which are excited by the fields and re-scatter light themselves. The values of ${\bf p}_{j}(\omega)$ are not known a priori, since they depend on the polarizability and the total field at ${\bf r}_{j}$ [solving for ${\bf p}_{j}(\omega)$ will be discussed in following steps]. As the field at any given point in space is just the sum of the external or driving field and the field re-scattered by the dipoles, we find ${\bf E}({\bf r},\omega)={\bf E}_{\rm p}({\bf r},\omega)+\mu_{0}\omega^{2}\,\sum_{j=1}^{N}{\bf G}({\bf r},{\bf r}_{j},\omega)\cdot{\bf p}_{j}(\omega)$, where $\mu_{0}$ is the vacuum permeability.

The question is how to translate this classical equation into an equation for quantum operators. In fact, the quantum nature of the field is inherited from the quantum properties (e.g., correlations and fluctuations) of the sources, while the field propagation remains the same as both the quantum and classical fields obey Maxwell’s equations. Therefore, the above equation is valid for quantum fields, but replacing ${\bf p}_{j}(\omega)$ by the dipole moment operator $\hat{{\bf p}}_{j}(\omega)$, and ${\bf E}({\bf r},\omega)$ by the field operator $\hat{{\bf E}}^{+}({\bf r},\omega)$, where the superscript refers to the positive-frequency component. In the case that the quantum dipoles are atoms, one can make a further approximation, taking advantage of the fact that an atom only has a significant optical response in a narrow bandwidth around its resonance frequency $\omega_{0}$. Thus, one is able to approximate ${\bf G}({\bf r},{\bf r}_{j},\omega)$ by ${\bf G}({\bf r},{\bf r}_{j},\omega_{0})$, which allows the Fourier transform of the equation to become local in time. Then, one arrives to the generalized input-output equation in the time domain, which reads Dung et al. (2002); Caneva et al. (2015); Xu and Fan (2015)

[TABLE]

To obtain the above expression, we have made use of the fact that $\hat{{\bf p}}_{j}=\boldsymbol{\wp}^{*}\,\hat{\sigma}^{j}_{eg}+\boldsymbol{\wp}\,\hat{\sigma}^{j}_{ge}$, where $\hat{\sigma}_{eg}^{j}=\mathinner{|{e_{j}}\rangle}\mathinner{\langle{g_{j}}|}$ is the atomic coherence operator between the ground and excited states of atom $j$, and $\boldsymbol{\wp}$ is the dipole matrix element associated with that transition. This equation is valid in the Markovian regime, where the dispersion in the Green’s function can be neglected and the replacement of ${\bf G}({\bf r},{\bf r}_{j},\omega)$ by ${\bf G}({\bf r},{\bf r}_{j},\omega_{0})$ is well-founded. For this to be true, two conditions have to be fulfilled. First, the retardation arising from the physical distance between atoms can be ignored Chang et al. (2012); Shi et al. (2015). For atoms in free space, this means that they should sit much closer than the length of a spontaneously emitted photon ($\lesssim 1$ meter). Second, the electromagnetic environment itself should not have very narrow-bandwidth features (e.g., one must avoid the strong coupling regime of cavity QED Thompson et al. (1992)).

What remains now is to solve for the dipoles (in this case, $\hat{\sigma}_{ge}$) themselves. We do so by writing down Heisenberg-Langevin equations for the atomic internal degrees of freedom, starting from the full atom-field Hamiltonian. Intuitively, the atomic spin $\hat{\sigma}_{eg}^{i}$ will be driven by the quantum field at position ${\bf r}_{i}$. However, as the field itself only depends on other atoms via the input-output equation, the atomic dynamics can be fully derived from an equivalent master equation of the form $\dot{\hat{\rho}}_{\rm A}=-({\rm i}/\hbar)\,[\mathcal{H},\hat{\rho}_{\rm A}]+\mathcal{L}[\hat{\rho}_{\rm A}]$ Meystre and Sargent (2007), where $\hat{\rho}_{\rm A}$ is the atomic density matrix, and the Hamiltonian and Lindblad operators read

[TABLE]

In the above expressions, the rates for coherent and dissipative interactions between atoms $i$ and $j$ are respectively given by

[TABLE]

where the sign of $J^{ij}$ is taken to be opposite to that of Refs. Gruner and Welsch (1996); Dung et al. (2002); Buhmann and Welsch (2007); Buhmann (2012); Hood et al. (2016); Asenjo-Garcia et al. (2017). In the above Hamiltonian, we have neglected Casimir interactions between ground-state atoms (of the form $\hat{\sigma}_{gg}^{i}\hat{\sigma}_{gg}^{j}$), as their spatial decay is very fast ($\sim 1/d^{6}$ in free space, $d$ being the inter-atomic distance) Buhmann (2012).

The dynamics under the master equation can analogously be described in the quantum jump formalism of open systems Markel and Sarychev (2007). In this formalism, the atomic wave function evolves deterministically under an effective non-Hermitian Hamiltonian that reads $\mathcal{H}=\hbar\omega_{0}\sum_{i=1}^{N}\hat{\sigma}_{ee}^{i}+\mathcal{H}_{\rm eff}$, with

[TABLE]

along with stochastically applied “quantum jump” operators to account for the population recycling term ($\hat{\sigma}_{ge}^{j}\hat{\rho}_{\text{A}}\hat{\sigma}_{eg}^{i}$) of Eq. (3b). While $\mathcal{H}_{\rm eff}$ just describes the interaction of atoms through emission and re-absorption of photons, one can directly add other terms to the Hamiltonian to account for external driving fields.

To conclude, we point out that although the full formalism above has only been rigorously and generally developed in recent years, many aspects have long been used within atomic physics and quantum optics. For example, for a single atom or other quantum emitter, the spin model becomes trivial and just yields the total spontaneous emission rate. Thus, the calculation of enhancement of spontaneous emission near dielectric structures is standardly reduced to the calculation of the Green’s function Chance et al. (1978); Lodahl et al. (2004); Englund et al. (2005). Alternatively, such equations are often used to model the optical response of dense three-dimensional atomic gases Schilder et al. (2016); Jennewein et al. (2016a, b); Bromley et al. (2016); Guerin and Kaiser (2017).

III Free space: Subradiant states

We now apply the spin model we describe in the previous section to investigate the properties of subradiant states associated with ordered atomic arrays in free space. Recently, the peculiar linear optical properties of periodic atomic arrays have started to attract interest Bettles et al. (2015, 2016a, 2016b); E. Shahmoon and Yelin (2017); Kornovan et al. (2016); Sutherland and Robicheaux (2016); Zoubi and Ritsch (2011); Plankensteiner et al. (2015); Haakh et al. (2016); Jenkins and J.Ruostekoski (2012); Facchinetti et al. (2016). This includes the identification of guided modes supported by infinite arrays Zoubi and Ritsch (2010); E. Shahmoon and Yelin (2017); Sutherland and Robicheaux (2016), and states with very long lifetimes in finite arrays Bettles et al. (2015, 2016a, 2016b); Kornovan et al. (2016); Sutherland and Robicheaux (2016); Zoubi and Ritsch (2011); Plankensteiner et al. (2015); Jenkins and J.Ruostekoski (2012); Haakh et al. (2016); Facchinetti et al. (2016). Here, we provide a clear and intuitive connection between the existence of guided modes in infinite arrays and subradiant states in a finite system. We provide conditions for the lattice constants in 1D and 2D that enable single-excitation guided Bloch modes with zero decay rate to emerge, which are decoupled from free-space radiation due to wave vector mismatch. We then analyze a single excitation in a finite lattice, and show how the guided modes acquire a non-zero decay rate due to scattering into electromagnetic radiation at the system boundaries. We also analyze the scaling of the decay rates with system size and elucidate the spatial structure of subradiant states. Finally, we go beyond previous studies of single-excitation subradiance (where the atoms can equivalently be treated as classical dipoles) to the rich physics of the multi-excitation case. In particular, in one dimension, we show that multi-excitation subradiant states exist for any low density of excitations, and that their wave functions have fermionic character.

The atoms are assumed to be tightly trapped, so that we can treat the positions of the particles as classical points rather than dynamical variables. In this situation, we substitute in Eqs.(2-5) the free-space Green’s tensor ${\bf G}({\bf r}_{i},{\bf r}_{j},\omega_{0})={\bf G}_{0}({\bf r}_{ij},\omega_{0})$, with ${\bf r}_{ij}={\bf r}_{i}-{\bf r}_{j}$. Here ${\bf G}_{0}({\bf r},\omega_{0})$ is the solution to Eq.(1) when setting $\epsilon({\bf r},\omega)=1$, and can be written as:

[TABLE]

where $r=|{\bf r}|$ and $k_{0}=2\pi/\lambda_{0}=\omega_{0}/c$ is the wave number corresponding to the atomic transition energy. For a single atom, evaluating Eq. (4b) simply reproduces the well-known vacuum emission rate $\Gamma^{ii}=\Gamma_{0}$, where $\Gamma_{0}=\omega_{0}^{3}|\boldsymbol{\wp}|^{2}/3\pi\hbar\epsilon_{0}c^{3}$. The single-atom energy shift $J^{ii}$ in Eq.(4a) arising from ${\bf G}_{0}$ formally yields a divergence and will be set to zero in what follows, as it should be incorporated into a re-normalized resonance frequency $\omega_{0}$. In Sec. IV, the Green’s function of a nanofiber will be decomposed into a free-space and a scattered component, ${\bf G}={\bf G}_{0}+{\bf G}_{\text{sc}}$, where ${\bf G}_{\text{sc}}$ does produce a finite, observable contribution to $J^{ii}$.

For concreteness, we will restrict ourselves to the following lattice geometries: one-dimensional (1D) linear chains and closed circular rings, two-dimensional (2D) square and three-dimensional (3D) cubic lattices. However, it should become clear that the underlying principles should be general to other lattice structures as well. In the following, the number of atoms and lattice constant are denoted by $N$ and $d$, respectively.

III.1 Infinite Lattice (Single excitation)

Let us consider first a perfectly ordered infinite array of atoms. Despite the infinite lattice being unrealistic, it provides insight into the problem thanks to its mathematical simplicity. In this case, the system is perfectly translationally invariant by any lattice vector displacement, and thus, both atomic and electromagnetic eigenmodes must obey Bloch’s theorem.

For a single excitation stored in the system, the eigenstates of the effective atomic Hamiltonian of Eq.(5) are spin waves, with well-defined quasi-momentum k, which can always be chosen to be within the first Brillouin zone. For such states, whose creation operators can be written as $S^{\dagger}_{\textbf{k}}=N^{-1/2}\sum_{j}e^{{\rm i}\textbf{k}\cdot{\bf r}_{j}}\hat{\sigma}_{eg}^{j}$, the single excitation is delocalized and shared in a coherent way among all the atoms. Classically, these states are analogous to oscillating dipoles where the phase of dipole $j$ is given by $e^{{\rm i}\textbf{k}\cdot{\bf r}_{j}}$.

As the Bloch modes are eigenstates of the effective Hamiltonian, they satisfy $\mathcal{H}_{\rm eff}\,S^{\dagger}_{\textbf{k}}\mathinner{|{g}\rangle}^{\otimes N}=\hbar(J_{\textbf{k}}-{\rm i}\Gamma_{\textbf{k}}/2)\,S_{\textbf{k}}^{\dagger}\mathinner{|{g}\rangle}^{\otimes N}$. Here $J_{\textbf{k}}$ and $\Gamma_{\textbf{k}}$ are real quantities and can be identified as the frequency shift of mode k (relative to the bare atomic frequency $\omega_{0}$) and the decay rate, respectively. One can readily show that, in terms of the single-atom spontaneous emission rate $\Gamma_{0}$, they are given by:

[TABLE]

where $\tilde{{\bf G}}_{0}(\textbf{k})=\sum_{j}e^{-{\rm i}\textbf{k}\cdot{\bf r}_{j}}{\bf G}_{0}({\bf r}_{j})$ is the discrete Fourier transform of the free-space Green’s tensor.

We will now show that, when the atoms are placed at close enough distances, dipole-dipole interactions can dramatically modify the decay rates of collective states. As the simplest case, let us consider an infinite one-dimensional chain of atoms first, oriented along the $\hat{z}$ direction. In that case, the wave vector $k_{z}$ constitutes an index for the modes, and one can consider the dispersion relation of frequency $\omega(k_{z})=\omega_{0}+J_{k_{z}}$ versus $k_{z}$. For a periodic structure, regardless of the system details, one expects the dispersion relation to exhibit general characteristics [see Fig. 1(a)]. First, and as mentioned before, $k_{z}$ is only uniquely defined within the first Brillouin zone ($|k_{z}|\leq\pi/d$) and thus, it suffices to plot the dispersion relation in that region. Second, it is helpful to draw the “light line”, i.e., the dispersion relation $\omega=c|k_{z}|$ corresponding to light propagating in free space along the $\hat{z}$ direction [dashed line of Fig. 1(a)]. Physically, the light line is significant because it separates states of very different character, as we now describe.

To see this, let’s consider the field generated by a spin-wave excitation, which is given by Eq.(2) under the replacement $\hat{\sigma}_{ge}^{j}\rightarrow e^{{\rm i}k_{z}z_{j}}$ (it is sufficient to consider the limit of classical dipoles for this argument). One can always expand the field ${\bf E}({\bf r})$ in terms of plane wave components, ${\bf E}({\bf r})=\sum_{q_{z},{\bf q_{\perp}}}{\bf E}_{q_{z},{\bf q_{\perp}}}e^{{\rm i}q_{z}z+{\rm i}{\bf q_{\perp}}\cdot{\bf r_{\perp}}}$. The state is clearly of Bloch’s form, and thus, only a discrete set of wave vectors $q_{z}=k_{z}+g_{z}$ ($g_{z}$ being any reciprocal lattice vector) will contribute. At the same time, the wave equation requires that the axial and perpendicular components of the wave vector satisfy $(q_{z})^{2}+{\bf q}_{\perp}\cdot{\bf q}_{\perp}=(\omega/c)^{2}$. Thus, one can readily verify that a spin wave outside the light line ($|k_{z}|>\omega/c$) has an associated electromagnetic field composed of axial wave vectors $|q_{z}|>\omega/c$. This in turn implies that ${\bf q_{\perp}}$ is imaginary, and the field is guided and decays evanescently away from the structure. Therefore, these guided modes are decoupled from all optical modes propagating in free space, and their inability to radiate away energy leads to perfect subradiance (exactly zero decay rate). Conversely, modes within the light line are generally unguided and can radiate energy out to infinity.

The concepts outlined above regarding the separation of the dispersion relation into guided and radiative regions are actually quite general, and well-known in the context of periodically modulated dielectric waveguides (“photonic crystals” Joannopoulos et al. (2008)). An atomic chain might appear quite different physically, but mathematically the same set of principles apply. Furthermore, while it is difficult to prove independent of lattice geometry and atomic level configuration, one would generally expect that for atoms any guided modes would occur within a narrow bandwidth (on the order of the atomic transition linewidth $\Gamma_{0}/2\pi\lesssim 10{\rm~{}MHz}$) around the resonance frequency ($\omega_{0}/2\pi\sim 300{\rm~{}THz}$), where the atoms have a significant optical response. Thus, in Fig. 1(a) the band structure will appear rather flat. Then, a sufficient condition for guided modes to exist in an atomic chain is essentially that the light line intersects the edge of the Brillouin zone $k_{z}=\pi/d$ at a frequency $\omega(k_{z})$ greater than the atomic resonance. This condition can be rewritten as a condition on the lattice constant $d<\lambda_{0}/2$ required to support guided modes.

Equipped with this general intuition, we now quantitatively investigate the dispersion relation for the 1D infinite chain of two-level atoms with polarization parallel or transverse to the array. The collective frequency shifts are derived in greater detail in Appendix A, and read:

[TABLE]

where $\textrm{Li}_{n}(x)$ is the PolyLogarithm of order $n$. These expressions are plotted in Fig. 1(b), for the particular value of $d/\lambda_{0}=0.2$. Here, the light line is indicated as before by a dashed line, but since $\Gamma_{0}/\omega_{0}\sim 10^{-8}$, it appears essentially as a vertical line.

As anticipated, we can see in this figure that the bands occupy only a narrow bandwidth around the resonance frequency, except close to the light line for transverse polarization. The exact shape of the bands depends on the value of $d/\lambda_{0}$ and the polarization direction, and for instance, the effective mass at the zone edge ($|k_{z}|=\pi/d$) is negative (positive) for parallel (transverse) polarization. Exactly at the light line the expression for $J_{k_{z}}^{\perp}$ ($J_{k_{z}}^{||}$) becomes non-analytic and diverges (has a derivative that diverges).

The collective decay rates can also be analytically derived (see Appendix A):

[TABLE]

These summations run over reciprocal lattice vectors that satisfy $|g_{z}+k_{z}|\leq k_{0}$. That is, only the diffracted waves enclosed within the light line will contribute to the decay rate. When $|k_{z}|>k_{0}$, there are no values of $g_{z}$ satisfying the above condition. Thus the decay rates are zero and we mathematically recover the result previously anticipated – modes beyond the light line are perfectly guided without radiative losses. The decay rates are plotted in Fig. 1(c). As we can see from the expressions above, at the light line the state can be subradiant or radiant depending on the polarization direction. This results in a discontinuity at the light line for transverse polarization.

A similar set of results can be obtained for a 2D array. Considering a square lattice in the $\hat{y}$-$\hat{z}$ plane [Fig. 2(a)], the corresponding first Brillouin zone for Bloch wave vectors extends over the region $|k_{y}|,|k_{z}|\leq\pi/d$, as shown in Fig. 2(b). The set of electromagnetic fields propagating in the plane at the atomic frequency $\omega_{0}$ have a wave vector of magnitude $k_{0}$ which defines a circle centered around the origin in k-space, as illustrated in Fig. 2(b). Similar to 1D, a sufficient condition for spin-wave excitations to be guided is that the wave vector lies outside of this circle. It should be noted that the longest “distance” in the first Brillouin zone from the origin extends along the diagonal, and has magnitude $k^{\text{max}}=\sqrt{2}\pi/d$. Thus, in 2D, guided modes exist as long as $k_{0}<k^{\text{max}}$, which translates into a maximum allowed lattice constant $d/\lambda_{0}=1/\sqrt{2}$.

Analogous to the 1D case, we can obtain closed mathematical expressions for the decay rates in the 2D lattice. They are given by:

[TABLE]

from which we recover again the important result that Bloch states with $|\textbf{k}|>k_{0}$ do not radiate out to infinity. For these states, the electromagnetic field is now confined within the plane, and evanescently decays away from the lattice in the transverse direction. The decay rates are plotted in Fig. 2(c)-(d) for atomic polarizations along $\hat{z}$ and $\hat{x}$ directions, for the particular value of $d/\lambda_{0}=0.2$, which defines the light line as the circle $k_{0}=0.4\pi/d$, beyond which the decay rate is exactly zero.

We would like to remark that the previous considerations are valid regardless of the specific atomic structure, provided that the atom in question only contains a single ground state (see Sec. V for a discussion of the subtleties associated with a ground-state manifold).

The previous analysis for the 2D lattice provides another interesting result. As Fig. 2(d) shows, for transverse atomic polarization in a 2D square lattice, subradiance can emerge not only outside the light line, but also at the center of the Brillouin zone, that is, for Bloch states with quasi-momentum $\textbf{k}\sim(0,0)$. Physically, the origin of this effect can be understood as follows. On one hand, and as previously discussed, for $d/\lambda_{0}<1$ the field created by such a state is generally evanescent at all diffraction orders except for the component ${\bf g}=0$ [c.f. Eq.(10b)], which corresponds to a plane wave propagating perpendicularly to the atomic plane. On the other hand, the state $\textbf{k}=(0,0)$ corresponds to an array of dipoles that are in phase. However, dipoles oscillating in phase and perpendicularly to the atomic plane are forbidden to radiate energy in the perpendicular direction, and thus, the state must be subradiant. In contrast, as soon as $d/\lambda_{0}>1$, there will be other ${\bf g}$ components that are not evanescent, yielding a radiative state. We note that, although only the state $\textbf{k}=(0,0)$ has exactly zero decay rate, other modes around this point will also show a strong suppression in the emission rate relative to $\Gamma_{0}$.

III.2 Finite Lattice (Single excitation)

In this section, we analyze the decay rates and spatial properties of single-excitation eigenstates, for a lattice of finite size. We show that all eigenstates now acquire a non-zero decay rate, and subradiant states can be identified as those for which the rate is suppressed with increasing system size. The small value of the decay can be interpreted as arising from the finite system boundaries, which scatter a mode propagating in the bulk into free space.

**1D Linear Chain.

**In the following, we consider a finite chain of atoms along $\hat{z}$, with a linear polarization along the chain (unless otherwise stated). However, a similar set of conclusions is obtained for the transverse polarization case.

Scaling of the most subradiant decay rates with system size.– The effective atomic Hamiltonian of Eq.(5) conserves the excitation number in the system, and thus, it can be diagonalized in blocks with fixed excitation number. Before proceeding futher, we discuss a technical but important point. Since the effective Hamiltonian is non-Hermitian, in general the eigenstates will not be orthogonal in the standard quantum mechanical sense (i.e., two eigenstates $\mathinner{|{\psi_{i}}\rangle}$ and $\mathinner{|{\psi_{j}}\rangle}$ will not satisfy $\langle{\psi_{i}}|{\psi_{j}}\rangle=\delta_{ij}$) Asenjo-Garcia et al. (2017). The infinite lattice case presented an exception, as Bloch’s theorem is still enforced. While this implies that general quantum mechanical rules, such as for eigenstate decompositions of states and observables, do not apply, we will nonetheless investigate the properties of the eigenstates further. This is physically motivated as they still represent non-evolving states under the Hamiltonian (aside from an overall phase and amplitude); thus, for example, they might be expected to shed light on how a general state behaves at long times.

We consider the case of a 1D chain of $N$ atoms with lattice constant $d$, for which numerical diagonalization of $\mathcal{H}_{\textrm{eff}}$ in the one-excitation manifold produces $N$ eigenstates (denoted by $\mathinner{|{\psi_{\xi}}\rangle}$, $1\leq\xi\leq N$) and complex eigenvalues. As in the infinite case, the eigenvalues can be written in the form $J_{\xi}-i\Gamma_{\xi}/2$, with $J_{\xi}$ and $\Gamma_{\xi}$ representing the frequency shift relative to $\omega_{0}$ and decay rate, respectively. For concreteness, the eigenstates will be ordered in increasing decay rate, such that $\xi=1$ represents the most subradiant state and $\xi=N$ the most radiant one. To start understanding the properties of this system we fix the atomic number $N=50$ and change the lattice constant $d$. For each value of $d$, we diagonalize $\mathcal{H}_{\textrm{eff}}$, and obtain the $N$ different values for the decay rates and frequency shifts associated with each of the collective modes. Figure 3(a) shows the resulting single-excitation decay rates $\Gamma_{\xi}$ for each collective mode, normalized by the free-space single-atom emission rate $\Gamma_{0}$, in the case of atomic polarization parallel to the chain. In this plot, a vertical cut at a fixed value of $d/\lambda_{0}$ contains the $N$ different values of $\Gamma_{\xi}$. As expected, for large interparticle distances, the collective decay rates tend to the spontaneous emission rate of a single atom. As the distance decreases, they are periodically modulated, showing for $d/\lambda_{0}<1/2$ a qualitatively distinct behavior. In this region, the decay rates of some of the modes are dramatically suppressed ($\Gamma_{\xi}/\Gamma_{0}\ll 1$), in accordance with the condition for the emergence of modes with zero decay rate derived in the infinite lattice case.

The subradiant modes in the finite chain are closely related to those derived in the infinite chain. First, having established that subradiant states in the infinite chain correspond to guided modes, the nonzero decay rates in the finite chain can be interpreted as emerging from scattering of these guided modes from the ends of the system. This can be seen in Fig. 3(b), where we have plotted the field intensity in the plane $x=5d$ generated by the most subradiant state when the atomic polarization is parallel to the chain (we choose a distance $x$ offset from the $x=0$ plane containing the atomic chain in order to avoid seeing the divergent near-fields associated with each atom). Clearly, this figure shows that the field vanishes when moving away from the chain transversally, while it is very intense at the tips of the chain, where the spin wave scatters into an outgoing photon. The field intensity is computed from Eq.(2), by taking $\mathinner{\langle{\psi_{1}}|}\hat{{\bf E}}^{-}({\bf r})\hat{{\bf E}}^{+}({\bf r})\mathinner{|{\psi_{1}}\rangle}$. As the input field is vacuum, the intensity only involves calculating two-body correlations $\hat{\sigma}_{eg}^{i}\hat{\sigma}_{ge}^{j}$ of the eigenstate.

While the wave vector $k_{z}$ is strictly a good index for the modes only in the case of the infinite chain, in practice one can also unambiguously associate a distinct, dominant wave vector $k$ with each of the modes $\xi$ in the finite case. Specifically, the discrete Fourier transform of the coefficients that define each mode is peaked around a different value $k_{\xi}$, which can be used to label the state. In particular, let us consider a general single-excitation state, which can be written as $\mathinner{|{\psi_{\xi}}\rangle}=\sum_{j}c_{\xi}^{j}\mathinner{|{e_{j}}\rangle}$, where $\mathinner{|{e_{j}}\rangle}\equiv\hat{\sigma}_{eg}^{j}\mathinner{|{g}\rangle}^{\otimes{N}}$ is defined as the state where atom $j$ is excited while all others are in their ground states. Then, we define the discrete Fourier transform of the associated coefficients $\tilde{c}_{\xi}^{k}\equiv N^{-1/2}\sum_{j}e^{{\rm i}kjd}c_{\xi}^{j}$, for discrete values of $k=2\pi m/Nd$ ($1\leq m\leq N$). For each value of $\xi$, the function $\tilde{c}_{\xi}^{k}$ shows a well defined peak at a distinct value of $k=k_{\xi}$. In Figs. 1(b),(c) (circles) we plot the decay rates $\Gamma_{\xi}$ and energy shifts $J_{\xi}$ of each mode as indexed by the dominant wave vector, for $N=50$ atoms and both transverse and parallel polarizations, overlaid with the infinite lattice result. There is good agreement between them. For the decay rates of the finite chain, the points also correspond to those along a vertical cut in Fig. 3(a), at the fixed value of $d/\lambda_{0}=0.2$.

The exact behavior of $\Gamma_{\xi}$ depends on the microscopic details, such as the polarization of the atoms. For instance, for two-level atoms, the smallest decay rate decreases monotonically as $d/\lambda_{0}\rightarrow 0$, while for transverse polarization it oscillates. Regardless of these details, however, the scaling with $N$ of the few lowest decay rates seems to show a universal behavior, going like $\Gamma_{\xi}/\Gamma_{0}\sim\xi^{2}/N^{3}$. In Fig. 3(c), we show the $1/N^{3}$ scaling for the three lowest eigenstates as a function of $N$, while in Fig. 3(d) we show the $\xi^{2}$ scaling for fixed $N=20$. The scaling with $\xi$ is satisfied for all $\Gamma_{\xi}/\Gamma_{0}\ll 1$. For transverse polarization, there is a particular value of lattice constant (that tends to $d/\lambda_{0}\sim 0.25$ as the atom number increases), for which the decay rates do not follow exactly the scaling with $\xi$. We believe that this is related to the fact that for transverse polarization and $d/\lambda_{0}=0.25$ the band structure becomes flat at the edge of the Brillouin zone. Nevertheless, and as we discuss in Appendix B, the scaling $\Gamma_{\xi}\sim\xi^{2}/N^{3}$ seems to appear rather generically for finite-size, one-dimensional photonic crystal structures.

Ansatz for single-excitation collective modes.– If the chain is finite the single-excitation collective modes are not spin waves with pure wave vector $k_{z}$, and contrary to the infinite lattice case, they are not orthonormal in general. Nevertheless, as discussed earlier, the eigenstates can be characterized by a dominant wave vector $k$ that connects well with the infinite case. Furthermore, we find that the states far from the light line (including the most subradiant modes as well as those where $\textbf{k}\sim 0$) are almost orthonormal, and display a relatively simple spatial structure. This motivates us to find an orthonormal set of functions that approximates well these modes. For an even number of sites $N$, the wave-function coefficients $c_{\xi}^{j}$ of the exact collective modes are close to the orthonormal set of functions defined by wave vector $k_{n}$:

[TABLE]

Here $k_{n}d=\pi n/(N+1)$, $n=1,2,...,N$ and the atomic positions $x_{j}=jd-x_{0}$ ($1\leq j\leq N$) and $x_{0}=(N+1)d/2$. Figure 4(a)-(c) shows the exact coefficients $c_{\xi}^{j}$ for the most subradiant state ($\xi=1$), a state with dominant wave vector near the light line, and the most radiant state ($\xi=N$), for $N=20$ atoms, together with the corresponding ansatz coefficients. The error between the exact wave function and ansatz can be quantified by considering the mismatch in overlap between the two states, $\varepsilon=1-|\langle{\psi_{\text{ans}}}|{\psi_{\xi}}\rangle|^{2}$. In Fig. 4(d) this is plotted as a function of the wave-number $k$ associated to each of the modes. Generally, the error is negligible, except for states close to the light line. In Fig. 4(e) we show that for the most subradiant state this error vanishes with chain length $N$ as $\varepsilon\sim N^{-2}$.

Moreover, this ansatz not only approaches the spatial pattern of the wave function, but its decay rate defined as $\Gamma_{\rm ans}=-(2/\hbar)\textrm{Im}\mathinner{\langle{\psi_{\text{ans}}}|}\mathcal{H}_{\textrm{eff}}\mathinner{|{\psi_{\text{ans}}}\rangle}$ captures the same scaling with the index $\xi$ and $N$: $\Gamma_{\textrm{ans}}\propto\xi^{2}/N^{3}$. The overall proportionality constant varies depending on the microscopic details (such as the polarization or the value of $d/\lambda_{0}$). For instance, for $d/\lambda_{0}=0.3$, $\Gamma_{\text{ans}}/\Gamma_{\xi=1}\approx 3/2$ (parallel polarization) and $\Gamma_{\text{ans}}/\Gamma_{\xi=1}\approx 8$ (transverse polarization). The comparison between $\Gamma_{\xi}$ and $\Gamma_{\text{ans}}$ is shown in Fig. 3(c)-(d) (solid circles correspond to the ansatz).

**2D Square Array.

**The previous results are not specific to the one-dimensional chain and can be extended to other lattice geometries. As an example, let us consider a finite square array of $N\times N$ atoms spanning the $\hat{y}$-$\hat{z}$ plane. Just like in the linear chain, we can diagonalize the block Hamiltonian with a single excitation and find the decay rates associated with the $N^{2}$ collective modes. We can also define an ansatz wave function with coefficients $c_{{\rm ans,}\textbf{k}}^{\textrm{\bf j}}\propto c_{{\rm ans,}k_{y}}^{j_{y}}c_{{\rm ans,}k_{z}}^{j_{z}}$, where $c_{{\rm ans,}k}^{j}$ are the coefficients for the one-dimensional ansatz Eq. (11). We can then associate with each of the collective modes a pair of values $\left(k_{y},k_{z}\right)$, which lies within the first Brillouin zone, and for which the corresponding ansatz produces the highest overlap with the exact state.

In Fig. 5 we have plotted the decay rates as a function of $\left(k_{y},k_{z}\right)$ after diagonalizing the Hamiltonian for an array of $50\times 50$ atoms, for different values of $d/\lambda_{0}=\left\{1,0.55,0.3\right\}$. Figures 5(a)-(c) depict the case of polarization parallel to the array, and show the emergence of subradiant states (corresponding to wave vectors beyond the light line) for $d/\lambda_{0}<1/\sqrt{2}$. As it can be seen in this figure, the most subradiant modes correspond to those at the edges of the Brillouin zone, i.e., $\left(|k_{y}|,|k_{z}|\right)\sim(\pi/d,\pi/d)$. The wave-function amplitude of this mode is a generalization of the one shown in Fig.˜4(a) for the 1D chain, where now the alternating plus and minus sign in the amplitude exhibits a checkerboard pattern. Figures 5(d)-(f) depict the decay rates for the case of transverse polarization. Here, one sees a set of subradiant states emerges beyond the light line for $d/\lambda_{0}<1/\sqrt{2}$ as before, and also a set of subradiant modes around $(k_{x},k_{y})=(0,0)$ for $d/\lambda_{0}<1$, in agreement with the infinite lattice analysis.

While we can expect that in general the decay rate of the most subradiant modes will be suppressed with the system size, the scaling is more complex than for the linear chain. Nevertheless, we have numerically verified that for collective modes at the edge of the Brilllouin zone, and if $d/\lambda_{0}$ is small enough, the decay rate will scale as $\Gamma\sim N^{-\alpha}$. In particular, we find that for $\left(k_{y},k_{z}\right)=(\pi/d,\pi/d)$ (most subradiant state), $\alpha=6$ for $d/\lambda_{0}<1/2$ and $\alpha=3$ for $1/2\leq d/\lambda_{0}<1/\sqrt{2}$, while for $\left(k_{y},k_{z}\right)=(\pi/d,0)$, $\alpha=3$ for $d/\lambda_{0}<1/2$. For ranges of $d/\lambda_{0}$ not included above, the decay rates are not suppressed with increasing system size, since in that case the wave vector $\left(k_{x},k_{y}\right)$ of the two states lies within the light line. These scalings are shown for the two states in Fig. 5(g).

**3D Cubic Array.

**While the extension from 1D chains to 2D arrays is conceptually straightforward, it appears that three-dimensional lattices are governed by different physics. In particular, in infinite 1D and 2D arrays, Bloch modes diagonalize the system, and subradiant modes can be characterized as “guided” as the associated electromagnetic fields are evanescent in the spatial directions transverse to the array. In contrast, while Bloch modes still diagonalize the system in 3D, the associated fields are necessarily extended over space in all directions. Thus, for a finite-size system, it does not appear that subradiant states can be identified based on the infinite-system results, as was possible in 1D and 2D.

Nonetheless, for completeness, we can still numerically investigate the decay rates for a single excitation in a 3D finite-size lattice. In Fig. 6, we plot the decay rates $\Gamma_{\xi}$ for the $N=10^{3}$ eigenstates associated with a $10\times 10\times 10$ lattice of two-level atoms, in the case that the polarization of the transition is aligned with one of the cubic axes. It can be seen that while decay rates are still most prominently suppressed for lattice constants $d\ll\lambda_{0}$, the effect of subradiance can survive even for lattice constants $d>\lambda_{0}$. It would be interesting to further explore the nature of subradiance in 3D systems in future work, and identify conceptual similarities it has to arrays in lower dimensions, if any.

**Atoms in a ring configuration.

**The result that we have found for 1D linear chains, indicating that subradiant modes are guided and that radiation leakage is primarily from the system ends, motivates us to study the decay rates when the atoms form a closed configuration, since this might lead to a stronger suppression in the decay. In particular, we consider now that the atoms are sitting on a circular ring separated by an equal distance $d$ (see sketch in Fig. 7) and with linear polarization transverse to the plane of the ring.

In Fig. 7, for a distance of $d/\lambda_{0}=0.3$, we numerically diagonalize the single-excitation block Hamiltonian, and plot the decay rate $\Gamma_{\xi=1}$ of the most subradiant state versus atom number $N$. It can be seen that an exponential suppression emerges, $\Gamma_{\xi=1}\sim\exp(-N)$. For the chosen parameters, the minimum decay rate for a ring drops below that of an open chain for $N\gtrsim 20$ atoms. The subradiant modes of the ring can be interpreted as “whispering gallery modes”, which weakly radiate into free space only via the finite radius of curvature. The exponential suppression with ring radius is analogous to the scaling of radiation losses in a conventional whispering gallery resonator Buck and Kimble (2003).

**Localized resonance in an atomic chain.

**Here, we also show how to achieve a spatially confined mode in a linear 1D chain of atoms, which also exhibits an exponential suppression of decay rate with atom number. This can be achieved by introducing a smooth, local variation in the lattice constant, in analogy to the principles that govern the design of a conventional photonic crystal cavity Joannopoulos et al. (2008).

To illustrate this, we consider the geometry schematically depicted in Fig. 8(a). The atomic chain along $\hat{z}$ has been divided into three regions: in the left and right regions the lattice constant $d$ is uniform and equal to $d_{\rm max}$, while in the middle it changes slowly (from $d_{\rm max}$ to $d_{\rm min}=0.75d_{\rm max}$) following the red line and creating a defect. The lattice constant in the middle is chosen to follow a sinusoidal variation, $d(z_{i})=d_{\rm max}+(d_{\rm min}-d_{\rm max})\sin^{2}(3\pi z_{i}/N)$, where $z_{i}$ ($i=N/3,...,2N/3$) denotes the atom position. In the same figure, the band structures for an infinite lattice with constant $d_{\rm max}=0.4\lambda_{0}$ and $d_{\text{min}}=0.3\lambda_{0}$ are plotted for the case of atomic polarization parallel to the chain. It can be seen that the smaller lattice constant $d_{\text{min}}$ supports propagating modes over a range of frequencies (pink shaded region) that lies within the band gap of lattice $d_{\rm max}$. Thus, for the system with slowly varying lattice constant, a set of localized resonances can appear in the middle region, with frequencies situated in the gap of lattice $d_{\rm max}$, and unable to propagate into the left and right regions.

The atomic excited state population associated with the fundamental localized mode (i.e., the mode with no nodes in the population) is illustrated in Fig. 8(b), for a representative case of $N=90$ atoms. For a smooth variation in the lattice constant (here occurring over a region of size $N/3$), one expects that the fundamental mode will have a Fourier transform with an exponentially small weight of wave vectors lying within the light line. Likewise, the leakage of this mode through the left and right regions to the ends of the chain will be exponentially suppressed, leading to an overall exponential suppression in decay rate with increasing atom number $N$. The numerically calculated decay rate of this mode is plotted in Fig. 8(c) as a function of $N$, and clearly confirms the expected behavior.

III.3 Multi-excitation modes

We now turn to the problem of multiple excitations stored in an atomic 1D lattice. As in the previous sub-section, we consider that the atoms are linearly polarized along the chain axis. However, similar results are found for the transverse polarization case. First, we note that if the Hamiltonian of Eq.(5) were composed of bosonic particles instead of spins (i.e., $\hat{\sigma}_{ge},\hat{\sigma}_{eg}\rightarrow a,a^{\dagger}$) the multiple excitation case would be trivial. In particular, the resulting Hamiltonian would be quadratic in the creation and annihilation operators, and a Fock state of $n$ excitations in a given mode would simply have a decay rate $n$ times that of a single excitation. The fact that we are dealing with spins, where a single spin cannot be excited twice ($\hat{\sigma}_{eg}^{2}=0$), leads to highly non-trivial properties of multiply excited states. In this section, we will analyze in detail the spatial properties and scaling of decay rates of multi-excitation subradiant states. While we will not explicitly utilize these states in later sections, these findings might help to provide some initial insight into how many-body physics can be encoded into subradiant manifolds.

Let us first consider the two-excitation manifold. A general state within this manifold can be written as $\mathinner{|{\psi^{(2)}}\rangle}=\sum_{i<j}c^{ij}\mathinner{|{e_{i}e_{j}}\rangle}$ where now $\mathinner{|{e_{i}e_{j}}\rangle}=\hat{\sigma}_{eg}^{i}\hat{\sigma}_{eg}^{j}\mathinner{|{g}\rangle}^{\otimes N}$ corresponds to the state whith atoms $i$ and $j$ excited while the rest remain in the ground state. Although it is necessary only to specify $c^{ij}$ for $i<j$ to define the wave function, in the following plots and for visual appeal we also assign values to $c^{ij}$ for $j<i$, by simply defining $c^{ij}=c^{ji}$. To illustrate that the spin system behaves differently than a bosonic system, we begin by considering the two-excitation state formed by occupying the same single-excitation mode twice. In particular, we construct the two-excitation state given by $\mathinner{|{\psi^{(2)}_{\rm b}}\rangle}\propto(S_{\xi=1}^{\dagger})^{2}\mathinner{|{g}\rangle}^{\otimes N}$ (properly normalized). Here, $S_{\xi=1}^{\dagger}=N^{-1/2}\sum_{j}c^{j}_{\xi=1}\hat{\sigma}_{eg}^{j}$ is the collective operator that creates the most subradiant single excitation in a chain of $N$ atoms, when applied to the ground state. In Fig. 9(a) we plot the corresponding probability density $|c_{\rm b}^{ij}|^{2}$ for the case $N=50$. The two-excitation wave function appears relatively smooth, except for a sharp cut $c_{\rm b}^{ii}=0$ along the diagonal, owing to the fact that a single spin cannot be excited twice. As a result of this feature, the Fourier transform of this state will be relatively broad, and in particular, will contain many components that lie within the light line and can subsequently radiate. Its decay rate, defined as $\Gamma^{(2)}_{\text{b}}=-(2/\hbar)\textrm{Im}\mathinner{\langle{\psi^{(2)}_{\textrm{b}}}|}\mathcal{H}_{\textrm{eff}}\mathinner{|{\psi^{(2)}_{\textrm{b}}}\rangle}$ is only suppressed with the length of the chain as $\Gamma^{(2)}_{\text{b}}\sim N^{-1}$, in stark contrast to the single-excitation case. This scaling is shown in Fig. 9(d) (orange circles).

Numerically, we now exactly diagonalize the Hamiltonian in the two-excitation manifold, and identify the most subradiant state. The scaling of its decay rate with $N$ is plotted in Fig. 9(d) (blue open circles), and is seen to preserve the scaling $\Gamma\sim N^{-3}$ present in the single-excitation manifold. The probability density $|c^{ij}_{\xi=1}|^{2}$ for the case of $N=50$ is plotted in Fig. 9(b). The wave function appears distinctly different than that of Fig. 9(a), and in particular, it appears that the two excitations are smoothly repelled from one another.

From the previous considerations, it is apparent that the most subradiant states should simultaneously satisfy that they are composed predominantly of wave vectors beyond the light line, and that the excitations are smoothly repelled from one another in order to avoid sharp kinks in the wave function. This inspires us to try an antisymmetric (or fermionized) ansatz $\mathinner{|{\psi_{\rm ans}^{(2)}}\rangle}$ for the wave function of the form $c_{{\rm ans,}k_{1}k_{2}}^{ij}\propto c_{\text{ans},k_{1}}^{i}c_{\text{ans},k_{2}}^{j}-c_{\text{ans},k_{2}}^{i}c_{\text{ans},k_{1}}^{j}$ (properly normalized). Here $c_{\text{ans,}k_{n}}^{i}$ denote the coefficients of the single-excitation orthonormal ansatz Eq.(11) associated with the wave vector $k_{n}$. Such an ansatz naturally constructs a state that incorporates “Pauli exclusion”, and a smooth separation of excitations in space. Taking $k_{1}d=\pi N/(N+1)$ and $k_{2}d=\pi(N-1)/(N+1)$, i.e., building a two-excitation state from the two most subradiant single-excitation states, yields a wave function $c_{\rm{ans}}^{ij}$ that agrees well with the exact one, as seen in Fig. 9(c) where the probability density is plotted. Moreover, the decay rate associated with this state scales again with the particle number as $\Gamma^{(2)}_{\textrm{ans}}\sim N^{-3}$, as it is shown in Fig. 9(d) (blue solid circles).

We can then associate with each of the exact two-excitation collective states a pair of quasi-momentum values $\left\{k_{1},k_{2}\right\}$ for which the wave-function overlap with the ansatz is maximum. The exact decay rates as a function of these values are plotted in Fig. 9(e). This figure shows that when both $k_{1},k_{2}>k_{0}$ the decay rates are strongly suppressed, and we can identify this region as the one containing the subradiant states. In fact, the sum of decay rates of the single-excitation modes used to construct the ansatz, i.e., $\Gamma^{(2)}_{\textrm{sum}}=\Gamma_{\textrm{ans},k_{1}}+\Gamma_{\textrm{ans},k_{2}}$, is not far from the exact value. This is quantified by the relative error $\delta\Gamma\equiv|\Gamma^{(2)}-\Gamma^{(2)}_{\textrm{sum}}|/\Gamma^{(2)}$ and it is plotted in Fig. 9(f). For completeness, we also show in Fig. 9(g) the error in overlap between each two-excitation eigenstate and the best-matched ansatz state, $\varepsilon=1-|\langle{\psi}|{\psi^{(2)}_{\text{ans}}}\rangle|^{2}$. This error is very small in the subradiant region.

In the more general case of $n$ excitations, the most subradiant mode and its decay rate $\Gamma^{(n)}_{\xi=1}$ can be found as in the previous cases by exactly diagonalizing the corresponding block Hamiltonian. For a low density of excitations, the scaling of the decay rate with the chain length is still as in the single and two-excitation manifolds, i.e., $\Gamma^{(n)}_{\xi=1}\sim N^{-3}$, as shown in Fig. 10(a) (open blue symbols). For comparison, we also show in the same figure (orange symbols) the decay rate of the state with $n$ excitations in the most subradiant mode, i.e., $\mathinner{|{\psi^{(n)}_{\rm b}}\rangle}\propto(S_{\xi=1}^{\dagger})^{n}\mathinner{|{g}\rangle}^{\otimes N}$, which scales as in the two-excitation case, $\Gamma_{\rm b}^{(n)}\sim N^{-1}$.

One can also numerically evaluate the error $\varepsilon$ in overlap between the most subradiant state and an ansatz state, $\mathinner{|{\psi_{\rm ans}^{(n)}}\rangle}=\sum_{i_{1}<i_{2}<...<i_{N}}c_{\textrm{ans},k_{1}k_{2}...k_{N}}^{i_{1}i_{2}...i_{N}}$. Here the wave-function amplitudes $c_{\textrm{ans},k_{1}k_{2}...k_{N}}^{i_{1}i_{2}...i_{N}}$ are generalized from the two-excitation case, and constructed from the Slater determinant of $n$ single-excitation wave-function ansatz coefficients $c_{\textrm{ans},k_{n}}^{i_{n}}$. For the most subradiant mode these correspond to the $n$ most subradiant single-excitation modes and in the large atom number limit the error is found to scale like $\varepsilon\sim N^{-2}$.

If the ansatz holds, then for $n$ excitations one expects that the decay rate for the most subradiant state scales like $\Gamma^{(n)}_{\xi=1}/\Gamma_{0}\sim\sum_{m=1}^{n}m^{2}/N^{3}\sim(n/N)^{3}$. In Fig. 10(b), we compare this predicted scaling with numerically calculated values of $\Gamma_{\xi=1}^{(n)}$, and find qualitatively good agreement for low excitation density $n/N\ll 1$. We note that the prediction of the ansatz also seems physically reasonable in that it can be extended to the thermodynamic limit, as it predicts a decay rate that only depends on the density $n/N$ of excitations.

IV Atoms coupled to a nanofiber: selectively-radiant states

In the previous section we elucidated the nature of subradiant states in atomic arrays, whose long-lived nature arises from weak coupling to all propagating electromagnetic modes. Subradiant manifolds themselves might be useful for many purposes, for example, to accumulate interactions without dissipation in order to realize strongly correlated states. However, to realize an efficient atom-light interface, one would instead prefer to utilize a set of atomic states that strongly radiate into a desired electromagnetic mode (or set of modes) through constructive interference, while destructive interference simultaneously suppresses the emission rate into all undesired modes. We term states that satisfy this property to be “selectively radiant,” as the overall emission rate might not be small, but the branching ratio into desired versus undesired channels could be extremely high. It should be noted that such a definition of “selectively radiant” is somewhat arbitrary – for example, even a single isolated atom emitting into a dipole radiation pattern is selectively radiant, if the preferred optical mode is defined to be the dipole pattern itself. In practice, however, the collection efficiency of a dipole pattern with realistic optics is quite small van Enk and Kimble (2000, 2001); Darquié et al. (2005); Tey et al. (2008); Wrigge et al. (2008); Hetet et al. (2011), and a functionally useful definition should involve a mode (e.g., focused Gaussian beam or guided mode of a dielectric structure) that is generally accepted to be efficient to match to.

Here, we show that one natural way to realize and utilize selectively radiant states is to couple one-dimensional atomic chains to the guided modes of an optical nano-structure (such as an nanofiber). Qualitatively, for sufficiently small lattice constants $d<\lambda_{0}/2$, a set of spin-wave excitations with associated wave vector $|k_{z}|>k_{0}$ emerge, which inefficiently radiate into free space as the wave vector lies beyond the light line. However, as an optical mode guided by a high-index dielectric itself has a wave vector $|k_{z}|>k_{0}$, we show that it is possible that a set of spin-wave excitations simultaneously experiences an enhanced emission rate into the guided modes while being subradiant to free space. We will provide an explicit construction of a protocol where selectively radiant states are exploited, involving a quantum memory or photon storage. We find in particular that these states enable an exponential improvement in the error probability versus atom number, over previously known bounds.

This section is organized as follows. Section IV A describes the nanofiber and provides the Hamiltonian that governs the interactions between the atoms located in the vicinity of the nanostructure. We introduce the “collective emission” model, which accounts for atom-atom interactions both through the guided and non-guided modes of the fiber. We also present the “independent emission” model, in which atoms interact through the guided modes but coupling via free-space modes is neglected. This thus represents the “standard model” of atom-light interactions specifically applied to nanofibers. In particular, it reproduces previously accepted bounds for fidelities of photon storage, against which the “collective emission” model can be compared. Section IV B describes linear optical processes (i.e., single-photon transmission and reflection) for two-level atoms coupled to the fiber. We show how the conventional figure of merit, the optical depth, is not sufficient to characterize optical transport through the array when the collective emission into non-guided modes is taken into account. In Section IV C we study how selective radiance influences electromagnetically induced transparency Harris et al. (1990); Boller et al. (1991); Fleischhauer and Lukin (2000); Lukin (2003); Fleischhauer et al. (2005), a phenomenon that is commonly used in photon storage protocols. In particular, we show that the bandwidth-delay product, which quantifies the number of photons that can be stored in the atomic medium, scales linearly with the number of atoms. The linear scaling is characteristic of ideal, non-lossy systems, and is in contrast to the independent emission model, which predicts a scaling that goes with the square root of the optical depth. In this sub-section, we also provide the first glimpse of improvement in photon storage beyond traditional bounds. Finally, in Section IV D we demonstrate how to achieve an exponential suppression with the atom number on the infidelity of a quantum memory.

IV.1 Description of the nanofiber

The possibility of enhancing atom-light interactions through selective radiance should exist for any nanophotonic structure where atoms can be periodically trapped, including in nanofibers Vetsch et al. (2010); Goban et al. (2012); Sørensen et al. (2016); Corzo et al. (2016) and 1D and 2D photonic crystal waveguides Douglas et al. (2015); González Tudela et al. (2015); Goban et al. (2015); Hood et al. (2016). For complicated structures, however, the Green’s function cannot be obtained analytically. Furthermore, while the Green’s function can be calculated numerically Sullivan (2000), to do so with sufficient accuracy appears quite challenging (in particular, it must be calculated with enough accuracy so that diagonalization correctly captures subradiant emission rates that scale like large inverse powers of $N$). Motivated by this observation, here we focus on a special geometry where the Green’s function can be exactly obtained, which consists of a chain of atoms coupled to guided modes of an infinite, cylindrical nanofiber.

We consider that the chain of atoms lies parallel to the axis of a dielectric nanofiber oriented along $\hat{z}$, with radius $r$ and relative permittivity $\epsilon$ (or corresponding refractive index $n_{\text{fiber}}=\sqrt{\epsilon}$). As shown in Fig. 11, the distance between the atoms and the center of the nanofiber is $\rho_{\text{a}}$, and the orientation of their dipole transition is directed along $\hat{\rho}$, perpendicular to the axis of the nanofiber (i.e., $\boldsymbol{\wp}=\wp\hat{\rho}$). The Green’s function for such a nanofiber can be found analytically Stratton (1941); Yeh and Shimabukuro (2008); Le Kien et al. (2005); Dzsotjan et al. (2010); Le Kien and Rauschenbeutel (2014); Qi et al. (2016); Le Kien and Rauschenbeutel (2017). In particular, we follow the work of Klimov and Ducloy  Klimov and Ducloy (2004). In the following, we provide a qualitative description of the derivation, while details are given in Appendix C. The first step in the derivation is to separate the Green’s function into two terms, i.e., $\mathbf{G}({\bf r}_{i},{\bf r}_{j},\omega_{0})=\mathbf{G}_{0}({\bf r}_{i},{\bf r}_{j},\omega_{0})+\mathbf{G}_{\text{sc}}({\bf r}_{i},{\bf r}_{j},\omega_{0})$. Here $\mathbf{G}_{0}$ is the already-known vacuum Green’s function given by Eq. (6), which corresponds to the field emitted by a dipole in free space, and $\mathbf{G}_{\text{sc}}$ is a general solution to the sourceless wave equation, which will physically correspond to the (thus far unknown) field scattered by the nanofiber. Exploiting the cylindrical symmetry of the problem, one can employ separation of variables and expand the vacuum and scattered Green’s functions using a set of functions $f_{m,k_{\parallel}}(\rho)e^{ik_{\parallel}z+im\phi}$. Here $k_{\parallel}$ is the longitudinal wave vector and $m$ denotes angular momentum. The coefficients in the expansion of $\mathbf{G}_{\text{sc}}$ associated with each value of $k_{\parallel}$ and $m$ inside and outside the fiber are a priori unknown, but can be solved for through equations that enforce electromagnetic field continuity relations at the surface of the fiber.

The fiber supports a set of guided modes, i.e., electromagnetic modes that propagate along the nanostructure and are confined in the transversal direction. These modes are denoted by their angular momentum $m$, and their associated wave vectors $k_{m}(\omega_{0})$ always satisfy $|k_{m}(\omega_{0})|>\omega_{0}/c$, as the guiding mechanism is by total internal reflection [here we have dropped the “$\parallel$” subscript associated with the guided mode wave vector $k_{m}(\omega_{0})$, for notational simplicity]. In other words, these modes are evanescent, and their dispersion relations are situated beyond the light line. The number of guided modes is determined by the fiber radius and dielectric constant. We will restrict ourselves to a single-mode fiber (with $m=\pm 1$), which occurs for a sufficiently small fiber radius. Instead of working with $\mathbf{G}_{0}$ and $\mathbf{G}_{\text{sc}}$, for our purposes it is convenient to isolate the guided mode contribution and separate the Green’s function $\mathbf{G}({\bf r}_{i},{\bf r}_{j},\omega_{0})=\mathbf{G}_{\text{1D}}({\bf r}_{i},{\bf r}_{j},\omega_{0})+\mathbf{G}^{\prime}({\bf r}_{i},{\bf r}_{j},\omega_{0})$ into two terms: one that characterizes the excitation of the guided mode of the fiber, $\mathbf{G}_{\text{1D}}({\bf r}_{i},{\bf r}_{j},\omega_{0})$, and another that describes the non-guided electromagnetic modes, $\mathbf{G}^{\prime}({\bf r}_{i},{\bf r}_{j},\omega_{0})$. In particular, the guided Green’s function takes the form $\mathbf{G}_{\text{1D}}({\bf r}_{i},{\bf r}_{j},\omega_{0})=\mathbf{g}(\boldsymbol{\rho}_{\text{a}})\,e^{{\rm i}k_{\text{1D}}|z_{i}-z_{j}|}$, where $\mathbf{g}(\boldsymbol{\rho}_{\text{a}})$ is a tensor that only depends on the radial and azimuthal position of the atoms (assumed to be identical), and $k_{\rm 1D}=|k_{\pm 1}(\omega_{0})|$.

The dynamics of the atoms is governed by the non-Hermitian Hamiltonian $\mathcal{H}_{\rm eff}$ of Eq. (5), which can be similarly split, i.e., $\mathcal{H}_{\text{eff}}=\mathcal{H}_{\rm 1D}+\mathcal{H}^{\prime}$. From the form of $\mathbf{G}_{\text{1D}}$ given above, it follows that

[TABLE]

where $\Gamma_{\rm 1D}=(2\mu_{0}\,\omega_{0}^{2}\wp^{2}/\hbar)\,\text{Im}\,G^{\rm 1D}_{\rho\rho}({\bf r}_{i},{\bf r}_{i},\omega_{0})$ is the spontaneous emission rate of a single atom into the fiber guided mode. The plane-wave dependence reflects the fact that the guided photon propagates without diffraction between two atoms and thus produces an infinite-range interaction.

The non-guided term

[TABLE]

accounts for the interaction through the remaining non-guided electromagnetic modes. Already for just a single atom, the self term of the non-guided Green’s function $G^{\prime}_{\rho\rho}({\bf r}_{i},{\bf r}_{i},\omega_{0})$ gives rise to both a frequency shift and a decay rate that we will denote as $J^{\prime}-{\rm i}\Gamma^{\prime}/2=-(\mu_{0}\,\omega_{0}^{2}\wp^{2}/\hbar)\,G^{\prime}_{\rho\rho}({\bf r}_{i},{\bf r}_{i},\omega_{0})$. This self-term reflects the fact that the modification of electromagnetic modes by the nanofiber causes a single atom to have a resonance frequency $\omega_{0}+J^{\prime}$ shifted from its vacuum value, and a decay rate into radiative modes $\Gamma^{\prime}$ different than $\Gamma_{0}$. For many atoms, the above Hamiltonian accounts for collective emission into non-guided modes, as it takes into account atom-atom interactions that are not mediated by the guided mode. Unlike $\mathbf{G}_{\text{1D}}$, $\mathbf{G}^{\prime}$ does not admit a simple form, and in what follows it will be evaluated numerically using the prescription detailed in Appendix C. Throughout this manuscript, we will refer to the dynamics generated by $\mathcal{H}_{\text{1D}}+\mathcal{H}^{\prime}$ as the “collective emission” model.

Whether in free space or a nanofiber (or other guided structures), exact collective effects involving modes that are not directly of interest (such as those captured by $\mathcal{H}^{\prime}$) are typically difficult to treat in the context of applications of atomic ensembles. It is usually heuristically argued that photon-mediated interactions through these modes are not relevant, particularly for disordered or dilute atomic gases, and the “standard” model within quantum optics is to ignore such terms Hammerer et al. (2010); Gorshkov et al. (2007a); Chang et al. (2012); Qi et al. (2016). Specifically, the terms of $\mathbf{G}^{\prime}({\bf r}_{i},{\bf r}_{j},\omega_{0})$ involving two different atoms (${\bf r}_{i}\neq{\bf r}_{j}$) are assumed to be zero, and the Hamiltonian accounting for emission into non-guided modes reduces to

[TABLE]

In this approximation, the non-guided modes of the fiber introduce a modified Lamb shift due to the presence of the fiber surface, and more importantly, provide independent baths for each atom to emit into, at a rate $\Gamma^{\prime}$. We will refer to the dynamics generated by $\mathcal{H}_{\text{1D}}+\mathcal{H}^{\prime}_{\text{indep}}$ as the “independent emission” model.

Before proceeding further, we digress to clarify the different usages of the terms super/sub-radiance in literature. Within the independent emission model, the concept of superradiance and subradiance has also been discussed, since $\mathcal{H}_{\text{1D}}$ alone yields a set of collective atomic states that radiate strongly or weakly into the waveguide Chang et al. (2012); van Loo et al. (2013); Asenjo-Garcia et al. (2017). Similar effects have also been pointed out in cavities (with collective states emitting strongly or weakly into the cavity mode) Bohnet et al. (2012); Reimann et al. (2015); Casabone et al. (2015). Protocols for photon generation and storage and other quantum information tasks have been built around the manipulation of these states Chang et al. (2012); González-Tudela et al. (2015). However, as these models still assume independent emission into free space, these protocols do not surpass conventional error bounds.

Throughout this manuscript, the nanofiber radius is taken to be $k_{0}r=1.2$, the distance between the atoms and the center of the nanofiber is $\rho_{\text{a}}=1.5r$, and the dielectric constant of the fiber is $\epsilon=4$ (as that of silicon nitride). As an illustration, for the D1 line of Cesium (of resonance frequency $\omega_{0}=2\pi\times 335.1$ THz), the radius of the fiber would be $r\simeq 170$ nm, and the distance between the atoms and the fiber surface would be approximately $85$ nm. The wave vector of the photonic guided mode is found to be $k_{\rm 1D}\simeq 1.3k_{0}$, larger than any wave vector within the light line, as the guided mode is confined. The single-atom decay rates are calculated to be $\Gamma_{\rm 1D}\simeq 0.4\Gamma_{0}$ to the guided mode, and $\Gamma^{\prime}\simeq 1.3\Gamma_{0}$ to the non-guided modes. The modified Lamb shift due to the fiber is $J^{\prime}\simeq-0.5\Gamma_{0}.$

In the following sub-sections, we will utilize the formalism introduced above to identify novel phenomena that emerge when collective emission is exactly accounted for, which cannot be predicted from the independent emission approximation, and will show how collective emission enables exponential improvement in performance for quantum memories of light.

IV.2 Linear optics for two-level atoms

We begin by studying the transmission and reflection properties of a chain of atoms coupled to the fiber within the “independent emission” model. The effective Hamiltonian that describes the atomic dynamics under a coherent-state guided-mode probe field of frequency $\omega_{\text{p}}$ reads $\mathcal{H}_{\text{tot}}=\mathcal{H}_{\text{drive}}+\mathcal{H}^{\prime}_{\text{indep}}+\mathcal{H}_{\text{1D}}$. The Hamiltonians $\mathcal{H}_{\text{1D}}$ and $\mathcal{H}^{\prime}_{\text{indep}}$ are defined in Eqs. (12), and (14), respectively, and the driving term is given by

[TABLE]

In the above expression, $\Delta=\omega_{\rm p}-\omega_{0}$ is the detuning between the probe field frequency and the atomic resonance frequency. We have also defined the Rabi frequency of the guided-mode probe field as $\Omega=\boldsymbol{\wp}^{*}\cdot{\bf E}_{\rm p}^{+}(\rho_{\text{a}})/\hbar$, where ${\bf E}_{\rm p}^{+}({\bf r})\equiv\mathinner{\langle{\hat{{\bf E}}_{\rm p}^{+}(\rho_{\text{a}})}\rangle}$ is the amplitude of a coherent probe field that implicitly contains the radial position $\rho_{\text{a}}$ of the atoms. For the remainder of this subsection, we consider that the probe field is weak and does not saturate the atoms. Therefore, all the calculations can be performed in the single-excitation manifold, the realm of classical linear optics.

In the single excitation subspace, the wave function of the atomic ensemble is written as the superposition $\mathinner{|{\psi(t)}\rangle}=c_{g}(t)\mathinner{|{g}\rangle}^{\otimes N}+\sum_{j=1}^{N}c_{e}^{\,j}(t)\mathinner{|{e_{j}}\rangle}$, where $\mathinner{|{e_{j}}\rangle}=\hat{\sigma}_{eg}^{j}\mathinner{|{g}\rangle}^{\otimes N}$. In the low saturation regime, with $c_{g}\simeq 1$, the evolution equations for the amplitude of the $\mathinner{|{e}\rangle}$ states are found to be

[TABLE]

The generalized input-output expression of Eq. (2) allows us to calculate the guided-mode field at any point of the fiber, which reads

[TABLE]

It is important to notice that the Green’s function appearing in the field equation is not the total one, but just that of the guided mode, as it describes the propagation of the photonic guided field along the nanostructure. Due to the cylindrical symmetry of the fiber, the guided modes with angular momenta $m=\pm 1$ are degenerate. One can alternatively take superpositions of these to obtain quasilinearly-polarized $H$ and $V$ modes Le Kien and Rauschenbeutel (2014); Qi et al. (2016). The polarization basis of the fiber can always be set so that the $H$ mode at the atomic positions has polarization components along the $\hat{\rho}$ and $\hat{z}$ directions. We will consider the case where the probe field is $H$-polarized, in which case the atoms scatter solely back into $H$, and the $V$-polarized mode de-couples from the problem.

We can thus project the input-output equation into 1D equations for the $H$-modes, and further separate the guided fields into left- and right-propagating components. The resulting equations are given by

[TABLE]

where $\hat{E}_{\text{in,R(L)}}^{+}(z)$ are the right(left)-going vacuum fluctuation fields, and $\Theta$ is the Heaviside function. The vacuum fluctuations do not contribute to any of our observables of interest. For convenience, we have re-scaled the fields so that the atomic parameters $\Omega$ and $\Gamma_{\rm 1D}$ directly appear.

In the quasistatic limit ($\dot{c}_{e}^{\,j}=0$), the solutions of Eq. (16) for the $\mathinner{|{e}\rangle}$-state amplitudes are directly proportional to the probe field Rabi frequency $\Omega$. Together with Eqs. (18a,18b), this allows us to find the reflection and transmission coefficients for the guided field. For example, the transmittance is found by evaluating $T=\mathinner{\langle{\psi|\hat{E}_{\text{1D,R}}^{-}(z)\hat{E}_{\text{1D,R}}^{+}(z)|\psi}\rangle}/\Omega^{2}$, where $z$ is a position immediately after the right-end of the atomic chain, and $\hat{E}_{\text{1D,R}}^{-}(z)$ is the Hermitian conjugate of $\hat{E}_{\text{1D,R}}^{+}(z)$. A similar expression can be found for the reflectance $R$. One can also calculate the loss probability due to scattering into free space, which is given by $\kappa=1-T-R$.

We choose the distance between the atoms to be $d=\lambda_{\text{1D}}/4$, with $\lambda_{\text{1D}}=2\pi/k_{\text{1D}}$ being the guided-mode wavelength. Any other separation except for the so-called mirror configuration, i.e., $d=\lambda_{\text{1D}}/2$ or integer multiples thereof Chang et al. (2012), would display qualitatively similar optical properties. In Appendix D, we analyze the linear optics of such a special configuration, which has been theoretically known and experimentally observed to become a very reflective mirror Deutsch et al. (1995); Schilke et al. (2011); Chang et al. (2012), around which powerful protocols for quantum information processing can be built Chang et al. (2012); González-Tudela et al. (2015); Paulisch et al. (2016). In the mirror configuration and within the independent emission model, there is only one atomic collective state that couples to the guided mode of the fiber, decaying superradiantly into it at a rate $N\Gamma_{\rm 1D}$.

In contrast, for any other separation, every atomic collective state is excited by the probe field, and contributes to light transmission and reflection. Therefore, the behavior of the atoms cannot be attributed to a single “super-atom” of enhanced decay rate, and the transmission spectrum – depicted by the red line of Fig. 12(a) – differs significantly from a Lorentzian Asenjo-Garcia et al. (2017). For large enough number of atoms, and for low single-atom coupling efficiency into the waveguide ($\Gamma_{\rm 1D}\lesssim\Gamma^{\prime}$), the transmittance approximately follows the expression

[TABLE]

in accordance with the result obtained for a free-space atomic gas Asenjo-Garcia et al. (2017). On resonance (when $\Delta-J^{\prime}=0$), the figure of merit that determines how much light is transmitted is the optical depth, $D=2N\Gamma_{\rm 1D}/\Gamma^{\prime}$. For a chain of $N=20$ atoms, the expression of Eq. (19) nicely reproduces the transmittance spectrum shown in Fig. 12(a). The corresponding reflectance spectrum is displayed by the red curve of Fig. 12(b), which shows a very small bump, as the distance $d=\lambda_{\text{1D}}/4$ minimizes reflection due to destructive interference Caneva et al. (2015). As both transmission and reflection are very small close to resonance, the dominant process is photon loss due to atom-mediated scattering into free space. The loss probability $\kappa$ is shown in Fig. 12(c).

If the atoms are closely packed, the above calculations are no longer valid due to the atomic interactions mediated by non-guided modes. Nevertheless, the previous techniques can be straightforwardly modified to calculate the new transmission and reflection coefficients. Within the “collective emission” model, the atoms evolve under the Hamiltonian $\mathcal{H}_{\text{tot}}=\mathcal{H}_{\text{drive}}+\mathcal{H}^{\prime}+\mathcal{H}_{\text{1D}}$, where $\mathcal{H}^{\prime}$ replaces $\mathcal{H}^{\prime}_{\text{indep}}$. In the low saturation limit, the evolution equations for the $\mathinner{|{e}\rangle}$ state amplitudes now read

[TABLE]

Once again, we evaluate Eqs. (18a,18b) using the steady state solution for the atomic wave function, in order to reconstruct the electromagnetic field along the nanofiber. We overlay our results for the transmission, reflection, and loss probability spectra in Figures 12(a-c). The transmittance spectrum displays many sharp peaks that are not observed within the independent emission model, similar to what was found in Ref. Kornovan et al. (2016). These peaks correspond to the interference between different collective atomic modes, whose response can be observed due to the diminished photon loss. Close to resonance, reflectance is significantly larger than that obtained within the independent emission model. As a matter of fact, accounting for cooperative emission into non-guided modes lowers significantly the probability $\kappa$ of photon scattering into free space, as can be observed in Fig. 12(c). However, this decrease in loss is not uniform for all detunings, and close to resonance this spectrum also showcases sharp peaks. Globally, the behavior is much more complex than that of a “standard” atomic ensemble, such as given by Eq. (19).

One interesting question is how the loss scales with the number of atoms. We find that, far from resonance, $\kappa/\kappa_{\text{indep}}\sim N^{-1}$, where $\kappa_{\text{indep}}$ is the photon loss probability when the collective emission into non-guided modes is neglected. However, in Section III we have found that for a chain of atoms in free space the decay rate of the most subradiant mode scales as $N^{-3}$. Taken together, Figs. 12(a-c) suggest a simple reason for this apparent discrepancy. Based on previous arguments of Sec. IIIA, both the infinite atomic chain and the fiber have sets of perfectly guided modes, which experience zero radiation into free space. However, the dispersion relation of the effective medium formed by the fiber and the atomic chain is different from that of the bare fiber alone (i.e., for a given guided-mode frequency, there is a different wave vector). This impedance mismatch leads to large scattering loss at the interface between the two different systems (bare fiber versus fiber with atoms), in close analogy to what occurs between different conventional waveguides Moerman et al. (1997); Almeida et al. (2003).

IV.3 Electromagnetically induced transparency

Having posited that scattering at the interface between the bare fiber and the atomic chain dominates the losses observed in the two-level case, we now attempt to reduce these losses by better matching the dispersion relations of the two regions, using three-level atoms under conditions of electromagnetically induced transparency (EIT) Harris et al. (1990); Fleischhauer and Lukin (2000); Lukin (2003); Fleischhauer et al. (2005).

The system under consideration is illustrated in Fig. 13. In addition to the $\mathinner{|{g}\rangle}$ to $\mathinner{|{e}\rangle}$ transition studied earlier, a third metastable level $\mathinner{|{s}\rangle}$ (of frequency $\omega_{s}$) is added. We assume that the $\mathinner{|{e}\rangle}$ to $\mathinner{|{s}\rangle}$ transition does not couple to the optical fiber (e.g., due to its dipole matrix element being orthogonal to the guided mode polarization), but can be addressed by an external classical control field of Rabi frequency $\Omega_{c}$ that propagates through free space. Through a two-photon interference effect mediated by the control field, a near-resonant guided photon interacting with an atom originally in state $\mathinner{|{g}\rangle}$ can be coherently mapped to state $\mathinner{|{s}\rangle}$, with minimal excitation of state $\mathinner{|{e}\rangle}$. The lack of population in $\mathinner{|{e}\rangle}$ and associated photon scattering causes the otherwise optically opaque medium to become transparent, and thus EIT nominally preserves the effective refractive index of the guided mode.

Again, we consider the case where atoms are separated by a distance $d=\lambda_{\text{1D}}/4$ ($k_{\rm 1D}d=\pi/2$), to guarantee minimal reflection. Any other distance of the form $d=n\lambda_{\text{1D}}/4$ (with $n$ being odd) also strongly suppresses reflection and should suffice, as long as $d$ fulfils the subradiance condition. In particular, as atoms nominally do not alter the effective index under EIT, the guided wave vector $k_{\rm 1D}$ itself should lie outside the light line. Without atoms this is clearly always true, as the fiber mode is guided. With atoms, however, one must ensure that $k_{\rm 1D}$ lies outside the light line when folded back into the first Brillouin zone. If we set the distance between the atoms to be such that $k_{\rm 1D}d=\pi/2$, then $k_{\rm 1D}$ lies within the first Brillouin zone and automatically satisfies this constraint, $k_{0}<k_{\rm 1D}\leq\pi/d$. However, if $k_{\rm 1D}d=3\pi/2$, the condition on the guided mode wave vector, $k_{\rm 1D}>3k_{0}$, becomes much more stringent. In fact, for the radius and dielectric constant of the fiber here considered, the subradiance condition is not met for $d=3\lambda_{\text{1D}}/4$.

We begin by solving for the characteristics of EIT under the independent emission model, which we find to reproduce previously derived and well-known results in free-space atomic ensembles. In particular, we consider the system evolving under the effective Hamiltonian $\mathcal{H}_{\text{tot}}=\mathcal{H}_{\text{1D}}+\mathcal{H}^{\prime}_{\text{indep}}+\mathcal{H}_{\text{drive}}+\mathcal{H}_{\text{c}}$. The Hamiltonians $\mathcal{H}_{\text{1D}}$, $\mathcal{H}^{\prime}_{\text{indep}}$, and $\mathcal{H}_{\text{drive}}$ are defined in Eqs. (12), (14), and (15), respectively. $\mathcal{H}_{\text{c}}$ captures the interaction of the atoms with the control field, and is given by

[TABLE]

where $\Delta_{\rm s}=\omega_{\rm p}+\omega_{c}-\omega_{\rm s}$ is the two-photon detuning. We take the control field Rabi frequency $\Omega_{c}$ to be real, and allow for a possible spatial dependence. We have also assumed that $\mathinner{|{e}\rangle}$ has a negligible decay rate into $\mathinner{|{s}\rangle}$, as in the case of a dipole-forbidden transition or ladder system. For EIT within the independent emission model, this assumption is not necessary, as the emission rate from $\mathinner{|{e}\rangle}$ to $\mathinner{|{s}\rangle}$ can be incorporated into $\Gamma^{\prime}$ and simply leads to a moderate decrease of optical depth $D$. Such a condition, however, becomes important when considering the collective emission case (see more detailed discussion about multi-level structure in Sec. V).

Within the single excitation manifold, the wave function of the atomic ensemble is $\mathinner{|{\psi(t)}\rangle}=c_{g}(t)\mathinner{|{g}\rangle}^{\otimes N}+\sum_{j=1}^{N}c_{e}^{\,j}(t)\mathinner{|{e_{j}}\rangle}+\sum_{j=1}^{N}c_{s}^{\,j}(t)\mathinner{|{s_{j}}\rangle}$, with $\mathinner{|{s_{j}}\rangle}=\hat{\sigma}_{sg}^{j}\mathinner{|{g}\rangle}^{\otimes N}$. For a uniform control field [$\Omega_{\text{c}}(z_{i})\equiv\Omega_{\text{c}}$] and in the low saturation limit [$c_{g}(t)\simeq 1$], the equations for the evolution of the amplitudes of the $\mathinner{|{e}\rangle}$ and $\mathinner{|{s}\rangle}$ states read

[TABLE]

We solve these equations in the steady state and readily find $c_{s}^{\,i}=-(\Omega_{c}/\Delta_{s})\,c_{e}^{\,i}$, and

[TABLE]

Here, we have chosen $\Delta_{\rm s}=\Delta-J^{\prime}$, which assures total transparency when the probe field is resonant with the (shifted) $\mathinner{|{e}\rangle}-\mathinner{|{g}\rangle}$ transition ($\Delta-J^{\prime}=0$). Having found the steady state solution of the spin wave function, it is now possible to calculate the transmitted guided-mode field by means of the input-output expression of Eq. (17), and, therefore, the transmission coefficient of the array under EIT, $t_{\rm EIT}$.

The transmission coefficient gives us enough information to calculate two key quantities that describe the EIT medium: the group velocity of the polariton, and the bandwidth-delay product, a parameter that quantifies how many spatially separate photons can be stored in the atomic ensemble Chang et al. (2011). After propagating along the atomic chain, the guided mode field acquires a phase $t_{\rm EIT}\equiv e^{{\rm i}k_{\rm eff}Nd}$, where $k_{\rm eff}$ is a complex effective wave vector that encodes both light absorption and dispersion. Up to second order in the atom-probe detuning, the effective wave vector reads Chang et al. (2011); Asenjo-Garcia et al. (2017)

[TABLE]

where $\eta=0(1)$ for an even (odd) number of atoms. It can be seen that when $\Delta=J^{\prime}$, the effective wave vector perfectly matches that of the bare fiber, $k_{\text{eff}}=k_{\text{1D}}$, as originally desired. From the above expression, the group velocity at the center of the transparency window is found to be $v_{\text{g}}=\left(\partial k_{\rm eff}/\partial\Delta\right)^{-1}=2\Omega_{c}^{2}d/\Gamma_{\rm 1D}$. The delay time, i.e., the time it takes for this slow polariton to traverse the medium is $\tau=Nd/v_{\text{g}}=N\Gamma_{\rm 1D}/2\Omega_{c}^{2}$. The bandwidth of the transparency window, which dictates how spectrally narrow a photon has to be to propagate with high transparency, is defined as

[TABLE]

where $\delta$ is the detuning for which $|t_{\text{EIT}}|^{2}=1/e$. Therefore, the bandwidth-delay product, $\mathcal{P}=\tau\Delta_{\text{EIT}}=\sqrt{2N\Gamma_{\rm 1D}/(\Gamma^{\prime}+\eta\Gamma_{\rm 1D}/2N)}\simeq\sqrt{D}$, scales with the square root of the optical depth $D=2N\Gamma_{\rm 1D}/\Gamma^{\prime}$, for realistic values of $\Gamma^{\prime}$. This is the same scaling that is predicted in free-space atomic ensembles, when atoms are assumed to emit independently Lukin (2003); Fleischhauer et al. (2005). In contrast, for some idealized system without loss (i.e., $\Gamma^{\prime}=0$), the bandwidth-delay product scales simply with the number of atoms, $\mathcal{P}\sim N$ Chang et al. (2011); Scheuer et al. (2005). This result does not follow from the perturbative expansion of Eq. (24). Rather, one can perform an exact calculation of the optical band structure and the bandwidth, $\Delta_{\text{EIT}}\sim\Omega_{c}^{2}/\Gamma_{\rm 1D}$, all of which is usable in the absence of loss Chang et al. (2011).

Figure 14(a) depicts a representative transmittance spectrum for a chain of $N=20$ atoms (red curve). For $\Delta-J^{\prime}=0$, i.e., when the probe field is in resonance with the (shifted) $\mathinner{|{e}\rangle}-\mathinner{|{g}\rangle}$ transition, the transmittance is perfect. However, total transparency is only exactly achieved at this precise frequency, decreasing with the detuning from resonance. The medium can be considered roughly transparent within a small window of bandwidth $\sim\Delta_{\text{EIT}}$, for which the transmittance $T>1/e$. The scaling of the bandwidth-delay product with the number of atoms is shown in Figure 14(b). The numerical results (red dots) are obtained by solving Eq. (23), then calculating the transmission as a function of the atom-probe detuning, and finally numerically finding the values of $\Delta$ where the transmittance drops to $1/e$. The calculations follow perfectly the simple scaling $\mathcal{P}=\sqrt{D}$ derived above (continuous red line). We should note that the usual definition of $D$ –in terms of exponential reduction of transmittance on resonance– does not apply any more to EIT. However, we maintain the definition $D=2N\Gamma_{\rm 1D}/\Gamma^{\prime}$, as it represents a physical resource.

As the final step in our summary, we turn our attention to the problem of the efficiency of an EIT-based quantum memory. Qualitatively, the large bandwidth-delay product associated with an optically dense ensemble enables an incident pulse to become spatially compressed and localized completely within the ensemble, while propagating with a reduced group velocity $v_{\text{g}}\ll c$. The slow group velocity is associated with the photon mixing strongly with a collective spin excitation $\hat{\sigma}_{sg}$, to form a “dark-state” polariton. Once the pulse is completely inside, the pump field can be adiabatically decreased to zero ($\Omega_{c}=0$), in which case $v_{\text{g}}\rightarrow 0$ and the pulse becomes stored, while simultaneously the polariton becomes a pure spin excitation Lukin (2003); Fleischhauer et al. (2005). This process of photon mapping can be reversed by ramping up the control field intensity at a later time, which allows for an “on demand” retrieval of the stored photon. Gorshkov and co-workers demonstrated that, due to time reversal symmetry, the optimal efficiency of photon storage is identical to that of photon retrieval Gorshkov et al. (2007a). Therefore, our discussion will focus on the latter. Neglecting collective emission into non-guided modes, Gorshkov et al. predicted that any smooth spin wave fitting inside the atomic medium should be retrieved with error $\varepsilon\sim 1/D$ Gorshkov et al. (2007a, b). The reason for such scaling is that the optical depth sets the branching ratio between emitting the photon into the desired channel (the guided mode of the fiber) and into the undesired reservoir (free space).

In order to demonstrate that our calculations match the previously known results, we initialize a single-excitation spin wave of the form

[TABLE]

where $\mathcal{N}$ is a normalization constant, and the phase $e^{{\rm i}k_{\rm 1D}z_{j}}$ guarantees retrieval of this excitation as a photon in the right-propagating guided mode, Eq. (18a). This peculiar-shaped spin wave (in particular, the relative population of atom $j$ grows like $j^{2}$) presents a balance between the pulse being sufficiently smooth (such that its wave vector components fit within the transparency window), and the majority of the pulse sitting at the forward end of the medium, such that it does not accumulate propagation losses over a large distance. In the limit of large optical depth, such a polariton is predicted to be of the optimal shape to yield maximal retrieval efficiency (in particular, $\varepsilon\simeq 5.8/D$ Gorshkov et al. (2007b)). At $t=0$ we switch on the control field and let the atomic wave function evolve under the effective Hamiltonian $\mathcal{H}_{\text{1D}}+\mathcal{H}^{\prime}_{\text{indep}}+\mathcal{H}_{\text{c}}$ until no excitation is left in the atomic chain (having been emitted into the waveguide or free space). We calculate the infidelity in the photon retrieval in two different manners, which yield identical results. The first method consists in integrating over time the radiative emission into non-guided photonic modes. The error is thus $\varepsilon=\int_{0}^{\infty}\,dt\,\kappa^{\prime}_{\text{indep}}(t)$, where $\kappa^{\prime}_{\text{indep}}(t)=-(2/\hbar)\,\text{Im}\mathinner{\langle{\mathcal{H}^{\prime}_{\text{indep}}}\rangle}=\Gamma^{\prime}\mathinner{\langle{\sum_{j=1}^{N}\hat{\sigma}_{ee}^{j}}\rangle}$. The second way to calculate the infidelity is to realize that a successful retrieval occurs whenever the photon is emitted to the guided mode of the fiber. Then, the error is $\varepsilon=1-\int_{0}^{\infty}\,dt\,\kappa_{\text{1D}}(t)$, where the time-dependent decay rate into the guided mode is $\kappa_{\text{1D}}(t)=-(2/\hbar)\,\text{Im}\mathinner{\langle{\mathcal{H}_{\text{1D}}}\rangle}$. Technically, the efficiency should be calculated only accounting for emission into the preferred (right-going) direction of the fiber, using Eq. (18a). We have checked that this gives nearly an identical answer, as emission in the left-going direction is negligible. The scaling of the retrieval infidelity with the number of atoms is shown by red circles in Fig. 14(c). The numerical results agree very well with the expected scaling ($\varepsilon\simeq 5.8/D$ Gorshkov et al. (2007b) 111In Ref. Gorshkov et al. (2007b) the scaling is $\varepsilon\simeq 2.9/D$, due to a factor of 2 difference in their definition of $D$., red line) for large number of atoms. In principle, the shape of the outgoing photon can be further tailored via a time-dependent control field, but we will not treat that case here.

Now that we have reviewed the basic parameters characterizing an EIT medium as well as the fidelity of a quantum memory, we analyze how collective emission into non-guided modes modifies the relevant figures of merit. In this case, the system evolves under the effective Hamiltonian $\mathcal{H}_{\text{tot}}=\mathcal{H}_{\text{1D}}+\mathcal{H}^{\prime}+\mathcal{H}_{\text{drive}}+\mathcal{H}_{\text{c}}$, where now collective emission is taken into account through the $\mathcal{H}^{\prime}$ term (instead of the previous $\mathcal{H}^{\prime}_{\text{indep}}$). Before proceeding to the calculation of the optical properties, we would like to discuss the decay rates of the eigenstates of the system without guided-mode driving, i.e., of $\mathcal{H}_{\text{1D}}+\mathcal{H}^{\prime}+\mathcal{H}_{\text{c}}$. Due to the presence of the $s$-states, the number of eigenstates in the single excitation subspace is $2N$. If the population of the $s$-states of a given eigenstate is larger than that of the $e$-states, we say that this eigenstate belongs to the “$s$-branch”, and vice versa. For any finite control field, there is mixing between the $e$- and $s$-branches, meaning that the eigenstates do not purely consist of $\mathinner{|{e}\rangle}$ or $\mathinner{|{s}\rangle}$ states.

In Fig. 15(a) we show the guided [$\Gamma_{\rm 1D}(k_{z})=-(2/\hbar)\,\text{Im}\mathinner{\langle{\mathcal{H}_{\text{1D}}}\rangle}$] and non-guided [$\Gamma^{\prime}(k_{z})=-(2/\hbar)\,\text{Im}\mathinner{\langle{\mathcal{H}^{\prime}}\rangle}$] decay rates of the numerically calculated eigenstates that belong to the $s$-branch, for a fixed number of atoms $N=200$. As in Section III A, we have performed a finite Fourier transform to associate an effective wave vector $k_{z}$ to each of the atomic spin eigenstates. As expected, the non-guided decay rates are negligible when the dominant wave vector $k_{z}$ lies beyond the light line. On the contrary, the guided decay rates peak strongly outside the light line, at $k_{z}=\pm k_{\rm 1D}$. It can also be seen that these same states experience a decay rate into free space of $\Gamma^{\prime}(k_{z})/\Gamma_{0}\ll 1$, and are thus the “selectively radiant” states that we previously anticipated. Some of the eigenstates with $|k_{z}|<k_{0}$ have a non-zero $\Gamma_{\rm 1D}(k_{z})$ decay rate into the guided mode. This occurs because the eigenstates are not purely Bloch waves with a perfectly-determined $k_{z}$, but instead can have some finite contributions from all $k_{z}$. As a technical note, we remark that only when the chain of atoms is infinite can $\mathcal{H}_{\text{1D}}$ and $\mathcal{H}^{\prime}$ be simultaneously diagonalized. For any finite number of atoms, the eigenstates of $\mathcal{H}_{\text{1D}}+\mathcal{H}^{\prime}+\mathcal{H}_{c}$ are not simultaneously eigenstates of its guided and non-guided parts.

One can also consider the behavior of the selectively radiant states, as a function of atom number. In particular, of interest is the maximum possible branching ratio $\Gamma_{\rm 1D}(k_{z})/\Gamma^{\prime}(k_{z})$ of all the eigenstates, as a function of $N$. We plot this quantity in Fig. 15(b), where we find an approximate scaling of max$\{\Gamma_{\rm 1D}(k_{z})/\Gamma^{\prime}(k_{z})\}\propto N^{2}$. We find that this scaling is in fact independent of the magnitude of the control field, and is in contrast to the $\sim N$ scaling in the case of the independent emission model. We later show that this same scaling manifests itself in the photon storage/retrieval error probabilities.

Let’s now calculate EIT transmittance spectra. Under the same conditions as before (low saturation, uniform control field), the evolution equations for the state amplitudes are found to be

[TABLE]

While analytical approximations are not as readily obtained, the numerical procedures follow exactly as presented for the case of independent emission. The blue curve in Fig. 14(a) shows how the transmittance spectrum is modified by selective radiance. The first noticeable consequence of collective suppression of the emission into non-guided modes is that the transparency window becomes wider, as expected if the loss becomes smaller. This is further confirmed in Fig. 14(b), which displays a linear scaling of the bandwidth-delay product with the atom number, in contrast to the conventional square root dependence. As mentioned before, such a scaling is characteristic of a system without losses Scheuer et al. (2005). This scaling, along with the conclusion in Sec. III that suppression of emission into free space can occur for low densities of excitations, suggests that it might be possible to store a number of photons in an atomic medium that scales linearly with atom number (in contrast to the $\sim\sqrt{D}$ scaling within the independent emission model).

Finally, Fig. 14(c) shows the improvement in the infidelity of retrieval of the spin wave given by Eq. (26). The error is now calculated including collective emission as $\varepsilon=\int_{0}^{\infty}\,dt\,\kappa^{\prime}(t)$, where

[TABLE]

Again, this error matches the one calculated by taking into account the component of the photon that is released into the guided mode, i.e., $\varepsilon=1-\int_{0}^{\infty}\,dt\,\kappa_{\text{1D}}(t)$, with $\kappa_{\text{1D}}(t)=-(2/\hbar)\,\text{Im}\mathinner{\langle{\mathcal{H}_{\text{1D}}}\rangle}$. As before, we have checked that emission in the left-going direction is negligible. We find that by exploiting collective emission, the error decreases with atom number like $\varepsilon\propto 1/N^{2}$. This result is consistent with the scaling of branching ratios for the most selectively radiant eigenstates, previously plotted in Fig. 15(b). Moreover, by varying the radial positions of the atoms over a limited range, thus modifying the ratio $\Gamma_{\rm 1D}/\Gamma^{\prime}$, we are able to separate the contributions of the number of atoms and the optical depth to the infidelity. We obtain $\varepsilon\simeq 15/(ND)$, where the numerical prefactor is not necessarily universal, as it probably depends on the fiber properties.

An interesting question is why the error of photon storage improves ‘only’ by a factor of $N$ (from $1/N$ to $1/N^{2}$). In particular, given that single excitations in a free-space chain can experience a suppression in the emission rate of up to $1/N^{3}$, one might have expected a greater suppression of errors of up to $1/N^{4}$ in photon storage. An initial – but somewhat erroneous – guess would be to attribute this “bad scaling” to an unfavorable spatial profile of the initial spin wave. Perhaps surprisingly, although EIT nominally matches the effective guided mode indices of the bare fiber and the composite system of fiber and atomic chain, we show in the next subsection that the slight impedance mismatch away from perfect resonance $\Delta=J^{\prime}$ is still responsible for the majority of scattering losses into free space. We thus present an improved impedance matching scheme, which allows for exponential improvement with $N$ of the quantum memory fidelity.

IV.4 Quantum memory with exponential fidelity

The importance of residual impedance mismatch can be seen in a simple example, where one considers an initial Gaussian spin-wave profile,

[TABLE]

and investigates the dynamics of the retrieval process more carefully. In the above expression, $\mathcal{N}$ is a normalization constant, $z_{c}=(N-1)d/2$ is the center of the atom chain, and $\sigma=\sqrt{N}d$ is the standard deviation of the Gaussian spin wave. Figure 16(a) shows the evolution of the spin excitation at different times $\tau$, for a chain of $N=200$ atoms. Here, to aid in visualization, we have defined a re-scaled dimensionless time $\tau\in[0,1]$, where $\tau$ represents the total amount of atomic population that has decayed (i.e., at $\tau$=1 the spin wave has fully decayed and the photon has been completely released). We have plotted not only the population in the $\mathinner{|{s}\rangle}$-state through the ensemble, which essentially matches population of the dark-state polariton, but also the excited state population, which is ultimately responsible for any emission into free space.

For a spatially uniform control field, and for a theory of EIT within a uniform medium (i.e., where the atomic density is treated as smooth rather than discrete points), it can readily be shown Fleischhauer and Lukin (2002) that the excited state population is proportional to $|\partial_{z}\hat{\sigma}_{gs}(z)|^{2}$ (also see Appendix E), a result that also agrees well with our numerical results. This excited state population is necessarily associated with a pulse of finite extent or bandwidth, and in complementary ways reflects the fact that perfect transparency in EIT only occurs at a single frequency, or that there are non-adiabatic corrections to the formation of a dark-state polariton. At certain times, such as at $\tau=0.07$ or $\tau=0.9$, the excited-state spin wave has a large amplitude at the edge of the atomic chain. Most of the error on retrieval occurs at these times, as can be seen in Fig. 16(b), which shows the instantaneous loss $\kappa^{\prime}(\tau)$. Here, $\kappa^{\prime}(\tau)$ is re-scaled as well, so that its integral provides the total loss, $\int_{0}^{1}d\tau\kappa^{\prime}(\tau)=\varepsilon$.

A plausible cause of this behavior is the discontinuity of the excited state population at the system’s end, which is associated with a large number of wave vector components $k_{z}$ that lie within the light line and couple to free space. To further confirm this intuition, we have developed a model for the time-dependent loss based upon the Fourier decomposition of the spatial profile of the excited state amplitude. At every time of the evolution, we calculate $c_{e}(k_{z},\tau)$ by doing a finite Fourier transform of the excited state amplitudes, $c_{e}(z_{j},\tau)\equiv c_{e}^{\,j}(\tau)$. Then, we find the Fourier-based instantaneous loss as $\tilde{\kappa}^{\prime}(\tau)=(d/2\pi)\int_{-k_{0}}^{k_{0}}dk_{z}\,\Gamma^{\prime}(k_{z})|c_{e}(k_{z},\tau)|^{2}$, where $\Gamma^{\prime}(k_{z})$ is obtained from classifying the decay rates of the eigenstates of $\mathcal{H}^{\prime}$ according to their dominant wave vector. This calculation, represented by the dashed line in Fig. 16(b), shows good agreement with the numerics. Figure 16(c) depicts the components of $|c_{e}(k_{z})|^{2}$ inside the light line for different times [corresponding to the snapshots of Fig. 16(a)]. For initial times (purple curve), the wave function has minimal population within the light line, suggesting it propagates with little loss down the atomic chain. The population drastically increases as the pulse hits the end of the chain (brown curve). A large population of wave vector components within the light line is correlated with increased instantaneous loss and sharp features in the profile of the excited state population at the system edge. Finally, Fig. 16(d) displays the scaling of the infidelity in photon retrieval with the atom number. The scaling is poor ($\varepsilon=4.1/N$), as large losses occur when the polariton hits the ends of the atomic chain.

We now describe how to smoothly reduce the excited state population at the end of the chain, by introducing a spatially-dependent control field. Heuristically, the idea is to increase the control field at the ends of the chain, as shown in Fig. 17. As the EIT bandwidth is proportional to the control field intensity [see Eq. (25)], the atomic medium becomes more transparent at the edges, where the excited state population is reduced. One can develop an effective continuum wave equation to predict the evolution of the populations in $\mathinner{|{s}\rangle}$ and $\mathinner{|{e}\rangle}$ during the retrieval process, in the presence of a spatially dependent control field (see Appendix E) Fleischhauer and Lukin (2002). Similar to the case of a uniform control field presented earlier, in principle the scattering loss can then be estimated and minimized from these populations. This optimization process seems quite challenging in practice, however, as it depends on the initial spin wave, the control field profile, and on the integral of momentum components over the entire history of evolution. We will not do such an optimization here, but instead show that rather simple choices can already lead to significant improvements over the infidelity scalings that we found in the previous subsection.

For the initial spin wave of Eq. (29), Fig. 18(a) shows the evolution of the e- and s-states populations of a chain of $N=200$ atoms for the spatially-dependent control field $\Omega_{c}(z_{j})=\Omega_{c}^{(0)}\sqrt{N/(N+1-j)}$, which exhibits a rapid increase at the right edge of the chain [see Fig. 17]. Since the control field is switched on suddenly in our simulations, its magnitude is taken to be very small ($\Omega_{c}^{(0)}=0.005\Gamma_{0}$) to minimize rapid, non-adiabatic accelerations of the spin wave that would artificially increase the losses at initial times. Both the $\mathinner{|{e}\rangle}$- and $\mathinner{|{s}\rangle}$-state populations exhibit smooth profiles at every time of the evolution, and in particular, one sees that the excited state population smoothly vanishes at the edge of the system. The dashed lines in the plots show the results from the analytical model developed in Appendix E, which agree well with the numerics. Both the instantaneous loss and the amount of spin-wave population lying within the light line are several orders of magnitude smaller than in the case of a uniform control field, as can be seen in Figs. 18(b) and (c), respectively. Moreover, the excited state population within the light line does not significantly increase from its initial values. Finally, Fig. 18(d) shows the exponential decrease of the retrieval infidelity $\varepsilon$ as a function of the atom number $N$ for this given profile. We anticipate that optimized initial spin waves and control field profiles will result in a much steeper exponential scaling. Nonetheless, we have demonstrated that a fairly trivial selection of those settings already exponentially improves previously-known bounds for photon storage.

V Physical implementations and possible challenges

Having analyzed the physics of both subradiance and selective radiance, we devote this section to discussing suitable experimental platforms as well as the challenges that might be encountered to observe this physics.

V.1 Physical implementations

To potentially observe the physics that we described in Sec. III requires atoms to be regularly trapped, forming an ordered lattice. Moreover, among all the possible collective atomic states, one should be able to access the subradiant manifold. We shall start our discussion with possible physical platforms. As we have demonstrated, the minimal distance at which subradiant states appear depends on the dimensionality of the atomic array ($\lambda_{0}$ in 2D, and $\lambda_{0}/2$ in 1D). Standard free-space optical lattices Bloch (2005) can achieve lattice constants of $d\sim\lambda_{0}/2$, and quantum gas microscopes Bakr et al. (2009) are able to generate single 2D arrays. In such systems, both bosonic Bakr et al. (2010) and fermionic Greif et al. (2016) Mott insulator phases – where the number of atoms per site can be limited to one – have been realized.

Very recently, several experimental groups have built almost defect-free 1D Endres et al. (2016) and 2D Lester et al. (2015); Barredo et al. (2016) lattices in an atom-by-atom manner using optical tweezer arrays. While the inter-atomic distance achieved so far is still larger than the free space wavelength, due to the problem of interference between the tweezers at close distances, it could be possible that further improvements enable sub-wavelength distances to be reached. It might also be possible to employ a transition with a shorter wavelength for the trapping scheme, and use another of longer wavelength to explore subradiant phenomena. Finally, periodically patterned 1D or 2D dielectric structures can readily yield sub-wavelength trapping potentials, with the periodicity of the structure itself Hung et al. (2013); González Tudela et al. (2015). While cold atoms can now routinely be trapped near dielectric structures Vetsch et al. (2010); Thompson et al. (2013); Goban et al. (2014, 2015); Hood et al. (2016), the filling fractions remain quite low and new approaches (such as integration with tweezer arrays) must be developed to achieve near-perfect filling.

Overcoming the second requirement, that of exciting the subradiant manifold efficiently, is not trivial. As subradiant states are characterized by a wave vector that lies beyond the light line, they do not naturally couple to a laser beam that propagates through free space. There are several options to overcome this hurdle. An already-suggested possibility Plankensteiner et al. (2015) is to map superradiant states, which are easy to excite, to subradiant ones via magnetic field gradients. Specifically, a laser can efficiently excite a spin wave $\sim\sum_{j}e^{{\rm i}k_{z}z_{j}}\mathinner{|{e_{j}}\rangle}$ whose wave vector $k_{z}$ lies within the light line. In the case that the excited state is magnetic-field sensitive, a field gradient would then imprint a spatially dependent phase shift $\mathinner{|{e_{j}}\rangle}\rightarrow e^{{\rm i}(\beta t)z_{j}}\mathinner{|{e_{j}}\rangle}$ in time, which then could allow the wave vector $k_{z}\rightarrow k_{z}+\beta t$ to be mapped outside of the light line. In 1D chains in free space, one might exploit the fact that the emission of subradiant states occurs primarily from the ends, and in a pattern that can be collected reasonably well with conventional optics. Note that the question of efficient excitation does not come up with selectively radiant states, as by definition they are well-coupled to a mode of interest.

Finally, we would like to stress that the exploration of both subradiance and selective radiance is not restricted to atoms. Molecules Faez et al. (2014) and solid-state emitters should also exhibit these properties, although they pose a different set of challenges. The ability of deterministically placing quantum dots Lund-Hansen et al. (2008); Young et al. (2015); Lodahl et al. (2015), rare earth ions Kornher et al. (2016); Jacob et al. (2016), and color centers Sipahigil et al. (2016) has significantly improved in the past years, thus putting ordered arrays within reach. However, one of the main appeals of atoms is that they are identical to each other and that their decay is purely radiative. An open question is how to translate these features into the domain of solid-state emitters, as it would require high homogeneity among them as well as a large emission into the zero-phonon line (minimizing non-radiative losses).

V.2 Atomic level structure

There are a number of potential imperfections that could limit the observation of subradiance and selective radiance, and the performance of protocols that exploit them. A number of these imperfections are conceptually clear (if not necessarily straightforward to analyze), such as disorder in atomic positions, classical and quantum motion, dephasing, and imperfect site filling. A more thorough investigation of these effects will be left to future work.

Here, we would like to discuss a more subtle issue, related to the complications associated with multi-level atomic structure. For most of this manuscript, we have assumed that atoms are two-level systems, with a single ground state and excited state. For atoms with hyperfine structure (and thus a ground state manifold), an effective two-level system is often achieved by exploiting a cycling transition Metcalf and van der Straten (1999), where an excited state of maximum angular momentum can only decay back into a single ground state, also of maximum angular momentum. Such a transition only responds to pure circularly polarized light. In the case of multiple atoms, an important issue is that light scattered from one atom does not display circular polarization globally in space. Thus, for example, the resulting dipole-dipole interactions can potentially drive other atoms to excited states outside of the cycling transition. This can be avoided in the specific case of a 1D chain (where the re-scattered field has the same polarization along the axis of the chain), but not in general.

Another possibility to avoid the full complexity of hyperfine structure is to use atoms without nuclear spin, such as bosonic Ytterbium or Strontium Xu et al. (2003); SSK . In this case, there is a single ground state but three excited states with orthogonal dipole matrix elements (giving an isotropic optical response to the atoms). Then, one can exploit the fact that in 1D arrays, dipole-dipole interactions involving different excited states decouple from each other, to effectively yield two-level physics (similar to the case of circular polarization described above). In 2D arrays, the transitions involving a dipole matrix element out of the plane decouple, while the two in-plane transitions can hybridize, and calculating the band structure for an infinite system involves diagonalizing a $2\times 2$ matrix associated with the Fourier transform of the in-plane components of the Green’s function, $\tilde{G}_{0,\alpha\beta}(\textbf{k})$, with $\{\alpha,\beta\}=\{y,z\}$. This solution qualitatively maintains the same properties as the two-level case analyzed in Sec. III [for example, see the discussion surrounding Eq. (10)].

In the presence of hyperfine structure, and excluding the special case of a 1D array described previously, the complication with regard to subradiance can be understood with the following simple example. Suppose that atoms are initialized in a single ground state $\mathinner{|{g}\rangle}$, from which a single-excitation spin wave of the form $\mathinner{|{\psi}\rangle}\sim\sum_{j}e^{{\rm i}kz_{j}}\mathinner{|{e_{j}}\rangle}$ is somehow generated. If the wave vector $k$ is beyond the light line, then as argued in Sec. III, collective dissipative interactions [such as those encoded in the $\hat{\sigma}_{eg}^{i}\hat{\sigma}_{ge}^{j}$ term of Eq. (3a)] will suppress emission of an excited state back into $\mathinner{|{g}\rangle}$ through destructive interference. However, dipole-dipole interactions will also generally exist between that excited state and any other ground state $\mathinner{|{s}\rangle}$ connected by a dipole-allowed transition, e.g., of the form $\hat{\sigma}_{es}^{i}\hat{\sigma}_{se}^{j}$. Since the initially prepared spin wave $\mathinner{|{\psi}\rangle}$ does not contain any population in $\mathinner{|{s}\rangle}$, there is no interference that prevents decay via this channel, and thus the spin wave would experience a decay rate into $\mathinner{|{s}\rangle}$ equal to that of a single, isolated atom excited to $\mathinner{|{e}\rangle}$. Interesting recent work suggests that it is possible to encode subradiance in a more complex initial state beyond the simple product state $\mathinner{|{g}\rangle}^{\otimes N}$ Hebenstreit et al. (2017), which would be worth exploring further.

Within the context of the enhancement of EIT-based storage protocols through collective emission, studied in Sec. IV, this implies that the state $\mathinner{|{s}\rangle}$ cannot be another state in the ground-state manifold that is directly connected to $\mathinner{|{e}\rangle}$ by a dipole-allowed transition. Various possibilities to implement EIT and retain the desired collective interference effects include the use of a ladder scheme, with the state $\mathinner{|{s}\rangle}$ being a long-lived excited state (e.g., a Rydberg level), or to use a state $\mathinner{|{s}\rangle}$ in the ground-state manifold that is connected only through a two-photon transition van Enk et al. (1997); Porras and Cirac (2008).

VI Summary and outlook

In summary, we have shown that subradiant states acquire an elegant interpretation in 1D and 2D atomic arrays, in terms of optically guided modes whose decay rate is only limited by the system boundaries. We have provided a first glimpse into the nature of subradiance for multiple excitations, and introduced a new concept of selective radiance that should enable the construction of more efficient atom-light interfaces. As a concrete example, we have constructed a protocol for quantum memories for light using selectively radiant states in an optical nanofiber, whose infidelity decreases with atom number at a rate exponentially better than previously known bounds.

Even though memories are a very relevant quantum technology, the improvement in their performance is just an example of the bountiful possibilities spawned by subradiance and selective radiance. We anticipate that exploiting these phenomena could yield new error bounds and protocols for many applications of interest, ranging from nonlinear optics to metrology. At the same time, the nature of subradiance for multiple excitations or internal states could itself constitute a rich new many-body problem.

Acknowledgments.– We are grateful to C. Regal, A. M. Rey, A. Gorshkov, E. Polzik, S. Yelin, M. Lukin, H. Ritsch, E. Shahmoon, and J. Muniz for stimulating discussions. H.J.K. funding is provided by the AFOSR Quantum Many-Body Physics with Photons and QuMPASS MURI, NSF Grant PHY-1205729, the Office of Naval Research (ONR) Award N00014-16-1-2399; the ONR QOMAND MURI; and the IQIM, an NSF Physics Frontiers Center. A.A.-G. was supported by the IQIM Postdoctoral Fellowship and the Global Marie Curie Fellowship LANTERN (655701). D.E.C. acknowledges support from Fundació Privada Cellex, Marie Curie CIG ATOMNANO, Spanish MINECO Severo Ochoa Programme SEV-2015-0522, MINECO Plan Nacional Grant CANS, CERCA Programme/Generalitat de Catalunya, and ERC Starting Grant FOQAL.

Appendix A Analytical expressions for the infinite lattice case

In this appendix we provide a formal derivation for the collective frequency shifts and decay rates of an infinite 1D linear and 2D square lattice of two-level atoms. As discussed in the main text, in an infinite lattice the effective Hamiltonain $\mathcal{H}_{\textrm{eff}}$ has discrete translational invariance and it can be written as:

[TABLE]

Here $S^{(\dagger)}_{\textbf{k}}=N^{-1/2}\sum_{j}e^{-(+){\rm i}\textbf{k}\cdot{\bf r}_{j}}\hat{\sigma}^{j}_{ge(eg)}$ represents the annihilation (creation) operator of one of the collective modes, which in this case corresponds to an atomic spin wave with momentum k (defined within the first Brillouin zone). The quantites $J_{\textbf{k}}$ and $\Gamma_{\textbf{k}}$ are real and correspond to the collective frequency shift and decay rate of the mode k, calculated by taking the discrete Fourier transform of Eq.(4), that is

[TABLE]

A.1 Collective frequency shifts

Here we evaluate the collective frequency shifts for an infinite 1D chain of atoms along the $\hat{z}$ direction. The collective frequency shifts when the atoms are polarized along or transversally to the chain (denoted by $J^{||}_{k_{z}}$ and $J^{\perp}_{k_{z}}$, respectively) are given by:

[TABLE]

Using the representation of the PolyLogarithmic function as an infinite sum: $\textrm{Li}_{s}(z)=\sum_{\ell=1}^{\infty}z^{\ell}\,\ell^{-s}$ [for $z\in\mathbb{C}$, and $\textrm{Li}_{s=1}(z)=-\textrm{ln}(1-z)$], the previous expressions reduce to Eqs.(8a,8b) in the main text.

A.2 Collective decay rates

In order to compute the Fourier transform $\tilde{{\bf G}}_{0}(\textbf{k})$ it is useful to express the free space Green’s tensor in terms of the spherical wave function:

[TABLE]

and make use of the spherical wave decomposition into plane waves:

[TABLE]

with $Q_{x}=\sqrt{k_{0}^{2}-\textbf{Q}^{2}_{||}}$. Here, we have chosen the axis $\hat{x}$ to be perpendicular to the array, while the components ${\bf r}_{||}$ are parallel to a plane that contains the atomic chain (in 1D) or the atomic array (in 2D). Then, the free space Green’s tensor can be expressed as

[TABLE]

where we have defined $\bar{\textbf{Q}}\equiv(Q_{x}\textrm{sign}(x),Q_{y},Q_{z})$. For 1D and 2D lattice geometries we can choose that the corresponding line or plane of atoms sits at $x=0$. Then, since we are only interested in evaluating ${\bf G}_{0}$ at the atomic positions, we can set $x\rightarrow 0$ in the above expressions. Making use of the Dirac delta representation in D dimensions,

[TABLE]

it is possible to express $\tilde{{\bf G}}_{0}(\textbf{k})$ as a sum over reciprocal lattice vectors ${\bf g}$. For a 1D linear chain along $\hat{z}$,

[TABLE]

where we have defined $\textbf{q}\equiv(q_{x}\textrm{sign}(x),Q_{y},k_{z}+g_{z})$, with $q_{x}=\sqrt{k_{0}^{2}-(k_{z}+g_{z})^{2}-Q_{y}^{2}}$. In the case of a 2D square lattice in the $\hat{y}$-$\hat{z}$ plane this reads:

[TABLE]

with $\textbf{q}\equiv(q_{x}\textrm{sign}(x),k_{y}+g_{y},k_{z}+g_{z})$ and $q_{x}=\sqrt{k_{0}^{2}-|\textbf{k}+{\bf g}|^{2}}$.

Eqs.(10a) and (10b) are ill-defined for the crossed components $xy$ and $xz$ of $\tilde{{\bf G}}_{0}(\bf k)$. However, we note that $G_{0,x\alpha}({\bf r}_{\perp},x\rightarrow 0)=0$ ($\alpha=y,z$), since the electromagnetic field emitted by a dipole is always parallel to the dipole itself, at any point of its normal plane. Thus, also $\tilde{G}_{0,x\alpha}(\textbf{k})=0$ ($\alpha=y,z$). The fact that this crossed term vanishes is relevant, since it implies that the modes with transverse and in-plane polarization will not be mixed when dealing with multi-level atoms. This is true specifically for 1D and 2D lattices.

From Eq.(31b) and Eqs.(37,38) the collective decay rates can be evaluated. We first note that for a wave vector beyond the light line (i.e., $|\textbf{k}|>k_{0}$), $q_{x}$ is purely imaginary for any reciprocal lattice vector ${\bf g}$. Thus, the imaginary part of all diagonal tensor components vanishes, and the decay rates are exactly zero. This mathematically demonstrates that any state beyond the light line is necessarily subradiant. In order to have states satisfying this condition, the maximum magnitude of the wave vectors defined in the first Brillouin zone must be larger than the one defining the light line, i.e., $k^{\textrm{max}}>k_{0}$. In a linear chain in 1D, one has $k^{\textrm{max}}=\pi/d$, and this sets the condition $d/\lambda_{0}<1/2$ for the existence of these states. In a 2D square lattice, for which the first Brillouin zone is a square, $k^{\textrm{max}}=\sqrt{2}\pi/d$, yielding instead the condition $d/\lambda_{0}<1/\sqrt{2}$.

For states with wave vector $|\textbf{k}|\leq k_{0}$ the only contribution in the decay rate is from reciprocal lattice vectors satisfying $|\textbf{k}+{\bf g}|\leq k_{0}$. In the 1D case we obtain, for parallel and transverse polarization, and after performing the integral in Eq.(37):

[TABLE]

For the 2D square lattice one gets:

[TABLE]

Appendix B Transfer Matrix Formalism

In Sec. III.2, we found that a linear chain of $N$ atoms has a set of subradiant single-excitation modes, whose decay rates scale like $\Gamma_{\xi}\sim\xi^{2}/N^{3}$. Here $\xi=1,2,3,...$ serves as an index for the subradiant modes.

Here we present a simple one-dimensional model of light interacting with a periodic system of scatterers. It is important to note that one cannot establish a formal mapping from the original system to this one. Heuristically, however, one might hope that the simple model is sufficient to capture the salient features of a generic pseudo-1D system. In particular, we will find that the simple model also produces a set of resonances, whose decay rates scale like $\Gamma_{\xi}\sim\xi^{2}/N^{3}$.

One-dimensional scattering through several optical elements (such as an array of scatterers) can be efficiently described using the Transfer Matrix formalism Bendickson et al. (1996); Joannopoulos et al. (2008). This method takes advantage of the fact that in a one-dimensional scattering model there are only two propagation directions (left and right). The transfer matrix $\mathcal{M}_{\textrm{sc}}$ (see Fig. 19) relates the fields on one side ($E_{R}^{-}$, $E_{L}^{-}$) and on the other side ($E_{R}^{+}$, $E_{L}^{+}$) of a point scatterer:

[TABLE]

Propagation through a unit cell can also be described via a transfer matrix, which itself is a product of transfer matrices describing interaction with the point scatterer (described by reflection and transmission coefficients $r$ and $t$), and the one-dimensional free-space propagation at frequency $\omega=ck$ over a distance $d$. That is, $\mathcal{M}=\mathcal{M}_{\text{sc}}\cdot\mathcal{M}_{\text{free}}$, with

[TABLE]

and $\zeta=-{\rm i}r/t$. The relation of $\mathcal{M}_{\textrm{sc}}$ to $\zeta$ is determined by the additional constraint that $1+r=t$, which states that the field should be continuous on each side of the point scatterer.

It is useful to decompose the matrix $\mathcal{M}$ as:

[TABLE]

with

[TABLE]

since $\textrm{Tr}\mathcal{A}=0$ and $\mathcal{A}^{2}=\mathbb{1}$. Here $v_{g}$ is the group velocity for an infinite lattice, and it will be given by Eq.(53).

Dispersion relation.– The two eigenvalues that diagonalize the transfer matrix $\mathcal{M}$ are necessarily of the form $\lambda_{\pm}=e^{\pm{\rm i}qd}$, since $\det\mathcal{M}=1$. As we will see, $\pm q$ corresponds to the Bloch index or quasi-momentum. Moreover, since the trace of $\mathcal{M}$ is independent of the basis, one can obtain the dispersion relation (the relationship between $\omega$ and $q$):

[TABLE]

Given the quasi-momentum value $|q|$, there are two possible solutions (or branches) that fulfill this relation. In particular,

[TABLE]

Moreover, from Eq.(51) the group velocity for an infinite system can be derived:

[TABLE]

Group velocity close to the band edge.– At the band edge of the Brillouin zone ($q=\pi/d$) the dispersion relation exhibits a band gap, with the lower ($\omega_{0-}$) and upper ($\omega_{0+}$) frequencies of the gap given by:

[TABLE]

Near the band edge the dispersion relation can be approximated by

[TABLE]

and the group velocity can be identified as $v_{g}=\mp c(\pi-qd)/\zeta$.

Finite array. Transmission coefficient and resonances.– In an ordered array of $N$ point scatterers (separated by the distance $d$) the total transfer matrix is simply the product of matrices $\mathcal{M}_{N}=\mathcal{M}^{N}$. Thus, the eigenvectors of $\mathcal{M}$ are also eigenvectors of $\mathcal{M}_{N}$. Eq.(47) is very useful to compute $\mathcal{M}_{N}=e^{{\rm i}Nqd\mathcal{A}}=\cos(Nqd)\mathbb{1}+{\rm i}\sin{(Nqd)}\mathcal{A}$ Deutsch et al. (1995).

As any transfer matrix, $\mathcal{M}_{N}$ can be written as:

[TABLE]

where now $t_{N}$ and $r_{N}$ represent the reflection and transmission coefficients throughout the whole array. Thus, one can obtain the transmission and reflection coefficients from the elements $\mathcal{M}_{N}^{22}$ and $\mathcal{M}_{N}^{12}$:

[TABLE]

The previous expressions allow us to identify resonances of the finite array. Indeed, for particular values of the Bloch index, namely $q_{\xi}d=\pi(N-\xi)/N$ (where $\xi$ is an integer number), the transmission coefficient $t_{N}(q_{\xi})\rightarrow(-1)^{N-\xi}$. That is, the transmission probability through the array is maximum and equal to one.

It is also interesting to analyze the physics close to those resonances. As discussed, the transmission spectra exhibits peaks at each value $q_{\xi}$. Around each peak $\xi$, it can readily be shown that the transmission spectrum behaves approximately as a Lorentzian, $t_{N}\approx(-1)^{(N-\xi)}({\rm i}\Gamma_{\xi}/2)/({\rm i}\Gamma_{\xi}/2-\delta\omega_{\xi})$, where $\delta\omega_{\xi}$ denotes the detuning from the resonance frequency of mode $\xi$ and $\Gamma_{\xi}$ its linewidth. For small values of $\xi$, one finds approximately that

[TABLE]

Appendix C Green’s function of a nanofiber

In this section we provide the expressions for the radial components of the Green’s function of an infinite nanofiber of radius $r$ directed along $\hat{z}$, following Klimov and Ducloy Klimov and Ducloy (2004). The total field produced by a dipole near a fiber can be expressed in terms of a free field (solution in vacuum) and a field re-scattered by the fiber. To exploit separation of variables, the free field can be written in cylindrical coordinates as an expansion in longitudinal wave vector $k_{\parallel}$ and angular momentum $e^{im\phi}$. The full Green’s function $\textbf{G}({\bf r},{\bf r}^{\prime},\omega_{0})$ is a rather complicated expression. Here, we are interested in the case where all atoms sit at identical distances from the fiber and a fixed azimuthal angle around the cylinder. Since we only need the Green’s function at the atomic positions themselves, we can construct a simplified version of the scattering Green’s function, only evaluated at the atomic positions $z_{j}$:

[TABLE]

where

[TABLE]

In the above expression, $\rho_{a}>r$ is the radial position of the atoms (assumed to be identical), $k_{\perp}=\sqrt{k_{0}^{2}-k_{\parallel}^{2}}$ is the perpendicular component of the wave vector outside the fiber, and $H^{(1)}_{m}$ is the Hankel function of the first kind. The coefficients $\,a_{m}(k_{\parallel})$ and $\,b_{m}(k_{\parallel})$ are found by matching boundary conditions for the electromagnetic field at the surface of the fiber, and can be taken from Eqs. (47)-(50) in Ref. Klimov and Ducloy (2004) (by choosing the value of the dipole moment equal to unity).

In the following, we detail how to perform the integral in Eq. (61). We only need to do it for the case where $z_{j}>z_{k}$, as $G_{\rho\rho}(z_{j},z_{k})=G_{\rho\rho}(z_{k},z_{j})$, due to Lorentz reciprocity. In the complex plane, the integrand $\tilde{G}_{m}(k_{\parallel})$ has two branch cuts (due to the square root form of $k_{\perp}$) and two poles, that correspond to the guided mode of the fiber (for a small enough nanofiber, all the guided modes with $|m|\neq 1$ are cut-off). This guided mode, the so-called HE11Yeh and Shimabukuro (2008), does not have a cut-off frequency and corresponds to the $m=\pm 1$ pole. In particular, the variation in the position of the pole with $\omega_{0}$ gives rise to the dispersion relation. In order to avoid the poles when integrating along the real line, we perform a contour integration and employ Cauchy’s theorem. The contour that we choose is shown in Fig. 20. Based on this image, we have $I=I_{\rm pole}-I_{r}-I_{+}-I_{-}-I_{c}-I_{cc}$, where $I$ is given by Eq. (61), and the other integrals are performed along the contours of Fig. 20. After performing the integrals, we find that the total Green’s function can be separated into guided, $G^{\rm 1D}$, and non-guided, $G^{\prime}$, contributions as

[TABLE]

In the above expressions, $k_{\rm 1D}$ is the wave vector of the guided mode, $G_{0,\rho\rho}(z_{j},z_{k})$ is the vacuum’s Green’s function, and $C_{pole}$ and $C_{cc}$ are the red and left-most green contours in Fig. 20, respectively. For the non-guided Green’s function, the integral $I_{r}$ produces both frequency shifts and decay, whereas the integrals $I_{\pm},I_{\text{c}},I_{\text{cc}}$ only contribute to the frequency shifts. The integrals $I_{\text{c}}+I_{\text{cc}}$ do not contribute for $z_{j}\neq z_{k}$. For the local Green’s function (i.e., $z_{j}=z_{k}$), both $I_{+}+I_{-}$ and $I_{\text{c}}+I_{\text{cc}}$ are divergent individually, but the infinity is cancelled when they are added.

Appendix D Linear optics for two-level atoms in the mirror configuration

Within the independent emission model, and when the atoms are spaced at distances such that $k_{\rm 1D}d=n\pi$, with $n$ being an integer, only a single collective atomic mode couples to the fiber. This case constitutes the so-called “mirror configuration”, as it has been shown that the ensemble behaves as a nearly perfect mirror with increasing atom number $N$. As we will show in this section, the resultant physics when suppression of emission into non-guided modes is accounted for becomes significantly more complicated. In particular, while one can dramatically enhance the reflectance of the atomic chain, we find that the atoms no longer respond as a single mode, and that the impedance mismatch between the atomic chain and the photonic guided mode of the fiber is a major issue.

For atoms spaced by a distance $k_{\rm 1D}d=2\pi$, the guided Hamiltonian of Eq. (12) simply reads

[TABLE]

In the single excitation manifold there is only one superradiant mode, with decay rate $N\Gamma_{\rm 1D}$, while all the others are completely dark (the same physics would be observed for any other separations of the form $k_{\rm 1D}d=n\pi$). In this configuration, the transmittance and reflectance coefficients can be found analytically, and read Chang et al. (2012)

[TABLE]

and the photon-loss probability is $\kappa_{\text{indep}}=1-T_{\text{indep}}-R_{\text{indep}}=2R_{\text{indep}}\Gamma^{\prime}/N\Gamma_{\rm 1D}$. The transmittance spectrum is a Lorentzian whose linewidth $N\Gamma_{\rm 1D}+\Gamma^{\prime}$ grows linearly with number of atoms, for sufficiently large $N$. On resonance (when $\Delta-J^{\prime}=0$), the atomic chain becomes a very good mirror, and the only relevant quantity that determines how much light is reflected is the ratio $D=2N\Gamma_{\rm 1D}/\Gamma^{\prime}$, which is in fact the optical depth of the system. In particular, on resonance and in the limit of large optical depth, the transmittance, reflectance, and loss become $T_{\text{indep}}=4/D^{2}$, $R_{\text{indep}}=1-4/D$, and $\kappa_{\text{indep}}=4/D$, respectively.

For atoms at close distances from each other, the independent emission model is not valid, and one cannot find analytical expressions for the transmission and reflection coefficients. Figures 21(a-c) show the transmission, reflection, and loss probability spectra of a chain of $N=20$ atoms coupled to the nanofiber, for both the independent emission model (red curves) and the collective emission calculation (blue curves). The differences are striking. For instance, close to resonance, the transmission and the loss decrease about four orders of magnitude when collective suppression into free-space is taken into consideration.

Repeating these calculations for chains with different number of atoms, we extract the scalings of the minimum transmittance and maximum reflectance within the collective emission model. In particular we find $T\sim 1/N^{8}$ and $1-R\sim 1/N^{6}$. Moreover, at the detuning that minimizes emission into free space, the loss scales as $\kappa\sim 1/N^{6}$. It thus seems apparent that collective emission cannot be captured by some trivial modification of the independent emission model [e.g., one cannot simply replace $\Gamma^{\prime}$ by some $\Gamma^{\prime}_{\text{eff}}(N)$ in Eqs. (65a,65b)]. In other words, the atom-light coupling can no longer be described by a single collective atomic eigenstate, but has become instead a multimode problem. This idea is further confirmed by a careful analysis of the lineshapes of Figs. 21(a-c). Although hard to appreciate in the figures, the transmittance and reflectance spectra are not smooth Lorentzians due to interference between different atomic eigenstates, all of which contribute to the optical response of the chain.

Beyond the added complexity, there is another important issue to notice: Fig. 21(c) shows that off-resonant losses are still quite large. We find that far away from resonance, the loss is independent of the atom number. Figure 21(d) sheds light on the reasons for such behavior, as it compares the excited state profile at the detuning where the loss is minimal [orange arrow and circles in (c) and (d), respectively] and far off resonance, at $\Delta-J^{\prime}=-20\Gamma_{0}$ [green arrow and circles in (c) and (d), respectively]. At the detuning where the loss is minimal, the atoms at the end of the chain are negligibly excited, and the atomic population smoothly increases toward the middle of the chain. In other words, the response of the atoms to the incoming field appears “impedance matched,” in that the smooth excitation profile reduces the amount of spin-wave components that sit within the light line and couple to free-space radiation. In contrast, far off resonance, the prepared state mostly builds upon a single eigenstate with a very uniform spatial profile. This further supports the argument in the main text, that the response of light at the interface between the atoms and the bare fiber plays an important role in the observed scattering losses.

Appendix E Polariton model

In this appendix, we develop a model for the dynamics of the dark- and bright-state polaritons where losses into free space are completely neglected. We extend previous theory for EIT in continuous atomic media Fleischhauer and Lukin (2002), in order to study the effect of spatially-dependent control fields on the dynamics of the bright polariton. Instead of focusing on the spin-model equations only, and reconstructing fields using an input-output equation, we return to explicitly keeping track of the wave equation of the electric field as it propagates through the fiber. Employing continuous atomic operators (denoted by tildes), the equations of motion are:

[TABLE]

In the above equations, $n=1/d$ is the smoothed-out linear density associated with atoms spaced at distance $d$, and a phase $e^{{\rm i}k_{\rm 1D}z}$ has been incorporated into the field and atomic coherence operators to make them slowly varying in space. All operators depend on $z$ and $t$, and $\Omega_{c}=\Omega_{c}(z)$ is taken to be real. We now respectively define the dark and bright-state polaritons as

[TABLE]

where the mixing angle is given by $\tan\theta=\sqrt{cn\Gamma_{\rm 1D}/2\Omega_{c}^{2}}$. In the adiabatic and slow-light ($v_{g}\ll c$) limits, the equations of motion for these polaritons are

[TABLE]

We consider that the bright-state polariton only perturbatively affects the dynamics of the dark-state polariton Douglas et al. (2016). Therefore, as a first approximation, we can set $\Phi(t,z)=0$ and solve for $\Psi(t,z)$. Then, the equation of motion of the dark-state polariton reads

[TABLE]

Plugging the above expression into Eq. (72), we readily find that the bright-state polariton follows

[TABLE]

The formal solution of the equation for the dark polariton is $\Psi(t,z)=\sqrt{c/v_{g}(z)}\,\tilde{f}\left(t-\int_{0}^{z}v_{g}^{-1}(z^{\prime})dz^{\prime}\right)$, where $\tilde{f}$ is a function that fulfils the condition $\tilde{f}\left(-\int_{0}^{z}v_{g}^{-1}(z^{\prime})dz^{\prime}\right)=\sqrt{v_{g}(z)/c}\,\Psi(0,z)$. In the slow-light limit, the dark-state polariton is nearly a pure spin-wave excitation, and thus corresponds to the spatial profile of the s-state spin wave at t=0. In order to obtain $\tilde{f}$ we need to perform an inversion. How hard this function inversion is depends on the profile of $v_{g}(z)$ and on the initial dark-state polariton shape. Introducing this expression into the equation for the bright-state polariton, we find

[TABLE]