Quantum dynamics of propagating photons with strong interactions: a generalized input-output formalism

There has been rapid development of systems that yield strong interactions between freely propagating photons in one dimension via controlled coupling to quantum emitters. This raises interesting possibilities such as quantum information processing with photons or quantum many-body states of light, but treating such systems generally remains a difficult task theoretically. Here, we describe a novel technique in which the dynamics and correlations of a few photons can be exactly calculated, based upon knowledge of the initial photonic state and the solution of the reduced effective dynamics of the quantum emitters alone. We show that this generalized"input-output"formalism allows for a straightforward numerical implementation regardless of system details, such as emitter positions, external driving, and level structure. As a specific example, we apply our technique to show how atomic systems with infinite-range interactions and under conditions of electromagnetically induced transparency enable the selective transmission of correlated multi-photon states.

There has been rapid development of systems that yield strong interactions between freely propagating photons in one dimension via controlled coupling to quantum emitters. This raises interesting possibilities such as quantum information processing with photons or quantum many-body states of light, but treating such systems generally remains a difficult task theoretically. Here, we describe a novel technique in which the dynamics and correlations of a few photons can be exactly calculated, based upon knowledge of the initial photonic state and the solution of the reduced effective dynamics of the quantum emitters alone. We show that this generalized "input-output" formalism allows for a straightforward numerical implementation regardless of system details, such as emitter positions, external driving, and level structure. As a specific example, we apply our technique to show how atomic systems with infinite-range interactions and under conditions of electromagnetically induced transparency enable the selective transmission of correlated multi-photon states.

I. INTRODUCTION
Systems in which individual photons can interact strongly with each other constitute an exciting frontier for the fields of quantum and nonlinear optics [1]. Such systems enable the generation and manipulation of nonclassical light, which are crucial ingredients for quantum information processing and quantum networks [2]. At the many-body level, it has been predicted that these systems can produce phenomena such as quantum phase transitions of light [3,4] or photon crystallization [5,6]. Early examples of such systems where strong interactions between photons could be observed consisted of individual atoms coupled to single modes of high-finesse optical cavities, within the context of cavity quantum electrodynamics (QED) [7,8]. More recently, a number of systems have emerged that produce strong optical nonlinearities between freely propagating photons, including cold atomic gases coupled to guided modes of tapered fibers [9,10], photonic crystal fibers [11] or waveguides [12], cold Rydberg gases in free space [13], and superconducting qubits coupled to microwave waveguides [14,15].
The paradigm of cavity QED admits exact analytical or numerical solutions at the few-photon (and/or fewatom) level through the elegant "input-output" formalism [16], as illustrated conceptually in Fig. 1(a). In this formalism, all of the properties of the field exiting the system (the output) can be determined based upon knowledge of the input field and the dynamics of the atomcavity system alone. For a few excitations, the latter can be solved due to the small Hilbert space associated with a small number of excitations of a discrete cavity mode. In contrast, despite rapid development on the experimental front, theoretical techniques to treat the dynamics of freely propagating photons interacting with spatially distributed emitters are generally lacking. At first glance, the challenge compared to the cavity case arises from the fact that a two-level system is a nonlinear frequency mixer, which is capable of generating a continuum of new frequencies from an initial pulse, as schematically depicted in Fig. 1(b). A priori, keeping track of this continuum as it propagates and re-scatters from other emitters appears to be a difficult task. An exception is the weak excitation limit, in which atoms can be treated as linear scatterers and the powerful transfer matrix method of linear optics can be employed [17,18].
The full quantum case has been solved exactly in a limited number of situations in which nonlinear systems are coupled to 1D waveguides, such as up to three photons scattering on a single atom [19,20], one and two photons scattering on a cavity QED system [21], and many photons scattering on many atoms coupled to a chiral (monodirectional) waveguide [22]. The formalism employed in Ref. [20] is particularly elegant, because it establishes an input-output relation to determine the nonlinear scattering from a two-level atom. Here, we show that this technique can be efficiently generalized to many atoms, chiral or bi-directional waveguides, and arbitrary atomic configurations, providing a powerful tool to investigate nonlinear optical dynamics in all systems of interest.
This paper is organized in the following way: first, we present a generalized input-output formalism to treat few-photon propagation in waveguides coupled to many atoms. We show that the infinite degrees of freedom associated with the photonic modes can be effectively integrated out, yielding an open, interacting "spin" model that involves only the internal degrees of freedom of the atoms. This open system can be solved using a number of conventional, quantum optical techniques. Then, we show that the solution of the spin problem can be used to re-construct the optical fields. In particular, we provide a prescription to map spin correlations to S-matrix elements, which contain full information about the photon dynamics, and give explicit closed-form expressions for the one-and two-photon cases. Importantly, in analogy with the cavity QED case, our technique enables analytical solutions under some scenarios, but in general allows arXiv:1501.04427v1 [quant-ph] 19 Jan 2015 for simple numerical implementation under a wide variety of circumstances of interest, such as different level structures, external driving, atomic positions, atomic motion, etc. Finally, to illustrate the ease of usage, we apply our technique to the study of nonlinear field propagation through an optically dense ensemble of atoms with Rydberg-like interactions.

II. GENERALIZED INPUT-OUTPUT FORMALISM
In this section we consider a generic system composed of many atoms located at positions z i along a bidirectional waveguide. We assume that there is an optical transition between ground and excited-state levels |g and |e to which the waveguide couples, but otherwise we leave unspecified the atomic internal structure and the possible interactions between them (e.g., Rydberg interactions), as such terms do not affect the derivation presented here. The bare Hamiltonian of the system is composed of a term describing the energy levels of the atoms H at , and a waveguide part where k is the wavevector and v = ± is an index for the direction of propagation, with the plus (minus) denoting propagation towards the right (left) direction. We assume that within the bandwidth of modes to which the atoms significantly couple, the dispersion relation for the guided modes can be linearized as ω k = c|k|. The interaction between atoms and photons is given by which describes the process where excited atoms can emit photons into the waveguide, or ground-state atoms can become excited by absorbing a photon. The coupling amplitude g is assumed to be identical for all atoms, while the coupling phase depends on the atomic position (e ivkzi ). Here, we will explicitly treat the more complicated bidirectional case, although all of the results readily generalize to the case of a single direction of propagation. Our immediate goal is to generalize the well-known input-output formalism of cavity QED [16] and the more recent formalism for a single atom coupled to a waveguide [20], to the present situation of many spatially distributed atoms coupled to a common waveguide. In short, we will eliminate the photonic degrees of freedom by formal integration, which reveals that the output field exiting the collection of atoms is completely describable in terms of the input field and atomic properties alone. This formal integration also provides a set of generalized Heisenberg-Langevin equations that governs the atomic evolution.
The Heisenberg equations of motion for σ i ge and b v,k can be readily obtained by calculating the commutators with H. To simplify the presentation of the resulting equations, we replace the spin operators σ i ge with the bosonic annihilation operator a i . The two-level nature a)

b)
FIG. 1. Cavity QED vs. many-atom waveguide (a) Schematic representation of a cavity QED system. The inputoutput formalism enables the outgoing field to be calculated based on knowledge of the input field and the internal dynamics of the cavity system. (b) System composed of many atoms coupled to a common waveguide. The nonlinearity of an atom enables the generation of a continuum of new frequencies upon scattering of an incoming state. This property, combined with multiple scattering from other atoms, appears to make this system more complicated than the cavity QED case.
of the atomic transition can be retained by introducing an interaction energy U 0 for multiple excitations, through in H at , and by taking the limit U 0 → ∞ at the end of calculations for observable quantities.
The Heisenberg equations for b v,k can be formally integrated and Fourier transformed, to obtain the real-space wave equation Here b v,in is the homogeneous solution, physically corresponding to the freely propagating field in the waveguide, while the second term on the right consists of the part of the field emitted by the atoms. Inserting Eq. (2) into the equation for a i , we obtaiṅ In realistic systems, time retardation can be neglected, resulting in the Markov approximation a j (t − |z i − z j |/c) ≈ a j (t)e iωin|zi−zj |/c . Here, ω in is a central frequency around which the atomic dynamics is centered (typically the atomic resonance frequency). This approximation is valid when the difference in free-space propagation phases ∆ωL/c 1 is small across the characteristic system size L and over the bandwidth of photons ∆ω involved in the dynamics. As a simple example, the characteristic bandwidth of an atomic system is given by its spontaneous emission rate, corresponding to a few MHz, which results in a significant free-space phase difference only over lengths L 1 m much longer than realistic atomic ensembles.
We have thus obtained the generalized Heisenberg-Langevin equatioṅ where we have identified Γ 1D = 4πg 2 /c as the singleatom spontaneous emission rate into the waveguide modes. If we keep separated the terms proportional to a j coming from the right and left-going photonic fields, we can find easily that the Lindblad jump operators corresponding to the decay of the atoms into the waveguide are O ± = Γ 1D /4 j a j e ∓iωinzj /c , in terms of which we can write the master equation for the atomic density matriẋ We also see that we can derive Eq. (4) from a non-Hermitian effective Hamiltonian which can be used for a quantum jump description of the atomic dynamics. The resulting infinite-range interaction between a pair of atoms i,j intuitively results from the propagation of a mediating photon between that pair, with a phase factor proportional to the separation distance. Within the same approximations employed above to derive the Heisenberg-Langevin equations we can obtain a generalized input-output rela- for right(left)-going fields. However, since the right-going output field propagates freely after z R , it is convenient to simply define b +,out (t) = b +,out (z R + , t) as the field immediately past the right-most atom (where is an infinitesimal positive number), and similarly for the left-going output. The derived relation shows that the out-going field properties are obtainable from those of the atoms alone.
The emergence of infinite-range interactions between emitters mediated by guided photons, and input-output relationships between these emitters and the outgoing field, have been discussed before in a number of contexts [18,20,23], but the idea that such concepts could be used to study quantum interactions of photons in extended systems has not been fully appreciated. In the remaining sections, we will demonstrate the effectiveness of this approach to quantum nonlinear optics. In particular, the infinite-dimensional continuum of the photons is effectively reduced to a Hilbert space of dimension where n is the maximum number of atomic excitations (for n = N we have dim[H] = 2 N ). The atomic dynamics, on the other hand, having been reduced to standard Heisenberg-Langevin equations, quantum jump, or master equations, are solvable by conventional prescriptions [24]. It is also possible to derive generalized master equations describing the atomic dynamics in response to arbitrary (e.g., non-classical) incident states of light, as detailed in Appendix C.

III. RELATION TO S-MATRIX ELEMENTS
The S-matrix characterizes how an incoming state of monochromatic photons evolves via interaction with atoms into a superposition of outgoing monochromatic photons. Because monochromatic photons form a complete basis, the S-matrix thus contains all information about photon dynamics. Here, we will show that Smatrix elements can be calculated from the dynamics of the atoms evolving under the effective spin model.
Formally, the S-matrix for the interaction of an nphoton state with an arbitrary system is defined as where the input and output creation operators respectively create freely propagating incoming and outgoing photonic states. The vectors p and k denote the outgoing and incoming frequencies of the n photons. The input and output operators can be any combination of + and -propagation directions (we have omitted this index here for simplicity). In the last line we have used a global Fourier transformation , to express the S-matrix in time.
The S-matrix has been calculated exactly before in a limited number of situations [19][20][21][22]. Here, we provide a general prescription to numerically obtain the Smatrix starting from the input-output formalism. First, it should be noted that the operators b in , b out correspond to the input and output operators defined in the previous section [20]. On the other hand, the input-output relation enables the correlator of Eq. (6) to be written purely in terms of atomic operators. For notational simplicity, we give a derivation for a single spin and a monodirectional waveguide, but its generalization to the bidirectional waveguide and many atoms is straightforward. For our purpose it is enough to have an input-output relation of the form where a is in our case the spin operator σ ge . We summarize the main idea of the derivation here, while the details can be found in Appendix A.1. We begin by noting that Eq. (7) enables one to replace output operators by a combination of system and input operators, or input operators by system and output operators. Selectively using these substitutions, one can exploit favorable properties of either the input or output field, in order to gradually time order all of the system operators (where operators at later times appear to the left of those at earlier times), while removing input and output operators from the correlation. The favorable properties that can be used are that 1) system and input operators commute, [a(t), b in (t )] = 0, when t > t, and likewise [a(t), b out (t )] = 0 for t < t, 2) input annihilation operators at different times commute amongst themselves, as do output annihilation operators, and 3) the correlator can be reduced in size Through this procedure, the S-matrix can be expressed as a Fourier transform of a sum of terms involving only time-ordered atomic operators (indicated by the operator T ). These functions generally have the form with m ≤ n, multiplied by n − m delta functions in time. Moreover, using the general expression for the Heisenberg-Langevin equation and the quantum regression theorem, it can be proven that when external fields driving the system do not generate waveguide photons, the correlation function of Eq. (8) can be evaluated by evolving a(t) as e iH eff t a e −iH eff t (see Appendix A.2 for a formal derivation). Here, H eff is the effective Hamiltonian from Eq. (5) that contains only the spin operators. Although such a form for a(t) is not true in general due to quantum noise, these noise terms have no influence on the correlation. A similar procedure as above enables one to express other important observables of the field, such as the second-order correlation function g (2) (t), in terms of correlation functions involving atoms alone.
While the discussion has thus far been completely general, the case of S-matrix elements involving only one or two photons can be formally reduced to particularly simple expressions. For example, in Appendix B, we show that the transmission coefficient T k for the manyatom, bi-directional waveguide case is related to the Smatrix by S (1) Furthermore, it can be expressed in terms of a known ∼ N × N matrix corresponding to the single-excitation Green's function G 0 (whose form varies depending on the system details), Similarly, the two-photon S-matrix in transmission is generally given by where the first and second terms on the right describe the linear and nonlinear contributions, respectively. The latter term can be expressed in terms of matrices W related to the single-excitation Green's function, and a known ∼ N 2 × N 2 matrix T characterizing atomic nonlinearities and interactions.

IV. ELECTROMAGNETIC INDUCED TRANSPARENCY
In this section, we apply our formalism to a specific example involving three-level atoms under conditions of electromagnetically induced transparency (EIT) and with Rydberg-like interactions between atoms [13]. The linear susceptibility for a two-level atom with states |g and |e , in response to a weak probe field with detuning δ = ω p − ω eg from the atomic resonance, is shown in Fig. 2b. It can be seen that the response on resonance is primarily absorptive, as characterized by the imaginary part of the susceptibility (χ , red curve). In contrast, the response can become primarily dispersive near resonance if a third level |s is added, and if the transition |e − |s is driven by a control field (characterized by Rabi frequency Ω and single photon detuning δ L = ω L − ω es ). Specifically, via interference between the probe and control fields, the medium can become transparent to the probe field (χ = 0) when two-photon resonance is achieved, δ − δ L = 0, realizing EIT [25]. In this process, the incoming probe field strongly mixes with spin wave excitations σ sg to create "dark-state polaritons." The medium remains highly transparent within a characteristic bandwidth ∆ EIT around the two-photon resonance, which reduces to ∆ EIT ∼ 2Ω 2 /(Γ √ D), when δ L = 0. Here we have introduced the total single-atom linewidth Γ (see below) and the optical depth D, corresponding to the opacity of the medium in absence of EIT and defined in terms of input and output intensity as I out (Ω = 0)/I in (Ω = 0) = exp(−D). These polaritons propagate at a strongly reduced group velocity v g c, as indicated by the steep slope of the real part of the susceptibility χ in Fig. 2b, which is proportional to the control field intensity [25].
Taking s i to be the annihilation operator for the state |s i , EIT is described within our spin model by the effective spin Hamiltonian where the first line represents the explicit form of H at in Eq. (5) for the EIT three-level atomic structure. In addition to waveguide coupling, here we have added an independent atomic decay rate Γ into other channels (e.g., unguided modes), yielding a total single-atom linewidth of Γ = Γ + Γ 1D . Our theoretical model nearly ideally describes experiments in which atoms are coupled to one-dimensional waveguides such as nanofibers [9] or photonic crystals [10], in which the interaction probability of a single photon and atom is Γ 1D /Γ ∼ 0.1 is significant and up to N ∼ 10 3 atoms are trapped. One consequence of the large interaction probability is that even a single atom can yield a significant reflectance for single photons [10]. In free-space experiments, the atom-photon interaction probability is much smaller, while a much larger number of atoms is used to induce a significant optical response. While this large number cannot be implemented numerically with our model, our system can nonetheless be used to reproduce macroscopic observables, thus making our approach an excellent description for atom-light interfaces in general.
Intuitively (and as can be rigorously shown, see below), provided that free-space experiments only involve a single transverse mode of light, one expects that the system response to light only depends on the interaction probability and number of atoms via their product. This product in fact defines the optical depth D = N (2Γ 1D /Γ) [26]. As D 10 2 in free-space experiments, such systems can be evaluated using our spin model with several hundred atoms and suitably large Γ 1D /Γ. The only remaining issue is the potentially large reflection generated by atoms in the mathematical model, compared to free-space ensembles where reflection is negligible. Reflection can be suppressed in the model by giving the atoms a lattice constant d where k in d = (2m + 1)π/2, where m is a nonnegative integer, such that the reflection from atoms destructively interferes [29,30]. All subsequent numerical calculations are performed under this condition.
The spin model on a lattice describing EIT, Eq. (11), can be exactly solved in the linear regime using the transfer matrix formalism [29], which correctly reproduces the free-space result and dependence on optical depth for group velocity v g = 2Ω 2 n/Γ 1D and transparency window ∆ EIT , where n is the (linear) atomic density. The corresponding minimum spatial extent of a pulse that can propagate inside the medium with high transparency is given by σ EIT = v g /∆ EIT .
A single photon propagating inside an ensemble of atoms under EIT conditions is coherently mapped onto a single dark polariton, corresponding to a delocalized spin wave populating the single excitation subspace of the atomic ensemble. The polariton dynamics can be therefore visualized directly by monitoring the excitation probability σ j ss of the atoms in the ensemble. In Fig. 3, we initialize a single polariton inside the medium with an atomic wave function of the form |ψ = j f j σ j sg |g ⊗N , and we determine numerically the time evolution under H in Eq. (11) up to a final time t f . Choosing an initially Gaussian spin wave, f j = exp(ik in dj) exp(−(jd − µ) 2 /4σ 2 p )/(2πσ 2 p ) 1/4 , with spatial extent σ p (blue line), one sees that the wavepacket propagates a distance v g ·t f , and with little loss provided that σ p > σ EIT . Numerics (green line) show perfect agreement with theoretical predictions (red lines) obtained via the transfer matrix formalism [29].

V. INFINITE RANGE INTERACTION
The spin model formalism can be easily extended to include arbitrary atomic interactions, providing a powerful tool to study quantum nonlinear optical effects. As a concrete example, we consider a system in which atoms can interact directly over a long range, such as via Rydberg states [6,31,32] or photonic crystal bandgaps [33]. The total Hamiltonian is given by in which U ij represents a dispersive interaction between atoms i and j when they are simultaneously in state |s . As we are primarily interested in demonstrating the use of our technique, we take here a "toy model" where atoms experience a constant infinite-range interaction, U ij /2 ≡ C. Such a case enables the numerical results to be intuitively understood, although we note that other choices of U ij do not increase the numerical complexity.
In particular we are interested in studying the propagation of a constant weak coherent input field through the atomic ensemble. The corresponding driving then is given by Γ is the amplitude of the constant driving field, ∆ = δ − δ L the detuning from two photon resonance condition, and the initial state is given by the global atomic ground state |ψ i = |g ⊗N [27,28]. With infinite-range interaction, one spin flip to state |s j shifts the energies of all other states |s j by an amount C. A second photon should then be able to propagate with perfect transparency, provided it has a detuning compensating for the energy shift C, thus ensuring the two-photon resonance condition is satisfied. As we result, we expect to see a transparency window for two photons, whose central frequency shifts linearly with C.
This predicted behavior can be confirmed by plotting the transmitted intensity fraction, T 1 = I/I in = b † +,out (t)b +,out (t) /E 2 , and also the second-order correlation function which corresponds roughly to the two-photon transmission. Fig. 4 shows the single-photon transmission T 1 and two-photon transmission T 2 as a function of the interaction strength C and detuning from two photon resonance ∆. As expected, T 1 shows a peak at ∆ = 0 independently of the interaction intensity C; instead the peak in T 2 shifts towards ∆ = C/2 with increasing C. The decay of T 2 for increasing C can be intuitively understood by noting that we have a constant coherent state input, in which photons are randomly spaced, causing two photons to enter the medium at different times. Thus, until the second photon enters, the first photon propagates as a single polariton detuned by ∆ from the single-photon transparency condition, getting partially absorbed in the process. By increasing the interaction we increase the detuning for this single polariton and consequently its absorption, explaining the trend observed for T 2 in Fig. 4. A quantitative description of this phenomenon is given in Appendix D.
Field correlation functions like intensity I = b † +,out (t)b +,out (t) or g 2 (τ ) = b † +,out (t)b † +,out (t + τ )b +,out (t + τ )b +,out (t) /I 2 can be computed according to the following strategy. First we switch from Heisenberg representation to Schrödinger representation, so that for the intensity we get: The time evolved wave function, |ψ(t) , is determined by numerically evolving the initial spin state |ψ i under H for a time t. Then, the state immediately after detection of one photon, b +,out |ψ(t) = |φ , is evaluated by expressing b +,out in terms of spin operators using the input-output formalism: b +,out = Ee ikinz R − iΓ 1D /2 j σ ge e ikin(z R −zj ) . Finally we obtain the intensity by computing the probability of the one-photon detected state, I = φ|φ .
For g 2 (τ ) an extra step is needed. Its numerator describes the process of detecting two photons, the first at time t and the second at time t + τ : As before we can switch to Schrödinger picture, f (t + τ ) = ψ(t)|b † +,out e iHτ b † +,out b +,out e −iHτ b +,out |ψ(t) , and evaluate the state after dectection of the first photon, b +,out |ψ(t) = |φ . Then detection of a second photon after a time τ entails performing an extra evolution under H and annihilating a photon, that is: b +,out e −iHτ |φ = b +,out |φ(τ ) . Finally, we evaluate the quantity f (t + τ ) = φ(τ )|b † +,out b +,out |φ(τ ) , by again expressing b +,out in terms of spin operators. In Fig. 4c, we plot the numerically obtained result for g (2) (τ ), for the case where infinite-range interactions are turned on (C = 1) and for a weak coherent input state with detuning ∆ = 0. In such a situation, one expects for the singlephoton component of the coherent state to transmit perfectly, while the two-photon component is detuned from its transparency window and becomes absorbed. This nonlinear absorption intuitively yields the strong antibunching dip g (2) (τ = 0) < 1. We also evaluate this second-order correlation function using the analytical result for the two-photon S-matrix in Eq. (10), which shows perfect agreement as expected.

VI. CONCLUSION
We have shown that the dynamics of a few photons, propagating under strong interactions mediated by quantum emitters, is fully and efficiently characterized by the dynamics of an open quantum spin model. As the spin model is solvable by standard quantum optical techniques for open systems, our approach provides an easily implementable recipe for the exact numerical study of a large class of quantum nonlinear optical systems. As an example, it provides an attractive alternative to numerical simulations of propagation through cold Rydberg gases, where typically the continuous field is finely discretized and its degrees of freedom are explicitly kept track of [31]. Our technique can also be used to treat a number of strongly nonlinear systems based upon atomic saturation [5], where only approximate effective theories for field propagation were previously available. We anticipate that the availability of exact results will greatly aid in the development of effective theories for the generally challenging problem of quantum many-body states of light, and that our model will shed new insight on conceptual links between such systems and quantum spin systems. We present here the full derivation of the decomposition of the S-matrix elements in terms of time-ordered atomic correlation functions. The starting point is the definition of the S-matrix in Eq. (6). Since the output operators commute between themselves because of the indistinguishability of photons, they can be freely ordered by decreasing times. Introducing the time ordering operator T and also using Eq. (7), the operators in Eq. (6) can be written as It is natural to label the terms above by the number m of system operators a present in each term. Thanks to the fact that [a(t), b in (t )] = 0 for t > t, all the input operators can be moved to the right of the spin operators. Thus, the term of order m will be of the form and can be simplified using the commutation relations between input operators [b in (t), b † in (t )] = δ(t − t ). This manipulation results in a sum of (n!) 2 /(m!) 2 (n − m)! terms for each original term of order m. Each term of the sum consists of n − m delta functions multiplied by a correlation function of the form Since the a operators commute with the b † in operators at later times, the time ordering operator can be extended to all the operators in the correlation function. Using again Eq. (7) to express the input operators one gets However, since the operators b † out commute with all the operators a on the left, only remains, which reproduces Eq. (8).

Evolution under the effective Hamiltonian
In Eq. (8) the vacuum state |0 stands for |0 b |g ⊗N , i.e., the vacuum state of field modes and the ground state of all the atoms. The atomic operators are in the Heisenberg picture, a(t) = e iHt ae −iHt , and H is the Hamiltonian of the whole system without the driving field. By taking the term with t 1 > t 2 ... > t m > t 1 > t 2 ... > t m as an example (our argument holds for any time ordering), we show now that Eq. (8) can be evaluated by effectively evolving system operators as a(t) = e iH eff t ae −iH eff t . The quantum regression theorem is applied here to eliminate the bath or photonic degree of freedom, and results in where ρ(0) = |g g|⊗|0 b 0| and L is the Lindblad superoperator of the system, defined in the main text. L contains a deterministic part, which generates an evolution driven by H eff and which conserves the number of excitations, and a jump part, which reduces the number of atomic excitations. Because of the form of the correlators, which contain an equal number of atomic creation and annihilation operators, the jump part of the evolution of the operators gives a vanishing contribution to the correlation function, proving what was stated above.

Appendix B: Examples of S-matrix elements
In this section of the Appendix we provide explicit examples of how to use the results derived above to calculate specific matrix elements. In particular we present the case of the scattering of 1) n photons on a two-level atom coupled to a one-directional waveguide and 2) one and two photons scattering on a Rydberg-EIT system.

n photons scattering on a single two-level atom
The S-matrix element S (n) for the scattering of n photons can be in general decomposed in a part that is a series of products of lower-order elements and a part that cannot be expressed in such a way. The latter is called fully connected part of the matrix element and physically corresponds to a n-body interaction, and is denoted by iT (n) . Thus, knowing how to decompose the S-matrix, one has to calculate iT (j) with 1 ≤ j ≤ n to construct S (n) . Furthermore, it can shown that the fully connected part in an element of order n can be obtained by the system correlation function of order n. We show here how to calculate such elements using the results presented above, for the case of a single two-level atom coupled to a onedirectional waveguide. It will be possible to appreciate the simplicity of our formalism compared the more cumbersome method used in Ref. [22] to obtain the same result.
We start from the relation between the connected part of the matrix element and the correlation function of atomic operators iT (n) The time ordering in the correlator gives 2n! possible orderings, but it is easy to arrive to the conclusion that only orderings which start on the left with a σ − operator and alternate σ + and σ − give a non-zero contribution. It is also immediate to see that the number of this possible orderings is (n!) 2 . A possible ordering is for instance which, expressing explicitly the time dependence, is equal to Since H eff is diagonal in the |g , |e basis, we can insert identity operators between the σ operators in the form of |g g|+|e e| in order to evaluate the Hamiltonians in the exponents. We immediately end up with where the permutations are over the two sets of incoming and outgoing photon frequencies k i and p i . The integral where ∆ i = k i − p i , which coincides with the result of Ref. [22].

One and two-photon scattering on a Rydberg-EIT system
The effective Hamiltonian of the Rydberg-EIT system is given by where H 0 is given by Eq. (11), and the hardcore interaction is with U 0 → ∞ corresponding to the three-level atom. For the single incident right moving photon with momentum k, it follows from Eqs. (6) and (8) that the Smatrix element is with This element describes the amplitude of a single transmitted photon with momentum p, using the time ordered correlation function of the system operator. We set for notational simplicity the speed of light c = 1. The second term of S (1) p+;k+ , i.e., the Fourier transform of the correlator can be obtained from the Green's function G 0 (ω) = 1/(ω − H 1 ) with the single-particle Hamiltonian (B11) Here, we have expressed H 1 in the basis {|a i , |s i }, and [G 0 (k)] σσ ij denotes the element σ i | G 0 (k) σ j with σ, σ = a, s. The element S (1) p+;k+ ≡ T k δ pk gives rise to the transmission coefficient For two right-going incident photons with momenta k 1 and k 2 , the two-photon S-matrix element S (2) p1+,p2+;k1+,k2+ describes the amplitude to a final state with two transmitted photons of momenta p 1 and p 2 . By Eqs. (6) and (8) in the main text, we find that c denotes the connected four-point Green's function, which involves an interaction between the two excitations in the system.
The S-matrix also contains terms arising from disconnected correlation functions, e.g., describing the linear propagation of each excitation separately, which yields the term in (B14) proportional to T k1 T k2 .
The interaction between two excitations under the Hamiltonian H HC + ij U ij s † i s † j s j s i /2 can be represented by the ladder diagram shown in Fig. 6. This diagram is represented mathematically by the two-body T-matrix, T (E), which satisfies the Lippmann-Schwinger equation, Here, E = k 1 + k 2 is the total energy, the vacuum bubble Π 0 (E) = (E − H 2 ) −1 is given in terms of H 2 = H 1 ⊗ I 2N + I 2N ⊗ H 1 , and the interaction matrix U has the diagonal element U aa ij = U as ij = U sa ij = U 0 and U ss ij = U ij in the basis {|a i a j , |a i s j , |s i a j , |s i s j }. The solution of Eq. (B15) is T(E) = 1/(U −1 − Π 0 (E)). The S-matrix can then be written as (B17) The Fourier transform of S (2) p1+,p2+;k1+,k2+ results in the wavefunction ψ(x c , x) = dp 1 dp 2 2π S (2) p1p2;k1k2 e i(p1x1+p2x2) of two transmitted photons, where the relative momentum k = (k 1 − k 2 )/2, the center of mass coordinate x c = (x 1 +x 2 )/2, and the relative coordinate x = x 1 −x 2 . The symmetric function is defined by the eigenstates |χ l and |χ l of H 1 and H † 1 with the corresponding eigenenergies ε l and ε * l , σ i |χ l = χ l (σ i ) and σ i |χ l =χ l (σ i ). Knowledge of the wave function (B18) in the case where the two incident photons have the same momentum, k 1 = k 2 = E/2, enables one to calculate the second-order correlation function for the outgoing field, g (2) We compare the result of g (2) (x) from the scattering theory and that from numerically solving the effective spin model (12) with the weak driving field in Fig. 4c in the main text, which shows that they agree with each other perfectly.

Appendix C: Generalized master equation
In this section, we extend our derivation of the atomic master equation, in order to calculate the response of the system to an arbitrary few-photon input state (as opposed to a classical or coherent state). The initial state is generally written as |ψ in = |ϕ in ⊗ |χ in , where |χ in is the initial state of the system, and describes an arbitrary state of incident photons in the Fock state basis by the wavefunction ϕ in ({n k }). By the relation between the Fock state |n and the coherent state |J = n J n |n / √ n!, we rewrite |ϕ in = F |{J v,k } by the operator acting on the coherent state |{J k } . The evolution of the reduced density matrix ρ s (t) = T r b [U(t)ρ(0)U † (t)] is determined by where ρ(0) = |ψ in ψ in | is the density matrix of the initial state, and H(t) is the total Hamiltonian in the rotating frame with b v,k → b v,k e −ivkt . By the relation |ϕ in = F |{J v,k } , the reduced density matrix reads where ρ s (0) = |χ in χ in | is the initial density matrix of the system. The displacement transformation and the generating density matrix where the unitary transformation 6. An intuitive explanation of how the two-photon wave function evolves within the atomic medium (located in the region 0 < x, y < L) can be found by considering a larger space, where first the two photons are both outside the medium (x, y < 0) and where one photon enters the medium first.
We specifically consider the case where infinite-range interactions cause the two photons to align with the EIT transparency condition when they are both inside the medium, as discussed further in the main text. The two-photon transmission is proportional to the value of the two-photon wave function at x, y = L.
is given by where H J (s) describes the system under the driving fields J v,k . In the density matrix ρ J (t), the initial state becomes the vacuum state due to the displacement transformation, thus the evolution of ρ J (t) satisfies the master equation where L is the Lindblad super-operator of the system without H J .
In conclusion, Eq. (C5) establishes the relation bewteen the evolution of the reduced density matrix ρ s (t) for a few incident photons (which could be a non-classical state) and the evolution of the reduced density matrix ρ J (t) for the system with classical (coherent state) driving fields J v,k . As a result, the few photon scattering problem can be understood as the perturbation expansion of the driving strength J v,k .
Appendix D: Linear approximation for two photon transmission T2 In this section of the Appendix we provide an intuitive explanation based on linear optics to the behavior of the two-photon transmission T 2 depicted in Fig. 4a,b. By indicating with x and y the coordinates of the first and second photons, we can divide the space into four regions according to their positions: both photons outside the atomic medium, one photon inside and one outside, and both photons inside the medium, see Fig. 6. As we are studying an infinite-range interaction, within each region the evolution is effectively linear. The different dispersion relations in each region, however, leads to non-trivial boundary conditions at their edges. Relevant to our discussion is the value of the two-photon wave function at the boundaries 0 < x < L and 0 < y < L. In particular, evolution in the region 0 < x, y < L with these boundary conditions dictates the two-photon transmission, which is proportional to the value of the two-photon wave function at (x, y) = (L, L). We study the case where the two photons have detunings ∆ = C/2. In this case, infinite-range interactions cause these two photons to satisfy the EIT transparency condition when both photons are inside the medium. Then, there are two qualitatively different regimes for the boundary values of the two-photon wave function, depending on the detuning from two-photon resonance δ (for simplicity we assume δ L = 0, so that ∆ = δ−δ L = δ): the small detuning regime, δ < ∆ EIT , in which a single photon can travel with high transmission through the medium and reach the detection point; and the large detuning regime, δ > ∆ EIT , in which a single photon is absorbed and cannot reach the detection point.
Small detuning regime.-Within the first regime δ = C/2 < ∆ EIT , a single polariton still fits within the EIT transparency window and exhibits high transmission through the medium. A representative probability distribution for a single excitation in this regime is illustrated in Fig. 7a. Now, we consider the two-excitation manifold. After detection of the first outgoing photon, the second photon can be at any position inside the ensemble with nearly uniform probability, and on average has to travel over a distance N d/2 before reaching the detection point. Therefore we expect T 2 = | exp(−ik(δ)N d/2)| 2 = √ T 1 . The predicted behaviour is in good agreement with full simulations results as shown in Fig. 8 for small detuning.
Large detuning regime.-In this regime a single photon is strongly absorbed. This is illustrated in Fig. 7b, where one observes a strong decay of the single-polariton probability and field intensity inside the medium. Within the context of Fig. 6, this localized single excitation serves as the boundary condition on the segments 0 < x < L and 0 < y < L. Furthermore, the evolution in the region 0 < x, y < L is effectively linear. Thus, the effect of this boundary condition at the detection point (a twoparticle problem) can be mapped onto a simpler problem, wherein one studies how a single excitation, initialized in the shape of the localized photon, transmits through the medium. In this regime, we therefore can estimate that T 2 is twice the value of the maximum transmission associated with this single localized photon, whose evolution we calculate numerically.
The good agreement of our simplified estimates with full simulations is demonstrated in Fig. 8, in which T 2 as a function of C = 2δ is compared for the different cases.