Effective Field Theory for Rydberg Polaritons
M. J. Gullans, J. D. Thompson, Y. Wang, Q.-Y. Liang, V. Vuletic, M. D. Lukin, A. V. Gorshkov
EEffective Field Theory for Rydberg Polaritons
M. J. Gullans, J. D. Thompson, Y. Wang, Q.-Y. Liang, V. Vuleti´c, M. D. Lukin, and A. V. Gorshkov Joint Quantum Institute and Joint Center for Quantum Information and Computer Science,National Institute of Standards and Technology and University of Maryland, College Park, Maryland 20742, USA Physics Department, Massachusetts Institute of Technology, Cambridge, Massachusetts 02138, USA Physics Department, Harvard University, Cambridge, Massachusetts 02138, USA
We develop an effective field theory (EFT) to describe the few- and many-body propagation ofone dimensional Rydberg polaritons. We show that the photonic transmission through the Rydbergmedium can be found by mapping the propagation problem to a non-equilibrium quench, where therole of time and space are reversed. We include effective range corrections in the EFT and showthat they dominate the dynamics near scattering resonances in the presence of deep bound states.Finally, we show how the long-range nature of the Rydberg-Rydberg interactions induces strongeffective N -body interactions between Rydberg polaritons. These results pave the way towardsstudying non-perturbative effects in quantum field theories using Rydberg polaritons. Photons can be made to strongly interact by dress-ing them with atomic Rydberg states under conditionsof electromagnetic induced transparency (EIT) [1–3].Probing such Rydberg polaritons in the few-body limit,recent experiments were able to observe non-perturbativeeffects including the formation of bound states [4], single-photon blockade [5–7] and transistors [8–10], and two-photon phase gates [11]. Theoretical work on quantumnonlinear optics with Rydberg polaritons has focused ontwo-body effects or dilute systems [2–5, 12–17]; however,these theoretical methods often fail in dense systems withmore than two photons.Effective field theory (EFT) aims to describe low en-ergy physics without resorting to a microscopic model atshort distances or high energies [18]. In few-body sys-tems, it is a useful approach to describe particle scat-tering and bound states when the momentum k involvedis much less than the inverse range of the interactions[18, 19] At the two-body level, the EFT depends onlyon the scattering length a . For scattering at momenta ka (cid:28)
1, one can solve the EFT perturbatively [19, 20].However, describing unitarity ( a → ±∞ ) or bound statesrequires inclusion of all orders in perturbation theory,which can be re-summed, provided the EFT parametersare properly renormalized [21, 22].In this Letter, we develop an EFT to describe thefew- and many-body transmission of photons through adispersive, one dimensional Rydberg polariton medium.We first consider the renormalized theory, which dependsonly on the local two-body scattering length, the effec-tive mass, and the group velocity of the Rydberg po-laritons. By switching the role of time and space inthe Lagrangian, we map the transmission problem toa non-equilibrium quench, which greatly simplifies thedescription of the dynamics. We then consider correc-tions to the EFT arising from the long range natureof the Rydberg interactions and the corrections to themassive dispersion. We evaluate the so-called “effectiverange corrections” to the EFT and show that they domi-nate the dynamics near unitarity in the presence of deep bound states. We then find the non-perturbative solutionfor the many-body Rydberg polariton problem at largemomenta. Integrating out this momentum scale leadsto strong N -body interactions, which appear as contactforces in the EFT.A schematic of a Rydberg polariton transmission ex-periment is shown in Fig. 1(a) [14]. A spatially inhomo-geneous atomic cloud is probed with a classical controlfield, with frequency ω c and Rabi frequency Ω, and afew-photon probe beam focused into a 1d channel. Thecontrol and probe beams are configured for EIT on a two-photon resonance from the ground state | g (cid:105) to a Rydbergstate | s (cid:105) via an intermediate state | p (cid:105) . We analyze thedispersive limit where the detuning of the control field Scattering Length ! … ! N ! N ! … Rydberg mediumControl field µ m) V (r) / h ( M H z ) gp s ⌦ (a)(b) Interaction potential V ( r ) r b r ( µ m) V ( r ) / h ( M H z ) ˆ z z/r b (c) n ( z ) /n (0)
10 5 . . . . ⌫ = ! + . . . + ! N ⌫ = ! + . . . + ! N r b /a ( z ) FIG. 1: (a) Rydberg polariton transmission experiment: anatomic cloud is probed with a few-photon beam, focused intoa 1d channel, and a classical control field. Under dispersiveconditions, the total energy (cid:126) ν and number N of probe pho-tons are conserved. (b) Interaction potential V ( r ) = C /r for two 100S / Rydberg states in Rb. The range is givenby the blockade radius r b . (Inset) Level diagram of an in-teracting atom for different r . (c) Dimensionless density n ( z ) /n (0) = exp( − z / σ ax ) and inverse scattering length r b /a ( z ) for σ ax = 36 µ m, r b = 18 µ m, OD = 25, Ω / π =5 MHz, and ∆ / π = 20 MHz. a r X i v : . [ phy s i c s . a t o m - ph ] A ug ∆ = ω ps − ω c is much greater than the p -state halfwidth γ ; here ω ab is the atomic transition frequency from | a (cid:105) to | b (cid:105) . For large enough atomic density n ( z ), the probephotons transform into Rydberg polaritons upon enteringthe medium because the collective, single-photon Rabifrequency of the probe g c ( z ) = [6 πγ c n ( z ) /ω gp ] / ismuch greater than Ω [23]. We use the dimensionless mea-sure of the density given by the resonant optical depthOD = (cid:82) dz [ g c ( z )] / γc .Consider two Rydberg atoms interacting through thevan der Waals potential V ( r ) = C /r . This interactionis strong enough that a single Rydberg excitation modi-fies the optical response over a region large compared tothe optical wavelength (see Fig. 1(b) and inset). The sizeof this region is given by the blockade radius r b , definedby the condition that V ( r b ) is equal to the off-resonantEIT linewidth 2Ω / | ∆ | [14, 24].To see how these effects lead to strong photon-photoninteractions, one can use a gedanken experiment whereone photon (polariton) is held at fixed position z , thenany photon that passes by will pick up a nonlinear phaseshift ϕ ( z ) ≈ [ g c ( z )] r b /c ∆; here we assume g c ( z ) variesslowly over r b . For atomic densities achievable with laser-cooled atoms ( n (cid:38) cm − ), this nonlinear phaseshift can be a sizable fraction of π [4]. An alterna-tive metric is the two-body scattering length a , whichwas mapped out for Rydberg polaritons in a uniformmedium in Ref. [12]. For an inhomogeneous medium, wecan similarly define a local scattering length a ( z ). Forsmall ϕ ( z ), these two metrics are closely related because a ( z ) ≈ (3 /π ) r b / [ ϕ ( z )] . We show in the supplementalmaterial that a ( z ) is well defined when the density variesslowly over r b [25]. Figure 1(c) shows a ( z ) calculatedfor a Gaussian density profile with parameters similar torecent experiments [4].In the absence of interactions, the propagation of Ry-dberg polaritons is captured by the local EIT dispersionrelation [23, 25] q ( ω, z ) = ωc (cid:18) g c ( z )] Ω − ω (∆ + ω ) (cid:19) , (1)where ω = ω (cid:96) − ω is the detuning of the probe fre-quency ω (cid:96) from the two-photon resonance ω = ω gs − ω c .The electric field of the probe evolves as E ( ω, z ) = E ( ω, z ) exp (cid:2) iω ( z − z ) /c + i (cid:82) zz dz (cid:48) q ( ω, z (cid:48) ) (cid:3) . For a suf-ficiently slowly varying density, we can define a localgroup velocity v g ( z ) = c/ (1 + [ g c ( z )] / Ω ) and mass m ( z ) = − (cid:126) Ω / v g ( z )] by solving Eq. (1) for ω andexpanding near q = 0: ω ≈ v g ( z ) q + (cid:126) q / m ( z ) [23, 26].For non-relativistic bosons in 1d, the only interactionterm that is relevant under renormalization is the two-body contact interaction [21]. As a result, the renormal-ized Lagrangian density for Rydberg polaritons is L = ˆ ψ † (cid:34) i (cid:126) ∂ t − i (cid:126) v g ( z ) ∂ z − (cid:126) ∂ z m ( z ) (cid:35) ˆ ψ − (cid:126) ˆ ψ † ˆ ψ m ( z ) a ( z ) , (2) where [ ˆ ψ ( t, z ) , ˆ ψ † ( t, z (cid:48) )] = δ ( z − z (cid:48) ) and ˆ ψ is a singlecomponent field because there is only a single polaritonbranch near the two-photon resonance [25]. Outside themedium, ˆ ψ is the quantum field for the probe photons,while inside it corresponds to the Rydberg polariton field.The scaling of the contact interaction as 1 /a is the uni-versal behavior for bosons in 1d, in contrast to higherdimensions where it scales as a [27].Despite its relative simplicity compared to the micro-scopic model [25], the theory is still difficult to solve be-cause it has z -dependent parameters combined with sec-ond derivatives in z . To overcome this we define a newEFT with time and space exchanged via the local trans-formation ( t, z ) → ( z/v g ( z ) , tv g ( z )). Similar transforma-tions have been used to study propagation of quantumlight in nonlinear optical fibers [28, 29]. For the steadystate transmission with a uniform density, this transfor-mation is equivalent to the rotated boundary conditionsused in Ref. [4]. The resulting EFT is L = ˆ ψ † (cid:34) i (cid:126) v g ( z ) ∂ z − i (cid:126) ∂ t − (cid:126) ∂ t m ( z )[ v g ( z )] (cid:35) ˆ ψ − (cid:126) ˆ ψ † ˆ ψ m ( z ) a ( z ) v g ( z ) , (3)where [ ˆ ψ ( z, t ) , ˆ ψ † ( z, t (cid:48) )] = δ ( t − t (cid:48) ). Up to higher orderderivatives in t (which can be neglected under renormal-ization), Eq. (3) is equivalent to Eq. (2); however, the sec-ond derivative is now in t rather than z , which makes iteasier to account for the z -dependence of the parameters.In particular, Eq. (3) gives rise to propagation equationsakin to a time-dependent Schr¨odinger equation − i (cid:126) v g ( z ) ∂ z ˆ ψ ( z, t ) = (cid:90) dt (cid:48) [ H ( z, t (cid:48) ) , ˆ ψ ( z, t )] , (4)where H ( z, t ) is given by the last three terms in Eq. (3).In the dispersive regime, this propagation equation con-serves the total photon number N , which simplifies thetransmission problem. Benchmarking the EFT.—
We now compare the pre-dictions of the renormalized EFT for the two-photontransmission through a finite Rydberg medium with nu-merical simulations [5] of the exact wavefunction prop-agation. We decompose the two-photon wavefunctionat the exit of the medium ( z = L ) as ψ ( L, t , t ) = (cid:112) g (2) ( τ ) e iφ ( τ )+ i φ , where t are the time coordinatesof the two photons and τ = t − t is the relative time.The probability density g (2) ( τ ) can be measured in two-photon coincidence measurements of the output light fora weak coherent state input [5], while the nonlinear phase φ ( τ ) is defined relative to phase of the non-interactingmedium with a single-photon phase shift φ [4].The results are shown in Fig. 2 for a representativeset of parameters similar to Ref. [4]. We take a steady-state probe input on two-photon resonance ( ω = 0) with − − τ / τ d P h a s e ( R a d ) SimulationEFT: Non − uniform densityEFT: Uniform density ( ⌧ ) Simulation Non-uniform densityUniform density ⌧ /⌧ d . . . ( ⌧ ) A m p lit ud e . . . g ( ) ( ⌧ ) ⌧ /⌧ d (b)(a) FIG. 2: (a) Photon correlation function g (2) ( τ ) and (b) phase φ ( τ ) (in radians) of transmitted two-photon state calculatedusing EFT (solid lines) and numerical simulations (circles).We took Ω / π = 5 MHz, ∆ / π = 20 MHz, r b = 10 µ m,OD = 10, σ ax = 36 µ m, L = 4 σ ax , and L (cid:48) = 2 . σ ax . To aidcomparison we neglect decay from the | p (cid:105) and | s (cid:105) states. a Gaussian density profile n ( z ) ∝ exp[ − ( z − L/ / σ ax ]with a cutoff at the entrance to ( z = 0) and exit from( z = L ) the medium. We compare g (2) ( τ ) and φ ( τ )found with three different methods: numerical simula-tions, EFT with no free parameters, and EFT with a uni-form density with g c a free parameter and medium length L (cid:48) chosen to match the time delay τ d = (cid:82) L dz [1 /v g ( z )].For an intermediate time window, we see that bothEFT results capture many of the qualitative features ofthe simulations, but the inhomogeneous EFT capturesmore features and obtains better quantitative agreement.We can understand the deviations at long and short timesas follows. The long-time deviations arise because theEFT has a low momentum cutoff associated with spatialvariations in the density profile [25]. For a Gaussian oruniform density profile, this scale is given by 1 /L , withthe associated low-frequency cutoff 1 /τ d . The short-timedeviations arise from corrections to the EFT associatedwith: our use of a massive polariton dispersion, the swapof time and space, and the finite interaction range. Thefirst two effects contribute on timescales shorter than τ m ≈ max(∆ / Ω , / ∆), while the effect of the finite inter-action range appears on timescales less than r b /v g . Forthe parameters in Fig. (2), r b /v g , τ m (cid:28) τ d , which is con-sistent with the good agreement we find at intermediatetimes r b /v g , τ m (cid:46) τ (cid:46) τ d .In related work, we have shown that this renormalizedEFT also gives good agreement with numerical simula-tions of the three-photon transmission [30]. Yet, for in-creasing N , simulations of the full transmission becomeintractable and it is natural to ask: what are the lead-ing corrections to the theory? In the framework of EFT,these corrections can be found systematically by evaluat-ing higher order corrections in kr b . We show below thatthe terms in this expansion arise from two intertwinedeffects: (i) the finite range of the interactions and (ii) de-viations of the dispersion from that of a massive particle. Effective range corrections.—
A standard approach toinclude finite range effects for massive particles is throughthe effective range expansion. In this treatment, higher order corrections to the scattering phase shift δ ( k ) aretaken into account [19]. For bosons in 1d, the expansiontakes the form [31] k tan δ ( k ) = 1 a + r k . . . , (5)where r is the so-called “effective range” parameter.These corrections can be included in the EFT by addingterms to the Lagrangian that contain higher derivativesin ˆ ψ , e.g., (after switching time and space) L → L + C ˆ ψ † ( ∂ t ˆ ψ † )( ∂ t ˆ ψ ) ˆ ψ, (6)where C = (cid:126) r / mv g is fixed by Eq. (5) [32]. Includingthese terms extends the validity of the EFT to higher po-lariton densities. Most notably, this approach allows oneto study unitarity in the presence of deep bound states,which occur when ϕ (cid:29)
1. In this regime, we can solvefor a and r analytically [33, 34], and we find that thetwo-body contact vanishes near a scattering resonance,but r ≈ . √ ϕ r b remains finite [25].Scattering resonances associated with the appearanceof additional two-body bound states can be achieved forRydberg polaritons at sufficiently high atomic density[12]. Current experiments, however, are limited to den-sities such that only a single two-body bound state ispresent. In this case, we find r ≈ (2 / r b /a and thesecorrections are suppressed. We now show that the dom-inant corrections to the theory in this regime arise fromeffective 3-body interactions. N -body interactions.— The strong long-range Rydberginteractions that result in blockade are also expected toinduce large effective N -body interactions [35, 36]. Thisis illustrated at the three-body level because, when twopolaritons are less than r b from a Rydberg atom, theydo not interact with each other. As a result, one expectsa three-body force of the same magnitude and oppositesign as the two-body force.More formally, effective N -body interactions emergefrom integrating out virtual processes with high energyor large momenta. In the case of Rydberg polaritons,this can be done in a surprisingly straightforward man-ner because the theory dramatically simplifies at largemomenta. In particular, the single-body propagator forthe Rydberg polaritons projected onto the s -states g ss saturates to a constant (see supplemental material [25])lim q →∞ (cid:126) g ss ( q, ν ) = (cid:126) χ ( ν ) = (∆ + ν )(∆ + ν ) ν − Ω , (7)for momentum q (cid:29) /v g τ m . The physical origin ofthis can be seen in Eq. (1), where the local momen-tum q ( ω, z ) diverges at the Raman resonance conditionsΩ − ω (∆ + ω ) = 0 . These Raman excitations have afrequency close to two-photon resonance and effectivelyinfinite mass; therefore, near ν = 0, they dominate vir-tual processes with large internal momentum. In the con-text of EFT, these effects can be included by adding two N - B ody I n t e r ac ti on (e)(a) Connected Disconnected (c) E E E E (b)(d) r/r b V ( r, , , V ( r, , V ( r, FIG. 3: Examples of a (a) connected and (b) disconnectedscattering diagram for N = 5. The diagram in (a) con-tributes to V N eff , while (b) does not. (c,d) Graph represen-tations of (a,b), respectively, neglecting the ordering of scat-tering events. The graph in (c) is a tree graph, which impliesthat (a) is a lowest order diagram for V N eff . (e) Cut of the non-perturbative solution for V N eff in units of χ − N up to N = 4. fictitious, infinitely massive particles to the theory associ-ated with the Raman resonances [25]. Due to their high-energy, the Rydberg interactions can only excite these“particles” virtually. Integrating them out of the theoryresults in effective N -body interactions for the ˆ ψ field.The associated N -body interaction potential V N eff canbe found by accounting for all of the virtual processeswhere N of these fictitious particles exchange momen-tum. These contributions to the scattering amplitudesare represented by connected diagrams of the type shownin Fig. 3(a), where the particles cannot be broken intodisjoint clusters. Particles are connected by two-body in-teractions (curly lines), with the insertion of the N -bodypropagator in between (vertical lines). Figure 3(b) showsan example of a disconnected diagram, which is separa-ble into two disjoint clusters (12) and (345) and does notcontribute to V N eff .Integral equations for the connected contributions tomulti-particle scattering amplitudes were first formulatedby Weinberg [37] and Rosenberg [38]. The full integralequations have only been solved for N ≤ V N eff : V N eff ( z ; ν ) ≈ ( N − χ N ( ν )] N − (cid:88) T ( N,E ) V E ...V E N − , (8)where E k = ( i k , j k ) denotes a particle pair, V E k = V ( z i k − z j k ), and the sum is over all labeled treegraphs T ( N, E ) with N vertices and N − E = { E , ... , E N − } . Here the N -body propagator is (cid:126) χ N ( ν ) = (cid:126) (cid:82) dωχ ( ω ) χ N − ( ν − ω ) ≈ ( ν − N Ω / ∆) − (for Ω (cid:28) ∆). If r (cid:28) r b , then | V ( r ) χ N | >
1, and theperturbative approach of Eq. (8) is no longer valid. Wederive the non-perturbative solution for V N eff in the sup-plemental material [25]. Figure 3(e) shows a cut of thissolution up to N = 4. Consistent with the blockade ef-fects described above, we see that V has the oppositesign from V . More generally, we find V N eff alternateswith N between attraction and repulsion [25].During low-momenta processes kr (cid:28)
1, the polari-tons hardly probe the blockaded region of the potential.In this case, we can replace V ( r ) with the renormalizedinteraction U ( r ) = − (2 (cid:126) /ma ) δ ( r ) and apply the pertur-bative result from Eq. (8) to find the N -body interactions[40]. After switching time and space, the resulting EFTis governed by Eq. (4) with the Hamiltonian density H = ˆ ψ † (cid:20) − i (cid:126) ∂ t − (cid:126) ∂ t mv g (cid:35) ˆ ψ + (cid:88) N h N ˆ ψ † N ˆ ψ N , (9) h N = ( − N − N (cid:18) (cid:126) mav g (cid:19) N − ( N χ N ) N − , (10)where the two-body interaction h is the same as inEq. (3) and we used Cayley’s tree formula N N − for thenumber of labeled tree graphs with N vertices in eval-uating Eq. (8) [41]. Using approximate expressions for a near a scattering resonance [25], we find the genericscaling h N ∼ ( r b /a ) N − .To determine the importance of the N -body interac-tions for non-perturbative effects in the EFT, h N shouldbe compared with the effective range corrections at themomentum scale k ∼ /a . For large ϕ , r ∼ r b and,from Eq. (6), we see that the effective range correctionscontribute at the same order as h ∼ r b /a . On theother hand, for ϕ (cid:28) r is suppressed by an additionalpower of r b /a and the effective range corrections scaleas h ∼ r b /a . Thus, for weak interactions, we find thesurprising result that the 3-body force dominates the cor-rections to the theory for all momentum scales k (cid:46) /a .The nature of these corrections has important impli-cations for the propagation dynamics of 1d Rydberg po-laritons. The largest corrections to the theory will de-termine the deviations from the universal predictions forthe shallow bound state clusters when a > a < Conclusion.—
We developed an EFT to describe thefew- and many-body propagation of 1d Rydberg polari-tons. The broad applicability of EFT to describe thesesystems opens a new perspective on the design of ex-periments aimed at probing non-perturbative effects inquantum field theories using Rydberg polaritons. In par-ticular, Rydberg polariton experiments can involve com-plex geometries [46] or more Rydberg levels [11], dimen-sions [47], and external control fields [48]. The theoret-ical methods developed here can be naturally extendedto these more complex configurations. For example, us-ing additional control fields or atomic levels to modify χ N ( ν ) would allow precise control over the range andstrength of the N -body potentials. This could be usedto realize exotic situations where, e.g., a single M -bodyforce dominates over all N -body forces with N (cid:54) = M . Asanother example, accounting for light diffraction intro-duces 3d effects, where EFT predicts the emergence ofan Efimov effect in the vicinity of a scattering resonance[19, 49]. Further extending these theoretical methods toinclude dissipative interactions of the type demonstratedin Ref. [5] may uncover new universality classes for few-body physics, as well as new phases of non-equilibrium,strongly-correlated light and matter. Note added.—
During completion of this work, we be-came aware of related work on the three-body problemfor Rydberg polaritons [50].
Acknowledgements.—
We thank I. Carusotto, S. Diehl,P. Julienne, B. Ruzic, O. Firstenberg, M. Maghrebi, andR. Qi for helpful discussions. This research was sup-ported in part by the Kavli Institute for TheoreticalPhysics through the NSF under Grant No. NSF PHY11-25915, the NSF PFC at the JQI, the Harvard-MIT CUA,ARL CDQI, NSF QIS, AFOSR, and ARO. [1] J. D. Pritchard, D. Maxwell, A. Gauguet, K. J. Weath-erill, M. P. A. Jones, and C. S. Adams, Phys. Rev. Lett. , 193603 (2010).[2] D. Petrosyan, J. Otterbach, and M. Fleischhauer, Phys.Rev. Lett. , 213601 (2011).[3] A. V. Gorshkov, J. Otterbach, M. Fleischhauer, T. Pohl,and M. D. Lukin, Phys. Rev. Lett. , 133602 (2011).[4] O. Firstenberg, T. Peyronel, Q.-Y. Liang, A. V. Gor-shkov, M. D. Lukin, and V. Vuleti´c, Nature , 71(2013).[5] T. Peyronel, O. Firstenberg, Q.-Y. Liang, S. Hofferberth,A. V. Gorshkov, T. Pohl, M. D. Lukin, and V. Vuletic,Nature , 57 (2012).[6] Y. O. Dudin and A. Kuzmich, Science , 887 (2012).[7] D. Maxwell, D. J. Szwer, D. Paredes-Barato, H. Busche,J. D. Pritchard, A. Gauguet, K. J. Weatherill, M. P. A.Jones, and C. S. Adams, Phys. Rev. Lett. , 103001(2013).[8] H. Gorniaczyk, C. Tresp, J. Schmidt, H. Fedder, andS. Hofferberth, Phys. Rev. Lett. , 053601 (2014).[9] D. Tiarks, S. Baur, K. Schneider, S. D¨urr, and G. Rempe,Phys. Rev. Lett. , 053602 (2014).[10] H. Gorniaczyk, C. Tresp, P. Bienias, A. Paris-Mandoki,W. Li, I. Mirgorodskiy, H. P. B¨uchler, I. Lesanovsky, andS. Hofferberth, arXiv:1511.09445 (2015).[11] D. Tiarks, S. Schmidt, G. Rempe, and S. D¨urr, Sci. Adv. , e1600036 (2016).[12] P. Bienias, S. Choi, O. Firstenberg, M. F. Maghrebi,M. Gullans, M. D. Lukin, A. V. Gorshkov, and H. P.B¨uchler, Phys. Rev. A , 053804 (2014).[13] M. F. Maghrebi, M. J. Gullans, P. Bienias, S. Choi, I. Martin, O. Firstenberg, M. D. Lukin, H. P. B¨uchler,and A. V. Gorshkov, Phys. Rev. Lett. , 123601(2015).[14] O. Firstenberg, C. S. Adams, and S. Hofferberth,arXiv:1602.06117.[15] J. Otterbach, M. Moos, D. Muth, and M. Fleischhauer,Phys. Rev. Lett. , 113001 (2013).[16] M. Moos, M. H¨oning, R. Unanyan, and M. Fleischhauer,Phys. Rev. A , 053846 (2015).[17] A. Grankin, E. Brion, E. Bimbard, R. Boddeda, I. Us-mani, A. Ourjoumtsev, and P. Grangier, Phys. Rev. A , 043841 (2015).[18] See, e.g., H. Georgi, Ann. Rev. Part. Sci. , 209(1994).; D. B. Kaplan, arXiv:nucl-th/0510023 (2005).; C.P. Burgress, Ann. Rev. Part. Sci., , 329 (2007).; H. W.Hammer, J. Phys. G: Nucl. Part. Phys. 31, S1253 (2005).[19] E. Braaten and H. W. Hammer, Phys. Rep. , 259(2006).[20] E. Braaten and A. Nieto, Phys. Rev. B , 14 745 (1997);Phys. Rev. B , 8090 (1997).[21] S. K. Adhikari, T. Frederico, and I. D. Goldman, Phys.Rev. Lett. , 487 (1995).[22] P. F. Bedaque, H.-W. Hammer, and U. van Kolck, Phys.Rev. Lett. , 463 (1999).[23] M. Fleischhauer, A. Imamoglu, and J. P. Marangos, Rev.Mod. Phys. , 633 (2005).[24] More precisely, we define r b via V ( r b ) = 1 / | χ (0) | =2 (cid:126) Ω | ∆ | / | Ω − ∆ | [12]. For our parameters, this defini-tion is nearly equivalent to the one in the main text.[25] See Supplemental Material for a derivation of the spa-tially varying scattering length, effective range correc-tions, EFT including the Raman resonances, and thenon-perturbative solution for V N eff .[26] In the definition of m ( z ) below Eq. (1), we neglected asmall correction of order Ω / [ g c ( z )] .[27] M. Olshanii, Phys. Rev. Lett. , 938 (1998).[28] Y. Lai and H. A. Haus, Phys. Rev. A , 844 (1989).[29] P.-´E. Larr´e and I. Carusotto, Phys. Rev. A , 043802(2015).[30] M. J. Gullans et al. in preparation.[31] V. E. Barlette, M. M. Leite, and S. K. Adhikari, Eur. J.Phys. , 435 (2000).[32] U. van Kolck, Nucl. Phys. A , 273 (1999).[33] G. F. Gribakin and V. V. Flambaum, Phys. Rev. A ,546 (1993).[34] V. V. Flambaum, G. F. Gribakin, and C. Harabati, Phys.Rev. A , 1998 (1999).[35] H. P. B¨uchler, A. Micheli, and P. Zoller, Nat Phys , 726(2007).[36] H. Weimer, M. M¨uller, I. Lesanovsky, P. Zoller, andH. B¨uchler, Nature Phys. , 382 (2010).[37] S. Weinberg, Phys. Rev. , B232 (1964).[38] L. Rosenberg, Phys. Rev. , B217 (1965).[39] S. K. Adhikari and K. L. Kowalski, Dynamical CollisionTheory and Its Applications (Academic Press, San Diego,CA, 1991).[40] Here we are assuming the delta-function is regularized bya high-momentum cutoff Λ satisfying 1 /r b (cid:28) Λ (cid:28) k .[41] R. P. Stanley, Enumerative Combinatorics, Volume 1 (Cambridge Univ. Press, Cambridge, 2011), 2nd ed.[42] J. B. McGuire, J. Math. Phys. , 622 (1964).[43] E. H. Lieb and W. Liniger, Phys. Rev. , 1605 (1963).[44] L. Pricoupenko, Phys. Rev. Lett. , 170404 (2008).45] A. Imambekov, A. A. Lukyanov, L. I. Glazman, andV. Gritsev, Phys. Rev. Lett. , 040402 (2010).[46] A. Sommer, H. P. B¨uchler, and J. Simon,arXiv:1506.00341.[47] S. Sevin¸cli, N. Henkel, C. Ates, and T. Pohl, Phys. Rev.Let. , 153001 (2011).[48] D. Maxwell, D. J. Szwer, D. Paredes-Barato, H. Busche, J. D. Pritchard, A. Gauguet, M. P. A. Jones, and C. S.Adams, Phys. Rev. A , 043827 (2014).[49] V. N. Efimov, Sov. J. Nucl. Phys. , 589 (1971); Phys.Rev. C , 1876 (1993).[50] K. Jachymski, P. Bienias, and H. P. B¨uchler, Phys. Rev.Lett. , 053601 (2016). upplemental Material to the Manuscript: “Effective Field Theory for RydbergPolaritons” Contents
I. Microscopic Model II. Effective Range Corrections
21. Weak Attractive Interactions 22. Strong Attractive Interactions 3
III. EFT Including Raman Resonances IV. Non-Perturbative N -body InteractionPotential References I. MICROSCOPIC MODEL
In this section, we derive the two-body Lippmann-Schwinger equation for an inhomogeneous density andshow that the spatially varying scattering length a ( z ) iswell defined when the density varies slowly compared tothe blockade radius.The effective Hamiltonian, including decay, that de-scribes the Rydberg polariton system is ( ~ = 1) H = − ic Z dz ˆ E † ( z ) ∂ z ˆ E ( z ) (S1) − Z dz g c ( z ) (cid:2) ˆ P ( z ) ˆ E † ( z ) + h.c. (cid:3) + H p + H int ,H p = − Z dz (∆ + iγ ) ˆ P † ( z ) ˆ P ( z ) (S2) − iγ s ˆ S † ( z ) ˆ S ( z ) − Ω[ ˆ P † ( z ) ˆ S ( z ) + h.c. ] ,H int = Z dzdz V ( z − z ) ˆ S † ( z ) ˆ S † ( z ) ˆ S ( z ) ˆ S ( z ) , (S3)where ˆ E ( z ), ˆ P ( z ), and ˆ S ( z ) are bosonic annihilation op-erators for a photon, excited atom, and Rydberg state atposition z . They satisfy [ ˆ E ( z ) , ˆ E † ( z )] = [ ˆ P ( z ) , ˆ P † ( z )] =[ ˆ S ( z ) , ˆ S † ( z )] = δ ( z − z ). The parameters ∆, Ω, γ , and g c ( z ) are defined in the main text and γ s is the halfwidthof the s -state.For an inhomogeneous medium it is convenient to solvethe scattering problem in real space. For a single polari-ton we can find the propagator at frequency ω from theequations of motion − iωE ( z ) = − c∂ z E + ig c ( z ) P ( z ) , (S4) − iωP ( z ) = − ( γ − i ∆) P ( z ) + ig c ( z ) E ( z ) + i Ω S ( z ) , (S5) − iωS ( z ) = − γ s S ( z ) + i Ω P ( z ) . (S6) The solution is given by E ω ( z, z ) = 1 N ω ( z ) exp (cid:20) i Z zz dz q ( ω, z ) (cid:21) , (S7) q ( ω, z ) = ωc (cid:18) − [ g c ( z )] ˜∆˜ δ (cid:19) , (S8) P ω ( z, z ) = − g c ( z )˜∆ (cid:18) ˜∆˜ δ (cid:19) E ω ( z, z ) , (S9) S ω ( z, z ) = g c ( z )Ω˜∆˜ δ E ω ( z, z ) , (S10)where ˜∆ = ∆ + ω + iγ , ˜ δ = − Ω / ˜∆ + ω + iγ s , N ω ( z )is a normalization constant chosen to satisfy | E ω ( z ) | + | P ω ( z ) | + | S ω ( z ) | = 1, and q ( ω, z ) is defined in Eq. (1)in the main text for γ = γ s = 0.To each of these wavefunctions, we associate the oper-ator ˆ ψ † ω ( z ) = Z dz p ρ ( ω, z ) (cid:2) E ω ( z, z ) ˆ E † ( z )+ P ω ( z, z ) ˆ P † ( z ) + S ω ( z, z ) ˆ S † ( z ) (cid:3) . (S11)where ρ ( ω, z ) = dq ( ω, z ) /dω is the local density ofstates. For an infinite, homogeneous medium this opera-tor creates a dark-state polariton [S1]. For a sufficientlyslowly varying density, this field becomes approximatelybosonic with the commutation relation[ ˆ ψ ω ( z ) , ˆ ψ † ω ( z )] = Z dzρ ( ω, z ) e i R zz dz [ q ( ω,z ) − q ( ω ,z )] ≈ ρ ( ω, z ) δ [ q ( ω, z ) − q ( ω , z )] = δ ( ω − ω ) . (S12)When ω is on the Raman resonance ˜∆˜ δ = 0, thesesolutions need to be treated with care because E ω → δ = ( ω + iγ s )(∆ + ω + iγ ) − Ω = 0 (S13)has the solutions ω ± = − ∆ + i ( γ + γ s )2 ± r [∆ + i ( γ + γ s )] + ( γ − i ∆) γ s . (S14)We work in the limit of large ∆ and small γ s , where theseresonances can be treated as if they were on the real axis.For these two eigenstates, ˆ ψ †± ( z ) = P ± ˆ P † ( z )+ S ± ˆ S † ( z )with P ± = ω ± q Ω + ω ± , (S15) S ± = Ω q Ω + ω ± . (S16) = + z z z z ⌫ ! ! z z z z z z z z z z z z T ( ⌫ ) T ( ⌫ ) FIG. S1: Diagrammatic representation of the two-bodyLippmann-Schwinger equation in real space, ν is the totalfrequency of the two polaritons, dots are V ( z − z ), and thelines are the single-polariton propagator g ss ( z, z , ν ). To solve the interacting problem, Eq. (S3) implies thatwe only need the propagator projected onto the Rydbergstates. We can use the eigenstates ψ † ω ( z ) | i to find thepropagator in a vicinity of z g ss ( z, z , ν ) = h ˆ S ( z ) 1 ν − H + i + ˆ S † ( z ) i (S17)= Z dω h ˆ S ( z ) ψ † ω ( z ) ih ψ ω ( z ) ˆ S † ( z ) i ν − ω + i + = ρ ( ν, z ) S ν ( z, z ) S ∗ ν ( z , z ) + χ ( ν ) δ ( z − z ) ,χ ( ν ) = X s = ± | S s | ν − ω s = ∆ + ν (∆ + ν ) ν − Ω . (S18)As discussed in the main text, the χ ( ν ) contributionarises from the two Raman resonances and accounts forthe saturation of the propagator at large momentum.With this representation of the single particle prop-agator, we can now write down the explicit Lippmann-Schwinger equation for the transition matrix for Rydbergpolaritons in an inhomogeneous medium, represented di-agrammatically in Fig. S1, T ( z , z , ν ) = V ( z − z ) (cid:2) δ ( z − z ) (S19)+ Z d z g ss,ss ( z , z , ν ) T ( z , z , ν ) (cid:3) ,g ss,ss ( z , z , ν ) = Z dωg ss ( z , z , ω ) g ss ( z , z , ν − ω ) . (S20)To find the EFT parameters, we replace g ss ( z , z , ω )with the first term in the last line of Eq. (S17) andapproximate the dispersion by the formula given be-low Eq. (1) in the main text ω = v g ( z ) q + q / m ( z ).We then replace V ( z − z ) with the pseuodopotential − ~ δ ( r ) /m ( z ) a ( z ), where z = ( z + z ) / r = z − z and derive an effective solution for T . This can thenbe matched to the asymptotic solution for the exact T .When the propagator varies slowly with z compared tothe potential (i.e., the density changes slowly on the scaleof r b ), we can find T using the anlysis of Ref. [S2] fora uniform density. In the next section, we discuss sev-eral regimes where a can be found analytically with thisapproach. II. EFFECTIVE RANGE CORRECTIONS
In this section we give approximate formulas for thescattering length a and effective range parameter r ateach scattering resonance for weak and strong attractiveinteractions.These parameters can be found by solving the mi-croscopic two-body problem using Eq. (S19). Assum-ing the density is slowly varying, we can make the localdensity approximation discussed in the main text andeliminate the center of mass momentum from Eq. (S19).In Ref. [S2], it was shown that the solution to the re-sulting integral equation can be found by solving a 1dSchr¨odinger-like equation − m ∂ r ψ + V ( r ) ψ = νψ, (S21)where r = z − z is the relative distance between thetwo polaritons and V ( r ) = V ( r ) / [1 − χ V ( r )] is shownin Fig. 3(e) in the main text. For χ C < r < r b , is approximately flat, while for r > r b itdecays as 1 /r . The relative importance of the core ver-sus the tail of the potential can be determined by com-paring the associated length scales: the blockade radius r b versus the van der Waals length r vdw = ( mC ) / , re-spectively. Here we focus on the regime Ω (cid:28) | ∆ | , where r vdw can be expressed in terms of the interaction param-eter ϕ = g r b /c ∆ as r vdw ≈ r b √ ϕ . The scaling of r vdw with ϕ indicates that, for weak interactions ϕ (cid:28)
1, thelow-energy scattering will be dominated by the core ofthe potential, while, for strong interactions ϕ (cid:29)
1, thelow-energy scattering will be dominated by the van derWaals tail. We now give approximate expressions for a and r in these two regimes for attractive interactions( m/χ >
1. Weak Attractive Interactions
In the regime ϕ (cid:28)
1, the low-energy scattering is dom-inated by the core of the potential, which can be wellapproximated by a square well of width 2 r b . We parame-terize the depth as − β r b /χ , where β is a free parameterchosen to match the observed scattering resonances. Inthis case, the scattering states can easily be found ana-lytically and a and r take the form [S3] a = r b + r b β ϕ tan( β ϕ ) , (S22) r r b = 2 − r b a + 23 r b a − (cid:18) − r b a (cid:19) (cid:18) tan βϕβϕ + 1cos βϕ (cid:19) . (S23)In this approximation, the scattering resonances occurwhen ϕ crosses nπ/β . Expanding near each resonancegives n = 0 , a = r b β ϕ , r r b = 23 β ϕ , (S24) n > , a = r b βnπ δϕ , r r b = 1 − βδϕnπ , (S25)where δϕ = ϕ − nπ/β is assumed to be small.We fix β by comparing Eq. (S24) to the asymptoticresult for ϕ →
0. In this limit, one can replace theeffective potential by a delta function 2 v δ ( r ), with v = R ∞ drV ( r ) = − ( π/ r b /χ . The scatteringlength takes the form a ≈ (3 /π ) r b /ϕ [S2], which fixes β = p π/
3. With this choice of β , we find Eq. (S22) is ingood agreement with the n = 0 and n = 1 scattering res-onances characterized in Ref. [S2], but begins to deviateat the n = 2 scattering resonance.
2. Strong Attractive Interactions
When ϕ (cid:29)
1, the effective potential for the polari-tons has many features in common with the potentialsconsidered in models of atomic scattering [S4]. In thesemodels, the atomic potential U ( r ) is treated as having adeep attractive core, while for large rU ( r ) ≈ − C n r n , (S26)where n = 6 for van der Waals forces. For s-wave scatter-ing, a and r can be found from the zero energy solutionto the radial Schr¨odinger equation ∂ r ψ + [ p ( r )] ψ = 0 , (S27)where p ( r ) = p − mU ( r ) and the boundary condition fors-wave scattering is ψ (0) = 0.For the 1d Rydberg polariton problem consideredhere, Eq. (S27) is equivalent to Eq. (S21) with p ( r ) = ϕ/ p r b + r /r b , but with the boundary condition ∂ r ψ | r =0 = 0. Due to the similarity in the equations,for ϕ (cid:29)
1, we can follow Ref. [S4, S5] to find analyti-cal solutions for a and r . In particular, we can solvefor ψ for small r using a WKB approximation, whichcan then be matched to the known asymptotic solutionfor ψ at large r . The WKB solution is valid in the re-gion r (cid:28) √ ϕ r b , while the asymptotic solution is validwhen r (cid:29) r b [S4]. Thus, the existence of an intermediateregime r b (cid:28) r (cid:28) √ ϕ r b is equivalent to the requirementof strong interactions ϕ (cid:29) ψ ( x ) = C √ p ( r ) cos (cid:2)R r dr p ( r ) (cid:3) , r (cid:28) √ ϕ r b √ x (cid:20) AJ / (cid:18) ϕ r b r (cid:19) − BN / (cid:18) ϕ r b r (cid:19)(cid:21) , r (cid:29) r b (S28) where A, B, and C are unkown coefficients which haveto be determined by matching the two solutions in theintermediate region and J α ( N α ) are Bessel functions ofthe first (second) kind. The WKB solution is chosen tosatisfy the boundary condition that ψ has zero derivativeat the origin. Following a similar analysis to Ref. [S4, S5]we solve for the coefficients A , B , and C , which determine a and r as a = ¯ a h − tan (cid:0) Φ + π/ (cid:1)i , (S29)¯ a = Γ(3 / / p ϕ r b ≈ . √ ϕ r b , (S30)Φ = Z ∞ drp ( r ) = Γ(1 / / √ π ϕ ≈ . ϕ, (S31) r = 1 . √ ϕ r b − . ϕ r b a + 0 . ϕ / r b a (S32)where ¯ a is the scattering length averaged over Φ (exclud-ing the resonances) and Γ( · ) is the gamma function. Thescattering resonances occur when Φ = Φ n = nπ + 3 π/ n th resonance gives a ≈ ¯ a Φ − Φ n , (S33)while the effective range becomes r ≈ . r vdw . III. EFT INCLUDING RAMAN RESONANCES
In this section, we write down an EFT that describesthe coupling between the polariton field ˆ ψ and the Ramanresonance excitations.We account for the presence of the constant term in thepropagator g ss ( q, ν ) by adding a fictitious pair of parti-cles d ± to the EFT H = − m ˆ ψ † ∂ z ˆ ψ + X s = ± ω s d † s d s (S34)+ Z dz Ψ † ( z )Ψ † ( z ) V ( z − z )Ψ( z )Ψ( z ) , Ψ( z ) = α ˆ ψ ( z ) + S + d + ( z ) + S − d − ( z ) , (S35)where ω ± is given by Eq. (S14) and the interaction termaccounts for all of the allowed interactions between thefictitious particles. The terms α = g c / p Ω + g c ≈ S ± (given by Eq. (S16)) account for the overlap of theseparticles with the | s i state. Integrating out the fields d ± gives rise to the N -body interactions discussed in themain text and in the section below. IV. NON-PERTURBATIVE N -BODYINTERACTION POTENTIAL When | χ N V ( r ) | >
1, the perturbative solution for V N eff given in Eq. (8) of the main text breaks down. We nowshow how to find the non-perturbative solution to V N eff recursively using the Rosenberg integral equations forthe connected transition matrix [S6]. We explicitly solvethese equations for N = 2 , , and 4. Finally, we showthat, inside the blockade radius, V N eff oscillates betweenattraction and repulsion with every increase in N .Before writing the Rosenberg equations, we first intro-duce some basic concepts needed to describe N -particlescattering [S7]. The notion of connected scattering dia-grams, illustrated in Fig. 3 in the main text, leads to thecluster decomposition for the N -body transition matrix T N = X α ∈ P [ T N ] α , (S36)where P is the set of all partitions of N particles into dis-joint clusters. For example, for three particles there arefive such partitions (1)(2)(3), (12)(3), (13)(2), (23)(1),and (123). We define n α as the number of clusters withinthe partition α . We denote the particles representedin each cluster as i , . . . , i m n , where 1 ≤ n ≤ n α and m n is the length of the n th cluster. For the partition(12)(3) n α = 2, m = 2 with i = 1 , i = 2 and m =1with i = 3. We also introduce the notion of an order-ing ≺ of the clusters: α ≺ β if every cluster in α is asubset of the elements in clusters of β . For example,(12)(3)(45) ≺ (123)(456), but (12)(3)(45) ⊀ (1234)(56).To find each term in Eq. (S36), one needs to evalu-ate the sum of all scattering diagrams where each clus-ter in the partition is disconnected from the others, butfully connected internally. For example, the diagramin Fig. 3(b) of the main text contributes to [ T N ] α with α = (12)(345). In Eq. (8) of the main text we derived aperturbative solution for the fully connected contributionto T N , V N eff = [ T N ] α with α = (1 , . . . , N ). However, thisequation breaks down for r (cid:28) r b where | χ N V ( r ) | > V N eff has to be found non-perturbatively.The key insight into the connected N -body scatteringdiagrams is that they can each be written as the product M ‘α χ N L ‘ , where L ‘ is a diagram that ends with interac-tion V ‘ on the left and M ‘α can be broken into two disjointclusters α such that ‘ ⊀ α [S8, S6]. For the connected dia-gram in Fig. 3(a) in the main text M ‘α = ( χ N ) V V V and L ‘ = V , with α = (1234)(5) and ‘ = (45). TheRosenberg integral equations (which are algebraic equa-tions for a constant propagator) take advantage of thisstructure to recursively define [S6, S7] T cN‘ = X ‘ X α (cid:31) ‘,n α =2 [ T ‘ ] α ¯∆ α‘ χ N T ‘ , (S37) T ‘ = V ‘ + V ‘ χ N T N = V ‘ − χ N P ‘ V ‘ , (S38) V N eff ( z ; ν ) = X ‘ T cN‘ ( z ; ν ) , (S39)where T N = P ‘ T ‘ , ‘ = ( ij ) denotes a particle pair andranges over all N ( N − / T ‘ groups all dia-grams contributing to T N that end with the interaction V ‘ on the left, and ν is the total frequency of the incom-ing photons. The sum in Eq. (S37) is over all partitions α with two clusters, which contain the pair ‘ . The ma-trix ¯∆ α‘ = 1 if ‘ ⊀ α and zero otherwise. ¯∆ α‘ reflectsthe structure of the connected diagrams described aboveand enforces that all the terms in Eq. (S37) are fully con-nected. Using the results from Sec. I, we can also give anexplicit expression for the N -body propagator χ N ( ν ) = X { ( s ,...,s N ) ,s i = ±} ν − P Ni =1 ω s i N Y i =1 | S s i | . (S40)Equation (S37) is recursive because [ T ‘ ] α can be ex-pressed in terms of the connected transition matrices for1 ≤ k ≤ N − T ‘ ] α = T NcN − ,‘ ( z i , . . . , z i N − ) , ( k = 1) (S41)[ T ‘ ] α = (cid:18) N − k − (cid:19) χ N T NcN − k,‘ ( z i , . . . , z i N − k ) (S42) × X ‘ ≺ ( i ...i k ) T Nck‘ ( z i , . . . , z i k ) , ( k > N denotes that T Ncm‘ is found usingEq. (S37) for m particles, but with the propagator χ N replacing χ m . The binomial factor in front of Eq. (S42)counts the number of ways to arrange the scatteringevents between the two clusters, with the constraint thatthe pair ‘ always interact first.For N = 3, we find the non-perturbative solution T c = χ − χ P ‘ V ‘ V − χ V ( V + V ) , (S43) T c = χ − χ P ‘ V ‘ V − χ V ( V + V ) , (S44) T c = χ − χ P ‘ V ‘ V − χ V ( V + V ) , (S45) V ( z , z , z ; ν ) = X ‘ T c ‘ ( z , z , z ; ν ) , (S46)which agrees with the perturbative result from Eq. (8) inthe main text to lowest order in | χ V ‘ | . For N = 4, thefull expression involves considerably more terms: T c = χ − χ P ‘ V ‘ (cid:2) T c ( z , z , z )( V + V + V )+ T c ( z , z , z )( V + V + V ) (S47)+ χ T c T c ( V + V + V + V ) (cid:3) , and similarly for the other ‘ . Here T c ( z , z , z ) isgiven by T c ( z , z , z ) from Eq. (S43) with χ ( ν ) re-placed by χ ( ν ), T c = V / (1 − χ V ) and similarlyfor T c ( z , z , z ) and T c . The resulting expressionfor V , to lowest order in χ V ‘ , contains 96 tree diagrams(16 unique) and agrees with Eq. (8) in the main text.By construction, Eq. (S37) accounts for all connectedscattering diagrams. To check that these recursive for-mulas give the same number of terms as the perturbativeresult, we use Eq. (S37)-(S42) to find a recursive formulafor the number of terms contributing to T cN‘ t ‘ = 1 , (S48) t ‘N = t ‘N − ( N − N −
1) (S49)+ N − X k =2 t ‘N − k t ‘k k ( k − (cid:18) N − k − (cid:19)(cid:18) N − k (cid:19) ( N − k ) k, where (cid:0) N − k (cid:1) is the number two-cluster partitions con-taining ‘ with m = N − k and m = k and ( N − k ) k is the number of non-zero elements in ¯∆ α‘ for each two-cluster partition α . Based on Eq. (8) in the main textwe expect t ‘N = 2( N − N N − . This can be proved byinduction using Eq. (S49) combined with an applicationof the binomial formula [S9]( a + b ) n a = n X k =0 (cid:18) nk (cid:19) ( a − kc ) k − ( b + kc ) n − k , (S50) with a = 1 , c = − , b = N − , and n = N −
3. This resulthelps confirm that Eq. (8) of the main text is consistentwith our non-perturbative solution for V N eff .When all the photons are separated by much less thanthe blockade radius, we can approximate V ( r ) by ±∞ .In this limit, T Nck‘ saturates to a constant value that de-pends only on N and k . By adapting the counting ar-guments used to derive Eq. (S49), after some simplifica-tions, we then arrive at a similar recursive formula for V N eff ( z , . . . , z N ; ν ) in this regime V N eff ( z , . . . , z N ; ν ) ≈ ( − N − c N χ N , (S51) c N = 2 c N − + 2 N − X k =2 (cid:18) N − k − (cid:19)(cid:18) N − k − (cid:19) c N − k c k , (S52)where c = 1 and c = 2. Since c N is a positive integer forevery N , we find that, similar to the perturbative resultfrom Eq. (10) in the main text, V N eff alternates betweenattraction and repulsion for every increase in N . [S1] M. Fleischhauer, A. Imamoglu, and J. P. Marangos, Rev.Mod. Phys. , 633 (2005).[S2] P. Bienias, S. Choi, O. Firstenberg, M. F. Maghrebi,M. Gullans, M. D. Lukin, A. V. Gorshkov, and H. P.B¨uchler, Phys. Rev. A , 053804 (2014).[S3] V. E. Barlette, M. M. Leite, and S. K. Adhikari, Eur. J.Phys. , 435 (2000).[S4] G. F. Gribakin and V. V. Flambaum, Phys. Rev. A ,546 (1993).[S5] V. V. Flambaum, G. F. Gribakin, and C. Harabati, Phys. Rev. A , 1998 (1999).[S6] L. Rosenberg, Phys. Rev. , B217 (1965).[S7] S. K. Adhikari and K. L. Kowalski, Dynamical CollisionTheory and Its Applications (Academic Press, San Diego,CA, 1991).[S8] S. Weinberg, Phys. Rev. , B232 (1964).[S9] R. P. Stanley,