Cooperative resonances in light scattering from two-dimensional atomic arrays
Ephraim Shahmoon, Dominik S. Wild, Mikhail D. Lukin, Susanne F. Yelin
CCooperative resonances in light scattering from two-dimensional atomic arrays
Ephraim Shahmoon, Dominik S. Wild, Mikhail D. Lukin, and Susanne F. Yelin
1, 2 Department of Physics, Harvard University, Cambridge MA 02138, USA Department of Physics, University of Connecticut, Storrs, Connecticut 06269, USA (Dated: October 10, 2016)We consider light scattering off a two-dimensional (2D) dipolar array and show how it can betailored by properly choosing the lattice constant of the order of the incident wavelength. In par-ticular, we demonstrate that such arrays can operate as nearly perfect mirrors for a wide rangeof incident angles and frequencies close to the individual atomic resonance. These results can beunderstood in terms of the cooperative resonances of the surface modes supported by the 2D array.Experimental realizations are discussed, using ultracold arrays of trapped atoms and excitons in 2Dsemiconductor materials, as well as potential applications ranging from atomically thin metasurfacesto single photon nonlinear optics and nanomechanics.
Control over propagation and scattering of light fieldsplays a central role in optical science and photonics.In particular, it is well known that emitters exhibit astrongly modified linear and nonlinear optical responseon resonance. For example, enhanced optical scatter-ing in 2D arrays of linearly polarizable elements havebeen extensively studied in photonics [1–4]. Recently, ithas been shown that thin 2D metamaterials, known asmetasurfaces, whose constituent elements are optical an-tennas with varying resonances, can drastically alter thetransmitted field by enabling spatial control of its ampli-tude, phase and polarization [5, 6]. As a rule, these el-ements are micro-fabricated from macroscopic material,while the separation between the array elements is typ-ically much smaller than the operating wavelength. Atthe same time, resonant light can be completely reflectedby individual atoms when they are strongly coupled tonanophotonic devices with sub-wavelength localization oflight [7–11]. Intuitively, this originates from resonant en-hancement of the optical cross section of a polarizabledipole, which at resonance universally scales as λ , λ be-ing its resonant wavelength. Such single atom reflectorsyield extraordinary nonlinearities at the level of individ-ual photons [12–14].In this Letter we explore light scattering from a 2D or-dered and dilute array of atoms, with a lattice constantof the order of a wavelength, as can be realized, e.g. us-ing ultracold atoms loaded into optical lattices [15, 16].In such a case near-resonant operation can still lead tostrong scattering. Indeed, vanishing transmission at nor-mal incidence was recently discovered in a numericalstudy of a 2D atomic lattices for a specific frequency andlattice arrangement [17]. Due to resonant enhancement,one may na¨ıvely expect that a single layer of dipoles,even if they are as small as individual atoms, may “tile”the plane and thus act as a strong scatterer, providedthe density of dipoles exceeds 1 /λ (Fig. 1a). This rea-soning, though providing intuition for the possibility ofstrong scattering in dilute media, ignores the importanteffect of interactions between the dipoles [18–21].In what follows, we develop an analytical approach to (a) (b)(c) FIG. 1: (a) 2D array of atoms spanning the xy plane at z = 0, with inter-atomic spacing a on the order of the res-onant wavelength of the atoms, λ . For resonant light, theindividual atomic cross section is of order λ (dashed circles).(b) Light scattering off the array in the single diffraction orderregime: The incident field (yellow arrow) produces a forwardscattered field at z > z <
0. (c)Intensity transmission coefficient T and reflection coefficient R for a square-lattice at normal incident and resonant light( δ = 0) as a function of the lattice constant a . Strong scat-tering is observed with perfect reflection (vanishing transmis-sion) occurring at a/λ ≈ . , . the scattering problem and show that the near-unity re-flection is a generic phenomenon associated with the co-operative resonances of the dipolar array and its collec-tive surface-wave excitations. Strong scattering occurswhen the frequency of the incident light matches thatof the cooperative resonance of the array. The controlof scattering off the array can be achieved by adjustingthe lattice constant, which determines the cooperativeresonances. We demonstrate that the array can forma perfect mirror at almost all incident angles and po-larizations, where both the resonant frequency and thebandwidth can be tuned. a r X i v : . [ qu a n t - ph ] O c t The atomic array.—
We consider a 2D array of identi-cal point-like particles with a generic linear and isotropicpolarizability [22] α ( δ ) = − π (cid:15) λ a γ/ δ + i ( γ + γ nr ) / . (1)Here δ = ω − ω a , with | δ | (cid:28) ω a , is the detuning betweenthe frequency of the incident light ω = 2 πc/λ and that ofthe resonance of the particles ω a = 2 πc/λ a , and γ ( γ nr )is the radiative (non-radiative) width of this resonance.For a closed cycling transition in atoms we have γ nr = 0and the isotropic and linear response corresponds to a J = 0 to J = 1 transition far from saturation. Thearray is taken to be an infinite square lattice with latticeconstant a < λ , spanning the xy -plane at z = 0 (Fig. 1a). Scattering at normal incidence.—
We first focus on thesimplest case of a plane wave at normal incidence. Thecondition a < λ guarantees that only a single diffractionorder is present in the far field such that the scatteredfield on both sides of the array consists of plane wavespropagating in the z -direction (Fig. 1b for waves along z ). Figure 1c shows the transmission and reflection coef-ficients as a function of the lattice constant, computed forresonant light δ = 0 and in the absence of non-radiativelosses, γ nr = 0, using our analytical approach presentedbelow. We observe that the array scatters strongly overa wide range of lattice constants. In particular, com-plete reflection (zero transmission) is observed at latticeconstants a/λ ≈ .
2, 0 .
8. We note that the null trans-mission at a/λ ≈ . a < λ thetotal field can be written as E = (cid:104) e ikz + Se ik | z | (cid:105) E , (2)where E is the amplitude of the field polarized in the xy -plane, k = ω/c , and S is a scattering amplitude. For S = −
1, the transmitted field (at z >
0) vanishes and thecorresponding perfect reflection gives rise to a standingwave for z <
0. The scattering amplitude is determinedby the polarization p induced on the atoms by the inci-dent field, which is identical for all atoms in this case. Inturn, p is the result of multiple scattering of the incidentfield by all atoms in the array, and it can be character-ized by an effective polarizability of the atoms defined by p = α e ( δ ) E . A self-consistent solution of this multiple-scattering problem yields [23] S ( δ ) = iπ (cid:18) λa (cid:19) α e ( δ ) ε λ = − i ( γ + Γ) / δ − ∆ + i ( γ + γ nr + Γ) / . (3)By comparing the structure of this linear response to thatof an individual atom, Eq. (1), we infer that the dipolarinteraction between atoms in the array renormalize boththe width γ and the resonant frequency ω a . They are (a) (b) FIG. 2: (a) The cooperative shift ∆, Eq. (4), as a function ofthe lattice constant a (normal incidence). This plot is centralin the design of the scattering since the shift determines thecollective resonances of the array according to Eq. (3). Per-fect reflection occurs when the cooperative shift equals theincident detuning, δ = ∆. For example, ∆ = 0 at a/λ ≈ . . R as a function of lattice constant a and detuning δ . We note that the emerging contour of perfect reflection(bright yellow) coincides with the cooperative resonance plot-ted in (a) (marked here by the dashed black curve). now supplemented by their cooperative counterparts Γand ∆, respectively, given by∆ − i − γλ (cid:88) n (cid:54) =0 G (0 , r n ) , Γ = γ π (cid:18) λa (cid:19) − γ. (4)Here, G (0 , r n ) is the transverse component ( xx or yy )of the dyadic Green’s function of electrodynamics in freespace [24], evaluated between the central atom (“ n = 0”)at r = 0 and the atom n at r n . The explicit expressionfor Γ holds for a < λ and is in fact valid for any δ = ∆. Perfect reflection ( S = −
1) occursif additionally γ nr = 0. Therefore, the key ingredient thatdetermines the scattering properties of the array is the co-operative dipole-dipole shift ∆, given by the summation(readily evaluated numerically) of the dispersive dipole-dipole shift over all atoms, the real part of Eq. (4). Figure2a provides us with a central tool by which to understandand design the scattering off the array, as it presents thecooperative shift ∆ as a function of the lattice constant a [23]. For example, the vanishing cooperative shift ∆at a/λ ≈ . , . δ = 0. Moreover, Fig. 2a shows that scat-tering resonances exist for a wide range of incident fielddetunings δ near the individual-atom resonance. This isillustrated by Fig. 2b, in which the reflection coefficientis plotted as a function of both a and δ .For lossy particles, where γ nr (cid:54) = 0, the scattering ampli-tude (3) at resonance becomes S = − (Γ+ γ ) / (Γ+ γ + γ nr ).Therefore, high reflection requires that radiation damp-ing via scattering is dominant over all other dampingsources, γ + Γ (cid:29) γ nr . The scaling γ + Γ ∝ ( λ/a ) , origi-nating from cooperative enhancement, then implies thatthis can be achieved for a sufficiently small lattice con-stant even if the individual dipoles are poor radiators( γ < γ nr ). General angle of incidence.—
The foregoing analysiscan be generalized to all incident angles. We begin byconsidering a < λ/
2, which ensures a single diffractionorder for all incident plane waves, E , k (cid:107) e i k (cid:107) · r e ik z z , at anyangle. Here k (cid:107) = ( k x , k y ,
0) denotes the projection of theincident wave vector onto the xy plane and E , k (cid:107) canbe decomposed into the two possible transverse polariza-tions e + p,s ⊥ k , see Fig. 3a. The total field has the formof Eq. (2) where the scattering amplitude now becomesa 3 × e i k (cid:107) · r (cid:107) E , k (cid:107) replacing E . Thescattering amplitude is again determined by the polar-ization of the atoms, which is spatially modulated by thein-plane incident wavevector, according to Bloch’s theo-rem. The polarization of atom n can thus be written as p n = p ( k (cid:107) ) e i k (cid:107) · r n , where p ( k (cid:107) ) = α e ( k (cid:107) ) E , k (cid:107) (5)denotes the polarization in momentum space. Hence, theeffective polarizability is generally defined as the linearresponse of the polarization of the array in momentumspace, given by the tensor α e ( k (cid:107) ) = − π ε λ γ/ δ − ∆( k (cid:107) ) + i [ γ + γ nr + Γ( k (cid:107) )] / . (6)In analogy with Eqs. (1) and (3), ∆( k (cid:107) ) and Γ( k (cid:107) ) are thecooperative resonance and width tensors, respectively,given in terms of the dyadic Green’s function G by∆( k (cid:107) ) − i k (cid:107) ) = − γλ (cid:88) n (cid:54) =0 G (0 , r n ) e i k (cid:107) · r n . (7)An analytic expression can be obtained for Γ [23], while∆ has been evaluated numerically.The scattering amplitude is related to the effective po-larizability α e ( k (cid:107) ) by an expression similar to that inEq. (3), from which we can deduce the intensity reflec-tion and transmission coefficients, R µν and T µν , where µ and ν denote the polarization, either p or s , of thereflected/transmitted and the incident field, respectively[23]. Figures 3b and 3c display R ss and R pp as a func-tion of the incident angle (represented by k x , k y ) in thecase of a = 0 . λ and δ = 0, for which full reflection wasobtained at normal incidence. Remarkably, we observethat the reflection coefficient exceeds 0 .
99 at all incidentangles for s -polarized light. In the case of p -polarization,reflection exceeds 0 .
95 for a wide range of angles, far be-yond the paraxial approximation. Furthermore, mixingbetween polarizations is small T sp , T ps , R sp , R ps ∼ − [23], which demonstrates that the array operates as anexcellent mirror for almost all incident angles and both (d)(a)(c) (b) FIG. 3: Scattering at a general angle of incidence. (a)The polarizations of incident and forward scattered fieldswith wavevector k are spanned by the forward basis e + p,s ⊥ k ,whereas that of the reflected field is spanned by the backwardpolarization basis e − p,s ⊥ k [23]. (b) Scattering properties of anarray with lattice constant a = 0 . λ . Intensity reflection coef-ficient R ss for s -polarized incident and scattered fields at zerodetuning from the bare atomic resonance ( δ = 0) as a functionof the in plane components k x,y of the incident wavevector.(c) Same as (b) for p -polarized reflected and incident fields.(d) ss component of the cooperative shift matrix ∆. The vari-ation of the energy shift around the resonance ∆ ss = δ = 0is small compared to an atomic linewidth, which explains thehigh reflection R ss at all angles. incident polarizations. The fact that the scattering res-onance found at normal incidence persists for incidentangles well beyond the paraxial regime, implies that themirror should operate well for realistic finite size incidentbeams and arrays. This was verified for Gaussian beamswith a waist smaller than the array size, using a directnumerical approach [23].The high reflection at oblique angles may again be un-derstood in terms of cooperative resonances of the atomarray. For example, in Fig. 3d we plot the ss matrix ele-ment of ∆, which is seen to vary by less than an atomiclinewidth over all incident angles. This provides an ex-planation for the excellent reflection of s -polarized light.For p -polarized light, a similar explanation holds onlywithin the paraxial regime, beyond which the discussionis complicated by the polarization degree of freedom [23]. Beyond single diffraction order.—
An additionaldiffraction order can appear when the lattice constantexceeds λ/ ∼ π/a , to the incident wavevector dueto Bragg scattering by the lattice, may now result in apropagating wave. This situation can be analyzed by astraightforward extension of the above formalism. For (b)(a)(c) FIG. 4: (a) Scattering beyond a single diffraction order.Compared to Figs. 1b and 3a, an additional diffraction orderappears on both sides of the array. The figure presents re-sults of a numerical calculation [23] for an s -polarized Gaus-sian beam of waist 3 λ and δ = ∆, incident at a 30 ◦ anglerelative to the z axis. The beam is scattered off an arrayof 26 ×
26 atoms with a = 0 . λ . All field intensities arein units of | E | . (b) Retro-reflector effect for the same ar-ray and incident Gaussian beam, this time p -polarized at anincident angle of 45 ◦ . Scattering is strongly suppressed intodirections for which the p -polarization is orthogonal to the p -polarization of the incident beam. As a result, most of thelight is back scattered into the first diffraction order, whichis parallel to the incident beam. (c) Band structure of thecollective surface modes of the atom array for a = 0 . λ . Thethree bands correspond to the three eigenvalues of the coop-erative shift ∆( k (cid:107) ). The modes in the yellow (dashed) bandare polarized along the z direction, while the polarizations ofthe other two bands lie in the xy plane. The inset shows thelocation of the special points Γ, X , M in the first Brillouinzone. The dashed vertical lines (main figure) and the dashedcircle (inset) indicate the light cone | k (cid:107) | = 2 π/λ . example, in the case of p -polarization, where the propa-gation along certain directions is suppressed, new possi-bilities arise, such as the retro-reflection effect in Fig. 4b. Surface dipole excitations.—
More insight into thephysics of the array is gained by noting that the cooper-ative shift ∆( k (cid:107) ) in fact describes the dispersion relationof collective surface dipole excitations. The nature ofthese surface modes is revealed by Eq. (5) as the nor-mal modes of the atomic dipoles on the surface, p ( k (cid:107) ).The resonant frequencies of the modes and their corre-sponding polarizations can be deduced from their linearresponse α e ( k (cid:107) ) in Eq. (6) as the three eigenvalues andeigenvectors of ∆( k (cid:107) ). This interpretation also followsfrom the quantum master equation governing the dynam-ics of the atoms [23]. In that fully equivalent treatment,the eigenvalues of ∆( k (cid:107) ) arise naturally as the energies ofthe Bloch modes of atomic excitations, while Γ( k (cid:107) ) givestheir respective decay rate. When considering the scattering problem, k x and k y are restricted to be less than 2 π/λ . On the other hand,the discussion of collective surface excitation requires nosuch assumption and we may extend k x and k y to theentire Brillouin zone [ − π/a, π/a ]. By diagonalizing ∆ foreach k (cid:107) we hence obtain the band structure of the surfacemodes shown in Fig. 4c. The modes near the center ofthe Brillouin zone (Γ), between the vertical dashed lines,carry a crystal momentum smaller than 2 π/λ and maythus couple to far-field radiation. These are preciselythe modes that give rise to high reflection when drivenresonantly. Beyond the vertical dashed lines in Fig. 4c,where | k (cid:107) | > π/λ , the surface can no longer couple tofar-field radiation and these modes are thus protectedagainst decay, satisfying Γ + γ = 0. As such, these modescan only be excited in the near field and could potentiallyguide light along the surface. Experimental considerations.—
Three promising ap-proaches for experimental realization of the 2D array canbe considered. First, an array of ultracold atoms may betrapped in an optical lattice. For the single scattering or-der regime, a < λ/
2, a blue-detuned trapping laser can beused [25]. In a second implementation, the dipole arraycan be realized using plasmonic nano-particles [2, 4]. Ourresults should readily apply to this situation, providedthe collective radiative decay exceeds the individual non-radiative losses, as discussed above. Finally, 2D semicon-ductors such as monolayers of transition metal dichalco-genides [26] can be used, in which excitons and trionswith near lifetime limited linewidths have recently beenobserved [27]. A lattice structure for the excitons (trions)can be created by periodically modulating strain [28, 29],varying the dielectric environment [30], or by applying anexternal electrostatic potential. In addition to giving riseto the unusual dispersion relation discussed above, con-fining the excitons (trions) to sub-wavelength scales mayimprove their optical properties by overcoming defect-induced localization [31].Any of these realizations may also be subject to dis-order which may affect the cooperative resonances andscattering of the array. Interestingly, we show in [23]that the cooperative resonances are robust to fluctua-tions in the atomic positions when the fluctuations aremuch smaller than the lattice period. In addition, wefind that for small values of a/λ , scattering is also robustto defects such as a missing atom.
Discussion and prospects.—
The above considerationsdemonstrate that the scattering properties of light at anyincident polarization and detuning δ off a 2D atomic ar-ray can be controlled by choosing the lattice spacing ofthe array a , enabling, for instance, its operation as anearly perfect mirror.These results suggest the potential use of such 2D ar-rays as powerful platforms for classical and quantum op-tics. While our results are obtained for a uniform atomicarray, it should be possible to generalize our approachto non-homogeneous arrays, resulting in the realizationof “atomic-scale metasurfaces” with desired properties.One particularly intriguing possibility is to engineer thesurfaces such that a given incident mode is effectivelyfocused on a single impurity atom. While even for theuniform array we find optical nonlinearities at the level of ∼
10 photons [23], such engineered arrays will likely dis-play nonlinearities down to a single photon level similarto 1D systems [7]. Alternatively, exceptional nonlinear-ities can be realized by combing the present approachwith Rydberg blockade and electromagnetically-inducedtransparency [32–34].This work also opens up new prospects in the field ofoptomechanics. Since the atoms are very light but at thesame time collectively exhibit nearly perfect reflection,they form a highly mechanically-susceptible mirror, po-tentially very useful for the exploration of optomechanicsat the quantum level [35].Finally, we stress the universality of our approachfor any physical system of waves and dipole-like scat-terers. In particular, our analysis of cooperative reso-nances is not limited to electric dipoles, and the elec-tromagnetic Green’s function can be readily replaced bya Green’s function describing an entirely different wavephenomenon at any wavelength.We thank J´anos Perczel for insightful comments con-cerning the quantum description and dispersion relationof the array. We also acknowledge valuable discussionswith Vladimir Shalaev, Markus Greiner, Peter Zoller,Darrick Chang, Hongkun Park, Alex High and Kristi-aan de Greve, and financial support from NSF and theMIT-Harvard Center for Ultracold Atoms. [1] Ebbesen, T. W., H. J. Lezec, H. F. Ghaemi, T. Thio, andP. A. Wolff, Nature (London) , 667 (1998).[2] F. J. Garc´ıa de Abajo, Rev. Mod. Phys. , 1267 (2007).[3] S. Xiao and N. A. Mortensen, Optics Lett. , 37 (2011).[4] S. Thongrattanasiri, F. H. L. Koppens and F. J. Garc´ıade Abajo, Phys. Rev. Lett. , 047401 (2012).[5] N. Yu and F. Capasso, Nat. Mater. , 139 (2014).[6] A. V. Kildishev, A. Boltasseva and V. M. Shalaev, Sci-ence , 1289 (2013).[7] D. E. Chang, A. S. Sørensen, E. A. Demler, andM. D. Lukin, Nat. Phys. , 807 (2007).[8] J.-T. Shen and S. Fan, Phys. Rev. Lett. ,153003(2007).[9] A. Asenjo-Garcia, J. D. Hood, D. E. Chang and H. J.Kimble, arXiv:1606.04977.[10] N. V. Corzo, B. Gouraud, A. Chandra, A. Goban, A. S.Sheremet, D. V. Kupriyanov, and J. Laurat, Phys. Rev.Lett. , 133603 (2016).[11] H. L. Sørensen, J.-B. B´eguin, K. W. Kluge, I. Iakoupov,A. S. Sørensen, J. H. M¨uller, E. S. Polzik and J. Appel,Phys. Rev. Lett. , 133604 (2016). [12] T. G. Tiecke, J. D. Thompson, N. P. de Leon, L. R. Liu,V. Vuleti´c and M. D. Lukin, Nature , 241 (2014).[13] I. Shomroni, S. Rosenblum, Y. Lovsky, O. Bechler, G.Guendelman and B. Dayan, Science , 903 (2014).[14] A. Sipahigil, R. E. Evans, D. D. Sukachev, M. J. Bu-rek, J. Borregaard, M. K. Bhaskar, C. T. Nguyen, J. L.Pacheco, H. A. Atikian, C. Meuwly, R. M. Camacho, F.Jelezko, E. Bielejec, H. Park, M. Lon˘car, M. D. Lukin,arXiv:1608.05147.[15] I. Bloch, Nat. Phys. , 23 (2005).[16] H. Labuhn, D. Barredo, S. Ravets, S. de L´es´eleuc, T.Macr`ı, T. Lahaye and A. Browaeys, Nature , 667(2016).[17] R. J. Bettles, S. A. Gardiner and C. S. Adams, Phys.Rev. Lett. , 103602 (2016).[18] J. Pellegrino, R. Bourgain, S. Jennewein, Y. R. P. Sortais,A. Browaeys, S. D. Jenkins, and J. Ruostekoski, Phys.Rev. Lett. , 133602 (2014).[19] K. Kemp, S. J. Roof, M. D. Havey, I. M. Sokolov, andD. V. Kupriyanov, arXiv:1410.2497.[20] M. O. Scully and A. A. Svidzinsky, Science , 1239(2010).[21] R. T. Sutherland and F. Robicheaux, Phys. Rev. A ,013847 (2016).[22] P. Lambropoulos and D. Petrosyan, Fundamentals ofQuantum Optics and Quantum Information (Springer,2007).[23] see Supplemental Material for more details.[24] L. Novotny and B. Hecht,
Principles of Nano-Optics (Cambridge University Press, 2006).[25] T. Kinoshita, T. Wenger and D. S. Weiss, Science ,1125 (2004).[26] Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Colemanand M. S. Strano, Nat. Nanotech. , 699 (2012).[27] C. Robert, D. Lagarde, F. Cadiz, G. Wang, B. Lassagne,T. Amand, A. Balocchi, P. Renucci, S. Tongay, B. Ur-baszek, and X. Marie, Phys. Rev. B , 205423 (2016).[28] H. Li, A. W. Contryman, X. Qian, S. Moeini Ardakani,Y. Gong, X. Wang, J. M. Weisse, C. H. Lee, J. Zhao, P.M. Ajayan, J. Li, H. C. Manoharan and X. Zheng, Nat.Commun. , 7381 (2015).[29] C. Palacios-Berraquero, D. M. Kara, A. R.-P. Mont-blanch, M. Barbone, P. Latawiec, D. Yoon, A. K. Ott,M. Loncar, A. C. Ferrari, M. Atat¨ure, arXiv preprintarXiv:1609.04244 (2016).[30] Y. Lin, X. Ling, L. Yu, S. Huang, A. L. Hsu, Y.-H. Lee,J. Kong, M. S. Dresselhaus, and T. Palacios, Nano Lett. , 5569 (2014).[31] A. Srivastava, M. Sidler, A. V. Allain, D. S. Lembke, A.Kis and A. Imamo˘glu, Nat. Nanotech. , 491 (2015).[32] M. D. Lukin, M. Fleischhauer, R. Cote, L. M. Duan, D.Jaksch, J. I. Cirac, and P. Zoller, Phys. Rev. Lett. ,037901 (2001).[33] O. Firstenberg, C. S. Adams and S. Hofferberth, J. Phys.B: At. Mol. Opt. Phys. ,163601 (2007).[35] M. Aspelmeyer, T. J. Kippenberg and F. Marquardt,Rev. Mod. Phys. , 1391 (2014). Supplemental Materials: Cooperative resonances in light scattering fromtwo-dimensional atomic arrays
This Supplementary document is organized as follows. Sec. 1 reviews the general approach we use for the scatteringproblem by an array of dipole scatterers. Sec. 2 is dedicated for the derivation of our analytical formulation ofscattering by a 2D lattice, including its quantum version. Sec. 3-5 discuss in more depth some of the results andconsequences that emerge from the our analytical treatment. Finally, Sec. 6 presents the direct numerical approachwe employ to verify our theoretical results in a realistic finite size scenario.
1. THE SCATTERING PROBLEM: GENERAL APPROACH
We consider the scattering of a general incident electric field by a collection of atoms modeled as point-like dipolescatterers characterized by their linear polarizability. This description naturally applies for the scattering due to theelectric linear response of any collection of scatterers which are much smaller than the incident wavelength, such asatoms in the optical frequency domain, and has been used before by many authors in the photonics community, e.g.for the treatment of light scattering by arrays of nano-particles [S1].
Assuming an incident field E ( r ) and atoms n = 1 , .., N at positions r n , the goal is to find the electric field at anygiven point r . As usual, we begin with Maxwell’s curl equations for fields at frequency ω = kc = k/ √ ε µ , obtainingthe wave equation ∇ × ∇ × E − k E = k ε P , (S1)whose formal solution is E i ( r ) = E ,i ( r ) + k ε (cid:88) j (cid:90) V d r (cid:48) G ij ( k, r , r (cid:48) ) P j ( r (cid:48) ) , (S2)with A i = e i · A the i = x, y, z component of a vector A and where G ij ( k, r , r (cid:48) ) is the Green’s function of Eq. (S1),namely the i component of electric field (in inverse length units) given a delta-function source at r (cid:48) polarized on the j axis; in free space this so-called dyadic Green’s function takes the form [S2] G ij ( k, r , r ) = e ikr πr (cid:20)(cid:18) ikr − k r (cid:19) δ ij + (cid:18) − − ikrk r (cid:19) r i r j r (cid:21) , (S3)with r = r − r , r = | r | and r i = e i · r .Considering now that polarization P exists only on the atoms and using their linear response α ( ω ) (generally atensor, but taken as scalar/isotropic and identical for all atoms), we have P ( r ) = N (cid:88) n =1 α E ( r ) δ ( r − r n ) , (S4)so that Eq. (S2) becomes E i ( r ) = E ,i ( r ) + 4 π αε λ (cid:88) j (cid:88) n λG ij ( k, r , r n ) E j ( r n ) , (S5)where λ = 2 π/k . Eq. (S5) has a Lippman-Schwinger form, and together with the Green’s function, Eq. (S3), it formsthe formal solution and the starting point of our scattering theory. The idea is to self-consistently evaluate the fields E j ( r n ), or more precisely the polarizations p j ( r n ) = αE j ( r n ), at the atomic positions r = r n , which in turn determinethe right-hand side of the equation, thus allowing for solution of the field at any point, E i ( r ) at the left-hand side ofthe equation. Writing Eq. (S5) at an atomic position r = r n , we obtain a self-consistent equation for the polarizationinduced on the atoms p ni = αE i ( r n ), p ni = p n ,i + (cid:88) m (cid:54) = n (cid:88) j π αε λ λG nmij p mj , (S6)with p n ,i = αE ,i ( r n ) and G nmij = G ij ( k, r n , r m ). Here the summation over atoms (index m ) skips the atom n atwhich we evaluate the field, in accordance with the convention we adopt for the atomic polarizability, as explainedbelow. The idea is that by inserting the solution of Eq. (S6) into Eq. (S5), we obtain the field at any point E i ( r ).The solution of Eq. (S6) and then of Eq. (S5) can be readily performed numerically for any general case of incidentfield and finite collection of atoms, as we show in Sec. below. Nevertheless, by exploiting the symmetry of an orderedarray, an analytical solution can be obtained, as shown in Sec. (similar to the approach taken in Ref. [S1]). Considering a J = 0 to J = 1 transition of an atom with a dipole matrix element d and frequency ω a the linearpolarizability is given by [S3, S4] α ( ω ) = 2 d ω a (cid:126) ω a − ω − iγ ω , (S7)where γ is a damping term which is addressed below. It is worth recalling here that a similar form of polaizability isfound for classical systems alike, so that our discussion is kept general for all point-like polarizable scatterers in theirlinear regime. For small detuning δ = ω − ω a with respect to ω a , | δ | (cid:28) ω a , and using the spontaneous emission rateof a single atom γ = d ω a / (3 πε (cid:126) c ), we obtain α ( ω ) = − π ε λ a γ/ δ + iγ / . (S8)Skipping the m = n term in the sum of Eq. (S6) implies that the radiative damping γ of an individual atom isincluded in the total damping γ . Namely, γ = γ + γ nr , (S9)where γ nr denotes the non-radiative loss rate. This also implies that ω a already includes the Lamb shift correction(or that we neglect it). We note, that one may chose to include the m = n term in Eq. (S6) by neglecting the realpart of G nnij (neglecting the Lamb shift) and taking γ = γ nr , yielding exactly the same results for p ni and E i ( r ).
2. ANALYTICAL APPROACH2.1 Polarization built on the atoms
We consider an infinite array of atoms spanning the xy plane at z = 0 and illuminated by an incident field representedby its decomposition to plane waves E ( r ) = (cid:80) k (cid:107) E , k (cid:107) e i k · r . The wavevector of each plane wave is characterized byits projections on the xy plane and z axis, k (cid:107) and k z = (cid:112) k − | k (cid:107) | , respectively, whereas its amplitude E , k (cid:107) isspanned by the two transverse polarizations e + p,s ⊥ k . At normal incidence ( k (cid:107) = 0) we have k z = k and the field ispolarized in the xy plane, e + p,s ∈ { e x , e y } .In order to exploit the discrete translational symmetry of the lattice (as in the Bloch theorem), we define the 2DFourier transform of a function f ( r ) sampled at the lattice sites n , f n = f ( r n ), as f ( k (cid:107) ) = 1 N N (cid:88) n =1 e − i k (cid:107) · r n f n , f n = N A (cid:90) BZ d k (cid:107) (2 π ) e i k (cid:107) · r n f ( k (cid:107) ) , (S10)where A is the area of the lattice unit cell and BZ means that the integral over k (cid:107) is performed within the firstBrillouin zone of the reciprocal lattice. For a square lattice, r n = r n x ,n y = a ( n x e x + n y e y ), we have A = a and (cid:82) BZ d k (cid:107) = (cid:82) π/a − π/a dk x (cid:82) π/a − π/a dk y . Applying this Fourier transformation on Eq. (S6) we find p ( k (cid:107) ) = α E , k (cid:107) + 4 π αε λ λg ( k (cid:107) ) p ( k (cid:107) ) (S11)for each k (cid:107) value contained in the incident field, where g ( k (cid:107) ) = (cid:88) n (cid:54) =0 G (0 , r n ) e − i k (cid:107) · r n , (S12)also noting the fact that the dyadic Green’s function in free space, G ( r , r (cid:48) ) depends only on the difference r − r (cid:48) , seeEq. (S3). The solution of (S11) is then p n = (cid:88) k (cid:107) p ( k (cid:107) ) e i k · r n , p ( k (cid:107) ) = α e ( k (cid:107) ) E , k (cid:107) , (S13)with the effective polarizability tensor α e ( k (cid:107) ) = − π ε λ a γ/ δ − ( λ a /λ ) ∆( k (cid:107) ) + i [ γ + γ nr + ( λ a /λ ) Γ( k (cid:107) )] / , (S14)and the collective radiative response from Eq. (7) of the main text, namely,∆( k (cid:107) ) = − γλ Re[ g ( k (cid:107) )] , Γ( k (cid:107) ) = 3 γλ Im[ g ( k (cid:107) )] . (S15)Finally, the effective polarizability from Eq. (6) of the main text is reached by taking λ a /λ ≈ ω ≈ ω a in quantum optics. We note that Eq. (S13) implies that thepolarization built on the atomic array due to an incident plane wave with in-plane wavevector k (cid:107) is not mixed withthat induced by a different wavevector, and that the collective response of the array polarization to a field with agiven wavenumber k (cid:107) is described by the effective polarizability α e ( k (cid:107) ).The result at normal incidence, Eq. (3) in the main text, is obtained by noting that the sum g (0) for k (cid:107) = 0 inEq. (S12) is symmetric in the x and y directions due to the structure of G from Eq. (S3) and the symmetry of thelattice. Therefore, the tensor g becomes diagonal g ij = δ ij g i with g x = g y . At normal incidence the field is polarizedin the xy plane so that g effectively appears as a scalar, g x = (cid:80) n (cid:54) =0 G xx (0 , r n ), and the collective response becomes∆ − i Γ / − (3 / γλg x . Both ∆ and Γ are obtained by a direct numerical summation over the atomic lattice points, as in Eq. (S15).Alternatively, one can begin by representing the sum g ( k (cid:107) ) in 2D reciprocal (wavevector) space by using the relation[S2] G ( x, y, z ) = e ik √ x + y + z π (cid:112) x + y + z = i π (cid:90) d k (cid:48)(cid:107) e − i k (cid:48)(cid:107) · r (cid:107) e ik (cid:48) z | z | k (cid:48) z , (S16)where r (cid:107) = ( x, y ) and k (cid:48) z = (cid:113) k − | k (cid:48)(cid:107) | . Then, writing g ij ( k (cid:107) ) = (cid:88) n (cid:90) d r (cid:107) G ij (0 , r (cid:107) ) e − i k (cid:107) · r (cid:107) δ ( r (cid:107) − r n ) − G ij (0 ,
0) (S17)(recalling that r n only exists on the xy plane), and using G ij = [ δ ij + (1 /k ) ∂ i ∂ j ] G , we find g ij ( k (cid:107) ) = (cid:40) i π (cid:82) d k (cid:48)(cid:107) k (cid:48) z (cid:16) δ ij − k (cid:48) i k (cid:48) j k (cid:17) (cid:80) n (cid:82) r (cid:107) δ ( r (cid:107) − r n ) e − i ( k (cid:48)(cid:107) + k (cid:107) ) · r (cid:107) − i λ δ ij , for i, j = { x, y } ∪ i = j = z, − i λ δ ij , otherwise . (S18)Here we used G ij (0 ,
0) = i λ δ ij obtained by taking the r → G ij ( r ,
0) and neglecting its real part (associatedwith the Lamb shift of the atomic resonance ω a ). We may now introduce the reciprocal lattice as the Fourier transformof the atomic lattice structure ρ ( k (cid:48)(cid:107) ) = (cid:90) d r (cid:107) e − i k (cid:48)(cid:107) · r (cid:107) (cid:88) n δ ( r (cid:107) − r n ) = (2 π ) A (cid:88) m δ ( k (cid:48)(cid:107) − q m ) , (S19)which forms a lattice in 2D wavevector space at lattice points q m (wavevectors) satisfying q m · r n = 2 πM with aninteger M and for any n and m . Using Eq. (S19) in (S18) we obtain g ij ( k (cid:107) ) = i A (cid:88) m (cid:112) k − | q m − k (cid:107) | (cid:20) δ ij − ( q m − k (cid:107) ) i ( q m − k (cid:107) ) j k (cid:21) − δ ij i λ , (S20)valid for i, j = { x, y } or i = j = z whereas g ij ( k (cid:107) ) = − δ ij i λ otherwise. For the square lattice we have q m x ,m y =(2 π/a )( m x e x + m x e y ) and the unit cell area is A = a so that Eq. (S20) becomes g ij ( k (cid:107) ) = i a ∞ (cid:88) m x = −∞ ∞ (cid:88) m y = −∞ (cid:0) δ ij − [(2 π/a ) m x − k x ][(2 π/a ) m y − k y ] /k (cid:1)(cid:112) k − [(2 π/a ) m x − k x ] − [(2 π/a ) m y − k y ] − δ ij i λ . (S21)The sums in Eqs. (S20,S21) can be used for a numerical evaluation of g (and hence of ∆ and Γ) in momentum space,in contrast to the sum (S12) that is performed in real space. This is achieved by using a proper regularization andrenormalization procedure [S5]. More importantly, this form of g can be used to find an analytical expression for Γin the single diffraction order case. Considering only the imaginary part of g ij in Eq. (S21), the square root in thedenominator has to be real. Since this square root originates in the wavenumber in the z direction, k (cid:48) z that appearsin the wave expansion of the Green’s function (Eqs. S16 and S18), this requirement means that radiation damping(imaginary part of g ij ) is related to waves that can propagate out of the atomic array. Taking for example k y = 0,we have | k x − m x π/a | < k and | k x | < k , leading to | m x | < a/λ (for k y (cid:54) = 0 this condition on m x becomes evenmore restrictive). Then, for the case a < λ/ k x , k y ) only the zeroth-order diffraction m x = m y = 0contributes to radiation and damping. At normal incidence k x = k y = 0, such single diffraction order exists even for a < λ as can be inferred from the condition on m x . Similar conclusions are reached for any lattice structure with ascale a of its unit cell by considering Eq. (S20). Then, from Eq. (S15) we obtainΓ ij = γ π (cid:18) λa (cid:19) kk z (cid:18) δ ij − k i k j k (cid:19) − γδ ij , (S22)again valid for i, j ∈ { x, y } or i = j = z , while Γ ij = − γδ ij otherwise. At normal incidence we have k x = k y = 0 and k z = k leading to Γ = Γ xx = γ π λ a − γ as in Eq. (4) of the main text. For larger values of a more diffraction ordersmay appear and analytical results for their contribution to Γ are obtained in the same way. Let us now find the scattered field due to a given plane wave component of the incident field E , k (cid:107) e i k · r with anin-plane wavevector k (cid:107) . Returning to Eq. (S5), we use E n = p n /α = α e ( k (cid:107) ) E , k (cid:107) e i k (cid:107) · r n /α in the right-hand side,and obtain E ( r ) = E , k (cid:107) e i k (cid:107) · r e ik z z + 4 π λg sc ( k (cid:107) , r ) α e ( k (cid:107) ) ε λ E , k (cid:107) , (S23)with g sc ( k (cid:107) , r ) = (cid:88) n λG ( r , r n ) e i k (cid:107) · r n . (S24)Using again the above methods [the expansion (S16) and the reciprocal lattice (S19)] we find g sc,ij ( k (cid:107) , r ) = i A (cid:88) m (cid:112) k − | k (cid:107) + q m | (cid:20) δ ij − ξ ij ( k (cid:107) + q m ) i ( k (cid:107) + q m ) j k (cid:21) e i ( k (cid:107) + q m ) · r (cid:107) e i √ k −| k (cid:107) + q m | | z | , (S25)with r (cid:107) and z the projections of r along the xy plane and the z axis, respectively, and where ξ ij = − i = z ∩ j (cid:54) = z or j = z ∩ i (cid:54) = z at z <
0, and ξ ij = 1 otherwise. For a square lattice we have as usual A = a , (cid:80) m → (cid:80) ∞ mx = −∞ (cid:80) ∞ my = −∞ and q m → (2 π/a )( m x e x + m y e y ). Therefore, the field scattered from all atoms to a point r , represented by the sum g sc ( k (cid:107) , r ) from Eq. (S24), is revealed by Eq. (S25) to be a sum of plane waves contributionsfrom all diffraction orders m with wavevectors k (cid:107) + q m and (cid:112) k − | k (cid:107) + q m | projected along the xy plane and z axis, respectively.Since we are mainly interested in fields that are scattered and which propagate away from the surface along the ± z directions, the z -projected wavenumber (cid:112) k − | k (cid:107) + q m | has to be real (as opposed to the case of surface modeswhich are evanescent along z ). This leads to the same conditions on the diffraction orders m as analyzed in theanalytical calculation of Γ, so that for the single diffraction order case, e.g. a < λ/ λg sc,ij ( k (cid:107) , r ) = i π (cid:18) λa (cid:19) kk z (cid:18) δ ij − ξ ij k i k j k (cid:19) e i k (cid:107) · r (cid:107) e ik z | z | . (S26)The total field at r = ( r (cid:107) , z ) with | z | (cid:38) a sufficiently far from the atomic array where the evanescent waves do notcontribute, is then found by inserting the above g sc ( k (cid:107) , r ) into Eq. (S23), yielding E ( r ) = (cid:104) e ik z z + S ± ( k (cid:107) ) e ik z | z | (cid:105) e i k (cid:107) · r (cid:107) E , k (cid:107) , (S27)with the scattering matrix S ± ( k (cid:107) ) = 4 π (cid:15) λ e − i k (cid:107) · r (cid:107) e − ik z | z | g sc ( k (cid:107) , r ) α e ( k (cid:107) ) . (S28)The subscript ± in S ± distinguishes between forward (+ sign for z >
0) and backward ( − sign for z <
0) scattering.The fact that S ± depends on the sign of z , without z explicitly appearing in it, is due to its dependence on ξ ij whichis sensitive to the sign of z (see explanations below Eq. S25).It is constructive to analyze the scattering within the field polarization basis which is most natural for propagatingplane waves, namely, that of the wavevector k = k e k and the transverse polarizations perpendicular to it. Eq. (S27)reveals that for an incoming plane wave with a wavevector k = ( k (cid:107) , k z ), an additional wavevector k (cid:48) = ( k (cid:107) , − k z )emerges for the backward scattered field (at z <
0) while the forward scattered field (at z >
0) possess the incidentwavenumber k . We can therefore define two relevant polarization basis, one for the incoming and transmitted fieldsand the other for reflected ones, { e k , e + p , e + s } and { e k (cid:48) , e − p , e − s } , respectively, with e + p,s ⊥ e k and e − p,s ⊥ e k (cid:48) (Fig. 1b inthe main text). In spherical coordinates these basis vectors read e k = (sin θ cos φ, sin θ sin φ, cos θ ) , e + p = (cos θ cos φ, cos θ sin φ, − sin θ ) , e + s = (sin φ, − cos φ, e k (cid:48) = (sin θ cos φ, sin θ sin φ, − cos θ ) , e − p = ( − cos θ cos φ, − cos θ sin φ, − sin θ ) , e − s = (sin φ, − cos φ, , (S29)where all column vectors here are written in cartesian basis, i.e. as ( A x , A y , A z ) for a vector A .Writing g sc from Eq. (S26) as λg sc ( k (cid:107) , r ) = i π (cid:18) λa (cid:19) kk z F ± e i k (cid:107) · r (cid:107) e ik z | z | , (S30)with F ± ,ij = δ ij − ξ ij k i k j /k (S31)being its tensor part (written in cartesian basis) and where the ± subscript is due to the dependence of ξ ij on the signof z , we notice that the { e k , e + p , e + s } and { e k (cid:48) , e − p , e − s } basis form the eigenvectors of F + and F − , respectively, witheigenvalues { , , } , respectively, for both ± cases. Considering also that the incident field must posses one of twoforward transverse polarizations e + p,s , we can now describe the total field in Eq. (S27) via its components E µ = e µ · E µ projected onto the transverse polarization basis µ = { p, s } , E µ ( r ) = (cid:88) ν = p,s (cid:104) δ µν e ik z z + S ± µν ( k (cid:107) ) e ik z | z | (cid:105) e i k (cid:107) · r (cid:107) E ν , k (cid:107) . (S32)Here the scattering matrix is represented by its 2 × S ± µν = e ±† µ S e + ν , which are obtained form Eqs.(S28,S30) and the matrix elements e ±† µ F ± e ± ν = δ µν as S ± µν ( k (cid:107) ) = iπ (cid:18) λa (cid:19) kk z ε λ e ±† µ α e ( k (cid:107) ) e + ν . (S33)The scattering properties of the array can then be further characterized by the intensity transmission and reflectionmatrices, T µν = | δ µν + S + µν | and R µν = | S − µν | , where δ µν is the Kronecker delta.At normal incidence we recall that α e is a scalar in the { x, y } basis (which forms the { p, s } -polarization basis).Then, S is also a scalar S µν = Sδ µν where S is given by Eq. (S33) with k z = k and with the scalar polarizability atnormal incidence (see text below Eq. S14), α e = − π ε λ a γ/ δ − ∆ + i ( γ + γ nr + Γ) / , (S34)equivalent to Eq. (3) in the main text. We now extend our discussion to treat both the atoms and the electromagnetic field quantum mechanically. Anisotropic atom must consist of at least four levels: a ground state | g (cid:105) and three degenerate excited states | i (cid:105) , where i = x, y, z . The excited states are labeled by the direction of the dipole moment d i = d e i associated with thetransition | g (cid:105) → | i (cid:105) , where e i is the unit vector along the i axis. We further define the atomic lowering operatorsˆ σ i = | g (cid:105)(cid:104) i | and the corresponding raising operators ˆ σ † i = | i (cid:105)(cid:104) g | . For an array of atoms, we require an additionalindex, ˆ σ mi , to label the site m of the atom. By tracing out over the photonic degrees of freedom and making theBorn–Markov approximation, one arrives at a quantum master equation describing the dynamics of the atoms [S7].The master equation can be written in terms of the reduced density operator ρ ( t ) asdd t ρ ( t ) = − i [ H S , ρ ( t )] + D [ ρ ( t )] , (S35)where H S = (cid:88) m,i ω a σ † mi σ mi + (cid:88) m (cid:54) = n,i,j ∆ ij ( r m , r n ) σ † mi σ nj (S36)captures the coherent dynamics, while D [ ρ ] = (cid:88) m,n,i,j Γ ij ( r m , r n ) (cid:18) σ † mi ρσ nj − (cid:110) σ † mi σ nj , ρ (cid:111)(cid:19) (S37)corresponds to incoherent, dissipative evolution. The parameters ∆ ij ( r m , r n ) specify the coherent dipole–dipoleinteraction between two atoms, whereas Γ ij ( r m , r n ) denote the strength of their collective decay. These parameterscan be expressed in terms of the dyadic Green’s function of free space as∆ ij ( r m , r n ) = − πγcω a Re [ G ij ( ω a , r m , r n )] , Γ ij ( r m , r n ) = 6 πγcω a Im [ G ij ( ω a , r m , r n )] . (S38)We note the sum involving ∆ ij ( r m , r n ) in Eq. (S36) excludes m = n since this term merely corresponds to a renor-malization of the transition frequency ω a [S8–S10].For an infinite atomic array with discrete translational symmetry, it is convenient to introduce the momentum spaceoperators σ k i = (cid:88) m σ mi e i k · r m , σ mi = A (cid:90) BZ d k (2 π ) σ k i e − i k · r m . (S39)where k is a 2D wavevector restricted to the first Brillouin zone of the reciprocal lattice and A is the area of a unitcell. The Hamiltonian H S and the dissipator D may then be written as H S = A (cid:90) BZ d k (2 π ) (cid:88) i,j [ ω a δ ij + ∆ ij ( k )] σ † k i σ k j (S40) D [ ρ ] = A (cid:90) BZ d k (2 π ) (cid:88) i,j [ γδ ij + Γ ij ( k )] (cid:18) σ † k i ρσ k j − (cid:110) σ † k i σ k j , ρ (cid:111)(cid:19) . (S41)Here, we defined ∆ ij ( k ) = (cid:88) m (cid:54) = n ∆ ij ( r m , r n ) e i k · ( r m − r n ) , Γ ij ( k ) = (cid:88) m (cid:54) = n Γ ij ( r m , r n ) e i k · ( r m − r n ) , (S42)independent of n since G ( ω a , r m , r n ) only depends on r m − r n . For a given momentum k , the evolution of the atomsis thus specified by the two 3 × k ) and Γ( k ). The former specifies the energies of the three Bloch modeswith momentum k , giving rise to a bandstructure, while the latter describes their decay.The expressions for ∆ and Γ derived here are identical to those from the classical treatment. Of course, this isexpected since the classical results should follow from the quantum treatment in the linear regime, where the atomictransitions are not saturated. Furthermore, the reflection and transmission coefficients computed from the linearresponse of the quantum system must agree with those obtained from the classical theory. One can show that this isindeed the case by introducing a driving term to the Hamiltonian, corresponding to an incident field, and computingthe field radiated by the atoms.
3. PROPERTIES OF THE COOPERATIVE RESONANCE AND DECAY3.1 Structure of the functions ∆( a ) and Γ( a ) at normal incidence Fig. 2a of the main text presents ∆ as a function of the lattice spacing a normalized to the incident wavelength λ .The general features of its dependence on a/λ can be elegantly explained by establishing its Kramers-Kronig relationwith the function Γ( a/λ ). We recall the expressions for ∆ and Γ at normal incidence (c.f. Eq. 4 in the main textor Eq. S15 here and text below it), which are in general functions of the incident frequency ω ( λ = 2 πc/ω ) and thelattice points ∆( ω, a ) = − γλ (cid:88) n (cid:54) =0 Re G xx ( ω, , r n ) , Γ( ω, a ) = 3 γλ (cid:88) n (cid:54) =0 Im G xx ( ω, , r n ) . (S43)Since the dyadic Green’s function is a linear response function (of the modes that form the electromagneticfield/vacum), it satisfies the Kramers-Kronig relation,Re G xx ( ω, r , r (cid:48) ) = 1 π (cid:90) ∞−∞ dω (cid:48) Im G xx ( ω (cid:48) , r , r (cid:48) ) ω − ω (cid:48) . (S44)Writing the position vector of an atom n as r n = r n x ,n y = a n with n = ( n x e x + n y e y ) and considering the dydadicGreen’s function in free space, Eq. (S3), we note that ω and a always appear together as 2 πωa/c = a/λ so that λG xx ( ω, , r n ) can be written as a function of only the dimensionless distance a/λ and the vector index n (independentof a and ω ). Therefore, the Kramers-Kronig relation, Eq. (S44), for G xx ( ω, , r n ) can be written asRe[ λG xx ( a/λ, n )] = 1 π (cid:90) ∞−∞ du Im[ λG xx ( u, n )] a/λ − u . (S45)Inserting this relation into the expression for ∆ in Eq. (S43) and performing the sum over n , the dependence on n disappears and we obtain a Kramers-Kronig relation between the real and imaginary part of the field response as afunction of the scaled lattice spacing a/λ ∆( a/λ ) = − π (cid:90) ∞−∞ du Γ( u ) / a/λ − u , (S46) λ − ∆ / γΓ / γ FIG. S1: Cooperative shift ∆ and width Γ as a function of a/λ at normal incidence (extension of Fig. 2a in the maintext; results obtained by a direct summation of Eq. (4) in the main text). The diffraction-orders-associated peaks and theKramers-Kronig structure of the relation between the ∆ and Γ as a function of the scaled lattice constant a/lambda are clearlyseen (see text, Sec. ). where both ∆( u ) and Γ( u ) are now understood to be functions of a single variable u = a/λ .Turning now to the structure of the function ∆( a/λ ) we begin with the region a/λ (cid:28)
1. In this limit the sum over n in Eq. (S43) can be converted into an integral and we analytically find ∆ ∝ /a which agrees with the divergencenear a/λ → + in Fig. 2a of the main text and Fig. S1 here (see also Ref. [S11]). This 1 /a scaling is easily understoodby recalling that the real part of the Green’s function Eq. (S3), amounting to dipole-dipole interaction, scales as 1 /r at the quasistatic, short distance limit, and that ∆ is a sum over such dipole-dipole interactions. In order to explainthe behavior of ∆( a/λ ) at all regions beyond the a/λ (cid:28) a/λ ). For a/λ < dissipation channel for the atomicarray at a/λ <
1, imposing a dissipation Γ which results only from the term m x , m y = 0 in Eq. (S21). For a/λ = 1,additional dissipation channels arise in the form of emission to the directions (diffraction orders) | m x | = 1 , m y = 0and | m y | = 1 , m x = 0. This physically explains the peak observed in Γ( a/λ ) for a/λ = 1 in Fig. S1 (mathematically,additional poles contribute to the imaginary part of the sum in Eq. S21). The meaning is that whenever a newdissipation channel appears, in the form of a new diffraction order in the far field, we expect a peak in Γ( a/λ ).Indeed, this can be seen in Fig. S1 by the additional peaks at a/λ = √ | m x | = | m y | = 1) and a/λ = 2 (diffraction orders | m x | = 2 , m y = 0 and | m y | = 2 , m x = 0). In turn, these peaks in Γ( a/λ ) then physicallyexplain the corresponding dispersive-like peaks in ∆( a/λ ) following the Kramers-Kronig relation Eq. (S46), nicelyseen in Fig. S1 (whereas Fig. 2a in the main text captures only the a/λ = 1 ”dispersive” peak). ∆ and Γ For a < λ/
2, we derived an explicit expression for Γ( k (cid:107) ) in Eq. (S22). Written in terms of spherical coordinates θ and φ , the expression becomesΓ( k (cid:107) ) = 3 γ π cos θ (cid:18) λa (cid:19) − sin θ cos φ − sin θ cos φ sin φ − sin θ cos φ sin φ − sin θ sin φ
00 0 sin θ − γ . (S47)The matrix can be straightforwardly diagonalized. One obtains the three eigenvaluesΓ = 3 γ cos θ π (cid:18) λa (cid:19) − γ, Γ = 3 γ π cos θ (cid:18) λa (cid:19) − γ, Γ = 3 γ sin θ π cos θ (cid:18) λa (cid:19) − γ, (S48)with the respective eigenvectors v = (cos φ, sin φ,
0) = k (cid:107) / | k (cid:107) | , v = ( − sin φ, cos φ,
0) = e s , v = (0 , ,
1) = e z . (S49)We observe that v points along the direction of s -polarized incident light, while v and v are superpositions of the p -polarization and a longitudinal component.The cooperative shift ∆( k (cid:107) ) has the same block-diagonal structure as Γ( k (cid:107) ), such that v = (0 , ,
1) is also aneigenvector of ∆( k (cid:107) ). The other two eigenvectors, which we denote by w , , are in general different. However, itturns out that they are approximately equal to v , as demonstrated in Fig. S2. As a result, the s -polarization is anapproximate eigenvector of both Γ and ∆ and hence of the scattering matrix. This further implies that the surfaceonly weakly mixes the s - and p -polarizations. Indeed, for a/λ = 0 .
2, the reflection coefficient R sp never exceeds4 × − as shown in Fig. S2c ( R sp = R ps by symmetry).The above discussion of the eigenvectors of ∆ and Γ also explains why the matrix element ∆ ss alone is sufficientto predict the reflection coefficient R ss to a good approximation (see main text). By contrast, no such simpleexplanation is available for the p -polarization since it forms a superposition of w and e z . Nevertheless, for smallenough incident angles within the paraxial regime where e ± p ≈ ± k (cid:107) / | k (cid:107) | ≈ ± w , the p -polarization becomes anapproximate eigenvector of Γ and ∆. The plot of ∆ pp = e −† p ∆ e + p in Fig. S2d then exhibits a resonance (∆ pp ≈ δ )matching that of R pp within the paraxial regime, but fails to reproduce many qualitative features of R pp beyond it(see Fig. 3c of the main text). (c)(a) (b)(d) FIG. S2: Overlap between the two in-plane eigenvectors w and w of ∆ and one of the in-plane eigenvectors v of Γ. (a)indicates that w is approximately orthogonal to v for all k (cid:107) inside the light cone. Similarly, it is evident from (b) that w is approximately parallel to v . (c) Off-diagonal reflection coefficient R sp (due to symmetry, R sp = R ps ). (d) Matrix element∆ pp ( k (cid:107) ) = e −† p ∆( k (cid:107) ) e + p of the cooperative shift ∆( k (cid:107) ). All plots were computed for lattice constant a/λ = 0 .
4. DISORDERED ARRAYS
Disorder in atomic positions can be described by writing r = r n + δ r n for the position of an atom n , where r n are the ordered lattice positions and δ r n are a small random fluctuations. The effect on the self-consistent scatteringequation, Eq. (S6), is a perturbation δG nm of the Green’s function, G nm = G ( r n , r m ) = G ( r n , r m ) + δG nm , δG nm ≈ (cid:88) a = n,m (cid:88) i = x,y,z ∂G ( r n , r m ) ∂r ia δr ia + ... (S50)0with δr in = e i · δ r in . We then derive a formal perturbation theory for Eq. (S6), finding that to lowest order the effecton the effective polarizability is a correction to the cooperative resonance, δ ∆( k (cid:107) ) = − γλ N (cid:88) n (cid:88) m (cid:54) = n Re (cid:104) δG nm e − i k (cid:107) · ( r n − r m ) (cid:105) . (S51)For a simple estimation of the effect of disorder, we take statistically independent and identically-distributed fluctu-ations δr in with zero mean and variance δr , (cid:104) δr in (cid:105) = 0 , (cid:104) δr in δr jm (cid:105) = δr δ nm δ ij . (S52)Considering the Taylor expansion of δG nm from Eq. (S50) up to second order (valid for δr (cid:28) a ), we find e.g., that atnormal incidence, the average of the disorder-induced shift is given by (cid:104) δ ∆ (cid:105) ≈ π ( δr/λ ) ∆.For disorder due to impurities, we consider the case of an absence or addition of an atom in the array. We findnumerically (see Sec. below) that for small values of a/λ (e.g. 0 .
2) the results are almost unchanged compared withthose of a perfect lattice, whereas for a/λ larger than 1 / .
8) they may change considerably. This may supportthe requirement for a high filling factor found in [S12] for the case a/λ = 0 . δ = 0).
5. OPTICAL SATURATION OF THE ARRAY
Here we address the optical nonlinearity of the array by estimating the incident power needed to saturate individualatoms. We first evaluate he fraction of incident power that is absorbed (and then scattered) in each individual atomof the array during the scattering process at the collective resonance ( δ = ∆). This allows us to asses the number ofatoms that actually participate in the interaction and scattering of the incident field, hinting at the potential degreeof nonlinearity of the array (main text).We consider a normal incident Gaussian beam, E ( r ) = E e L w w ( z ) e ikz e − iϕ ( z ) e − x y (cid:48) w z ) e ik x y R ( z ) , (S53)with the usual beam parameters, w ( z ) = w (cid:115) (cid:18) zz R (cid:19) , z R = πw λ , R ( z ) = z (cid:20) (cid:16) z R z (cid:17) (cid:21) , ϕ ( z ) = arctan (cid:18) z (cid:48) z R (cid:19) . (S54)and where w is the beam waist at its focal point and e L its polarization. The power carried by the beam is W = π w c(cid:15) | E | . (S55)Considering the paraxial character of the incident Gaussian beam, namely, that | k (cid:107) | /k (cid:28) k (cid:107) contained in E , k (cid:107) , we have Γ( k (cid:107) ) ≈ Γ(0) and ∆( k (cid:107) ) ≈ ∆(0) [see e.g. Eq. (S22) and Fig. S2d here or Fig. 3c of the main themain text], so that α e ( k (cid:107) ) ≈ α e (0) = α e δ ij and Eq. (S13) becomes p ( k (cid:107) ) ≈ α e E , k (cid:107) ⇒ p n ≈ α e E ( r n ) (S56)Therefore, within the paraxial approximation around normal incidence, it appears as if each atom individually possessan effective polarizability α e , such that the power dissipated in an atom n (solely by radiation/scattering for the case γ nr = 0) at cooperative resonance, δ = ∆, is given by W n = W n x ,n y = 12 Im[ α e ] ω | E ( r n ) | = (cid:15) c | E | a e − a w ( n x + n y ) , (S57)where Γ + γ = γ [3 / (4 π )]( λ/a ) was used in the expression for α e at normal incidence (Eq. S34). It can be verifiedthat at the cooperative resonance all of the incident power is absorbed (scattered) by the array, by summing over allatoms ( n x , n y ) and arriving at (cid:80) n x (cid:80) n y W n x ,n y ≈ (cid:82) ∞−∞ dn x (cid:82) ∞−∞ dn y W n x ,n y = W with W from Eq. (S55).1The fraction of power absorbed by an individual atom is then given by P n = W n W = a ( π/ w e − a w ( n x + n y ) . (S58)Therefore, the atoms that are within the beam waist, namely, those which significantly participate in the interactionwith light and its scattering, absorb a fraction of P n = a / [( π/ w ] of the incoming power (per atom). By noticingthat the cross section (area) of the Gaussian beam scales as ∼ ( π/ w , the result for P n can be interpreted as a ratiobetween atomic and light cross sections, with the atomic effective cross section given by ∼ a . This interpretationis also supported by comparing the polarizabilities at resonance of bare and dressed atoms: from Eq. (S8) and Eq.(S34) we arrive at α ( δ = 0) ∝ λ and α e ( δ = ∆) ∝ λa which demonstrates how the individual cross section λ of abare atom is replaced by the effective cross section a for the renormalized/dressed atom.Since atoms are highly nonlinear, namely, they become transparent to light upon their full excitation (saturation),then the mirror effect may vanish for sufficiently high incident power W . In order to estimate the power required forthe saturation of the atoms, we note that an atom n is likely to be excited upon absorbing N = 1 /P n photons ata time < (Γ + γ ) − before it decays, setting a saturation power of W ≈ N (cid:126) ω (Γ + γ ). Taking, e.g., a = 0 . λ and w = 1 . λ we obtain Γ + γ ≈ γ and N ≈
14 for atoms within the beam waist, so that W ≈ (cid:126) ωγ , only a single orderof magnitude larger than the saturation power of a single atom coupled to a 1d photonic channel.
6. DIRECT NUMERICAL APPROACH TO THE SCATTERING PROBLEM
The problem of scattering of an electromagnetic field by a finite collection of point dipoles, generally formulated byEqs. (S5) and (S6) can be solved numerically by a simple matrix inversion. Introducing the vector notation E for the3 N -dimensional vector E ni = p ni /α (3 first entries are the 3 vector components of the E -field at the position of atom1, next 3 are the field components at atom 2, etc.) and the notation G for the 3 N × N matrix G nmij , the solution ofEq. (S6) for the local fields on the atoms, E = E i ( r n ), is E = (cid:20) − π α(cid:15) λ λ G (cid:21) − E . (S59)Given the collection of atomic positions that form the array r n (e.g. square lattice), the incident field E ( r ) (seebelow) and the polarizability α ( ω ) [Eq. S8], we can numerically invert the matrix in Eq. (S59). Then, we insert thesolutions E = E ni = E i ( r n ) into the right-hand side of Eq. (S5) and obtain the solution for the field at any givenpoint r , E i ( r ).In order to compare the results of our analytical theory to those obtained numerically in a more realistic scenario,the infinite array is replaced by a finite square lattice of N atoms whereas all incident plane waves are replaced byGaussian beams, E ( x (cid:48) , y (cid:48) , z (cid:48) ) = E e L w w ( z (cid:48) ) e ikz (cid:48) e − iϕ ( z (cid:48) ) e − x (cid:48) y (cid:48) w z (cid:48) ) e ik x (cid:48) y R ( z (cid:48) ) , (S60)with the beam parameters from Eq. (S54). Here the coordinates ( x (cid:48) , y (cid:48) , z (cid:48) ) are written in the reference frame of thebeam, namely, where the beam propagates along the z (cid:48) direction, e k = e z (cid:48) . For normal incident Gaussian beam wehave x = x (cid:48) , y = y (cid:48) , z = z (cid:48) and a polarization e L ∈ { e x , e y } , whereas for a beam propagating along the xz -plane atan angle θ from the z axis, we find x (cid:48) = x cos θ − z sin θ, y (cid:48) = y, z (cid:48) = z cos θ + z sin θ, (S61)and e L ∈ { e p , e s } , with e p = cos θ e x − sin θ e z , e s = − e y . (S62)For a faithful comparison with the infinite array case assumed in the theory, the cross section of the Gaussian beamat the position of the array, z = 0, has to be sufficiently smaller than the area of the array, N a . At normal incidence( θ = 0) we then always take w ≤ . a √ N , whereas for a general angle θ the waist w has to be smaller by roughlycos θ .2 −8 −6 −4 −2 0 2 4 6 8−1−0.500.51 z/ λ R e [ E x ] (c) incidentscattered −8 −6 −4 −2 0 2 4 6 8−1−0.500.511.5 z/ λ I m [ E x ] (d) incidentscattered−5 0 5−2.5 2.501234 z/ λ | E / E | (e) a/ λ =0.5, δ =0a/ λ =0.5, δ = ∆ a/ λ =0.2, δ =0 0 0.2 0.4 0.6 0.800.20.40.60.81 a/ λ T = | E / E | (f) analyticnumerical FIG. S3: Numerical approach, normal incidence. (a) x -polarized incident Gaussian beam at normal incidence to an array of N = 26 ×
26 atoms with a lattice constant a = 0 . λ at z = 0. The incident beam is resonant with the bare atom, δ = 0, andits waist is w = 0 . √ Na = 1 . λ . In agreement with Fig. 1c of the main text, no transmission is observed, and the reflectedwave forms a standing wave with the incident field. All fields intensities are in units of | E | , the peak intensity of the incidentGaussian beam, Eq. (S60). (b) Same as (a) for a = 0 . w = 0 . √ N × . . λ . Here it appears that some of the fieldis transmitted since as in Fig. 1c (main text). (c) A closer look at the scattering and interference processes for the case plottedin (a). Here we plot the incident and scattered field along z for a constant arbitrary value of x and y ( x = y = 0 . a ). The realpart of the scattered field is exactly opposite in sign with respect to the real part of the incident field thus exactly cancellingit (all fields in units of E ). (d) Same as (c) for the imaginary part of the field. For z >
0, beginning at a short distance afterevanescent fields have decayed, there exists again exact cancellation, so that no transmitted (real and imaginary) field exists.For z < a = 0 . λ where at δ = 0 a transmission of about ∼ . a = 0 . λ by changing the frequency of the incident field to match the cooperativeresonance δ = ∆, which is found for a = 0 . λ using Fig. 2a in the main text. The result is presented by the red dotted linewhich almost exactly coincides with the blue curve apart from the somewhat longer evanescent tail at z > a ). (f) Transmitted field as a function of lattice spacing a/λ . Repeating the calculation from (d) for different a values atthe individual-atomic resonance, δ = 0, the transmission coefficient is extracted. For each a value the Gaussian beam waist istaken to be w = 0 . √ na , i.e. smaller than the size of the array (blue dots). This nevertheless limits a/λ to be above 0 . w > λ/ FIG. S4: Numerical approach, oblique incidence. (a) Reflection of a p -polarized Gaussian beam at an incident angle of θ = 30 .The beam waist is w = 1 . λ and the array contains 40 ×
40 atoms with separation a = 0 .
2. In accordance with the theoreticalpredictions (Fig. 3 in the main text), almost perfect reflection is observed. Quantitative agreement to the theory is foundby examining the amplitude of the numerically calculated field along the x -axis at the far field, e.g. at z = − λ (noting theinterference between the incident and reflected fields closer to the surface). (b) Same as (a) for s -polarized beam at an incidentangle θ = 60 o . Due to the larger angle of incidence, a larger array of 70 ×
70 is taken, so that diffraction from the edges becomesvery small and agreement to the theory (Fig. 3) is observed.
Figs. S3 and S4 illustrate the excellent qualitative and quantitative agreement between this numerical solution andthe theoretical results presented in the main text, both for the normal incident and oblique incident cases, respectively. [S1] F. J. Garc´ıa de Abajo, Rev. Mod. Phys. , 1267 (2007).[S2] L. Novotny and B. Hecht, Principles of Nano-Optics (Cambridge University Press, 2006).[S3] D. P. Craig and T. Thirunamachandran,
Molecular Quantum Electrodynamics (Academic, London, 1984).[S4] P. Lambropoulos and D. Petrosyan,
Fundamentals of Quantum Optics and Quantum Information (Springer, 2007).[S5] The sum diverge due to the contribution of the Lamb shift whose counter term, the real part of G ij (0 , g ij we restore the real part of G ij (0 , m . Then, by regularizing both the sum andthe integral with the same, e.g. exponential/Gaussian cutoff function, the difference of the two converges numerically andyields the same results obtained by the direct summation in real space, Eq. (S12).[S6] M. O. Scully and M. S. Zubairy, Quantum Optics (Cambridge University Press, 1997).[S7] R. H. Lehmberg, Phys. Rev. A 2, 883 (1970).[S8] P. de Vries, D. V. van Coevorden, and A. Lagendijk, Rev. Mod. Phys. , 447 (1998).[S9] J. A. Klugkist, M. Mostovoy, and J. Knoester, Phys. Rev. Lett. , 1 (2006).[S10] M. Antezza and Y. Castin, Phys. Rev. A , 13816 (2009).[S11] S. Thongrattanasiri, F. H. L. Koppens and F. J. Garc´ıa de Abajo, Phys. Rev. Lett. , 047401 (2012).[S12] R. J. Bettles, S. A. Gardiner and C. S. Adams, Phys. Rev. Lett.116