Quantum nonlinear metasurfaces from dual arrays of ultracold atoms

Atoms in a sub-wavelength lattices have remarkable optical properties that have become of high scientific and technological significance. Here, we show how the coupling of light to more than a single atomic array can expand these perspectives into the domain of quantum nonlinear optics. While a single array transmits and reflects light in a largely linear fashion, the combination of two arrays is found to induce strong photon-photon interactions that can convert an incoming classical beam into highly antibunched light. Such quantum metasurfaces open up new possibilities for coherently generating and manipulating nonclassical light, from optical quantum information processing to exploring quantum many-body phenomena in two-dimensional systems of strongly interacting photons.

Advances in controlling cold atomic ensembles at the single-particle level [1] have enabled the development of novel light-matter interfaces. Recent investigations  have explored various approaches towards strong and coherent coupling to propagating photons, using nanoscale optical waveguides, photonic crystals or regularly arranged atoms in free space. In particular, extended two-dimensional lattices of atoms suggest themselves as optical metasurfaces that can be designed and engineered on sub-wavelength scales at the level of individual quantum emitters. Experiments have demonstrated the strong coherent coupling to subradiant collective excitations of atoms in optical lattices [20], which can enable the near lossless interfacing with a single mode of freely propagating light fields [15][16][17][18]. This makes it possible to explore the functionalities of optical metasurfaces [24][25][26][27][28], which, e.g., offers new possibilities for highly coherent wavefront engineering [19]. While dense emitter arrangements with an overall sub-wavelength dimension can exhibit high nonlinearities akin to a single atom [29,30], the extended geometry of two-dimensional surfaces renders such systems intrinsically linear. In fact, the simultaneous photon interaction with a large number of atoms, necessary to achieve strong collective coupling, diminishes the otherwise strong optical nonlinearity of each individual quantum emitter, thereby, restricting their collective optical response to the domain of linear optics.
Here, we describe and analyse how strong optical nonlinearities can be obtained in metasurfaces composed of sub-wavelength atomic lattices. Specifically, we will discuss how combining two atomic arrays, which separately are only weakly nonlinear, can greatly enhance their combined optical nonlinearity to a degree that acts on the level of individual photons. The quantum optical nonlinearity in this system arises from narrow transmission resonances around which photons are strongly confined between the two arrays, forming a high finesse optical resonator. We show that such dual-array settings can effectively transform an incident classical beam into strongly correlated states of light, while the statistics of incident  5)]. (c) In the vicinity of such narrow transmission resonances, the optical response becomes highly nonlinear and generates effective photon-photon interactions that can transform an incident coherent beam into highly nonclassical light, as demonstrated by the depicted two-photon correlation function g (2) (t) of the transmitted light (solid line). In contrast, light reflected from a single identical array remains largely uncorrelated (dashed line). The chosen lattice spacing is a = 0.6λ. The transmission in (a) is obtained for infinitely extended arrays, while the calculations in (c) are performed for finite arrays with 9 × 9 atoms driven by a Gaussian laser beam with a waist of w = 1.5λ, ∆ = 0.472γ and L = 1.55λ.
photons is left virtually unchanged by each individual array (cf. Fig. 1c). The nonlinearity can be traced back 2 to the emergence of an effective photon-photon interaction in the two-dimensional plane of the surface, which is shown to generate in-plane collisions between individual photons. This suggests a promising approach to employing atomic lattices as quantum nonlinear metasurfaces for generating and manipulating nonclassical states of light, and exploring quantum many-body physics with photons [31,32].
Let us first consider a single two-dimensional square lattice of atoms with a lattice spacing a, and which we assume to be infinitely extended in the x-y plane. We focus the discussion on two-level systems that are resonantly driven by a coherent cw-field with an electric field amplitude E in and spatial mode function f (R) [33]. At the small lattice spacings considered here, exchange of photons across the array generates strong atomic interactions that can be efficiently described within input-output theory by integrating out the photonic degrees of freedom and using a Born-Markov approximation [34,35]. This yields a master equation ∂ tρ = −i[Ĥ,ρ] + L[ρ] (with = 1) for the density matrix,ρ, of the atomic lattice, where the Hamiltonian and Lindblad operator describe the exchange of excitations and corresponding collective decay processes due to the photon-mediated dipole-dipole interactions between the atoms [36]. Here, σ n = |g n e n | denotes the transition operator between the ground state, |g n , and excited state |e n of an atom at position R n in the lattice. The interaction coefficients J nm and decay rates Γ nm for two atoms at positions R n and R m are determined by the Green's function tensor of the free-space electromagnetic field [36,37]. The atomic transition is driven by the incident light with a frequency detuning ∆ and a single-atom Rabi frequency Ω n = dE in f (R n ), which is determined by the drivingfield amplitude and the transition dipole moment d of the two-level emitters. From the solution forρ, one can reconstruct the electromagnetic field generated by the driven atomic dipoles. Choosing f (R n ) as the detection mode for the transmitted light, one obtains [36,37] for the electric-field amplitude,Ê, of the detected photons, where k = 2π/λ is the wave number of the incident light, λ denotes its wavelength, γ = Γ nn is the decay rate of the individual atoms, and η = d 2 r|f (R)| 2 with r = (x, y).
For weak plane-wave driving with f ∼ e ikz , Eqs. (1a), (1b) and (2) yield simple expressions for the transmission and reflection spectra [17] that feature a Lorentzian resonance at the collective Lamb shift∆ = − n =0 J n0 with a widthΓ = 3πγ/k 2 a 2 . On resonance, the single atomic layer, thus, reflects incoming photons with unit efficiency and no losses from the incident mode, |t| 2 + |r| 2 = 1. Such high reflectivities are attainable already for remarkably small systems [15]. For example, a 9×9 atomic array with a = 0.6λ can reflect into the incident mode of a focused Gaussian beam with a waist of w = 1.5λ with a large reflection amplitude of |r| = 0.998. On the other hand, the nonlinear response is very small as the beam still covers a sizeable number of atoms, which substantially diminishes saturation effects. This is seen directly from the second order correlation function g (2) (t) of the reflected light, shown by the dashed line in Fig. 1b. Here, one finds only a marginal suppression of simultaneous two-photon reflection, indicating that the reflected light largely retains the classical coherent-state nature of the incident beam (g (2) ∼ 1).
This situation changes dramatically as we add a second atomic array. Fig. 1a shows the transmission coefficient |T | 2 for a dual-array configuration of two parallel atomic lattices as a function of the detuning ∆ and the distance L between the two arrays. The calculations reveal a series of narrow transmission resonances that extends towards large values of L and a pair of sharp reflection resonances at small array distances. Both regimes can be traced back to the photon-mediated interactions between the two arrays.
For small distances L, atoms in different arrays interact via a dipole-dipole coupling that scales as J L ≈ −3γ/2(kL) 3 [37]. We can, thus, define symmetric and antisymmetric superposition states, |± n , of a single excitation that is symmetrically (|+ n ) and antisymmetrically (|− n ) shared between two adjacent atoms at a given lattice site n. The atomic interaction shifts their respective energies by ±J L . Therefore, the two atomic dimer states become energetically isolated for small L and separately generate reflection resonances at the collective energies∆ ± ∼ ±L −3 , as indicated by the dashed lines in Fig. 1a [37]. Their respective widths are given byΓ ± =Γ[1 ± cos(kL)], such that one finds an ultranarrow reflection resonance withΓ − Γ , generated by an effective array of subradiant atomic dimer states as L decreases.
For larger values of L, the evanescent-field coupling vanishes and the interaction between the arrays is predominantly generated by propagating photons. The coupling strength J L ≈ 3γ cos(kL)/kL, therefore, acquires an oscillating behaviour from the propagation phase and leads to an energy difference∆ + −∆ − ∼ sin(kL) that varies periodically with the array distance L. As both collective dimer states are excited by the parity-breaking incident field, their interference leads to a series of narrow transmission resonances, akin to electromagnetically induced transparency in three-level systems [37,38]. One can obtain a simple expression for the dual-array transmission amplitude [37] in terms of the reflection and transmission amplitudes, r and t, of the individual arrays. Substituting their explicit form, Eq. (3), we obtain the following condition for perfect transmission (|T | = 1) that defines the series of transmission resonances. Along these resonances, transmitted photons can acquire a substantial group delay with a delay time [37] τ =2 that diverges for resonant detunings ∆ =∆. The width of the transmission resonances decreases as ∼ 1/τ around these values and, therefore, vanishes at kL = nπ for integer n > 0. In between the resonances (kL = (n + 1/2)π), the system features high reflection and behaves largely linear, which can be used to store several delocalized excitations across distant arrays at L λ [39]. Eq. (4) describes the transmission of a Fabry-Pérot resonator composed of two identical mirrors with respective reflection and transmission amplitudes r and t. Here, however, the single-particle saturation of each emitter within the atomic mirrors together with the narrow subradiant transmission resonances can enhance optical nonlinearities and generate exceedingly strong photonphoton interactions.
We have studied the signatures of such photon-photon interactions via quantum trajectory wave function simulations [40] of the atomic master equation with Eqs. (1a) and (1b) for finite arrays. Working with finite arrays and focused driving beams, generally entails photon losses that tend to broaden the otherwise ultra-narrow transmission resonances. These effects can be mitigated through a proper choice of the atomic lattices, matching the wavefront profile of the incident beam [37,39]. In fact, already rather small lattices of 9×9 atoms permit to generate narrow transmission resonances with linewidths of ∼ 10 −2 γ and high peak transmission of |T | ∼ 0.98. Fig. 1c shows the calculated second order correlation function The characteristic time dependence of the two-photon correlation function, g (2) (t), can be understood by the dynamics of the short-lived (|+ ) and long-lived (|− ) single atomic excitation that is symmetrically (|+ ) and antisymmetrically (|− ) delocalized between the two arrays. (b) Following the detection of a photon, the subsequent population dynamics, |c+| 2 (blue) and |c−| 2 (red), of these two states agrees with the characteristic time dependence of g (2) (t) shown in panel (c) (see text for more details). The parameters are the same as in Fig. 1c. of the transmitted light for an incident cw field with a Gaussian transverse beam profile whose waist is centred right in between the two arrays. We consider long times t → ∞, such that Eq. (7) yields the temporal photonphoton correlation in the steady state. Its dependence on the time delay, t, between consecutively detected photons indicates the generation of highly nonclassical light. Interestingly, one finds a rapid initial drop of the twophoton correlation function to small values g (2) ∼ 0, which extends over a broad range of delay times between two transmitted photons. These characteristic temporal correlations can be understood as follows. Let us denote the steady state of the two arrays as |ψ . Detection of a transmitted photon in the steady state then projects this state onto |ψ =Ê|ψ / ψ|Ê †Ê |ψ . The correlation function can thus be obtained as from the time evolved state |ψ(t +t) following detection min = mintg (2) (t). The color bars in panel (e) show the long timescale on which the correlation functions eventually approach unity, which matches the photon delay time, or photon confinement time, τ . Its asymptotic value for infinite arrays is given by Eq. (6), while the dashed line in panel (e) has been obtained numerically for the 9 × 9 arrays.
of a photon at time t . For weak driving this state is predominantly determined by the collective ground state, |0 , and the single-excitation manifold. It can hence be expressed as a superposition, |ψ = c 0 |0 + c + |+ + c − |− , of the collective superradiant (|+ ) and subradiant (|− ) states of a single atomic excitation that is shared (anti)symmetrically across the two arrays, as discussed above. Indeed, the time evolution of the populations |c ± | 2 resembles the temporal photon correlations, see Fig. 2. The initial drop of g (2) thus reflects the fast decay of the superradiant excitation on a short timescale τ + , leaving the dual array depleted of excitations. Its subsequent slow rise, on the other hand, can be traced back to the repopulation of the long-lived subradiant cavity state, |− , on a long timescale τ − given by the inverse width of the transmission resonance, which corresponds to the photon delay time, discussed above (cf. Fig. 3e and see [37]).
Consequently, we expect stronger effects of photonphoton interactions around more narrow transmission lines. This is demonstrated in Fig. 3, where we show the two-photon correlation function g (2) while scanning the frequency detuning, ∆, and the array distance, L, along the transmission maximum of one of the resonances, as illustrated in Fig. 3b. Indeed, we find that the minimum value, g min = min t g (2) (t), of the two-photon corre-lation function decreases as the resonance becomes more narrow and the delay time increases. Owing to the finite size of the array, the linear transmission decreases as we approach this regime by varying L (Fig. 3a). The associated losses tend to broaden the transmission lines and, as shown in Fig. 3e, lead to a maximum delay time of τ ∼ 1000/γ, instead of the otherwise diverging behaviour, discussed above for infinite arrays and plane wave driving [cf. Eq. (6)]. Larger arrays yield longer photon confinement times, τ , for a given transmission maximum, which, therefore, enhances both the temporal extend and the strength of photon-photon correlations.
Remarkably, however, one can reach large single-photon nonlinearities and highly nonclassical light with g (2) min ∼ 0 under conditions of high photon transmission already for moderate system sizes, which are achieved in ongoing optical lattice experiments [1,20].
We can gain further insights into the generated nonlinearity by considering the two-photon momentum densitỹ is the transverse Fourier transform of the electric field operator of the transmitted light [37]. Fig. 4 shows the in-plane momentum distribution of the two photons. In Fig. 4a and b, we have fixed the transverse momentum k 1 of one photon at the value indicated by the white star and show the equal-time momentum distribution of the other (t = 0). The distribution of k 2 is sharply peaked around −k 1 indicating that the optical nonlinearity can indeed be understood in terms of effective photon-photon collisions that preserve the total transverse momentum k 1 + k 2 0 of the incident beam. Note that this approximate conservation law already emerges for moderate system sizes of 9 × 9 atoms. Fig. 4c displays the equal-time momentum correlations between the y-components of both photons for k 1,x = k 2,x = 0. The sharp maximum around k 2,y = −k 1,y once more reflects the conservation of total momentum, while the variation along the anti-diagonal results from the momentum dependence of the scattering process, i.e., the momentum dependence of the effective photon-photon interaction. The signal in the corners are due to Bragg scattering of the interacting photons. Fig. 4d shows the same momentum density as in Fig. 4c, but for a finite delay time t = 10/γ. In this case we observe practically no correlations in the transverse momenta and the signal corresponds to the transverse mode of the incident beam. This behaviour can be readily understood from the conditioned dynamics of the collective single excitation, discussed above (cf. Fig. 2). At small delay times, the two-photon signal naturally stems from the superradiant excitation following detection of the first photon. While the narrow linear transmission line arises from long-lived subradiant excitations, the superradiant mode can only be populated by the interaction The large density near k1 + k2 0 shows the photons have scattered off each other with conserved momentum. (c) and (d):ρ(k1, k2, t) for k1,x = k2,x = 0 with a time t = 0 and t = 10/γ, respectively, between the detection of the two photons. The time t = 10/γ corresponds to g (2) (t) ∼ 0. Again, the anti-diagonal bar in panel (c) indicates a momentumconserving scattering among photons. After a time t = 10/γ this correlation has vanished (see text for details). The colour scheme is chosen such that each panel simply shows the qualitative structure. The system parameters are the same as in Fig. 1c. between multiple excitations. One therefore finds strong momentum correlations in the transverse two-photon signal, as shown in Fig. 4a-c for a vanishing delay time t = 0. For longer delay times beyond the superradiant lifetime, the superradiant state is depleted and the second photon predominantly stems from the subradiant mode that is repopulated by the continuous optical driving of the arrays. The two detected photons, thus, originate from cascaded excitation and emission processes such that we find negligible transverse-momentum correlations albeit strong temporal correlations between the two photons. Despite the local nature of the underlying nonlinearity of each individual atom, this mechanism makes it possible to generate strongly correlated photons without significant transverse mode mixing in the plane of the atomic arrays.
Sub-wavelength lattices of ultracold atoms provide a promising platform for the coherent manipulation of optical fields, and in this work we have described how these perspectives can be extended into the domain of quantum nonlinear optics by using two atomic arrays. While large optical nonlinearities can also be generated in atomic ensembles via strong interactions between highlying atomic Rydberg states [41-47], the present setting yields an alternative mechanism beyond the physics of excitation-blockaded superatoms. This is made possible by ultra-narrow transmission resonances that emerge from interference between collective superradiant and subradiant states of the dual array, bearing analogies to electromagnetically induced transparency [38] and the physics of Fano resonators [48]. We have demonstrated that the single-photon saturation of each individual atom can generate a strong and finite-ranged effective interaction between photons. Such emerging photon-photon interactions suggest a number of question for future work. We have identified a regime in which strong temporal photon correlations emerge under conditions of very low transverse-mode mixing, thus generating large nonlinearities for freely propagating single photonic modes at greatly suppressed losses. This motivates future explorations of applications as nonlinear quantum optical elements to generating and processing photonic quantum states [49-51], or to study the physics of propagating multi-photon quantum states [52-56] through a many of such nonlinear elements.
In this work, we have mainly focused on analysing temporal correlations of photons in single transverse modes, drawing analogies to waveguide QED settings [23]. In addition, however, the multi-mode physics of large planar arrays should yield an interesting framework for exploring the many-body physics of multiple photons in the two-dimensional plane of the dual-array resonator, and motivates future work on the potential formation and nonlinear dynamics of effective cavity polaritons [31]. Hereby, our results indicate that this should make it possible to reach the quantum regime of strong interactions between individual polaritons.
We thank Kristian Knakkergaard, Aurélien Dantan, and Peter Rabl for valuable discussions. This work was supported by the Carlsberg Foundation through the "Semper Ardens" Research Project QCooL, and by the DNRF through the Center of Excellence "CCQ" (Grant agreement no.: DNRF156

TRANSMISSION RESONANCES OF DUAL ATOMIC ARRAYS
The Greens function tensor G(R) of the free-space electromagnetic field has cartesian elements given by where R = |R| and R i (i = 1, 2, 3) denotes the cartesian component of the three-dimensional distance vector between a given point in space and a radiating dipole. The Greens function yields the coefficients J nm and Γ nm in Eq. (1) of the main text via J nm + iΓ nm = G nm , where and R n and R m denote the three-dimensional positions of two atoms in the arrays. Considering two parallel arrays, positioned in the x − y plane, we use the notation R n = (x n , y n , z n ) = (r n , z n ). As discussed in the main text, we focus on a single atomic transition driven by circularly polarized light, choosing e + = (1, i, 0) T / √ 2 as a specific polarization vector. Explicitly, one, therefore, obtains where R 2 nm = r 2 nm + 2 , and denotes the atomic separation along the z-axis. Hence, = 0 describes interactions between two atoms in the same array and > 0 corresponds to interactions between atoms in each of the two arrays. We employ Eq. (3) together with Eq. (1) of the main text for the numerical simulations of finite lattices with planar and non-planar geometries.
For the analytical estimates, discussed in the main text, we consider two parallel infinitely extended arrays, positioned at z = ±L, respectively. The Fourier transform of Eq. (3) in the x − y planẽ then defines contributions to the effective level shifts and decay rates from interactions between atoms in the same array ( = 0) and in different arrays ( = L). Here, the sum runs over the reciprocal lattice sites, q = (2π/a)m with m = (m x , m y ) and m x , m y ∈ Z. One can see from this expression that a < λ/2 ensures that the imaginary part ofG has only one contribution from the q = 0 term, since |k ⊥ | < k. This implies that the transverse momentum of the incident light is conserved by the interaction with the array, eliminating any losses from the incident mode due to photon scattering. If k ⊥ = 0, this condition is already met for a < λ. Divergences in the above expression originate from the single-atom Lamb shift and can be accounted for the in the detuning ∆. The expression can be evaluated numerically using a regularization procedure described in [1]. For plane wave driving at normal incidence (k ⊥ = 0), one obtains the simple result Γ =Γ cos(k ). (5b)

arXiv:2201.06544v3 [quant-ph] 30 Mar 2023
The intra-layer terms are given bỹ The inter-layer contribution,∆ L , to the collective shift consists of two terms. The first describes the energy shift due to exchange of photons that propagate between the two arrays and therefore oscillates at the wavelength, λ = 2π/k, of the transition resonance. The second term, decays rapidly with L and describes the direct near-field interaction between the two arrays. Either term dominates in the respective limits of large and small array distances, L, as discussed in the main text and as we shall describe in more detail below. To this end, let us denote the atomic transition operators for the first and second array with respect to the field incidence byσ (1) n andσ (2) n , respectively. One can then introduce (anti)symmetric superpositionsσ ±,n = (σ (1) n ± σ (2) n )/ √ 2 at a given lattice site r n of the dual array, which define the symmetric and antisymmetric dimer states |± n =σ † ±,n |0 , discussed in the main text. For plane wave driving with k ⊥ = 0 we can write the Rabi frequencies for each array as Ω where the interaction coefficients are obtained from J nm ( ) + iΓ nm ( ) = G nm ( ) and Eq. (3). The parity symmetry of the light-induced interactions implies that there is no direct coupling between symmetric and antisymmetric states, such that only the external driving field can transfer populations between the two parity eigenstates at a given lattice site. The Hamiltonian can be rewritten aŝ whereΩ + = √ 2Ω cos(kL/2),Ω − = √ 2Ω sin(kL/2),Ω = √ N Ω,∆ ± =∆ ±∆ L , andΓ ± =Γ ±Γ L =Γ[1 ± cos(kL)], as given in the main text. The operatorsσ ± = N −1/2 nσ ±,n generate the collective (anti)symmetric states discussed in the main text. Using a truncated basis with not more than a single excitation, the system can be described by a simple non-Hermitian Hamiltonian From the corresponding steady state of the driven atoms, one can determine the transmitted field from Eq. (2) of the main text asÊ (a) (b) to obtain the transmission amplitude of the dual array The transmission obtained from this expression is shown in Fig. 1a of the main text. Below we discuss the limiting behaviour of this simple expression for small and large values of L, respectively.

Small-L limit
For small values of L < a, the major contribution to the interaction stems from atomic pairs at identical sites in different arrays. We thus obtain from Eq. (3), that the energy shift∆ L ≈ 3γ/2(kL) 3 predominantly arises from the direct dipole-dipole interaction between adjacent atoms from each array. In Fig. 1, we compare this expression to the numerical Fourier transform Eq. (5a) and find good agreement for small L. This implies that the two resonances at∆ ± ≈∆±3γ/2(kL) 3 in Eq. (12) are energetically well separated, such that we find from Eq. (12) two individual transmission resonances as discussed in the main text and as shown in Fig. 1 of the article. As L decreases the width of the subradiant resonance vanishes asΓ − ≈Γ 2 (kL) 2 , while we find a superradiant rateΓ + ≈ 2Γ for the resonance at positive frequency detunings. These approximate results agree well with the numerical findings shown in Fig. 1a of the main text.

Large-L limit
For large distances between the arrays, we can neglect their evanescent field coupling, such that∆ L ≈Γ sin(kL), according to Eq. (5a), follows a simple sinusoidal oscillation (cf. Fig. 1). Since, the two reflection resonances in Eq. (12) now overlap significantly they both contribute to the linear response of the array. We can gain further insights by transforming to a new basis generated byB = cos(kL/2)σ + + i sin(kL/2)σ − andD = − sin(kL/2)σ + + i cos(kL/2)σ − . The state |b =B † |0 (15) corresponds to a bright state with maximum coupling to the incident driving field, while represents a dark state with respect to the external driving field that is only coupled to |b via the atomic interactions with a coupling strength∆ L sin(kL). From Eqs. (8a) and (8b) we obtain Here, we have usedΓ L =Γ cos(kL) and the large-L limit∆ L =Γ sin(kL). This describes an effective ladder-type three-level system that is illustrated in Fig. 2. While the transmitted field is determined by the bright state only, the interference with the dark-state coupling can nevertheless generate transmission resonances akin to electromagnetically induced transparency in three-level systems [2]. From the steady state one finds that the interference between the two excitation pathways leads to a simple expression if∆ − ∆ =Γ tan kL, which corresponds to Eq. (5) in the main text. Under this condition we have The interference between the two excitation pathways, thus, gives rise to a transmission resonance, where the medium becomes fully transparent, |T | 2 = 1. While this is indeed similar to electromagnetically induced transparency in threelevel systems with a ladder-type coupling scheme, the atomic polarization, ∼ B , only vanishes for kL = ±nπ/2 (for integer values of n). In the present case, we still find obtain perfect transparency for finite values of B such that the transmitted field acquires a tunable phase of π − 2kL at no photon losses. Moreover, the width of this transmission resonance depends on ∆ and vanishes upon approaching ∆ =∆. As we shall discuss below, this gives rise to a long confinement time of incident photons within the formed resonator and consequently strong photon-photon interactions. The transmission T = Ê /E in obtained from Eqs. (18) and (19) can be expressed as (23) in terms of the transmission and reflection amplitudes of each individual array, given by Eq. (3) of the main text. Expectedly, this expression coincides with the transmission coefficient of a Fabry-Pérot resonator composed of two identical mirrors with transmission and reflection amplitude t and r placed at a distance L.

< l a t e x i t s h a 1 _ b a s e 6 4 = " l 4 k Y N g h 7 J c v j 2 / t 5 3 x M Y C o 7 T n D U = " > A A A B 8 H i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 1 D J g Y x n B J E p y h L 3 N X L J k d + / Y 3 R N C z K + w s V D E 1 p 9 j 5 7 9 x k 1 y h i Q 8 G H u / N M D M v S g U 3 1 v e / v c L K 6 t r 6 R n G z t L W 9 s 7 t X 3 j 9 o m i T T D B s s E Y m + j 6 h B w R U 2 L L c C 7 1 O N V E Y C W 9 H w e u q 3 H l E b n q g 7 O 0 o x l L S v e M w Z t U 5 6 e P I 7 m q q + w G 6 5 4 l f 9 G c g y C X J S g R z 1 b v m r 0 0 t Y J l F Z J q g x 7 c B P b T i m 2 n I m c F L q Z A Z T y o a 0 j 2 1 H F Z V o w v H s 4 A k 5 c U q P x I l 2 p S y Z q b 8 n x l Q a M 5 K R 6 5 T U D s y i N x X / 8 9 q Z j a / C M V d p Z l G x + a I 4 E 8 Q m Z P o 9 6 X G N z I q R I 5 R p 7 m 4 l b E A 1 Z d Z l V H I h B I s v L 5 P m W T W 4 q A a 3 5 5 W a n 8 d R h C M 4 h l M I 4 B J q c A N 1 a A A D C c / w C m + e 9 l 6 8 d + 9 j 3 l r w 8 p l D + A P v 8 w e x t 5 B J < / l a t e x i t > |0i < l a t e x i t s h a 1 _ b a s e 6 4 = " k F C N n l e b W q U k B a M I + 0 v 2 U y 4 H w 6 A = " > A A A B 8 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B Z B E E o i o h 4 L X j x W s B / S h r L Z T t q l m 0 3 Y 3 Q g l 9 l d 4 8 a C I V 3 + O N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Q S K 4 N q 7 7 7 R R W V t f W N 4 q b p a 3 t n d 2 9 8 v 5 B U 8 e p Y t h g s Y h V O 6 A a B Z f Y M N w I b C c K a R Q I b A W j m 6 n f e k S l e S z v z T h B P 6 I D y U P O q L H S w 9 N Z V 1 E 5 E N g r V 9 y q O w N Z J l 5 O K p C j 3 i t / d f s x S y O U h g m q d c d z E + N n V B n O B E 5 K 3 V R j Q t m I D r B j q a Q R a j + b H T w h J 1 b p k z B W t q Q h M / X 3 R E Y j r c d R Y D s j a o Z 6 0 Z u K / 3 m d 1 I T X f s Z l k h q U b L 4 o T A U x M Z l + T / p c I T N i b A l l i t t b C R t S R Z m x G Z V s C N 7 i y 8 u k e V 7 1 L q v e 3 U W l 5 u Z x F O E I j u E U P L i C G t x C H R r A I I J n e I U 3 R z k v z r v z M W 8 t O P n M I f y B 8 / k D q g C Q R A = = < / l a t e x i t >
|+i < l a t e x i t s h a 1 _ b a s e 6 4 = " x 3 M o Q i b J l q 8 o o D G s c J K S a u w Y z n Y = " > A A A B 8 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B i y U R U Y 8 F L x 4 r 2 A 9 p Q 9 l s J + 3 S z S b s b o Q S + y u 8 e F D E q z / H m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g o r q 2 v r G 8 X N 0 t b 2 z u 5 e e f + g q e N U M W y w W M S q H V C N g k t s G G 4 E t h O F N A o E t o L R z d R v P a L S P J b 3 Z p y g H 9 G B 5 C F n 1 F j p 4 e m s q 6 g c C O y V K 2 7 V n Y E s E y 8 n F c h R 7 5 W / u v 2 Y p R F K w w T V u u O 5 i f E z q g x n A i e l b q o x o W x E B 9 i x V N I I t Z / N D p 6 Q E 6 v 0 S R g r W 9 K Q m f p 7 I q O R 1 u M o s J 0 R N U O 9 6 E 3 F / 7 x O a s J r P + M y S Q 1 K N l 8 U p o K Y m E y / J 3 2 u k B k x t o Q y x e 2 t h A 2 p o s z Y j E o 2 B G / x 5 W X S P K 9 6 l 1 X v 7 q J S c / M 4 i n A E x 3 A K H l x B D W 6 h D g 1 g E M E z v M K b o 5 w X 5 9 3 5 m L c W n H z m E P 7 A + f w B r R a Q R g = = < / l a t e x i t > | i < l a t e x i t s h a 1 _ b a s e 6 4 = " y T L R + F s M W f w l Q 4 T / Z m A m c a J B O j s = " > A A A B 8 H i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 1 D J g Y x n B J E p y h L 3 N X L J k d + / Y 3 R N C z K + w s V D E 1 p 9 j 5 7 9 x k 1 y h i Q 8 G H u / N M D M v S g U 3 1 v e / v c L K 6 t r 6 R n G z t L W 9 s 7 t X 3 j 9 o m i T T D B s s E Y m + j 6 h B w R U 2 L L c C 7 1 O N V E Y C W 9 H w e u q 3 H l E b n q g 7 O 0 o x l L S v e M w Z t U 5 6 e I o 6 m q q + w G 6 5 4 l f 9 G c g y C X J S g R z 1 b v m r 0 0 t Y J l F Z J q g x 7 c B P b T i m 2 n I m c F L q Z A Z T y o a 0 j 2 1 H F Z V o w v H s 4 A k 5 c U q P x I l 2 p S y Z q b 8 n x l Q a M 5 K R 6 5 T U D s y i N x X / 8 9 q Z j a / C M V d p Z l G x + a I 4 E 8 Q m Z P o 9 6 X G N z I q R I 5 R p 7 m 4 l b E A 1 Z d Z l V H I h B I s v L 5 P m W T W 4 q A a 3 5 5 W a n 8 d R h C M 4 h l M I 4 B J q c A N 1 a A A D C c / w C m + e 9 l 6 8 d + 9 j 3 l r w 8 p l D + A P v 8 w f + 3 Z B 7 < / l a t e x i t >

/ W L v P F x 3 R 8 T V 1 H u Q 6 v w I U = " > A A A B + X i c d V B N S 8 N A E N 3 U r 1 q / o h 6 9 L B Z B E E o i U j 0 W v H i z g v 2 A J o T N Z t I u 3 W z C 7 q Z Q Q v + J F w + K e P W f e P P f u G 0 j q O i D g c d 7 M 8 z M C z P O l H a c D 6 u y s r q 2 v l H d r G 1 t 7 + z u 2 f s H X Z X m k k K H p j y V / Z A o 4 E x A R z P N o Z 9 J I E n I o R e O r + d + b w J S s V T c 6
2 k G f k K G g s W M E m 2 k w L Y 9 z X g E h X e b w J D M g r P A r j s N Z w H 8 j T S X x C 2 V O i r R D u x 3 L 0 p p n o D Q l B O l B q 6 T a b 8 g U j P K Y V b z c g U Z o W M y h I G h g i S g / G J x + Q y f G C X C c S p N C Y 0 X 6 v e J g i R K T Z P Q d C Z E j 9 R v b y 7 + 5 Q 1 y H V / 5 B R N Z r k H Q 5 a I 4 5 1 i n e B 4 D j p g E q v n U E E I l M 7 d i O i K S U G 3 C q p k Q v j 7 F / 5 P u e c N t N t y 7 i 3 r L K e O o o i N 0 j E 6 R i y 5 R C 9 2 g N u o g i i b o A T 2 h Z 6 u w H q 0 X 6 3 X Z W r H K m U P 0

GROUP DELAY OF TRANSMITTED PHOTONS
The interaction of the incident light with the atoms in general leads to a group delay. One can determine the corresponding delay time, from the derivative of the transmission amplitude with respect to the frequency of the field according to [3] For a single array, this yields Expectedly, the delay time is simply given by the inverse decay rate on the reflection resonance. Equivalently, we can calculate the delay time for the dual-array setting. While the general expression is rather lengthy, one can find a simple form in the two limiting cases of small and large values of L. For small L we obtain on the two reflection resonances at ∆ = ∆ ± , which is consistent with the interpretation in terms of an effective single array composed of either subradiant or superradiant atomic dimers, as discussed in the main text. Here, |R| 2 = 1 − |T | 2 is the reflection coefficient of the dual array. For large array distances, where we can neglect the evanescent-field coupling between the arrays we obtain instead where δ = (∆ −∆)/Γ. Along the transmission maximum defined by Eq. (21), this yields the delay time given in Eq. (6) of the main text.

FINITE SIZE EFFECTS
In this section we provide additional details on the behaviour of finite arrays that will be relevant for future experiments and which have been used in our numerical simulations of the nonlinear optical response.

Gaussian driving beam
For our numerical simulations we choose the following normalized Gaussian mode for the in-coming field and the detection mode  where w(z) = w 1 + z zR 2 is the beam waist at z, z R = πw 2 /λ is the Rayleigh range, R(z) = z 1 + (z R /z) 2 is the z-dependent radius of the wavefront curvature of the beam. Moreover, ψ(z) = arctan (z/z R ) denotes the Gouy phase. In (k x , k y , z)-space (i.e. Fourier transformed with respect to x and y only) the mode takes on a more simple form where . This form of the field mode satisfies the paraxial wave equation. The paraxial approximation is valid if k z ≈ k 1 − k 2 ⊥ /2k 2 , i.e. if k 2 ⊥ /2k 2 = 1. In our simulations we use w = 1.5λ, for which this approximation is well justified, as illustrated in Fig. 3. In fact, for = 0.05, the mode function is already suppressed by an order of magnitude and the relative fraction p = ∞ k dk ⊥ k ⊥ |f (k ⊥ , 0)| 2 /2πw 2 = e −k 2 w 2 /2 of momenta that would give > 0.05 is p ∼ 0.01. Hence, the paraxial approximation gives reliable results for w ∼ 1.5λ and further improves rapidly for larger values of the beam waist.

Diffraction losses and curved arrays
Diffraction of the beams as photons propagate between the two mirrors causes small additional loss. As shown in [4], this can be mitigated by choosing slightly curved arrays that follow the profile of the optical wavefront. For our Gaussian beams, this gives a near-parabolic surface of the two arrays, such that the phase of the Gaussian field E in is constant across each array. Specifically, atoms at lattice sites r n = (x n , y n ) are now placed at z = ±z n determined by Such curved arrays improve the coherence properties for large array distances L, while diffraction losses generally become insignificant for smaller values of L ∼ λ, which makes it possible to realize narrow transmission resonances and strong photon-photon interactions with flat arrays, compatible with standard optical lattice geometries.

Photon scattering at finite transverse momenta
Further losses may also arise from the transverse intensity profile of the incident beam, which excites collective atomic excitations at finite transverse momenta with a Rabi frequencyΩ k⊥ . As a result, there is a finite population of collective excitations with varying shifts∆ k⊥ , such that the resonance condition for perfect transmission cannot be fulfilled for the entire beam. While this remains a minor effect for singlel atomic arrays [5], it can have significant consequences for the ultra-narrow resonances of the present dual-array setting. As shown in Fig. 4b, this effect can indeed alter the optical response in the resonant regime. A closer inspection of the transverse single-array dispersion ∆ k⊥ (cf. Eq. (4)), however, shows that the main broadening stems from sharp resonances of the collective frequency shift. These resonances originate from higher order Bragg scattering involving finite momenta q = (2π/a)m. Their effect can thus be suppressed by decreasing the lattice spacing a and eliminated entirely for a < λ/2. This is illustrated Fig. 4, which shows narrow transmission resonances for a beam waist of w = 1.5λ and a lattice spacing of a = 0.6λ, while they are diminished entirely for a = 0.8. The former parameters are used for all finite-array calculations discussed in the main text.

SUPER-AND SUBRADIANT SINGLE-EXCITATION STATES
As discussed in the main text, the characteristic time dependence of g (2) (t) is determined by the dynamics of the short lived and long-lived collective excitation. In the case of infinitely extended arrays and plane wave driving, these states correspond to the symmetric and antisymmetric superposition states |+ and |− , and the two characteristic timescales are set byΓ ± as introduced above. Following detection of a transmitted photon, the state of the system is projected onto |ψ = c 0 |0 + c + |+ + c − |− , and the subsequent time evolution of the coefficients can be determined from the non-Hermitian Hamiltonian which coincides with Eq. (10) above. The solution of the corresponding Schrödinger equation gives c ± (t) = c ± (0) + ig ± i∆ ± − γ ± e (i∆±−γ±)t − ig ± i∆ ± − γ ± (32) and yields the dynamics of the probabilities |c ± | 2 shown in Fig. 2a of the main text. In order to apply this picture to a finite system we numerically simulate its dynamics after photon detection in the absence of the driving field. At long times, this yields the long-lived state following the rapid decay of the superradiant states, which we then construct by orthogonality. Knowing the subradiant and superradiant states, we can then determine the dynamics of c ± (t) from the numerical solution of the exact master equation for the finite system, and obtain an excellent fit to Eq. (32). The extracted effective parameters of the subradiant and superradiant state are presented in Fig. 5 for the same parameters as in Fig. 3 of the main text.

DETERMINING THE LIGHT FIELD
Having solved for the dynamics of the atomic lattices, the photonic state can be readily reconstructed. The electric field operator of the light field can be obtained from the input-output relation [6] E(R, t) =Ê in (R, t) + 6πγ kd n G(R − R n ) · e +σn (t) (33) where the incident light is considered to be a coherent field, with a classical amplitude,Ê in = E in , given by E in = E in e + f (R). Choosing f (R) also as the detection mode, one can define the electric field operator E f (t) = η −1 d 2 r e † + ·Ê(R, t)f * (R) (34) of the transmitted photons (z L/2) in the detection mode f with d 2 r|f (R)| 2 = η. Assuming that all fields are paraxial, one then arrives at the simplified form of the input-output relation [5] given in Eq. (2) of the main text.
Finally, we consider the transverse Fourier transform of the field amplitude in the plane of the atomic arrays E(k ⊥ ) = d 2 r e † + ·Ê(R, t)e −ik⊥r e −ikzz , as used in Eq. (9) of the main text. The longitudinal wavenumber k z = k 2 − k 2 ⊥ is determined by the carrier frequency of the incident field, and we focus here on small transverse momenta k ⊥ < k for which k z remains real. The Fourier transformed field is readily obtained from Eq. (33) E(k ⊥ , t) = E in (t)g(k ⊥ ) + 6πγ kd n d 2 r e † + · G(R − R n ) · e + e −ikRσ n (t), = E in (t)g(k ⊥ ) + 6πγ kd e † + ·G(k ⊥ , z) · e +σk (t), whereG(k ⊥ , z) = d 2 r G(R)e −ikR andσ k = nσ n e −ikrn−ikzzn denote the transverse Fourier transforms of the Greens function tensor and the atomic transition operators, respectively. The Fourier transformed transverse mode functiong(k) of the incident field is given in Eq. (29). In the last line we have substituted the considered polarization vector e + = (1, i, 0) T / √ 2 and assumed z L/2, for which the transmitted field,Ẽ(k), is independent of z, due to the absence of evanescent field contributions. Eq. (36) has been used to determine the two-photon densities shown in Fig. 4 of the main text.