Exponential improvement in photon storage fidelities using subradiance and “selective radiance” in atomic arrays
Abstract
A central goal within quantum optics is to realize efficient, controlled interactions between photons and atomic media. A fundamental limit in nearly all applications based on such systems arises from spontaneous emission, in which photons are absorbed by atoms and then rescattered into undesired channels. In typical theoretical treatments of atomic ensembles, it is assumed that this rescattering 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 for light or photonic quantum gates. However, this assumption can in fact be dramatically violated. In particular, it has long been known that spontaneous emission of a collective atomic excitation can be significantly suppressed through strong interference in emission between atoms. While this concept of “subradiance” is not new, thus far the techniques to exploit the effect have not been wellunderstood. 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 elegant interpretation in terms of optical modes that are guided by the array, which only emit due to scattering from the ends of the finite system. We also go beyond the typically studied regime of a single atomic excitation, and elucidate the properties of subradiant states in the manyexcitation 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 atomlight 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 numerically that selectively radiant states allow for a photon storage error that scales exponentially better with number of atoms than previously known bounds.
pacs:
42.50.Ct, 42.50.NnI 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 proofofprinciple 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 quantumenhanced 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 manybody 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 rescattered into other inaccessible modes or channels. Within the context of quantum lightmatter 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” of the medium as or slower. The optical depth is given by , where is the atom number, and represents the interaction probability between a single atom and a single photon in the preferred optical mode ( being the wavelength associated with the atomic transition and the beam area). Intuitively, the (or ) 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 atomlight 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 and excited state 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 spinwave excitations that constitute these subradiant states have associated wave vectors that are mismatched from freespace radiation fields, which consequently prevents the decay of energy from these states.

In the case of a finite array, a set of singleexcitation collective atomic modes can exhibit decay rates which are polynomially suppressed with atom number . 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 singleexcitation manifold, and investigate the nature of multiexcitation 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 freespace radiation due to the “hard core” nature of spins. In 1D arrays, we find that a “fermionic” ansatz works well to describe multiexcitation subradiant modes, where these multiexcitation states are constructed from antisymmetric combinations of singleexcitation 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 freespace 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 atomlight 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 freespace radiation, it becomes possible for a set of atomic spin waves to efficiently couple to these guided modes, while retaining a suppressed coupling to freespace 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, . This scaling represents an exponential improvement over the best previously established error bound of 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 atomlight 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 multiexcitation 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 inputoutput 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 wellsuited to cases where a limited number of modes are assumed to be relevant (such as a highQ 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 is the fundamental solution of the electromagnetic wave equation, and obeys Novotny and Hecht (2006):
(1) 
where is the position and possibly frequencydependent relative permittivity of the medium. The Green’s function physically describes the field at point due to a normalized, oscillating dipole at . is a tensor quantity (), as and 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); AsenjoGarcia et al. (2017).
In the frequency domain, the analogous classical problem that one would like to solve is to find the total field at point , given a known input field and a collection of polarizable dipoles located at , which are excited by the fields and rescatter light themselves. The values of are not known a priori, since they depend on the polarizability and the total field at [solving for 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 rescattered by the dipoles, we find , where 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 by the dipole moment operator , and by the field operator , where the superscript refers to the positivefrequency 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 . Thus, one is able to approximate by , which allows the Fourier transform of the equation to become local in time. Then, one arrives to the generalized inputoutput equation in the time domain, which reads Dung et al. (2002); Caneva et al. (2015); Xu and Fan (2015)
(2) 
To obtain the above expression, we have made use of the fact that , where is the atomic coherence operator between the ground and excited states of atom , and 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 by is wellfounded. 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 ( meter). Second, the electromagnetic environment itself should not have very narrowbandwidth 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, ) themselves. We do so by writing down HeisenbergLangevin equations for the atomic internal degrees of freedom, starting from the full atomfield Hamiltonian. Intuitively, the atomic spin will be driven by the quantum field at position . However, as the field itself only depends on other atoms via the inputoutput equation, the atomic dynamics can be fully derived from an equivalent master equation of the form Meystre and Sargent (2007), where is the atomic density matrix, and the Hamiltonian and Lindblad operators read
(3a)  
(3b) 
In the above expressions, the rates for coherent and dissipative interactions between atoms and are respectively given by
(4a)  
(4b) 
where the sign of 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); AsenjoGarcia et al. (2017). In the above Hamiltonian, we have neglected Casimir interactions between groundstate atoms (of the form ), as their spatial decay is very fast ( in free space, being the interatomic 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 nonHermitian Hamiltonian that reads , with
(5) 
along with stochastically applied “quantum jump” operators to account for the population recycling term () of Eq. (3b). While just describes the interaction of atoms through emission and reabsorption 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 threedimensional 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 singleexcitation guided Bloch modes with zero decay rate to emerge, which are decoupled from freespace radiation due to wave vector mismatch. We then analyze a single excitation in a finite lattice, and show how the guided modes acquire a nonzero 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 singleexcitation subradiance (where the atoms can equivalently be treated as classical dipoles) to the rich physics of the multiexcitation case. In particular, in one dimension, we show that multiexcitation 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.(25) the freespace Green’s tensor , with . Here is the solution to Eq.(1) when setting , and can be written as:
(6) 
where and is the wave number corresponding to the atomic transition energy. For a single atom, evaluating Eq. (4b) simply reproduces the wellknown vacuum emission rate , where . The singleatom energy shift in Eq.(4a) arising from formally yields a divergence and will be set to zero in what follows, as it should be incorporated into a renormalized resonance frequency . In Sec. IV, the Green’s function of a nanofiber will be decomposed into a freespace and a scattered component, , where does produce a finite, observable contribution to .
For concreteness, we will restrict ourselves to the following lattice geometries: onedimensional (1D) linear chains and closed circular rings, twodimensional (2D) square and threedimensional (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 and , 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 welldefined quasimomentum k, which can always be chosen to be within the first Brillouin zone. For such states, whose creation operators can be written as , 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 is given by .
As the Bloch modes are eigenstates of the effective Hamiltonian, they satisfy . Here and are real quantities and can be identified as the frequency shift of mode k (relative to the bare atomic frequency ) and the decay rate, respectively. One can readily show that, in terms of the singleatom spontaneous emission rate , they are given by:
(7a)  
(7b) 
where is the discrete Fourier transform of the freespace Green’s tensor.
We will now show that, when the atoms are placed at close enough distances, dipoledipole interactions can dramatically modify the decay rates of collective states. As the simplest case, let us consider an infinite onedimensional chain of atoms first, oriented along the direction. In that case, the wave vector constitutes an index for the modes, and one can consider the dispersion relation of frequency versus . 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, is only uniquely defined within the first Brillouin zone () 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 corresponding to light propagating in free space along the 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 spinwave excitation, which is given by Eq.(2) under the replacement (it is sufficient to consider the limit of classical dipoles for this argument). One can always expand the field in terms of plane wave components, . The state is clearly of Bloch’s form, and thus, only a discrete set of wave vectors ( 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 . Thus, one can readily verify that a spin wave outside the light line () has an associated electromagnetic field composed of axial wave vectors . This in turn implies that 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 wellknown 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 ) around the resonance frequency (), 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 at a frequency greater than the atomic resonance. This condition can be rewritten as a condition on the lattice constant required to support guided modes.
Equipped with this general intuition, we now quantitatively investigate the dispersion relation for the 1D infinite chain of twolevel atoms with polarization parallel or transverse to the array. The collective frequency shifts are derived in greater detail in Appendix A, and read:
(8a)  
(8b) 
where is the PolyLogarithm of order . These expressions are plotted in Fig. 1(b), for the particular value of . Here, the light line is indicated as before by a dashed line, but since , 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 and the polarization direction, and for instance, the effective mass at the zone edge () is negative (positive) for parallel (transverse) polarization. Exactly at the light line the expression for () becomes nonanalytic and diverges (has a derivative that diverges).
The collective decay rates can also be analytically derived (see Appendix A):
(9a)  
(9b) 
These summations run over reciprocal lattice vectors that satisfy . That is, only the diffracted waves enclosed within the light line will contribute to the decay rate. When , there are no values of 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  plane [Fig. 2(a)], the corresponding first Brillouin zone for Bloch wave vectors extends over the region , as shown in Fig. 2(b). The set of electromagnetic fields propagating in the plane at the atomic frequency have a wave vector of magnitude which defines a circle centered around the origin in kspace, as illustrated in Fig. 2(b). Similar to 1D, a sufficient condition for spinwave 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 . Thus, in 2D, guided modes exist as long as , which translates into a maximum allowed lattice constant .
Analogous to the 1D case, we can obtain closed mathematical expressions for the decay rates in the 2D lattice. They are given by:
(10a)  
(10b) 
from which we recover again the important result that Bloch states with 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 and directions, for the particular value of , which defines the light line as the circle , 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 groundstate 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 quasimomentum . Physically, the origin of this effect can be understood as follows. On one hand, and as previously discussed, for the field created by such a state is generally evanescent at all diffraction orders except for the component [c.f. Eq.(10b)], which corresponds to a plane wave propagating perpendicularly to the atomic plane. On the other hand, the state 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 , there will be other components that are not evanescent, yielding a radiative state. We note that, although only the state has exactly zero decay rate, other modes around this point will also show a strong suppression in the emission rate relative to .
iii.2 Finite Lattice (Single excitation)
In this section, we analyze the decay rates and spatial properties of singleexcitation eigenstates, for a lattice of finite size. We show that all eigenstates now acquire a nonzero 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 , 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 nonHermitian, in general the eigenstates will not be orthogonal in the standard quantum mechanical sense (i.e., two eigenstates and will not satisfy ) AsenjoGarcia 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 nonevolving 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 atoms with lattice constant , for which numerical diagonalization of in the oneexcitation manifold produces eigenstates (denoted by , ) and complex eigenvalues. As in the infinite case, the eigenvalues can be written in the form , with and representing the frequency shift relative to and decay rate, respectively. For concreteness, the eigenstates will be ordered in increasing decay rate, such that represents the most subradiant state and the most radiant one. To start understanding the properties of this system we fix the atomic number and change the lattice constant . For each value of , we diagonalize , and obtain the different values for the decay rates and frequency shifts associated with each of the collective modes. Figure 3(a) shows the resulting singleexcitation decay rates for each collective mode, normalized by the freespace singleatom emission rate , in the case of atomic polarization parallel to the chain. In this plot, a vertical cut at a fixed value of contains the different values of . 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 a qualitatively distinct behavior. In this region, the decay rates of some of the modes are dramatically suppressed (), 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 generated by the most subradiant state when the atomic polarization is parallel to the chain (we choose a distance offset from the plane containing the atomic chain in order to avoid seeing the divergent nearfields 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 . As the input field is vacuum, the intensity only involves calculating twobody correlations of the eigenstate.
While the wave vector 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 with each of the modes in the finite case. Specifically, the discrete Fourier transform of the coefficients that define each mode is peaked around a different value , which can be used to label the state. In particular, let us consider a general singleexcitation state, which can be written as , where is defined as the state where atom is excited while all others are in their ground states. Then, we define the discrete Fourier transform of the associated coefficients , for discrete values of (). For each value of , the function shows a well defined peak at a distinct value of . In Figs. 1(b),(c) (circles) we plot the decay rates and energy shifts of each mode as indexed by the dominant wave vector, for 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 .
The exact behavior of depends on the microscopic details, such as the polarization of the atoms. For instance, for twolevel atoms, the smallest decay rate decreases monotonically as , while for transverse polarization it oscillates. Regardless of these details, however, the scaling with of the few lowest decay rates seems to show a universal behavior, going like . In Fig. 3(c), we show the scaling for the three lowest eigenstates as a function of , while in Fig. 3(d) we show the scaling for fixed . The scaling with is satisfied for all . For transverse polarization, there is a particular value of lattice constant (that tends to as the atom number increases), for which the decay rates do not follow exactly the scaling with . We believe that this is related to the fact that for transverse polarization and the band structure becomes flat at the edge of the Brillouin zone. Nevertheless, and as we discuss in Appendix B, the scaling seems to appear rather generically for finitesize, onedimensional photonic crystal structures.
Ansatz for singleexcitation collective modes.– If the chain is finite the singleexcitation collective modes are not spin waves with pure wave vector , 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 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 ) 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 , the wavefunction coefficients of the exact collective modes are close to the orthonormal set of functions defined by wave vector :
(11) 
Here , and the atomic positions () and . Figure 4(a)(c) shows the exact coefficients for the most subradiant state (), a state with dominant wave vector near the light line, and the most radiant state (), for 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, . In Fig. 4(d) this is plotted as a function of the wavenumber 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 as .
Moreover, this ansatz not only approaches the spatial pattern of the wave function, but its decay rate defined as captures the same scaling with the index and : . The overall proportionality constant varies depending on the microscopic details (such as the polarization or the value of ). For instance, for , (parallel polarization) and (transverse polarization). The comparison between and is shown in Fig. 3(c)(d) (solid circles correspond to the ansatz).
2D Square Array.
The previous results are not specific to the onedimensional chain and can be extended to other lattice geometries. As an example, let us consider a finite square array of atoms spanning the  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 collective modes. We can also define an ansatz wave function with coefficients , where are the coefficients for the onedimensional ansatz Eq. (11). We can then associate with each of the collective modes a pair of values , 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 after diagonalizing the Hamiltonian for an array of atoms, for different values of . 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 . As it can be seen in this figure, the most subradiant modes correspond to those at the edges of the Brillouin zone, i.e., . The wavefunction 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 as before, and also a set of subradiant modes around for , 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 is small enough, the decay rate will scale as . In particular, we find that for (most subradiant state), for and for , while for , for . For ranges of not included above, the decay rates are not suppressed with increasing system size, since in that case the wave vector 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 threedimensional 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 finitesize system, it does not appear that subradiant states can be identified based on the infinitesystem 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 finitesize lattice. In Fig. 6, we plot the decay rates for the eigenstates associated with a lattice of twolevel 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 , the effect of subradiance can survive even for lattice constants . 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 (see sketch in Fig. 7) and with linear polarization transverse to the plane of the ring.
In Fig. 7, for a distance of , we numerically diagonalize the singleexcitation block Hamiltonian, and plot the decay rate of the most subradiant state versus atom number . It can be seen that an exponential suppression emerges, . For the chosen parameters, the minimum decay rate for a ring drops below that of an open chain for 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 has been divided into three regions: in the left and right regions the lattice constant is uniform and equal to , while in the middle it changes slowly (from to ) following the red line and creating a defect. The lattice constant in the middle is chosen to follow a sinusoidal variation, , where () denotes the atom position. In the same figure, the band structures for an infinite lattice with constant and are plotted for the case of atomic polarization parallel to the chain. It can be seen that the smaller lattice constant supports propagating modes over a range of frequencies (pink shaded region) that lies within the band gap of lattice . 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 , 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 atoms. For a smooth variation in the lattice constant (here occurring over a region of size ), 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 . The numerically calculated decay rate of this mode is plotted in Fig. 8(c) as a function of , and clearly confirms the expected behavior.
iii.3 Multiexcitation modes
We now turn to the problem of multiple excitations stored in an atomic 1D lattice. As in the previous subsection, 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., ) 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 excitations in a given mode would simply have a decay rate times that of a single excitation. The fact that we are dealing with spins, where a single spin cannot be excited twice (), leads to highly nontrivial properties of multiply excited states. In this section, we will analyze in detail the spatial properties and scaling of decay rates of multiexcitation subradiant states. While we will not explicitly utilize these states in later sections, these findings might help to provide some initial insight into how manybody physics can be encoded into subradiant manifolds.
Let us first consider the twoexcitation manifold. A general state within this manifold can be written as where now corresponds to the state whith atoms and excited while the rest remain in the ground state. Although it is necessary only to specify for to define the wave function, in the following plots and for visual appeal we also assign values to for , by simply defining . To illustrate that the spin system behaves differently than a bosonic system, we begin by considering the twoexcitation state formed by occupying the same singleexcitation mode twice. In particular, we construct the twoexcitation state given by (properly normalized). Here, is the collective operator that creates the most subradiant single excitation in a chain of atoms, when applied to the ground state. In Fig. 9(a) we plot the corresponding probability density for the case . The twoexcitation wave function appears relatively smooth, except for a sharp cut 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 is only suppressed with the length of the chain as , in stark contrast to the singleexcitation case. This scaling is shown in Fig. 9(d) (orange circles).
Numerically, we now exactly diagonalize the Hamiltonian in the twoexcitation manifold, and identify the most subradiant state. The scaling of its decay rate with is plotted in Fig. 9(d) (blue open circles), and is seen to preserve the scaling present in the singleexcitation manifold. The probability density for the case of 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 for the wave function of the form (properly normalized). Here denote the coefficients of the singleexcitation orthonormal ansatz Eq.(11) associated with the wave vector . Such an ansatz naturally constructs a state that incorporates “Pauli exclusion”, and a smooth separation of excitations in space. Taking and , i.e., building a twoexcitation state from the two most subradiant singleexcitation states, yields a wave function 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 , as it is shown in Fig. 9(d) (blue solid circles).
We can then associate with each of the exact twoexcitation collective states a pair of quasimomentum values for which the wavefunction 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 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 singleexcitation modes used to construct the ansatz, i.e., , is not far from the exact value. This is quantified by the relative error and it is plotted in Fig. 9(f). For completeness, we also show in Fig. 9(g) the error in overlap between each twoexcitation eigenstate and the bestmatched ansatz state, . This error is very small in the subradiant region.
In the more general case of excitations, the most subradiant mode and its decay rate 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 twoexcitation manifolds, i.e., , 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 excitations in the most subradiant mode, i.e., , which scales as in the twoexcitation case, .
One can also numerically evaluate the error in overlap between the most subradiant state and an ansatz state, . Here the wavefunction amplitudes are generalized from the twoexcitation case, and constructed from the Slater determinant of singleexcitation wavefunction ansatz coefficients . For the most subradiant mode these correspond to the most subradiant singleexcitation modes and in the large atom number limit the error is found to scale like