Optimization of photon storage fidelity in ordered atomic arrays

A major application for atomic ensembles consists of a quantum memory for light, in which an optical state can be reversibly converted to a collective atomic excitation on demand. There exists a well-known fundamental bound on the storage error, when the ensemble is describable by a continuous medium governed by the Maxwell-Bloch equations. The validity of this model can break down, however, in systems such as dense, ordered atomic arrays, where strong interference in emission can give rise to phenomena such as subradiance and"selective"radiance. Here, we develop a general formalism that finds the maximum storage efficiency for a collection of atoms with discrete, known positions, and a given spatial mode in which an optical field is sent. As an example, we apply this technique to study a finite two-dimensional square array of atoms. We show that such a system enables a storage error that scales with atom number $N_\mathrm{a}$ like $\sim (\log N_\mathrm{a})^2/N_\mathrm{a}^2$, and that, remarkably, an array of just $4 \times 4$ atoms in principle allows for an efficiency comparable to a disordered ensemble with optical depth of around 600.

Atomic ensembles constitute an important platform for quantum light-matter interfaces [1], enabling applications from quantum memories [2,3,4,5] and few-photon nonlinear optics [6,7,8,9,10,11] to metrology [12,13,14,15]. In typical experiments, ensembles consist of disordered atomic clouds, with the propagation of light through them modeled phenomenologically by the Maxwell-Bloch equations [16,17]. Within this description, the atoms are treated as a smooth density and the discreteness of atomic positions is ignored. In addition, spatial interference that can arise from light arXiv:1710.06312v2 [quant-ph] 4 Sep 2018 scattering is neglected, and the emission into directions other than the mode of interest is treated as an independent atomic process. Within this formalism, one can derive standard limits of fidelity for applications of interest -for example, the storage error of a quantum memory scales inversely with the optical depth (D) of the ensemble [18].
Recently, novel experimental platforms have emerged where it is possible to produce small ordered arrays of atoms [19,20,21,22,23]. Intuitively, one expects that strong interference in light emission can emerge, which renders inoperable the typical theoretical approaches to modeling light-atom interfaces. Theoretically there has been growing interest in novel quantum optical effects in arrays, such as subradiance [24,25,26,27,28,29,30], topological effects [31,32], and complete reflection of light [33,34,35]. Indeed, it has already been shown numerically that an ordered onedimensional array of atoms coupled to a nanofiber allows for a storage error exponentially smaller than the previously known bound [29]. In this work, the exponential scaling was observed by considering a fixed, spatial waveform for the optical pulse. However, two interesting questions that arise are (i) whether it is possible to develop a theoretical technique to bound the error, which takes fully into account the atomic positions and the interference of emission in all directions, and (ii) whether an improved scaling is possible for atoms in free space, as opposed to coupled to a photonic structure. These questions are affirmatively answered in our work.
In particular, we provide a construction that enables the maximum storage efficiency to be found, given the atomic positions and the desired spatial mode of light. This procedure is based upon solving the dynamics of a "spin model", which encodes the multiple scattering and interference of light as it interacts with atoms, and then calculating the light emitted into the desired mode by an input-output equation. We show that the maximum efficiency is given by the maximum eigenvalue of a Hermitian matrix, whose elements are derived from the atomic positions and optical mode. While this technique is completely general, we apply it specifically to the case of a twodimensional square array of atoms. In particular, it has recently been shown that an infinite array can in principle form a 100% reflector for light [33,34,35], when the lattice constant d is smaller than the resonant wavelength λ 0 . While a mirror constitutes a "passive" optical system, it is natural to ask whether this implies a 100% success probability, if the system were functionalized into a quantum memory. For a finite array, we show that the minimum error decreases like ∼ (log N a ) 2 /N 2 a for storage from a Gaussian-like mode, and remarkably, that a 4 × 4 array in principle already enables an error below 1%.

The spin model
The full dynamics of light emission and re-scattering of an arbitrary collection of atoms in free space, specified only by their discrete, fixed positions r j , can be related to an effective model containing only the atomic degrees of freedom and the incident field [36,37,38,39,40,41]. We first review this formalism for two-level atoms with ground state |g and excited state |e , with the dipolar transition |g − |e coupled with free space optical modes. Within the Born-Markov approximation, these modes can be integrated out to yield effective dynamics for the atomic density matrixρ, which evolves asρ = −(i/ )[H,ρ] + L[ρ], where the Hamiltonian and Lindblad operators read [36,37,38,39,40,41,42] Here H in is associated with the input field that drives the atoms (which need not be specified for our purposes), d eg andd j are the dipole matrix element and unit atomic polarization vector associated with the transition, and σ βγ = |β γ| are atomic operators with {β, γ} ∈ {e, g}. G 0 (r j , r l , ω eg ) is the electromagnetic Green's function tensor in free space, which is the fundamental solution of the wave equation and fulfills where the curl is taken with respect to r. The Green's function explicitly takes the form [43] G 0 (r j , r l , ω eg ) = where R = |r j − r l | and k 0 = ω eg /c is the wavevector associated with the atomic transition frequency ω eg , with c being the speed of light. We note that the local term [i.e., G 0 (r j , r j , ω eg )] is divergent. This term is responsible for the Lamb shift and is incorporated into a renormalized resonance frequency ω eg . Physically, Eq. (1a) describes the coherent exchange of atomic excitations mediated by photons. On the other hand, Eq. (1b) describes the collective emission or dissipation of excited atoms, after integrating out the common reservoir of electromagnetic modes with which they interact (within the Born-Markov approximation). Instead of solving the density matrix evolution as governed by the master equation, one can equivalently work within the stochastic wave function or "quantum jump" formalism [44]. In that case, the system is described by a wave function, which deterministically evolves under an effective, non-Hermitian Hamiltonian This Hamiltonian captures both the coherent evolution of Eq. (1a) and the last two terms of the Lindblad operator in Eq. (1b). In addition, one must also stochastically apply quantum jump operators to the wave function, to capture the population recycling terms σ l ge ρσ j eg of Eq. (1b). Formally, the jump operators of our system will consist of superpositions of σ l ge , i.e. atomic lowering operators, which physically encode the emission of a photon. In the following, we will be interested in initial states with just a single excitation; thus, any jump operator trivially takes the system to the ground state |g ⊗N , where it cannot further evolve or contribute to observables of interest (e.g., the emission of a photon). Furthermore, the rate that jumps occur is exactly equal to the rate of population loss of the wave function evolving under H eff . Thus, in our case, jumps are effectively accounted for just by evolution under H eff alone. Any loss of population from the single-excitation manifold implies that a corresponding population is building up in the manifold |g ⊗N |1(r, t) , where all the atoms are in the ground state and a single photon is emitted in some spatial-temporal pattern. We next discuss how the photon-emission pattern and its overlap with a mode of interest can be calculated.
Given the evolution of the atomic state under H eff , any observables associated with the total field operatorÊ out (r) can be derived from the input-output relation [37,38,39,40,41] Formally, this equation states that the total field is a superposition of the incoming field and the fields emitted by the atoms, whose spatial pattern is contained in the Green's function. Equation (5) enables the field to be calculated at any point r, based upon the evaluation of an atomic correlation function ∼ G 0 (r, r j , ω eg ) ·d j σ ge j weighted by the Green's function. Evaluating the Green's function at each r and the corresponding atomic correlation function to construct the field everywhere can become tedious. However, in experiments one often cares about the projection of the field into a specific spatial mode, such as a Gaussian (see Fig. 1). It can be proven (see Appendix A) that this projection depends only on the amplitudes of the mode of the classical field E det (r) at the positions of the dipoles. We can thus define the quantum operator associated with the detector aŝ whereÊ det,in is the input field in the detection mode and F det = z=const d 2 r E * det (r) · E det (r) is a normalization factor. Here, the normalization is such that Ê † detÊ det represents the photon number per unit time emitted into the mode.
Before discussing the specifics of the retrieval efficiency, we would like to briefly discuss the validity of the Born-Markov approximation, which allows one to trace out the photonic degrees of freedom and arrive at an atomic master equation, as well as to write equations for the field operators that depend instantaneously on the atomic operators. This approximation is valid whenever (1) the photon bath correlations decay much faster than the atomic correlations and (2) retardation can be ignored. The Figure 1. Schematic of a quantum memory using a two-dimensional atomic array. An excitation initially stored in the |s -manifold is retrieved as a photon by turning on the classical control field Ω c (blue arrows), which then creates a Raman scattered photon from the |g − |e transition. The photon is detected in some given mode, illustrated here as a Gaussian beam.
first condition is obviously satisfied for atoms in free space, as the vacuum's Green's function has a frequency spectrum that is much broader than the atomic linewidth. Neglecting retardation in both the photon-mediated interactions between atoms and the field produced by the atoms requires the characteristic length L of the atomic system to be much smaller than that of a spontaneously-emitted photon, which is ∼ c/Γ 0 ≤ 1 m [45,46,47,48], where Γ 0 = µ 0 ω 3 eg d 2 eg /3π c is the single-atom spontaneous emission rate in vacuum. It should also be pointed out that for at most a single atomic excitation, the dynamics of atom-light interactions can readily be solved in an exact manner [46,48,49,50,51]. In this regime of linear optics, the dynamics can be analyzed for each frequency component in the Fourier domain, exploiting the fact that different frequency components do not couple to one another. However, the spin model presented above has a natural extension to the multi-excitation case (e.g., studying the storage of multiple photons and their subsequent nonlinear interaction [52,53,54]), whereas exact solutions are only available in a limited number of cases [55,47,56].

The retrieval efficiency
The typical quantum memory scheme consists of an ensemble of three-level atoms where an additional metastable state |s is coupled to the excited state |e by a classical control field with Rabi frequency Ω c (r, t) and detuning ∆ from the transition frequency ω se (see Fig. 1) [18]. While the state |s is typically associated with another state in the groundstate hyperfine manifold, in our case this would deleteriously reduce interference effects in emission. For example, in storage where all atoms begin in |g , there is no interference pathway to suppress spontaneous emission into |s once an incident photon excites an atom to |e . Thus, we assume that our atoms have no hyperfine structure and there is a unique ground state, as would be the case for bosonic Sr or Yb atoms, and that level |s is a long-lived, higher-lying excited state. Dipole-dipole interactions on the |e -|s transition have no effect, as they require at least two total excitations in the system. In the main text, we will furthermore take the conceptually simpler case where |e is the unique excited state coupled to |g (for concreteness, with polarizationd j =x). A more realistic model with three excited states |e x,y,z , providing an isotropic atomic response to light, is presented in Appendix C, but the results qualitatively remain the same.
Instead of storage, it is mathematically more convenient to optimize the retrieval problem, in which an initial collective spin excitation |ψ(t = 0) = j s j (t = 0)σ sg j |g ⊗Na is emitted as an outgoing photon on the |g − |e transition via a Raman process facilitated by the control field Ω c . The initial state then evolves under the total Hamiltonian H = H eff + H c , where the Hamiltonian associated with the control field is H c = j − ∆σ ee j + Ω j c (t)(σ es j + H.c.) and H in = 0 as there is no external field driving the |g − |e transition in retrieval. We take a spatially uniform, but , although it is straightforward to generalize the following discussion to the case of a spatially varying control field. Then, for a given detection mode and atomic spatial configuration, we want to find the initial spin amplitude s j (0) that maximizes the retrieval efficiency. By time-reversal symmetry, the storage efficiency for an incoming photon in the same mode and for the same atomic configuration is identical, when optimized over the temporal shapes of the incoming photon and control field [18]. Writing the general state in time as where the matrix M jl = 3πk −1 0d * j · G 0 (r j , r l , ω eg ) ·d l . While we explicitly consider the model above, we note that it is straightforward to add a number of other effects (e.g., decay of the |s state or dephasing) into the analysis.
From Eq. (6), we can evaluate the expected total photon number η = ∞ 0 dt Ê † det (t)Ê det (t) emitted into the detection mode. Assuming that the control field is turned on for long enough, it is guaranteed that one photon in total is emitted into all modes, and thus η also represents the retrieval efficiency. Evaluating the atomic operators in Eq. (6), we find that where we have defined the local scalar field E j = E det (r j ) ·d * j at the atom positions, and S λ 0 = (3/2π)λ 2 0 is the resonant atomic optical cross-section (λ 0 = 2π/k 0 being the resonant wavelength).
Equation (9) can be simplified by noting that M jl in Eq. (7) is a symmetric complex matrix. Thus, if M jl is diagonalizable (as we numerically verify in our cases of interest), its eigenvalues λ ξ are complex and its eigenmodes v ξ are non-orthogonal in the quantum mechanical sense, but obey the orthogonality and completeness conditions v T ξ · v ξ = δ ξξ and ξ v ξ ⊗ v T ξ = I [37]. Projecting the equations of motion into this basis results in N a decoupled pairs of equations: where e ξ = j v ξ,j e j , s ξ = j v ξ,j s j . Provided that the atomic excitation has left the system as t → ∞, one can derive that Inserting this equality into Eq. (9), we readily find where and E ξ = m v ξ,m E m . Importantly, K is an N a × N a Hermitian matrix which depends only on the positions of the atoms and the detection mode, but not on the specific time dependence of the control field (for example, one could apply a π pulse that transfers all of the excitation from state |s to |e at time t = 0). The maximum retrieval efficiency is thus given by the initial configuration corresponding to the eigenvector of K with the largest eigenvalue. We should note that while the efficiency η of retrieval is independent of the particular profile Ω c (t), the shape of the outgoing photon is completely determined by the control field. By time-reversal symmetry, if one wants to store an incoming photon with maximum efficiency, one must first consider its time-reversed shape (i.e., an outgoing photon), find the unique control field Ω c (t) that generates such a shape in retrieval, and then apply the time-reversed fieldΩ c (t) for storage. Before proceeding further, we briefly comment on the classical and quantum optical aspects of the calculation presented above. An equation analogous to Eq. (9) also applies if the atoms were replaced by classical oscillating dipoles with amplitudes e j (t). Such an equation corresponds to the projection of the total classical radiated field into a particular spatial mode. The equivalence between classical and quantum equations is not surprising, given that both the propagation of classical and quantum fields are given by Maxwell's equations. In our particular problem of interest, the quantum nature of the field manifests itself by considering field correlations. For example, using Eq. (6), one can calculate the second-order correlation function Ê †2 detÊ 2 det . As the atomic state that we consider contains at most one excitation, this correlation function is exactly zero, or perfectly "anti-bunched," reflecting the fact that only a single photon is emitted.

2D square array
While the formalism presented above is general to any ensemble of atoms with known positions, we now apply it to a 2D square array with lattice constant d. This case is particularly interesting, as an infinite array of two-level atoms can act as a perfect mirror for incoming light at normal incidence when d is smaller than the atomic resonant wavelength λ 0 [33,34,35]. Physically, the incoming field guarantees that all the induced atomic dipoles oscillate with the same phase. While such a configuration can in principle emit into various diffraction orders, for d < λ 0 , all of the orders except the one perpendicular to the plane become evanescent, and cannot radiate away energy. With only two channels of emission possible (forward and backward), the scattered field of the array perfectly interferes with an incident resonant photon in the forward direction, leading to complete reflection of light. Likewise, when an excitation is stored uniformly in the infinite array with d < λ 0 , it is "selectively radiant" [29], as interference guarantees that the retrieved photon is perfectly emitted into two plane waves normal to the array (we assume that this symmetric emission can be re-combined). While this simple argument hints that a finite array can also be very efficient, what remains is to quantify the error. We thus analyze the retrieval efficiency of an array made of N a = N × N atoms.
As far as the detection mode is concerned, a common mode to project into is a Gaussian beam. There is a technicality, however, since a Gaussian beam is only an approximate (paraxial) solution to Maxwell's equations. While such an approximation usually suffices, here we anticipate that one can achieve nearly perfect storage and retrieval efficiencies. Consequently, it is not obvious a priori that the small (actual) retrieval errors are not overwhelmed by the error of the paraxial approximation itself. Thus, we consider an exact mode solution for Maxwell's equations (see Appendix B for details), which approaches the Gaussian solution in the limit of large beam waist w 0 .
Before presenting the numerics, one can already intuitively argue the fundamental sources of error associated with a finite array by considering the reflectance problem. If the beam waist w 0 is too large with respect to the array dimensions, then part of the incoming light will not see the atoms and will be transmitted or scattered in other directions by the edges of the array. If w 0 is too small, the incoming mode contains a broad range of wavevectors with different propagation directions. Since different angles have maximum reflectance at different detunings relative to the bare transition frequency ω eg [35], the overall reflectance for a near-monochromatic photon will be reduced. For a given array, an optimal beam waist thus maximizes the reflectance of an incoming photon (at optimal detuning). The situation is analogous for the retrieval problem, where the optimization over the photon frequency is replaced by an optimization over the initial spatial distribution of the collective s-excitation.
To check this behavior, we numerically calculate the minimum retrieval error = 1 − η varying the beam waist w 0 , for several different atom numbers. In Fig. 2(a), the error is plotted as a function of the ratio between the array area S arr = d 2 N a and  w 2 0 . Here, we have taken the retrieval mode to consist of a symmetric superposition of Gaussian beams emitted in opposite directions from the array, with the view that these beams can in principle be recombined. For concreteness, we consider a lattice constant of d = 0.6λ 0 , although other choices d < λ 0 do not affect the general scalings. As S arr /w 2 0 grows, the error initially scales as ∼ 1 − Erf 2 (N d/ √ 2w 0 ) (illustrated by the dashed curve), where Erf(z) is the error function. Physically, this error corresponds to the fraction of the energy carried by the Gaussian beam beyond the array boundaries. In Fig. 2(b) we plot (in log-log scale) as a function of the ratio between w 0 and λ 0 (for values larger than one), again for different array sizes. Up to a point where the beam waist becomes comparable with the array dimension, the error scales roughly as ∼ (λ 0 /w 0 ) 4 (dashed line). This error physically arises from the range of wavevector components that make up the detection mode, which is inversely proportional to w 0 . An analysis of the reflectance of a beam of finite waist from an infinite array in fact shows the same scaling, when considering the fraction of light that is not reflected. Overall we have that the minimum error can be approximated by the expression The constant C can be obtained by fitting the error: for d = 0.6λ 0 we find C ≈ 2.4·10 −3 . One can use Eq. (15) to find the optimal beam waist. After optimizing w 0 we find that the leading term for the error is given by In Fig. 2(c) this approximate expression for the minimum retrieval error is compared with the value obtained by numerical optimization. The associated optimum beam waist for the retrieval mode is also plotted for completeness. Interestingly, even a 4 × 4 array of atoms can in principle already enable a storage/retrieval efficiency of above 99%. In comparison, an optical depth of nearly D ∼ 600 is needed to obtain the same error in a conventional ensemble [18]. In the case where the beam waist does not significantly diverge over the length of the ensemble, the optical depth is given by D ∼ S λ 0 N a /w 2 0 . For cold atoms, an atom number on the order of N a ∼ 10 6 -10 7 might be required to achieve a value of D ∼ 600.

Analysis of disorder
In this section, we analyze the effects of various types of disorder in the array. One useful attribute of our efficiency calculation is that it enables different spatial configurations to be studied. Thus, we can easily include imperfections such as the absence of atoms ("holes") in the array, or classical position disorder. We first examine the case of some number N def of holes in the array. Intuitively, one expects that the relative decrease in efficiency, (η − η def )/η, will be proportional to the ratio between the intensity of the detection mode hitting the empty sites, to the total intensity over the array. Here, η def and η denote the maximum retrieval efficiency with and without the holes, respectively, with the beam waist w 0 chosen to optimize η. In Fig. 3(a) we plot the relative loss as a function of j∈def |E j | 2 / l |E l | 2 , where the sums of the field intensities in the numerator and denominator run over sites of holes and all sites, respectively, sampling over 100 random configurations for different densities of holes (N def /N a up to 20%). One sees a clear statistical relation of the form The constant of proportionality α in Eq. (17)   configuration, which would be applicable if an experiment could resolve the positions of the holes in a single shot [21], we expect a similar scaling even if the positions of holes are unknown.
Classical disorder for the atomic positions consists in having the atoms displaced by random amounts δ j = (δ x,j , δ y,j ) from their position in the perfect lattice. It is shown in Ref. [35] for the case of reflectance of an infinite array that, when the δ's are extracted from a Gaussian distribution with standard deviation σ, then the decrease in reflectance introduced by the disorder scales as σ 2 /d 2 . We find numerically the same result for the retrieval error of the finite array. In particular, in Fig. 3(b) the error introduced by disorder is plotted as a function of σ for different array dimensions and fixed lattice constant. This error is defined as the difference between the optimized maximum retrieval efficiency η of a perfect lattice, and the mean retrieval efficiency η dis (sampled over many configurations) with the same initial atomic wave function and beam waist but with disorder in the atomic positions.

Finite detection time
When calculating the retrieval efficiency, given by Eq. (9), we have implicitly assumed that the detection time is infinite, such that all the energy emitted into the detection mode is collected. Practically, it might also be relevant to consider the retrieval efficiency given a finite time window 0 < t < T d for photon collection, such as if an experiment has other limiting time scales (i.e., atom trapping time, required fast readout, etc.).
The efficiency detected for an arbitrary detection time window T d is given by where e j (t) is obtained by integrating Eqs. (7)- (8). In general the temporal profile of the emitted field depends on the control field amplitude Ω c (t) and detuning ∆. If one wants to achieve a high efficiency in the shortest time, then the optimal strategy is to essentially use the control field to apply an instantaneous π-pulse at t = 0, thus instantly transferring the excitation stored in the metastable state |s to the rapidly emitting excited state |e . In an array, this collective excitation in |e will emit a photon at a rate ∼ Γ 0 comparable to the single-atom emission rate, ensuring that the errors due to finite time window T d become very small once T d is on the order of a few ∼ Γ −1 0 . In Fig. 4, we plot the relative error 1 − η T d /η due to the finite detection time, where η is the detection efficiency for an infinite time window, for an array of 10 × 10 atoms with d = 0.6λ 0 and optimal beam waist. We notice that for a detection time T d ∼ 10/Γ 0 the error is of the order of 10 −3 . The possibility of having a good retrieval efficiency even for a short detection time is a consequence of the fact that, while the array can support highly subradiant states [25,29,30,34,35], they form a negligible component of the optimized spin wave for storage and retrieval. This makes intuitive sense, as to interface with light efficiently, one should use radiant or "selectively radiant" atomic excitations rather than states that decouple from light.

Conclusions
In summary, we have introduced a prescription to calculate the maximum storage and retrieval efficiency of a quantum memory, which fully accounts for re-scattering and interference of light emission in all directions. Our approach is in principle applicable to any system where the positions of the emitters are known (or can be reasonably modelled, such as assigning random positions) and the spatial and spectral response of the dielectric environment (i.e., the Green's function) is also known [2,3,4,5,57,58,29,59,37,60,61,62,63]. As one particular application, we have shown an improved scaling of errors for atoms in free space, compared to the result predicted by the one-dimensional Maxwell-Bloch equations. We speculate that it is possible to obtain an exponential reduction of errors versus atom number in free space, by using arrays that are not completely periodic. The question of how to tailor the spatial positions will be left to future work.
More broadly, we expect that a significantly improved storage efficiency is possible whenever the excited state emission is largely radiative and coherent, which includes not only atoms but solid-state emitters with large zero-phonon line and Fourier-limited linewidths [63]. Techniques to reversibly map between photonic and atomic excitations in arrays should find a variety of exciting applications. For example, it would allow for photonic quantum gates, given some form of spin interactions in the array (such as between Rydberg levels [64]), or would allow for exotic spin states (like subradiant [24,25,26,27,28,29] or topological excitations [31,32]) to be detected optically. It would also be interesting to investigate whether the spin state itself could be engineered to produce a useful non-classical state of outgoing light. More broadly, the ability to formally map atom-light interactions to a long-range open spin model could provide new insights into quantum optical phenomena with atomic systems.

Appendix A. Green's function expansion in plane and evanescent waves
Here we derive Eq. (6) of the main text by using an expansion of the Green's function in terms of plane and evanescent waves. The Green's function Eq. (3) can be written in the angular spectrum representation, i.e. as an integral over k x and k y in Fourier space, as [43] G ± 0 (r, r , ω eg ) = and the ± denoting the sign of z − z . We can separate the integral in Eq. (A.1) into two separate integrals: for values of k x , k y lying inside and outside the disk defined by k 2 x + k 2 y = k 2 0 . This decomposition separates the plane waves from evanescent waves, i.e., we can write G ± (r, r , ω eg ) = G ± pl (r, r , ω eg ) + G ± ev (r, r , ω eg ), where The integral in the plane waves part can be rewritten in polar coordinates using k 0 = k 0 (sin θ cos φ, sin θ sin φ, cos θ), obtaining G ± pl (r, r , ω eg ) = i 8π 2 k 0 2π 0 dφ π/2 0 dθ sin θ Q ± e ik 0 (sin θ cos φ(x−x )+sin θ sin φ(y−y )±cos θ(z−z )) . orthogonal to k 0 and between them, G + pl (r, r , ω eg ) can be expressed as where we have defined a plane wave basis u k 0 ,θ,φ,α (r) =ê α k 0 e −ik 0 (sin θ cos φx+sin θ sin φy+cos θz) , (A.9) with the normalization z=const d 2 r u * k 0 ,θ,φ,α (r) · u k 0 ,θ ,φ ,β (r) = Similarly one has An analogous expression can be found for the evanescent wave part. Here it is convenient to define the vectork 0 = k 0 (cosh ξ cos φ, cosh ξ sin φ, i sinh ξ): orthogonal tok 0 and between them, one can indeed write Similarly one has Now let's consider a detection mode that does not contain evanescent components for simplicity, so that it can be expanded just in terms of monochromatic plane waves as E det (r) = 1 (2π) 2 α 2π 0 dφ π 0 dθ sin θ c k 0 ,θ,φ,α u k 0 ,θ,φ,α (r).
(A. 18) The overlap between this mode and the field generated by a dipole is where we have used Eq. (5) (without input field) to express the field generated by the dipole through the Green's function and Eqs. (A.8) and (A.11) for the Green's function decomposition. Adding the input field and normalizing the detection mode we finally obtain Eq. (6) of the main text.

Appendix B. Gaussian detection mode
Here we present the detection mode which we have chosen to study the retrieval efficiency of the 2D array. We choose a solution oscillating with frequency e −iωegt , and where the x-component of the electric field in wavevector space is given by is the Heaviside step function. That is, E x has a Gaussian distribution for k 2 x + k 2 y ≤ k 2 0 while it is zero for k 2 x + k 2 y > k 2 0 , such that the field does not contain evanescent components. In the y direction, we take the field to be identically zero. The value of the z-component is then determined by Maxwell's equations [65]. The real space profile of this mode can be obtained by Fourier transformation: where (ρ, z) are the cylindrical coordinates for r, while J 0 and J 1 are Bessel functions. If evanescent components were included, the field in real space would identically consist of a Gaussian in the z = 0 focal plane with beam waist w 0 . The step function in wavevector space enforces in real space a diffraction limit, and distorts the beam to prevent a focal spot smaller than ∼ λ 0 . For large w 0 the mode tends to the paraxial solution, i.e. E z det vanishes and E x det assumes the form of a fundamental Laguerre-Gauss mode [43].

Appendix C. Spin model for isotropic atoms
In the main text we have introduced a formalism to calculate the retrieval efficiency of an atomic ensemble of three-level atoms, with an excitation initially stored in a metastable state |s coupled to the excited state |e by a classical control field. Instead of a single excited state, a more realistic minimal model of an atom consists of three excited states |e α , where α = x, y, z denotes the three possible orientations of the dipole transitiond. The effective Hamiltonian (4) generalizes to H eff = H in − µ 0 d 2 eg ω 2 eg j,l αβ G αβ (r j , r l , ω eg ) σ eg α,j σ ge β,l , (C.1) where the sum over α and β are over x, y, z. Here, σ ge β,l = |g l e β | l is the lowering operator on atom l, which takes the excited state |e β to the ground state |g . It should be noted that in general, transitions with different orientations can mix together (e.g., one atom could decay from |e y and excite another atom from the ground state to |e x ), as a photon emitted from a given dipole orientation does not have the same global polarization everywhere in space.
In the case in which the state |s is coupled only to one of the three excited states, for concreteness |e x , it is straightforward to generalize the main result of the paper.  Figure C1. (a) Comparison between the optimal retrieval error opt (left axis, blue lines) and the corresponding optimal beam waistw 0 (right axis, red lines) for the case of a single excited state discussed in the main text (continuous lines) and the case of three excited states (dashed lines), as functions of the linear array dimension N . (b) Relative difference ( opt,iso − opt,TL )/ opt,TL between the retrieval errors of the isotropic and two-level atomic structures plotted in (a).
Eq. (13) indeed keeps the same form, but with the matix K generalized to where E ξ = m v ξ,m E x det (r m ) and the sum over the index ξ of the eigenvectors has 3N a values. In Fig. C1(a) we compare the minimum retrieval error for an N × N square array of atoms versus N , for the cases of a single excited state and for the three-fold degenerate excited states. We notice that, while the scaling of the error remains the same, a small reduction of the efficiency is observable in the isotropic case, a consequence of the fact that light polarized along y can be emitted from atoms in the state |e x with a reduction of the overlap between the output mode and detection mode. The increase of the error is better quantified in Fig. C1(b) where the relative difference is plotted. We observe that for the range of array sizes considered here the error increases between 50% and 90%. The value of the optimal beam waist is instead not particularly affected, as expected.