Photon induced atom recoil in collectively interacting planar arrays
aa r X i v : . [ qu a n t - ph ] F e b Photon induced atom recoil in collectively interacting planar arrays
Deepak A. Suresh and F. Robicheaux ∗ Department of Physics and Astronomy, Purdue University, West Lafayette, Indiana 47907, USA Purdue Quantum Science and Engineering Institute, Purdue University, West Lafayette, Indiana 47907, USA (Dated: March 1, 2021)The recoil of atoms in arrays due to the emission or absorption of photons is studied for sub-wavelengthinteratomic spacing. The atoms in the array interact with each other through collective dipole-dipole interactionsand with the incident laser field in the low intensity limit. Shining uniform light on the array gives rise to patternsof excitation and recoil in the array. These arise due to the interference of different eigenmodes of excitation.The relation between the recoil and the decay dynamics is studied when the array is in its excitation eigenstates.The recoil experienced by a subradiant collective decay is substantially larger than from independent atomdecay. A method to calculate the rate of recoil when steady state has been achieved with a constant influx ofphotons is also described.
I. INTRODUCTION
Collective dipole emissions have recently been used inmany novel and interesting applications in the control of theinteraction between light and matter. Reference [1] showedthat when an ensemble of atoms interact and radiate collec-tively, the dynamics are altered due to the interference of theoutgoing light waves. This has led to abundant research in-vestigating a wide variety of concepts including subradiance,superradiance, and collective Lamb shifts [2–12]. These con-cepts are also being used in establishing a link between atomsseparated by more than their resonant wavelength [13–15].Placing the atoms in uniformly spaced arrays will enhancethe co-operative response resulting in increased coupling be-tween the atoms and radiation field [16–24]. Recent experi-ments have achieved the realisation of arrays of atoms with ahigh level of filling efficiency [25–31]. Closely packed atomarrays, where the atom separation is less than the wavelengthof the light, have been found to have highly reflective prop-erties due to the cooperative dipole interactions [18, 19, 31].Such arrays have also been suggested to be efficient optionsfor photon storage [32–34].Collective dipole interactions have found many applica-tions in quantum information [13, 15, 28, 35, 36]. As such,quantum information processing requires an extremely highdegree of fidelity and coherence. Many proposals ignore theeffect of the recoil during the photon atom interactions. How-ever, the recoil on the atoms can cause the information in theinternal states to mix or entangle with the vibrational states ofthe atomic motion leading to a loss of coherence in the manyatom electronic states. The subsequent motion can also affectthe dynamics of the internal states of the system. Thus, theatom recoil can affect the overall robustness of the quantumsystem and introduce avenues for decoherence. The depth ofthe trapping potential required in experiments will be deter-mined by this recoil and having a depth too low could esca-late this problem. But, the role of the photon recoil may becounter intuitive. For example, Ref. [31] discusses that theirexperiments showed that increasing the depth beyond a certain ∗ [email protected] level resulted in a reduced cooperative response. This bringsabout the question of an optimal potential depth in experi-ments. Reference [37] have also proposed to utilize atomicarrays to drive opto-mechanical systems and the recoil energywould play a large role in such interactions. Hence, it is im-perative to accurately model the recoil force effects in suchcases.The recoil in atoms have been studied in various other con-texts. References [38–40] describe how the atom recoil can al-ter the nature of interaction with light, followed by experimen-tal observation of the same in Ref. [41]. Reference [42] havedescribed a theoretical framework for the recoil from sponta-neous emission in an atom near the interface of a topologicalmedium. Reference [43] have also described an approach toanalyze the recoil when light scatters from a Bose-Einsteincondensate of a dilute gas with weak interatomic interaction.The effect of quantized motion on the decay rates and shiftshave been calculated in Ref. [44].The calculations in the current paper continue the explo-ration of the role of recoil in atomic ensembles conducted inRef. [45]. In the current paper, we study the recoil energywhen photons interact with atomic ensembles, more specifi-cally, ordered planar atomic arrays with sub-wavelength inter-atomic spacing. See Fig. 1 for a schematic drawing. We de-scribe a method to determine the recoil momentum and energydeposited in such cases using density matrix master equations.We focus on the regime where the periods of the atomic vibra-tions are much larger than the lifetimes of the internal excitedstates and the duration of the light pulse. This implies thatthe photon interactions can be considered as impulsive forces.Hence, the individual atom’s trapping potential and the timevariation of the position of the atoms can be ignored. We alsowork in the low light intensity regime, limiting the system tohave one excited state at most, to reduce the computationalcomplexity.In Ref. [45], the focus was on the situation where the arraysize is considered near-infinite, thereby ignoring the effects ofthe edges and assuming that most of the bulk of the array ex-periences a recoil to the same degree. The calculations in Ref.[45] were done only for the central atom and scaled accordingto the total number of atoms. In the present work, we explorethe role of finite array size and situations where the excitationof an atom strongly depends on its position in the array. Cal-culations for the individual atoms of the array show that thefiniteness of the array has pronounced effects. When the ar-ray is uniformly driven by light, there are interesting patternsarising in the recoil distribution as well as in the distributionof the excitation itself. The patterns also vary according to thedetuning of the incident driving field. Reference [2] discussesthe interference of the different eigenmodes, each associatedwith a modified lifetime and shift, when interacting with in-cident light. The different eigenmodes interfere and give riseto the patterns observed. We study these eigenmodes of exci-tation which gives us an understanding of the patterns as wellas the relation between the decay lifetimes and the recoil en-ergies deposited in individual atoms.In an atomic array, the direction of recoil on each atom isdependent on the position of the atom in the array. Since thedirection of the photons facilitating the interaction betweenthe atoms will be along the plane of the array, the recoil ofan atom in these directions are dependent on the number andposition of the neighbouring atoms. For example, the cor-ner atom gets imparted a net momentum kick away from thecenter, while the atom at the center experiences no net mo-mentum due to symmetry. On the other hand, in the directionperpendicular to the plane, the dependence of the recoil on itsneighbours will be less direct. Hence, we study the directionalproperties of the recoil experienced by individual atoms withdifferent neighbourhoods.Alternately, when the atom array is driven by a constantlight field, the system reaches a steady state. Experiments us-ing the array as perfect mirrors or as mirrors in a cavity wouldsometimes have a constant influx of photons. There is thena need to calculate the rate of recoil imparted to the atoms.We describe a method to find the rate of recoil energy andmomentum and calculate the rate in a condition similar to theexperiment conducted in Ref. [31].The paper is structured as follows. Section II describes themethods used for the calculations. Section III presents resultsfrom a variety of situations. Section IV presents the conclu-sions followed by the Appendix. II. METHODS
The calculations are done for an ensemble of N atomsarranged in a planar array with sub-wavelength interatomicspacing, d . The atoms have two internal energy levels whichcouple to each other through collective dipole-dipole interac-tions and with an external light field. The array lies in the x-yplane (with minor deviations in the z-direction when usingcurved arrays) and the direction of propagation of the incidentlaser light is chosen to be the z-direction ( k = k ˆ z ) (Fig. 1).The light and the orientation of the excited dipoles is circu-larly polarized in the ˆ e + direction. The circular polarizationis chosen to form symmetric patterns in the array but the mainconclusions should remain valid for other polarizations. Boldfonts denote vectors.To describe the dynamics of the collective dipole emissions,we use the density matrix formalism where the density matrix y xz dk FIG. 1. A schematic of the system considered. The atoms are theblue spheres and are placed in an ordered array in the x-y plane withinteratomic separation d . The incoming light (red wavy arrow) is inthe ˆ z direction with wavevector k evolves according to d ˆ ρdt = − i ~ [ ˆ H, ˆ ρ ] + L (ˆ ρ ) (1)where ˆ ρ is the density matrix of the system, ˆ H is the effectiveHamiltonian and L (ˆ ρ ) is the Liouville super-operator whichdescribes the lossy collective dipole emissions. They shall bediscussed subsequently.To perform density matrix calculations for a large numberof atoms, we adopt a few simplifications. We work in the lowintensity limit where only the singly excited states are rele-vant, meaning only one atom can be excited at any time. Inthis limit, the density matrix dimensions scales more slowlywith increasing number of atoms. This simplification allowsus to perform calculations up to N = 300 for the effects stud-ied below.The ground state of the system which corresponds to all theatoms being in the electronic ground state, is represented by | g i = | g i ⊗ | g i ⊗ ... | g N i , where | g j i denotes the electronicground state of atom j . The raising and lowering operators of j th atom are represented by ˆ σ + j and ˆ σ − j respectively. The statewhere only the j th atom is excited is represented by | e j i =ˆ σ + j | g i . The density matrix describes both the internal degreesof freedom and the positional dependence of the atoms. Thepositional dependence is based on the atoms’ coordinates r j and the density matrix is represented as ˆ ρ = ρ gg ( r , r , ... r N ; r ′ , r ′ , ... r ′ N ) | g ih g | X i ρ e i g ( r , r , ... r N ; r ′ , r ′ , ... r ′ N ) | e i ih g | X j ρ ge j ( r , r , ... r N ; r ′ , r ′ , ... r ′ N ) | g ih e j | X i,j ρ e i e j ( r , r , ... r N ; r ′ , r ′ , ... r ′ N ) | e i ih e j | (2)where i, j go from 1 to N. The coefficients describe the spa-tial dependence of the density matrix and are functions of N positional coordinates.Since operators acting on the left and right sides of the den-sity operator will act on different indices, we define the fol-lowing convention r ij ≡ r i − r j ; r ′ ij ≡ r i − r ′ j ; r ′′ ij ≡ r ′ i − r ′ j (3)The effective Hamiltonian of the system will have (i) Thelaser interaction, (ii) the collective dipole interaction and (iii)the center of mass motional Hamiltonian consisting of the ki-netic energy and the trapping potential of each atom. For thelaser interaction, we work in the rotating wave approximationwith ˆ H L = ~ X j (cid:20) − δ ˆ σ + j ˆ σ − j + Ω2 ˆ σ + j e i k · r j + Ω ∗ σ − j e − i k · r j (cid:21) (4)which describes the Hamiltonian for a plane wave incidentlaser with Ω as the Rabi frequency, δ as the detuning and k = k ˆz as the initial wavevector of the incoming photons.The collective dipole emission also participates in the effec-tive Hamiltonian based on the imaginary part of the Green’sfunction ( Im { g ( r ij ) } ) and is given by ˆ H dd = ~ X i = j Im { g ( r ij ) } ˆ σ + i ˆ σ − j (5)where the g ( r ij ) is the free space dyadic Green’s functiongiven by Eq. (7). When calculating the commutator of the ef-fective Hamiltonian ( ˆ H = ˆ H L + ˆ H dd ) with the density matrix ˆ ρ in Eq. (1), indices of the position vector must be carefullyimplemented. When right-multiplying ˆ ρ with ˆ H , the corre-sponding primed indices must be used ( r ′ j for ˆ H L and r ′′ ij for ˆ H dd ).The decay effects of the collective dipole interaction will becaptured by the Lindblad term L (ˆ ρ ) = X i,j (cid:2) Re { g ( r ′ ij ) } ˆ σ − i ˆ ρ ˆ σ + j − Re { g ( r ij ) } ˆ σ + i ˆ σ − j ˆ ρ − Re { g ∗ ( r ′′ ij ) } ˆ ρ ˆ σ + i ˆ σ − j (cid:3) (6)where the i = j terms are the usual single atom decay part ofthe Lindblad operator. The Green’s function g ( r ij ) is given by g ( r ij ) = Γ2 (cid:20) r ij · ˆ q i )(ˆ r ij · ˆ q ∗ j ) − (ˆ q i · ˆ q ∗ j )2 h (1)2 ( kr ij )+ ( ˆ q i · ˆ q j ∗ ) h (1)0 ( kr ij ) (cid:21) (7)where ˆ q i is the dipole orientation of the i th atom, r ij = | r ij | is the norm of r ij , ˆ r ij = r ij /r ij is the unit vector along r ij , Γ is the decay rate of a single atom and h (1) l are the outgoingspherical Hankel function of angular momentum l . When i = j , that is when r ij = 0 , the imaginary part of the functionbecomes undefined, while the real part is defined. Hence, weredefine g ( r ij ) to be g ( r ij ) = g ( r ij ) f or i = j = Re { g ( r ij ) } f or i = j (8) A. Slow Oscillation approximation
In this paper, the calculations are primarily in the limitwhere the timescales of the atomic motion in the trap aremuch slower than the timescales for the evolution of the in-ternal states. This means that during the evolution of the in-ternal states, the position of the atom does not change. Thisis a sudden approximation and allows us to drop the motionalHamiltonian with the kinetic energy and trapping potential.Typically, the frequencies of the trapping potentials will bearound a few 10s of KHz and are much slower than the decayrates which are usually of the order of MHz. We also assumethat the spread of the initial center of mass position is smallcompared to the separations. This will allow using the perfectposition of the atoms in the calculations.In the calculations where the system settles in the groundstate, at infinite time, only the ρ gg term is left non-zero and itsspatial dependence will be of the form ρ gg ( r , r , ... ; r ′ , r ′ , ... ) = ρ ( r , r , ... ; r ′ , r ′ , ... ) × F ( r , r , ... ; r ′ , r ′ , ... ) (9)where ρ is the spatial dependence at the initial time and F de-scribes the evolution of the ground state spatial dependence.To calculate the momentum and kinetic energy change of theatoms, we take the expectation value of the corresponding op-erator using the density matrix. h p j i = ~ i Z [ ˆ ∇ j ρ gg ] δ ( r ′ ) δ ( r ′ ) ...d r d r ′ d r d r ′ ... = T r [ˆ p i ( ρ F )] = h p j i + ∆ p j (10)where ˆ ∇ j = ˆx ∂∂x j + ˆy ∂∂y j + ˆz ∂∂z j is the momentum operatorof the j th atom. The term h p j i denotes the initial expecta-tion value of the momentum and derives from ρ . Hence, thechange in momentum can be calculated from the function F using ∆ p j = ~ i Z [ ρ ( ˆ ∇ j F )] δ ( r ′ ) δ ( r ′ ) ...d r d r ′ d r d r ′ ... = T r [ ρ (ˆ p i F )] (11)For the Kinetic energy, we follow a similar reasoning anduse the KE operator K j = p j / (2 m ) to get ∆ K j = − ~ m Z [ ρ ( ˆ ∇ j F )] δ ( r ′ ) δ ( r ′ ) ...d r d r ′ d r d r ′ ... = Re { T r [ ρ ( ˆ K j F )] } (12)The density matrix is propagated in time using the Runge-Kutta second order integration until the system completely de-cays into the ground state. This ground state density matrixcoefficient is used to evaluate F ( r , r , ... ; r ′ , r ′ , ... ) . Thisis in turn used to calculate the ∆ p and ∆ K in three dimen-sions using Eq. (11) and Eq. (12) by using symmetric 2-pointand 3-point differentiation. To calculate the derivatives ∇ j ,the positions of either the primed ( r j ) or the unprimed coordi-nates ( r ′ j ) are shifted by a small distance δr and evaluated. Fora more detailed explanation of the methods followed, refer toSection II of Ref. [45]. B. Interatom distance
The process of subradiance can be thought of as the destruc-tive interference of the wavefunctions of the light emitted byspontaneous decay of excited atoms. For example, when thereare two atoms close together to the point of overlap ( d → ),the wavefunctions of the emitted light will cancel out if theyhave opposing phases. This results in a prolonged lifetime forthe excited state. For two atoms, the out of phase state remainsas the subradiant state until the separation of d ∼ λ/ . As theseparation goes beyond half a wavelength, the light acquiresan extra phase of π which causes the in-phase states to be sub-radiant until d ∼ λ . This behavior continues as we increase d and oscillates with a period of λ/ . The maximum effectof subradiance and superradiance possible decreases with in-creasing d and becomes small beyond d ∼ λ .This behavior carries over to the periodically placed atomsin arrays. For interatomic distances less than ∼ λ/ , the sub-radiant states have adjacent atoms out-of-phase, while the su-perradiant states have them in-phase. Between λ/ and λ , thesubradiant states have adjacent atoms in-phase. Since we wantto primarily focus on subradiant states, and exciting atomswith adjacent atoms being in-phase is easier to experimen-tally realise, we choose the range of interatomic distance tobe between . λ to . λ . III. PHOTON RECOIL ENERGY AND MOMENTUM
When a photon is absorbed or emitted by a single atom, thephoton imparts a momentum kick ~ k and recoil kinetic energy E r = ~ k / (2 m ) . But when there is an ensemble of atomsand collective dipole interactions take place, Ref. [45] showedthat the energy deposited is different and depends on collec-tive decay dynamics. We delve deeper into this topic and dis-cuss the directional properties of the kicks and its relation tothe decay properties of the collective ensemble. Since a sin-gle photon undergoing perfect reflection on an atom imparts ~ k momentum, the ∆ p z / (2 ~ k ) serves as a good measure ofreflectance.The two different factors that contribute to the kicks haveslightly different effects. In the out-of-plane direction, thecollective dipole emissions emit light symmetrically on bothsides and hence contribute to no net momentum kick in this di-rection, but will still contribute to the recoil energy deposited.In the in-plane direction, the atoms exchange photons amongeach other and, hence, the momentum and energy depositedwill depend on the position of the atom within the array. Thecontribution from the laser will only be in the out of plane di-rection and has non-zero contributions to both momentum andkinetic energy deposition. A. Eigenstates
Reference [2] discussed the interference of many eigen-modes of the system which contributes to the cooperativeemission. To get an understanding of the effect of the de-cay rate on the recoil energy, we analyze the photon recoilmomentum and energy deposited when the initial state is aneigenstate of the excitation. When there is no driving inter-action with the laser, the eigenstate of the excitation is theeigenstate of the dyadic Green’s function in free space. Thesestates decay exponentially into the ground state but are sta-tionary with respect to their excitation pattern, i.e., the atomsmaintain the relative phase and magnitude with respect to eachother. We define the Green’s tensor G ij = g ( r ij ) which is anN × N-dimensional matrix, and V α is an N-dimensional vec-tor. The eigenvalue equation is X j G ij V jα = G α V iα = ( γ α i ∆ α ) V iα (13)where i , j are atom indexes and α denotes the index for theeigenstate. G α is the eigenvalue and V α is the correspondingeigenvector. The rate of decay is given by the real part of theeigenvalue, γ α , and the shift in energy is given by the imag-inary part, ∆ α . Since the Green’s function is not Hermitian,the regular orthogonality conditions do not apply. Therefore,the vectors V α have to be normalized to satisfy, X i V iα V iα ′ = δ α,α ′ (14)This relation should be used with care when there are degen-erate eigenstates in the system. Each set of degenerate vec-tors must be orthogonalized to follow this condition. Theseeigenstates form interesting and symmetric patterns on the ar-ray and exhibit similarities to TEM modes of light. As notedearlier in Section II B, the adjacent atoms in the subradiantmodes tend to be in-phase while the superradiant modes havethe adjacent atoms out-of-phase in the specified range of d .The array is initialized to one of the eigenstates, β , lead-ing to the electronic part of the density matrix starting as ρ ( t = 0) = N β | V β ih V β | , where N β = ( P i | V iβ | ) − isa normalizing factor to have only one excitation at the initialtime; note the magnitude in the definition of N β . The systemis evolved in time until it completely reaches the ground state,at which point, the ∆ p and ∆ K of each atom are calculated.The total kinetic energy deposited in the array for eacheigenstate is inversely proportional to the decay rate of thatstate. This implies that highly subradiant states will have veryhigh photon recoil energies. Another trend observed was thatthe kick on each individual atom was roughly proportional tothe excitation probability of that atom. This, while being anexpected result, allows for an easier way to look at the distri-bution of the energy deposition over the array. Commonly, themost subradiant mode at this range of interatomic distances issimilar to a Gaussian distribution ( T EM like) on the array.This implies that the deposition of the recoil is concentratednear the center, while the atoms close to the edges have rela-tively low recoil.When the system is initialized to its eigenstates, it is possi-ble to analytically obtain an expression to calculate the recoilafter the system and decayed to the ground state. This expres-sion can then be used to find the recoil in any arbitrary initialstate configuration by decomposing it into its eigenstates. Thecoefficient of the ground state density matrix at infinite timewill be ρ gg ( ∞ ) = X α,α ′ X j,j ′ Re { g ( r ′ jj ′ ) } V jα U ∗ j ′ α ′ C α,α ′ (0) G α + G ′′∗ α ′ (15)where G ′′ α ′ and U jα ′ are the eigenvalue and eigenvector ofthe tensor G ′′ ij = g ( r ′′ ij ) . The contribution of each eigenstateis given by C α,α ′ (0) = P jj ′ V jα U ∗ j ′ α ′ ρ e j e j ′ (0) . If the ini-tial state is an eigenstate β , then it is given by C α,α ′ (0) = N α δ αβ P i V ∗ iβ U ∗ iα ′ . and the term where α = α ′ = β will bethe most dominant. The derivation can be found in AppendixA.When initialized to eigenstates, this expression can be usedto calculate the momentum and kinetic energy kicks usingEqs. (11) and (12). The calculation using Eq. (15) is mucheasier than solving for large density matrices over extendedtimes and it also provides intuition towards the trends ob-served. The denominator of the dominant term, G β + G ′′∗ β ≈ Re {G β } = γ β explains why the kicks are inversely propor-tional to the decay rates of the eigenstates. The term V jβ U ∗ j ′ β also explains why the kicks on each atom is proportional totheir excitation probabilities.An expression for kick in the out-of-plane direction can beobtained to express the ∆ K z depending only on the decayrate of the eigenvector ( γ α ). The derivation can be found inAppendix A 1. ∆ K z E r = 25 Γ γ α (16) ∆ K z / E r (γ α / Γ) −1 (a) ( ∆ K x + ∆ K y ) / E r ( γ α / Γ ) − (b) FIG. 2. The net kinetic energy kick and the linear fit for an array of × array, initially in an eigenstate, with respect to the inverse ofthe decay rate of the eigenvalue. Blue dots denote the calculated plotpoints and the red line denotes the linear fit of the data. (a) denotesthe kick in the out-of-plane direction which is proportional to thelifetime with slope 2/5, as expected from Eq. (16). (b) denotes thenet kick in the plane of the array. This was a surprising result because there is no photon me-diated atom-atom interaction in the z-direction and there isjust a single photon emitted. Hence one might expect the ∆ K z /E r to be independent of the lifetime of the excited state.However, the results from Eq. (16) show that when there isspontaneous emission, the decay lifetimes play an importantrole in recoil of the atoms in the array.Figure 2 shows the results for the net ∆ K deposited in theout-of-plane and in-plane directions with respect to the decaylifetimes of the eigenstates using full density matrix calcula-tions. While the out-of-plane kicks are highly proportional,the in-plane kicks have a more complicated dependence andare only roughly proportional.This section discussed the recoil momentum and energycalculation for the case when the system is in either an eigen-state or an arbitrary excited state and allowed to decay intothe ground state. This is a simplified case where there areno incoming photons and it is highly unlikely to be experi-mentally seen. Nevertheless, it provides insight into how thelifetimes of the state affects the recoil. We notice that sub-radiant states experience more recoil and superradiant statesexperience lesser recoil than expected from individual atoms. B. Cavity
After having seen the subradiant properties of atom arraysfrom Fig. 2, we can further increase the subradiance by trap-ping the light between two arrays and forming a cavity. Byslightly curving the arrays, the cavity formed can greatly im-prove light retention, in turn making the system more subra-diant. In this section, we calculate the recoil experienced bythe atoms in highly subradiant cavities.The system consists of N atoms making up two arrayseach with N atoms. The arrays are separated by a distance L .This is similar to the system describes in Fig. 1 of Ref. [13].To find the most efficient cavity mode, we find the eigenstatesof the Green’s function for the cavity and find the most subra-diant mode. We vary the separation L and the curvature of thearray until we find the optimal parameters for maximum sub-radiance. The imaginary part of the eigenvalue also providesthe detuning that the cavity mode should be driven at.The large lifetimes of the subradiant states prevent calcula-tion of the kick in all the atoms of the array within reasonablecomputation time using the density matrix method. Hence weuse the second method described in Sec. III A to calculate therecoil.To characterize the light retention capacity of the cavity,we calculate the finesse of the cavity. In this case, we definethe finesse using the intensity distribution of the cavity. Thecavity is driven using lasers of fixed detuning, while scanningthe separation of the arrays. The finesse, F , will be defined asthe free spectral range, λ/ , divided by the full width at halfmaximum in separation of the intensity: F = λ/ δL F W HM (17)This gives an estimate of the number of times a photonbounces in the cavity before decaying.In the following sub sections, we discuss the results cal-culated for a few typical cavities with different decay rates.The cavities are initialized into their most subradiant eigen-state with a single excitation and evolved until they reach theground state, similar to Sec. III A. It is also compared with thefinesse of the cavity. In each case, we take a × array andthe separation has been optimized around L = 20 λ .
1. Spherical Mirrors
Reference [13] utilizes a set of spherically curved atom ar-rays to form highly subradiant states. Since light trapped ina cavity will have a Gaussian intensity profile, sphericallycurved mirrors match the spherical wavefront of the lightand provides good energy retention. The mirrors are placedconfocally. Following similar parameters to Ref. [13], with d = 0 . λ , separation L = 19 . λ , we attain the most subra-diant mode to be extremely subradiant, with decay rates reach-ing ∼ − Γ . The finesse of this particular configuration wasfound to be 1250. The kick the center atom experiences wascalculated to be around E r in the out-of-plane direction and E r in the in-plane directions. This progressively re-duces as we go closer to the edge in the shape of a Gaussian(See Fig. 3). This trend is expected because the most subra-diant eigenmode is typically a T EM mode (as noted in thesupplementary of Ref. [13]). The atoms in the corner receiveda total of only . × − E r . The net kick on the array leadto an energy increase of ∼ E r .Unfortunately, this result is not completely valid becausethe decay rate is too small for the slow oscillation approxima-tion to be satisfied. If the duration of the force on the arrayis of the order of the oscillation period, the impulsive natureof the recoil force will not be valid. In this timescale, theexact nature of the way fields interact with the atoms is notprecisely known. A fully quantum mechanical treatment con-sidering the motion of the particles as quantum oscillators isrequired to get a better understanding of the dynamics. Whilesuch a calculation would be very difficult, the enormous sizeof the kicks from our simplified treatment suggest that an at-tempt should be made.
2. Parabolic Mirrors
Another typical type of mirrors used in cavities areparabolic mirrors. Confocal cavities are only marginally sta-ble and any small errors in the positions of the mirrors willcause destabilization. Hence, we use a cavity in the regionbetween confocal and concentric to provide some room forerrors. This configuration, with d = 0 . λ , parabolic focus f = 15 λ and separation L = 19 . λ , while not as subadiantas the previous case, has a decay rate of ∼ − Γ . This decayrate also is at the edge of the slow oscillation approximation.The finesse for this cavity was calculated to be 263. The centeratoms experienced E r recoil in the z-direction and . E r inthe x and y direction. A similar trend of decreasing kick in theedge atoms is seen and a total of only . × − E r was de-posited in the corner atom. The total energy deposited on thewhole array in this decay process is ∼ E r
3. Plane mirrors
Curving a plane array of atoms as required in the abovecases may prove to be complicated. Hence, to compare, wediscuss the case where a cavity is formed by two plane mir-rors. In this configuration, with d = 0 . λ and separation L = 20 . λ , the decay rates reach ∼ − Γ . The finessewas calculated to be approximately 14. The center atom re-ceived a total of . E r , while the corner atom received atotal of . E r . The array as a whole, received a kick of E r in the out-of-plane direction and E r in the in-planedirections. C. Pulsed Laser
To simulate the effects of a single photon interacting withthe array, a low intensity pulse of light is incident on it. A y / d x/d FIG. 3. The recoil energy deposited in the z-direction on one arrayof a cavity with two × curved arrays with d = 0 . λ and L = 19 . λ corresponding to the system in Sec. III B 1. Each celldenotes one atom in the array. The color scale is in units of the recoilenergy. Note the center atom has a kick of ∼ E r . planar array with N atoms is initially in the ground state and alaser pulse with a Gaussian time profile is incident. The atomsare then evolved until they reach ground state. The ∆ K and ∆ p of each atoms are calculated in this time frame and arecompared to the probabilities of excitation of each atom. TheGaussian time profile of the light pulse of the form Ω( t ) = Ω e − t /t w (18)with peak Rabi frequency Ω and pulse width t w . The Ω iskept low enough to not go beyond the single excitation limit.Similar to the trend in the previous section, we see that therecoil is proportional to the integral of the excitation probabil-ities of each atom. That is, the excitation patterns are similarto the recoil distribution pattern in the array.This is similar to the calculations done in Ref. [45], exceptthat Ref. [45] assumes that there is a nearly uniform excita-tion of the atoms in the array and only calculates the dynamicsof the central atom. By calculating the kicks in each atom ofthe array, we see patterns arise that vary with the detuning ofthe laser. Figure 4 shows the comparison between the exci-tation pattern and the recoil distribution in the array. Figure4(a) shows the time integrated excitation probability of eachatom, while Fig. 4(b) depicts the ∆ K z deposited on eachatom. There is significant variation in the amount of energydeposited on the atoms and the corner atoms experienced onlyhalf the recoil energy of the atoms with the maximum recoil.These patterns are a combination of the eigenstates of the ex-citation as seen in Section III A. The selection of the eigen-states depends on both the pattern as well as the detuning ofthe incoming light.Experimentally, achieving a filling fraction for the ar-ray is difficult. Hence, the effects of the array missing a singleatom have been studied. The scale of the disturbance causedby a missing atom depends on the contribution it would havemade if present. If an atom is missing where the excitation is naturally weak, it has less effect on the excitation patternand recoil (Fig. 5a). However, if the atom is missing wherethe excitation is large, there is a more substantial effect seenspanning a few nearby atoms (Fig. 5b). The recoil in the in-plane direction for the neighbouring atoms will also change.The recoil they experience in the direction towards the miss-ing atom will be increased.As explained in Section II B, we have chosen the range of d between . λ to . λ so that coupling subradiant modes,which have adjacent atoms in-phase, will be easier with uni-form illumination. For other experiments that prioritize work-ing with superradiant states, the interatomic distance could bechosen in the range d < . λ for easy coupling with uniformlight. D. Steady State
We have so far studied at the recoil received after all thedynamics of the system are completed and all atoms return tothe ground state. But in most experiments and applications,it is also important to understand the behavior of these kickswhen there is a steady state involved. Instead of finding thetotal momentum and energy deposited, a more useful quantitywill be the rate of energy deposited in each atom in steadystate.The methods described earlier as well as in Ref. [45] de-scribe the process in which the final state of the system is theground state. Hence we developed a different method to de-termine the rate of momentum and energy deposited. Whenthe system has reached steady state, we can approximate thatup to first order in time, all the density matrix elements exceptthe ground state coefficient reaches equilibrium. This meansthat the only linear time dependent term will be the groundstate density matrix coefficient ( ρ gg ). Taking the positionalderivatives on ˙ ρ gg similar to Eq. (11) and Eq. (12) yields therate of the momentum and energy kicks. The density matrixcan be evolved until steady state is reached and the ˙ ρ gg can becalculated. Alternatively, the master equation can be solvedanalytically to obtain an expression for ˙ ρ gg . ˙ ρ gg ( t ) = − i X j Ω ∗ ( r j )2 ρ e j g + i X j Ω( r ′ j )2 ρ ge j + X ij Re { g ( r ′ ij ) } ρ e i e j (19)where Ω( r j ) describes the spatial profile of the incoming laseras well as carrying a phase factor e − i k · r j . The derivation andthe procedure to calculate the density matrix terms at steadystate are discussed in Appendix B. For this approximation tohave good accuracy, it is important that we stay within the lowintensity limit and use low Rabi frequencies.The momentum and the energy deposited increase as afunction of Ω which is expected as it is proportional to thenumber of incoming photons. To calculate the number ofphotons incident on an atom, we can integrate the intensityarriving on the area corresponding to one atom ( d ) per unittime and divide by the energy of a single photon ( ~ c/λ ). The y / d x/d (cid:1) (cid:0) (cid:2) (cid:3) (cid:4) (cid:5) (cid:6) (cid:7) (cid:8) (a) y / d x/d (b) FIG. 4. The excitation and recoil distribution pattern on a × atom array with d = 0 . λ when excited by a laser pulse with peak Rabifrequency Ω = 0 . and pulse width t w = 16 / Γ at zero detuning. Each cell denotes one atom in the array. (a) shows the integral of theexcitation probability in time for each atom in arbitrary units and (b) shows the ∆ K z deposited on each atom of the array per photon incidenton the atom, as a function of E r . y / d x/d (cid:11) (cid:12) (cid:13) (cid:14) (cid:15) (cid:16) (cid:17) (cid:18) (cid:19) (a) y / d x/d (cid:20) (cid:21) (cid:22) (cid:23) (cid:24) (cid:25) (cid:26) (cid:27) (cid:28) (b) FIG. 5. The effects of an atom missing from an array of × atoms with d = 0 . λ and pulsed light illumination with no detuning. Theplots show the total integrated probability of excitation of each atom in arbitrary units. Compare with Fig. 4(a) which is the case where noatoms are missing. (a) has atom at coordinate (2,2) missing. It has a larger impact and affects up to second nearest neighbour atoms. (b) hasatom at coordinate (3,3) missing. It has a smaller impact and only affects the immediate neighbours. number of photons incident on each atom in one lifetime of asingle atom is ν Γ = 2 π (cid:18) dλ (cid:19) (cid:18) ΩΓ (cid:19) (20)Hence Ω can be calculated to correspond to different ratesof photon arrival. Using this, we can simulate the parametersof the recent experimental work of Ref. [31] to calculate therate of deposition. For a × array at peak resonance, with d = 0 . and an influx of . × − photons per lattice sitein one lifetime of a single atom, the calculations showed thatthe average total energy deposited is . E r per incoming pho-ton. Each individual atom experienced a total recoil rangingfrom . E r to . E r per photon incident on each atom. Thetotal average momentum imparted on each atom is . ~ k per incoming photon, which corresponds to a reflectance of .Reference [18] discusses a null transmission situation foran infinite plane array with d = 0 . λ . Figure 8 of Ref. [45]also shows that the reflectance is near unity for an infinite ar-ray and estimates the same for a × array using the ∆ p z of the center atom. The emergence of the excitation patternsimplies that this will be an overestimation for non-infinite ar-rays. The new calculations show that for a × array thereflectance only reaches when taking the average overall the atoms. Although, this is due to edge effects and recal-culating the average ignoring the edge atoms brings it up to .When the array is illuminated by the laser to excite par-ticular eigenstates, the momentum imparted by the laser wasinversely proportional to the eigenmode’s decay rate. This P r obab ili t y d / G FIG. 6. The decomposition of the state into its eigenmodes when ex-cited by uniform light for a × array with d = 0 . λ versus the de-tuning δ . The thick lines show the probability of each eigenmode inrescaled arbitrary units while varying the detuning. The thin verticallines show the line shift ∆ α of the eigenmode with the correspond-ing color and dash-type. Only the 5 states that have a significantcontribution are shown. is the result of the altered absorption rate of the incominglight. Since the system is at equilibrium, the absorption ratemust match the decay rate of that eigenmode. The energy de-posited due to the collective dipole interaction followed a sim-ilar trend to the previous discussions and was also inverselyproportional to the decay rate. This can also be seen by thepresence of the eigenvalue G α in the denominator of the dom-inant terms in the analytical solutions (Eq. (B3)).When driving with lasers, the detuning plays a role in de-ciding the contribution of the different eigenmodes as speci-fied in Section III C. The eigenmodes closest to the symmet-ric Dicke state will couple better with the uniform incominglight. Since each eigenstate is associated with an energy shift( Im {G α } = ∆ α ), matching the detuning also contributes toselecting the eigenstate. By using these factors, it is possibleto influence which eigenstates contributes to the coupling ofthe atoms to the laser, to a certain extent. By using the or-thogonality relations defined in Eq. (14), we can determinethe contribution from each eigenstate as we vary the detuning.As seen from Fig. 6, the maximum contribution from a par-ticular eigenstate is when the detuning matches with the shiftassociated with it. IV. CONCLUSION
We present the calculations for the recoil energy and mo-mentum deposited in an array of atoms interacting collectivelywith light. We studied the directional properties of the re-coil and explored the effects of the eigenmodes of excitationin the system. The recoil energies added to each atom areproportional to the decay lifetimes of the eigenmodes and themore subradiant states experience more recoil than might beexpected. This implies that recoil effects should be consid-ered more carefully in experiments involving highly subradi- ant states. We note that when driving the array with uniformlight, the excitation and recoil distribution in the array is notuniform but forms patterns which vary with the detuning. Thepatterns form due to the interference of the different eigen-modes coupling with the uniform light. These patterns alsostrongly contribute to the recoil of individual atoms. Calcu-lations were done to determine the rate of energy depositionwhen the system reaches a steady state due to a constant influxof light.We also calculate the recoil effects in cavities made ofcurved atomic arrays and found pronounced effects as the de-cay rate decreased. However, our calculations are not appli-cable when the decay lifetimes becomes comparable to thetimescale of the atomic vibration. Thus, our results for photoncavities (Section III B) can only be taken as a qualitative esti-mate. It is unclear how the momentum transfer occurs on suchtimescales. A fully quantum approach would require largeamounts of computational resources but exploring this situa-tion with a few simplifying assumptions would give a betteranswer to this question. Reference [37] explored the oppositeregime where the internal state lifetimes are far larger thanthe motional timescale. It would be interesting to explore thephysics at the interface of these two regimes.
ACKNOWLEDGMENTS
We thank Akilesh Venkatesh for critical reading an earlyversion of this paper. This work was supported by the Na-tional Science Foundation under Award No. 1804026-PHY.This research was supported in part through computational re-sources provided by Information Technology at Purdue Uni-versity, West Lafayette, Indiana.
Appendix A: Eigenstate decay analytical calculation
We derive an analytical method to calculate the kicks whenthe system starts off with a single excitation in any arbitraryexcitation pattern. There is no driving from the laser inter-action. We write the state in terms of the eigenmodes V α and U α which are the eigenvectors of G ij = g ( r ij ) and G ′′ ij = g ( r ′′ ij ) respectively. X j ′ G jj ′ V j ′ α = G α V jα X j ′ G ′′ jj ′ U j ′ α = G ′′ α U jα (A1)Since G and G ′′ are non-Hermitian, their eigenvectors are notorthonormal and V α and U α have to be normalized to follow X j V jα V jα ′ = δ αα ′ X j U jα U jα ′ = δ αα ′ (A2)Hence, we can write the g ( r ij )s as follows g ( r ij ) = X α V iα V jα G α g ∗ ( r ′′ ij ) = X α U ∗ iα U ∗ jα G ∗ α (A3)0In the case where there is no laser interaction, the terms ρ e j g and ρ ge j will vanish giving ˆ ρ ( t ) = ρ gg ( t ) | g ih g | + X jj ′ ρ e j e j ′ ( t )ˆ σ + j | g ih g | ˆ σ − j ′ (A4)where the positional dependence has been suppressed for no-tational convenience. We can write the density matrix equa-tion as d ˆ ρdt = X i,j (cid:2) − g ( r ij )ˆ σ + i ˆ σ − j ˆ ρ − g ∗ ( r ′′ ij )ˆ ρ ˆ σ + i ˆ σ − j + 2 Re { g ( r ′ ij ) } ˆ σ − j ˆ ρ ˆ σ + i (cid:3) (A5)Solving this we get, dρ e j e j ′ ( t ) dt = X k (cid:16) − g ( r jk ) ρ e k e j ′ − ρ e j e k g ∗ ( r ′′ kj ′ ) (cid:17) (A6)Decomposing the ρ e j e j ′ in terms of the eigenfunctions V α and U ∗ α ′ ρ e j e j ′ = X αα ′ V jα U ∗ j ′ α ′ C αα ′ (A7)Using this, we get dC αα ′ dt = − ( G α + G ′′∗ α ′ ) C αα ′ = ⇒ C αα ′ ( t ) = C αα ′ (0) e − ( G α + G ′′∗ α ′ ) t (A8)Since the ρ gg term only arises from the last term of Eq.(A5), we get dρ gg dt = X jj ′ Re { g ( r jj ′ ) } ρ e j e j ′ = X αα ′ X jj ′ Re { g ( r jj ′ ) } V jα U ∗ j ′ α ′ C αα ′ ( t ) (A9)Since C αα ′ ( t ) are exponentials, ρ gg ( ∞ ) = X α,α ′ X j,j ′ Re { g ( r ′ jj ′ ) } V jα U ∗ j ′ α ′ C α,α ′ (0) G α + G ′′∗ α ′ (A10)The expression for C αα ′ (0) can be obtained from Eq. (A7)using the orthogonality relations (Eq. (A2)) to be C α,α ′ = X jj ′ V jα U ∗ j ′ α ′ ρ e j e j ′ (A11)If the system is initialized to one of its eigenstates β with ρ ( t = 0) = N β | V β ih V β | , then C α,α ′ ( t = 0) = N β δ αβ X i V ∗ iβ U ∗ iα ′ (A12)where N β = ( P i | V iβ | ) − is the normalization factor.
1. Out-of-plane eigenstate recoil energy
When we analyze the kinetic energy kick imparted on theout-of-plane direction when the system is initialized to aneigenstate, we get a proportionality equation between the en-ergy deposited and the state lifetimes. This can be derivedby taking the second derivative from Eq. (12) only in the z-direction. When the system is initialized to an eigenstate β ,the dominant terms in the summation of Eq. (A10) will bethe terms with α = α ′ = β . Hence the summation over theeigenstates will vanish giving ρ gg ( ∞ ) = X j,j ′ Re { g ( r ′ jj ′ ) } V jβ U ∗ j ′ β C β,β (0) G β + G ′′∗ β (A13)To take derivatives in the z-direction, we shift the primedcoordinates, r ′ j = r j + δz . This would result an insignificantchange in G ′′ β = G β + 10 − O ( δz ) and a second order changein the eigenvectors, U jβ = V jβ + 10 − O ( δz ) . Hence whentaking up to second derivatives, we can assume G ′′ β ≈ G β and U jβ ≈ V jβ . This implies that G β + G ′′∗ β = 2 Re {G β } = γ β and V jβ U ∗ jβ = | V jβ | . Using Eq. (A12), C β,β (0) will alsobecome N β . This gives ρ gg ( ∞ ) = X j,j ′ Re { g ( r ′ jj ′ ) } V jβ U ∗ j ′ β N β γ β (A14)When taking the second derivative for the atom k , only theterm j = j ′ = k will be non vanishing. ρ gg ( ∞ ) = 2 Re { g ( r ′ kk ) }| V kβ | N β γ β (A15)At the limit of δz → , lim δz → ∂ Re { g ( r ′ kk ) } ∂z = − k Γ5 (A16)Hence the kinetic energy deposited on atom k will be ∆ K ( k ) z = ~ k m | V kβ | N β γ β (A17)The total kinetic energy deposited on the whole array in thez-direction will be (using N β = ( P i | V iβ | ) − ) ∆ K z E r = 25 Γ γ β (A18) Appendix B: Steady state analytical calculation
In this section, we derive analytically a method to find therecoil at steady state be decomposing the state into its eigen-states. For simplicity in this derivation, ρ e j g , ρ ge j and ρ e i e j w j , ˜ w j and ˜ ρ ij respectively and the po-sitional dependence will not be explicitly shown. That is, Eq.(2) will be represented as ˆ ρ = ρ gg | g ih g | + X i w i | e i ih g | + X j ˜ w j | g ih e j | + X i,j ˜ ρ ij | e i ih e j | (B1)A constant laser with detuning δ and Rabi frequency Ω( r j ) is incident on the j th atom of the array. By using Eq. (1), wecan obtain the rate of change of the coefficients of the densitymatrix. ˙ w j = − i Ω( r j )2 ρ gg + iδw j − X k g ( r jk ) w k (B2a) ˙˜ w j = i Ω ∗ ( r ′ j )2 ρ gg − iδ ˜ w j − X k g ∗ ( r ′′ jk ) ˜ w k (B2b) ˙˜ ρ ij = − i Ω( r i )2 ˜ w j + i Ω ∗ ( r ′ j )2 w i − X k g ( r ik )˜ ρ kj − X k g ∗ ( r ′′ kj )˜ ρ ik (B2c) ˙ ρ gg = X ij Re { g ( r ′ ij ) } ˜ ρ ij − i X j Ω ∗ ( r j )2 w j + i X j Ω( r ′ j )2 ˜ w j (B2d)where Ω( r j ) describes the spatial profile of the incominglaser as well as carrying a phase factor e − i k · r j . We can de- compose the coefficients w j , ˜ w j and ˜ ρ ij in terms of the eigen-states V jα and U jα (defined in Eq. (A1)) to get w ′ α , ˜ w ′ α and ˜ ρ ′ αα ′ using w ′ α = X j V jα w j w j = X α V jα w ′ α (B3a) ˜ w ′ α = X j U ∗ jα ˜ w j ˜ w j = X α U ∗ jα ˜ w ′ α (B3b) ˜ ρ ′ αα ′ = X ij V iα U ∗ jα ˜ ρ ij ˜ ρ ij = X αα ′ V iα U ∗ jα ˜ ρ ′ αα ′ (B3c)In the first order in time, the recoil will only build up onthe ground state and the change in the other coefficients willbe zero, i.e., ˙ w j , ˙˜ w j and ˙˜ ρ ij will be zero. Therefore, we cansolve for w ′ α , ˜ w ′ α and ˜ ρ ′ αα ′ using Eqs. (B2a), (B2b) and (B2c),which gives w ′ α = − iρ gg X i Ω( r i ) V iα G α − iδ (B4a) ˜ w ′ α = iρ gg X i Ω ∗ ( r ′ i ) U ∗ iα G ′′∗ α + iδ (B4b) ˜ ρ ′ αα ′ = ρ gg X i Ω( r i ) V iα G α − iδ X j Ω ∗ ( r ′ j ) U ∗ jα ′ G ′′∗ α ′ + iδ (B4c)The term ρ gg can be calculated using T r [ˆ ρ ] = 1 , but in thelow intensity limit, it will be close to . These equations canbe used in conjunction with Eqs. (B2d) and (B3) to calculate ˙ ρ gg to find the recoil kinetic energy and momentum. [1] R. H. Dicke, “Coherence in spontaneous radiation processes,”Phys. Rev. , 99–110 (1954).[2] Robert J. Bettles, Simon A. Gardiner, and Charles S.Adams, “Cooperative ordering in lattices of interacting two-level dipoles,” Phys. Rev. A , 063822 (2015).[3] Nicholas E. Rehler and Joseph H. Eberly, “Superradiance,”Phys. Rev. A , 1735–1751 (1971).[4] R. Friedberg, S.R. Hartmann, and J.T. Manassah, “Frequencyshifts in emission and absorption by resonant systems ot two-level atoms,” Physics Reports , 101 – 179 (1973).[5] M. Gross, C. Fabre, P. Pillet, and S. Haroche, “Observationof near-infrared dicke superradiance on cascading transitions inatomic sodium,” Phys. Rev. Lett. , 1035–1038 (1976).[6] Marlan O. Scully, “Collective lamb shift in single photon dickesuperradiance,” Phys. Rev. Lett. , 143601 (2009).[7] Z. Meir, O. Schwartz, E. Shahmoon, D. Oron, and R. Oz-eri, “Cooperative lamb shift in a mesoscopic atomic array,” Phys. Rev. Lett. , 193002 (2014).[8] J. Pellegrino, R. Bourgain, S. Jennewein, Y. R. P. Sor-tais, A. Browaeys, S. D. Jenkins, and J. Ruostekoski,“Observation of suppression of light scattering inducedby dipole-dipole interactions in a cold-atom ensemble,”Phys. Rev. Lett. , 133602 (2014).[9] Antoine Browaeys, Daniel Barredo, and Thierry La-haye, “Experimental investigations of dipole–dipoleinteractions between a few rydberg atoms,”Journal of Physics B: Atomic, Molecular and Optical Physics , 152001 (2016).[10] S. L. Bromley, B. Zhu, M. Bishof, X. Zhang, T. Bothwell,J. Schachenmayer, T. L. Nicholson, R. Kaiser, S. F. Yelin,M. D. Lukin, A. M. Rey, and J. Ye, “Collective atomicscattering and motional effects in a dense coherent medium,”Nature Communications , 11039 (2016).[11] David Plankensteiner, Christian Sommer, HelmutRitsch, and Claudiu Genes, “Cavity antiresonance spectroscopy of dipole coupled subradiant arrays,”Phys. Rev. Lett. , 093601 (2017).[12] Stephan Jennewein, Ludovic Brossard, Yvan R. P. Sortais, An-toine Browaeys, Patrick Cheinet, Jacques Robert, and PierrePillet, “Coherent scattering of near-resonant light by a dense,microscopic cloud of cold two-level atoms: Experiment versustheory,” Phys. Rev. A , 053816 (2018).[13] P.-O. Guimond, A. Grankin, D. V. Vasilyev, B. Vermersch,and P. Zoller, “Subradiant bell states in distant atomic arrays,”Phys. Rev. Lett. , 093601 (2019).[14] P. Solano, P. Barberis-Blostein, F. K. Fatemi, L. A.Orozco, and S. L. Rolston, “Super-radiance revealsinfinite-range dipole interactions through a nanofiber,”Nature Communications , 1857 (2017).[15] A. Grankin, P. O. Guimond, D. V. Vasilyev, B. Vermersch, andP. Zoller, “Free-space photonic quantum link and chiral quan-tum optics,” Phys. Rev. A , 043825 (2018).[16] L. F. Buchmann, K. Mølmer, and D. Petrosyan, “Creation andtransfer of nonclassical states of motion using rydberg dressingof atoms in a lattice,” Phys. Rev. A , 013403 (2017).[17] Yang Wang, Xianli Zhang, Theodore A. Corcovilos, Aish-warya Kumar, and David S. Weiss, “Coherent address-ing of individual neutral atoms in a 3d optical lattice,”Phys. Rev. Lett. , 043003 (2015).[18] Robert J. Bettles, Simon A. Gardiner, and Charles S.Adams, “Enhanced optical cross section via col-lective coupling of atomic dipoles in a 2d array,”Phys. Rev. Lett. , 103602 (2016).[19] Ephraim Shahmoon, Dominik S. Wild, Mikhail D.Lukin, and Susanne F. Yelin, “Cooperative resonancesin light scattering from two-dimensional atomic arrays,”Phys. Rev. Lett. , 113601 (2017).[20] J. P. Clemens, L. Horvath, B. C. Sanders, and H. J. Carmichael,“Collective spontaneous emission from a line of atoms,”Phys. Rev. A , 023809 (2003).[21] Bo Yan, Steven A. Moses, Bryce Gadway, Jacob P. Covey,Kaden R. A. Hazzard, Ana Maria Rey, Deborah S. Jin, andJun Ye, “Observation of dipolar spin-exchange interactions withlattice-confined polar molecules,” Nature , 521–525 (2013).[22] Vahagn Mkhitaryan, Lijun Meng, Andrea Marini,and F. Javier García de Abajo, “Lasing and am-plification from two-dimensional atom arrays,”Phys. Rev. Lett. , 163602 (2018).[23] Stewart D. Jenkins, Janne Ruostekoski, Nikitas Papasimakis,Salvatore Savo, and Nikolay I. Zheludev, “Many-body subradi-ant excitations in metamaterial arrays: Experiment and theory,”Phys. Rev. Lett. , 053901 (2017).[24] R. Bekenstein, I. Pikovski, H. Pichler, E. Shahmoon, S. F. Yelin,and M. D. Lukin, “Quantum metasurfaces with atom arrays,”Nature Physics , 676–681 (2020).[25] Daniel Barredo, Sylvain de Léséleuc, Vincent Lien-hard, Thierry Lahaye, and Antoine Browaeys, “Anatom-by-atom assembler of defect-free arbitrary two-dimensional atomic arrays,” Science , 1021–1023 (2016),https://science.sciencemag.org/content/354/6315/1021.full.pdf.[26] Manuel Endres, Hannes Bernien, Alexander Keesling, HarryLevine, Eric R. Anschuetz, Alexandre Krajenbrink, CrystalSenko, Vladan Vuletic, Markus Greiner, and Mikhail D. Lukin,“Atom-by-atom assembly of defect-free one-dimensionalcold atom arrays,” Science , 1024–1027 (2016),https://science.sciencemag.org/content/354/6315/1024.full.pdf.[27] Brian J. Lester, Niclas Luick, Adam M. Kaufman,Collin M. Reynolds, and Cindy A. Regal, “Rapid pro-duction of uniformly filled arrays of neutral atoms,” Phys. Rev. Lett. , 073003 (2015).[28] T. Xia, M. Lichtman, K. Maller, A. W. Carr, M. J. Piotrow-icz, L. Isenhower, and M. Saffman, “Randomized benchmark-ing of single-qubit gates in a 2d array of neutral-atom qubits,”Phys. Rev. Lett. , 100503 (2015).[29] F. Nogrette, H. Labuhn, S. Ravets, D. Barredo, L. Béguin,A. Vernier, T. Lahaye, and A. Browaeys, “Single-atom trap-ping in holographic 2d arrays of microtraps with arbitrary ge-ometries,” Phys. Rev. X , 021034 (2014).[30] A. González-Tudela, C.-L. Hung, D. E. Chang, J. I. Cirac,and H. J. Kimble, “Subwavelength vacuum lattices andatom–atom interactions in two-dimensional photonic crystals,”Nature Photonics , 320–325 (2015).[31] Jun Rui, David Wei, Antonio Rubio-Abadal, Simon Hollerith,Johannes Zeiher, Dan M. Stamper-Kurn, Christian Gross, andImmanuel Bloch, “A subradiant optical mirror formed by a sin-gle structured atomic layer,” Nature , 369–374 (2020).[32] G. Facchinetti, S. D. Jenkins, and J. Ruostekoski, “Stor-ing light with subradiant correlations in arrays of atoms,”Phys. Rev. Lett. , 243601 (2016).[33] M T Manzoni, M Moreno-Cardoner, A Asenjo-Garcia,J V Porto, A V Gorshkov, and D E Chang, “Optimiza-tion of photon storage fidelity in ordered atomic arrays,”New Journal of Physics , 083048 (2018).[34] A. Asenjo-Garcia, M. Moreno-Cardoner, A. Albrecht, H. J.Kimble, and D. E. Chang, “Exponential improvement in pho-ton storage fidelities using subradiance and “selective radiance”in atomic arrays,” Phys. Rev. X , 031024 (2017).[35] K. M. Maller, M. T. Lichtman, T. Xia, Y. Sun, M. J. Piotrowicz,A. W. Carr, L. Isenhower, and M. Saffman, “Rydberg-blockadecontrolled-not gate and entanglement in a two-dimensional ar-ray of neutral-atom qubits,” Phys. Rev. A , 022336 (2015).[36] Yang Wang, Aishwarya Kumar, Tsung-Yao Wu, and David S.Weiss, “Single-qubit gates based on targeted phase shiftsin a 3d neutral atom array,” Science , 1562–1565 (2016),https://science.sciencemag.org/content/352/6293/1562.full.pdf.[37] Ephraim Shahmoon, Mikhail D. Lukin, and Susanne F. Yelin,“Quantum optomechanics of a two-dimensional atomic array,”Phys. Rev. A , 063833 (2020).[38] M. M. Cola, D. Bigerni, and N. Piovella, “Recoil-induced subradiance in an ultracold atomic gas,”Phys. Rev. A , 053622 (2009).[39] J. Guo, P. R. Berman, B. Dubetsky, and G. Gryn-berg, “Recoil-induced resonances in nonlinear spectroscopy,”Phys. Rev. A , 1426–1437 (1992).[40] J. Javaloyes, M. Perrin, G. L. Lippi, and A. Politi, “Self-generated cooperative light emission induced by atomic recoil,”Phys. Rev. A , 023405 (2004).[41] P. Wolf, S. C. Schuster, D. Schmidt, S. Slama, and C. Zim-mermann, “Observation of subradiant atomic momentum stateswith bose-einstein condensates in a recoil resolving optical ringresonator,” Phys. Rev. Lett. , 173602 (2018).[42] S. Ali Hassani Gangaraj, George W. Hanson, Mauro An-tezza, and Mário G. Silveirinha, “Spontaneous lateralatomic recoil force close to a photonic topological material,”Phys. Rev. B , 201108 (2018).[43] Yu. A. Avetisyan, V. A. Malyshev, and E. D.Trifonov, “Photon recoil in light scattering bya bose–einstein condensate of a dilute gas,”Journal of Experimental and Theoretical Physics , 380–386 (2020).[44] François Damanet, Daniel Braun, and John Martin, “Masterequation for collective spontaneous emission with quantizedatomic motion,” Phys. Rev. A , 022124 (2016).[45] F. Robicheaux and Shihua Huang, “Atom recoil dur- ing coherent light scattering from many atoms,” Phys. Rev. A99