Photon molecules in atomic gases trapped near photonic crystal waveguides

Realizing systems that support robust, controlled interactions between individual photons is an exciting frontier of nonlinear optics. To this end, one approach that has emerged recently is to leverage atomic interactions to create strong and spatially non-local interactions between photons. In particular, effective interactions have been successfully created via interactions between atoms excited to Rydberg levels. Here, we investigate an alternative approach, in which atomic interactions arise via their common coupling to photonic crystal waveguides. This technique takes advantage of the ability to separately tailor the strength and range of interactions via the dispersion engineering of the structure itself, which can lead to qualitatively new types of phenomena. As an example, we discuss the formation of correlated transparency windows, in which photonic states of a certain number and shape selectively propagate through the system. Through this technique, we show in particular that one can create molecular-like potentials that lead to molecular bound states of photon pairs.

Realizing systems that support robust, controlled interactions between individual photons is an exciting frontier of nonlinear optics. To this end, one approach that has emerged recently is to leverage atomic interactions to create strong and spatially non-local interactions between photons. In particular, effective interactions have been successfully created via interactions between atoms excited to Rydberg levels. Here, we investigate an alternative approach, in which atomic interactions arise via their common coupling to photonic crystal waveguides. This technique takes advantage of the ability to separately tailor the strength and range of interactions via the dispersion engineering of the structure itself, which can lead to qualitatively new types of phenomena. As an example, we discuss the formation of correlated transparency windows, in which photonic states of a certain number and shape selectively propagate through the system. Through this technique, we show in particular that one can create molecular-like potentials that lead to molecular bound states of photon pairs.

I. INTRODUCTION
In recent years, there has been tremendous progress to realize systems that are capable of achieving strong interactions between individual photons [1]. One common approach has been to couple single atoms (or other quantum emitters) to high-finesse optical cavities, to take advantage of the intrinsic nonlinear nature of these twolevel systems [2][3][4][5][6][7][8]. More recently, gases of cold Rydberg atoms have been investigated [9][10][11][12][13][14][15]. In this case, the optical nonlinearities are effectively generated via strong atom-atom interactions, with a novel consequence being that the nonlinearity becomes spatially non-local in character [16]. It has been experimentally shown that this type of nonlinearity can give rise to exotic states such as a two-photon bound state [17], and a number of other few-and many-body states of light have been theoretically predicted [18][19][20].
Key to this approach is that atoms excited to a Rydberg state, in particular through the absorption of a photon, exert a strong dispersive effect on proximal atoms, shifting their transition frequencies to Rydberg states by an amount proportional to ∼ 1/r 6 , with r being the interatomic separation (see Fig. 1a) [21]. This in turn, leads to a shift in the optical susceptibility of the atoms, which can be interpreted as a change of refractive index that depends on the number and position of photons in the system. In practice, the r dependence of the Rydberg interaction presents somewhat of a constraint. Specifically, one must reach extremely high Rydberg states in order to induce a significant nonlocal effect (n ∼ 100 in Ref. [17] to achieve a blockade radius of 18µm). This in turn yields a large level shift (∼ 50GHz for atoms separated by 4µm in Ref. [17])); however, all of the interesting variation in the atomic refractive index occurs over much smaller * Corresponding author james.douglas@icfo.eu bandwidths of a few MHz (characteristic of the atomic linewidth). In this article, we demonstrate the new opportunities arising for nonlinear optics if the strength and range of the interaction can be independently adjusted. In particular, we show that the highly tunable nature of systems coupling atoms with nano-photonic devices can greatly extend the gamut of photon-photon interactions.
Our work is inspired by recent developments to interface cold atoms with photonic crystals [22][23][24][25][26]. It has been proposed that one can control the type and range of atomic interactions [27][28][29][30][31], and hence the interactions between photons [32], by engineering the underlying optical dispersion of these structures. We specifically study the dynamics of photons propagating through an atomic ensemble coupled to the waveguide under conditions of electromagnetically induced transparency (EIT), which occurs for three level atoms such as shown in Fig. 1a-b [33]. Without long-range interactions, EIT is an interference effect that allows an optical probe field to propagate without absorption due to the interference created by a second pump field when a "two-photon resonance" condition is satisfied, δ = 0 in Fig. 1b. As a linear optical effect, a probe pulse of any photon number and shape within a certain bandwidth propagates in this transparent manner.
Atomic interactions (that for example shift level |s in Fig. 1b) shift the two-photon resonance, altering the propagation of photons in the system and effective photon-photon interactions result. Provided these interactions are sufficiently long-range and smooth, we show that one is able to create "correlated transparent states," in which only pulses of a desired photon number and shape will propagate through. As an example, we provide an explicit construction to generate a two-photon molecule, in which two photons are bound at a fixed separation by an effective spring, and where "phonon" oscillations constitute the fundamental excitations of such a state.
The article is structured as follows. In Section II, we briefly review the physics governing the interaction between atoms and photonic crystal waveguides. In particular, we show how tunable long-range interactions between atoms emerge when the atomic transition frequency is situated within an optical "band gap" of the structure. In Section III, we show how EIT can be used to convert these atom-atom interactions into effective interactions between photons. In this regime we find an effective equation for photons propagating in the system that supports the presence of correlated transparent states.
In Section IV, we show how the photon interactions can be tuned to create molecular-like states of photons and, in Section V, we derive the conditions to observe such a state, demonstrating that the required parameters are within reach of state-of-the-art nano-photonic systems, before concluding in Section VI.

II. INTERACTIONS MEDIATED BY PHOTONIC CRYSTALS
Photonic crystals are dielectric structures in which the refractive index is modulated periodically [34]. At some frequencies, light input into the dielectric reflects constructively from these modulations and the system acts like a mirror, preventing propagation. This leads to the presence of band gaps in the dispersion relation for the photonic modes, as shown in Fig. 1(c). When an atom trapped nearby the photonic crystal is excited at a frequency in the band gap, it is forbidden from emitting a propagating photon into the dielectric. However, it can generate an exponentially decaying evanescent field that forms a localized photonic cloud around the atom [35]. It has been formally shown that this photonic cloud has the same functionality as a real cavity of the same size [30], with an effective vacuum Rabi splitting g c dictated by the mode volume and the atom-cavity detuning ∆ c = 2∆ b set by the separation ∆ b between the atomic frequency and band edge. Thus, exactly as in a real cavity in the far-detuned regime (|∆ c | > |g c |) [36], one atom can exchange its excitation via this virtual 'cavity' photon with another atom nearby with a characteristic strength g 2 c /∆ c . In this case the exchange is characterized by the dipole-dipole interaction Hamiltonian H = g 2 c ∆c j,l σ j ds σ l sd f (z j , z l ), where |s and |d are the ground and excited states coupled by the cavity mode, see Fig. 1(c) [27-30, 35, 37-39].
In the case where the atom is coupled predominantly to only a single band of a one-dimensional photonic crystal waveguide (achieved by tuning the atomic resonance close to the band edge), the interaction has spatial form for atoms at positions z j and z l , where E k b (z) is the Bloch function of the photonic crystal at the band edge (wavevector k b , frequency ω b ). Below, we perform numerical simulations with the atoms trapped on regularly spaced lattice with spacing z a = 2a, for photonic crystal unit cell length a, in which case E k b (z j ) = 1 at the trapping sites and the in- Here atom 1 is excited to Rydberg level R by two-photon excitation and the Rydberg interaction shifts the R level in other surrounding atoms by an amount much larger than the two-photon transition width, making the transition effectively two-level for a second photon entering the system. (b) Atoms with three internal levels support EIT. The input quantum probe field E couples to the transition |g − |e , while a second transition |s − |e couples to a classical drive field. By adding an additional level |d that couples to the band edge modes of a photonic crystal and driving the transition |s − |d off-resonantly an effective interaction between the |s levels can be generated. For simplicity we assume states |e and |d both have freespace spontaneous decay rate γ. (c) The periodic dielectric structure of a photonic crystal (right) leads to band structure (left) in the dispersion relation for photons propagating in the dielectric. A band gap occurs over frequencies in which propagation is forbidden due to strong interference from the underlying structure. For photonic crystals that support multiple modes, e.g., TE and TM (represented here by the solid and dotted lines), the band gaps of different modes may occur at different frequencies. The |s − |d transition of the atoms (green spheres) trapped near the photonic crystal couples to the photonic crystal bands. When the transition frequency ω sd is in the band gap (with detuning ∆ b from the band edge) atoms can couple with one-another via evanescent fields in the photonic crystal (illustrated in red), yielding a dipole-dipole interaction between the atoms.
teraction is purely in the relative coordinate |z j −z l |. The attenuation length of the interaction L = αω b /(k 2 b ∆ b ) corresponds to the length of the effective cavity mode and depends on the detuning ∆ b as well as the curvature of the band edge α. The resulting atom-atom interaction can be tuned in strength and length by adjusting these parameters, where a smaller detuning gives a longer range interaction and flatter bands lead to a shorter range interaction.
We wish to convert this dipole exchange interaction into a dispersive interaction between atoms in levels |s , as required to shift the EIT two-photon resonance. To achieve this, we off-resonantly drive the |s − |d transition of all atoms with Rabi frequency Ω s and detuning ∆ s (|∆ s | |Ω s |), as shown in Fig. 1(b). In this case, virtual excitations of state |d result in a Stark shift of level |s for the jth atom, which depends on the number of proximal atoms also in state |s . The interaction now has the general form given by the Hamiltonian s ∆c f (z j , z l ) for the single band edge coupling described above. For realistic experimental parameters, the coupling g c can be as large as 2π×10 GHz for a cavity that is just one wavelength long, L ∼ λ [30]. Importantly the strength of the interaction can be tuned to any smaller value by altering the ratio |Ω s /∆ s | for the driving laser, such that level shifts seen by the atoms can be on the order of the freespace linewidth γ, representing a significant distinction from the Rydberg case.
The photonic crystal mediated interactions are not without losses. In a realistic system imperfections in the dielectric medium cause loss of photons at rate κ. At the same time the photonic crystals used in current experiments are not three-dimensional and the excited state |d can spontaneously emit into free space at rate Γ comparable to γ. To understand the effect of these loss channels we take advantage of the mapping of the photonic cloud surrounding the atom to an effective cavity mode. For an optimal choice of the detuning ∆ b from the band edge and the curvature α of the band, the two-atom dispersive interaction strength can be made larger than the residual dissipation by a factor √ C/2, where C = g 2 c /(κγ d ) is the cooperativity of the cavity [30]. In state of the art systems the cooperativity can be a large as C λ ∼ 10 4 for a photonic cloud with an attenuation length L equal to the wavelength λ. For longer range couplings with decay length L, the cooperativity scales as C λ λ/L.

III. EIT WITH INTERACTIONS.
The basis of EIT in an atomic medium is atoms with three internal levels |g , |e and |s , as shown in Fig. 1(b). We consider a system where the atoms are coupled by the |g − |e transition to a quantum probe beam E propagating in a one-dimensional waveguide, while the |e −|s transition is coupled to an external control laser field (Rabi frequency Ω) and the photonic-crystal mediated interaction acts on level |s as described above. An experimental system that allows for this configuration has been demonstrated in Refs. [23,24,26], where atoms are trapped near a photonic crystal waveguide. In this system the band gaps for the transverse electric (TE) and transverse magnetic (TM) modes of the waveguide occur over different frequency ranges, see Fig. 1(c). This enables the probe field E to be guided in the TM mode and couple resonantly to the |g − |e transition, while the

FIG. 2. EIT in atomic medium with uniform interactions
Hss = U j =l σ j ss σ l ss leads to number correlated transparency windows. (a) For a single photon, EIT occurs at the normal two-photon resonance δ = 0. (b) For two photons, EIT occurs at a different two-photon detuning δ = U , in which case transfer to the |ss two-excitation state is resonant.
Ei 2 for an EIT system as a function of the strength of the interaction U and the two-photon detuning δ. The system has optical depth D = 400 and is continuously pumped with |Ei| 2 = 10 −4 Γ and ∆ = 0. (Note that δ = U can be achieved by adjusting the probe or control field frequencies, here we have adjusted the control as in this case this results in higher two-photon coincidence counts for |U | > 0.5). For weak input intensity, I1 corresponds to the transmittance of single photons, as the contribution from higher photon number states is negligible. The maximum transmittance occurs when twophoton resonance is satisfied, and is independent of the twobody interaction strength U . The coincidence rate I2 shows a maximum when δ ∼ U , indicating a maximum transmittance of two photons under this condition. Other simulation parameters used are Ω = 2Γ , Γ1D = 2Γ and Na = 100.
long-range atom-atom interactions are mediated on the |s −|d transition by the TE mode that has a bandgap at the same frequency. At the same time the control beam Ω may illuminate the atoms from free space.
With the control laser switched off (Ω = 0), and atoms initialized in state |g , probe photons entering the atomic ensemble encounter a two level medium. When the probe beam frequency ω p is on resonance with the atomic transition (∆ = ω p − ω ge = 0), the probe beam transmission is attenuated by a factor of exp(−D) after propagating through N a atoms, where D = 2N a Γ 1D /Γ is the optical depth (see Appendix A). Here Γ 1D denotes the emission rate of state |e into the guided mode of the photonic crystal, while Γ is the emission into all other modes (which in current experiments is comparable to the free-space emission rate γ) and Γ = Γ 1D + Γ is the total single-atom linewidth [40].
Instead when the control is switched on, the twophoton process of absorbing a probe photon and stimulated emission into the control field at frequency ω L transfers atomic population to the internal |s state from state |g . When this two photon process is resonant, δ = ω p − ω L − ω gs = 0 (in the absence of atom-atom interactions), interference between excitation from |g to |e and |s to |e leads to zero population in the excited state and the medium becomes transparent as loss due to spontaneous emission is no longer possible [33]. In this case the input probe photons are mapped to dark state polaritons -superpositions of electromagnetic fields and spin coherence between the |g and |s states [41]. The presence of the atomic spin component allows the speed of propagation to be controlled by the external laser field Ω, where the group velocity becomes v g = 2|Ω| 2 /(Γ 1D n) for an atomic medium with linear density n in the direction of propagation [42]. For realistic parameters, in the slow-light regime when v g c, the polariton excitations are almost completely spin wave in character, i.e., each polariton excitation effectively results in one atom being flipped to state |s . When atom-atom interactions are not present in the system, this transparency is a linear optical effect and applies for probe pulses of any photon number.
The situation changes dramatically when the photonic crystal interaction H ss = − j =l σ j ss σ l ss V (z j − z l ) is added (note that we have excluded the diagonal level shift from this Hamiltonian, which is assumed to be absorbed into the definition of our control laser frequency ω L ). In particular, one polariton excitation results in a spatially dependent shift of level |s for all surrounding atoms. This implies that a second proximal polariton would generally not match the two-photon resonance condition and its propagation is suppressed. Such an effect has already been exploited using Rydberg interactions to produce strong single-photon level nonlinearities [9][10][11][12][13][14][15]. However, in our case, when the strength and range of interaction associated with H ss are chosen appropriately, a novel phenomenon can appear that has not been previously identified. In particular, it is possible for a second polariton to propagate through, provided that its frequency and shape are specially matched to an altered two-photon resonance condition. In other words, here we will investigate further the possibility of forming correlated transparency windows that depend both on the photon number and wavepacket shape.
We first provide some intuitive examples of such behavior, before considering the physics in greater detail. We begin with the limit where the interaction is uniform over the atomic sample, which could be achieved by tuning the interaction length L to be much larger than the size of the atomic ensemble so that V (z j − z l ) ∼ −U . In this case the second photon sees a constant shift U of the level |s , independent of where the first polariton is [42]. This second photon could satisfy a modified two-photon resonance condition, and thus propagate in a transpar-ent fashion, if its frequency were also compensated by an amount U (for self-consistency, the first photon should also be shifted by frequency U so that both are transparent). More generally, we have a number correlated transparency window, where a single photon propagating by itself will be maximally transmitted through the atomic medium when δ = 0, while N p photons will be transmitted maximally when δ = (N p − 1)U , see Fig. 2 In Fig. 2(c)-(d), we plot measures of single and twophoton output from an EIT medium given an input of continuous coherent light over a range of two-photon detuning δ and interaction strength U . The output probe field E o is produced numerically for a particular input E i by simulating the full EIT system using the spin-model from Ref. [42], which we further describe below. We plot the normalized intensity which for weak input probe intensity corresponds to the transmission of single photons through the system. I 1 peaks at the normal EIT resonance δ = 0 and is independent of the interaction strength. To visualize the two photon behavior we plot the normalized two-photon coincidence which is proportional to the rate at which two photons are detected leaving the system at the same time. I 2 peaks when the two-photon transmission is resonant, which occurs when δ ∼ U and confirms the intuition presented earlier. The deviation from unit two-photon transmission for larger values of U (|U | Γ ) is due to the fact that for a continuous wave input, the two photons are entering and leaving at random times. Thus, each photon partially experiences the loss associated with the single-photon dispersion relation at a frequency displaced from the transparency condition [42].
An extension of this simple example is to consider a step potential, where the interaction is a constant U up to a certain atom separation r s and then zero at larger distances, as shown in Fig. 3(a). This situation leads to the presence of spatially correlated transparent states of photons, where two photons may propagate together with separation r < r s for an input such that δ ∼ U , or separated with r > r s when δ ∼ 0. Since the dark-state polariton excitations essentially consist of atoms excited to state |s , we can visualize this behavior by plotting the probability |ψ j,l ss | 2 of having two atoms (j, l) both excited to |s . In Fig. 3(b) we plot this quantity for the case where r > r s and δ = 0, and in steady state when the system is driven by a weak continuous probe input. We can then see that when δ = 0 the step potential leads to minimal population of polaritons separated by less than r s . On the other hand when δ ∼ U we expect the polaritons to propagate together within the step separation r s . In fact, one expects that if the wave-function of the relative coordinate of a two-photon state is prepared as a harmonic eigenmode of the square well potential ( Fig. 3(a)), it can propagate maintaining its pulse shape in that direction. This is illustrated in Fig. 3(c) for the case where the relative wave-function is initially prepared in the third harmonic. The relative coordinate maintains its shape well, while the wave-function in the center-of-mass coordinate diffuses out.
Motivated by these basic intuitive and numerical observations, we seek to better understand the dynamics of photons within spatial and frequency correlated transparency windows. To this end, a more rigorous description of the EIT system can be developed by finding the effective Hamiltonian for the polariton dynamics in the presence of the atom-atom interactions. By perturbatively solving the Heisenberg equations of motion around points where such correlated windows are created (as described in Appendix B), we find for the dark-state polariton field-operator Ψ = cos θE − √ n sin θσ gs e −ikpz . This effective Hamiltonian highlights a number of features of polariton evolution, the first being the reduced group velocity, v g = c cos 2 θ c, of the polariton in the medium that is observed in slow-light experiments, seen here as the first term in the brackets. This term is followed by one describing the dispersion of the group velocity, which is equivalent to the polariton having an effective mass in the EIT medium m = − The effective mass is a complex parameter with real part depending on a renormalized probe detuning ∆ M and the imaginary part on the spontaneous decay rate into other modes Γ . The renomalization of the detuning, ∆ M = ∆ + δ, results from the different bare two-photon detuning δ required to achieve EIT transparency for different number of photons or spatial configurations in the interacting system. For example, for the case of uniform interactions above, we require δ = U (N p − 1) for N p photon transparency and the mass becomes dependent on the number of polaritons in the system [32] (see Appendix B for further details).
The final term in the brackets describes the interaction between the polaritons that results from the levelshift caused by the atom-atom interactions. The spininteraction is also renormalized by the input bare twophoton detuning, where the interaction between the polaritons is described by V R (z − z ) = δ + V (z − z ). It is worth noting that in this case the interaction V (z − z ) does not need to be small or correspond to a level shift within the normal EIT transparency bandwidth (e.g., as assumed in Ref. [18]). Instead we only require that the renormalized V R (z − z ) fits within the EIT transparency within the spatial region where there is significant polariton population. Indeed Eq. (1) can describe different self-consistent transparency windows, for different photon number and spatial configurations, and these need not overlap.
For two photons, Eq. (1) separates into center of mass and relative coordinates, where each coordinate has a mass/diffusion term exactly like in non-interacting EIT. In particular, a well-localized wave-function experiences a similar loss as in normal EIT. However, the relative coordinate also sees the effective potential, which could, for example, support bound states as we describe in the next section. In general we can map a multi-polariton system to a collection of massive particles in a potential moving in a reference frame traveling at v g .
The effective Hamiltonian of Eq. 1 is valid around the two-photon resonance where we expect it to provide an approximate description of the full system dynamics. To confirm these dynamics we numerically simulate the full system using the spin model introduced in Ref. [42]. The spin model describes the behaviour of many atoms coupled to the photonic modes of a one-dimensional waveguide, generalizing the powerful input-output formalism of cavity QED. The key insight of the model is that the exchange of photons (wavevector k p ) between atoms via the waveguide reduces to an effective interaction between the atoms, given by The dynamics of the electromagnetic field as it propagates through the system are now completely captured in the dynamics of the interacting spin chain with Hamiltonian H spin = H wg + H atom + H pump + H ss , where H atom describes the atomic Hamiltonian in the absence of the waveguide mediated interactions and H pump describes the coupling of the atoms to input light. Below we consider the case where the input probe light is coherent, in which case the coupling to the atoms can be accounted for by using the Mollow transformation [43], and we have . Here E i (t) is a classical field that may either be constant, corresponding to continuous driving of the system, pulsed or zero when we initialize the system with a spin wave.
The output field operator is determined from the coherent input, with classical part E i and vacuum fluctations ε i , and the atomic coherence σ j ge found by solving the spin evolution. The numerical description of the system is then reduced to following the dynamics of the discrete atomic Hilbert space under H spin instead of the in-principle continuous electromagnetic field. This allows us to fully model our EIT system including arbitrary interactions between the atoms in a way that is independent of the assumptions made in deriving the effective polariton dynamics in Eq. (1), providing numerical verification of the intuitive picture of polaritons as massive interacting particles.

IV. PHOTONIC MOLECULES
We now use our effective propagation equation along with the spin model to describe a regime where photons can bind to one another in the EIT medium. Here we rely on the tunability of the interactions in the photonic crystal setting to achieve a molecular-like interaction potential between polaritons. In particular, we consider the case where our photonic crystal has two band edges close to one another (see Fig. 1(b)) to which the |s − |d transition couples. In this case the resulting effective interaction V (z j − z l ) = − |Ωs| 2 ∆ 2 s g 2 u ∆u exp(−|z j − z l |/L u ) + g 2 ∆ exp(−|z j − z l |/L ) will be the sum of the interaction due the upper (u) and lower ( ) band edges (see Appendix C). Crucially, the sign of each contribution depends on the detuning from the band edge in question, and the contributions hence have opposite signs. Furthermore, by engineering the detunings and band curvatures correctly the contributions will also have different interaction length scales, leading the total interaction to have a minimum (or maximum) at some atomic separation r 0 , as plotted in Fig. 4(a) for the case where the interaction strengths have equal magnitude G ≡ |Ωs| 2 For the case of two polaritons, with coordinates z 1 and z 2 , Eq. 1 separates into a Hamiltonian H cm = − 2 4m ∂ 2 ∂R 2 − i v g ∂ ∂R for the center-of-mass coordinate R = (z 1 +z 2 )/2, whose dispersion corresponds to that of a free massive particle, and a Hamiltonian H rel = − 2 m ∂ 2 ∂r 2 − 2 V (r) for the relative coordinate r = z 1 − z 2 . From this we can find the bound states of pairs of polaritons in the sys-tem, as shown in Fig. 4(a). Initializing the two-polariton wavefunction in the motional ground state of the potential then leads to the two polaritons propagating through the system with fixed relative position. We may also perturb the state by off-setting the initial relative position, in which case we observe phononic molecular vibrations about the ground state separation (slightly offset from r 0 due to the asymmetry of the potential). In Fig. 4(b) we plot the decaying oscillation of the relative position about the ground state separation for a full spin model simulation of a two polariton wavefunction that has an initial offset and see the good agreement with the oscillation predicted by the simple model. In Fig. 4(c) the polariton wavepacket in the two excitation manifold is plotted at the extrema of the oscillations, where we see that the relative coordinate remains tightly bound, while the wave function in the center-of-mass coordinate disperses as a free massive particle.
In an experiment this type of oscillation may be observable by spin-dependent imaging of the cold atoms; however, the most straight forward measurements are of photons output from the system into the waveguide. Using the spin model, we can simulate an experiment where a weak coherent pulse is input into the system at a frequency resonant with the two photon molecule (δ ∼ V (r 0 )) and the output fields are recorded after the photons are allowed to propagate through the system. In Fig. 5(a)-(b) we plot the single and two-photon parts of the input and output fields for a weak pulse with a Gaussian spatial envelop. Here, we see that the correlated transparency window damps most of the twophoton component of the initial Gaussian input, and what remains at the output largely coincides with the well-separated two-photon molecule. This leads to a peak in the second order correlation function Fig. 5(c), which shows large bunching due to the presence of the bound photon pair.
We can also observe the phononic oscillations of the pair of photons in the photon output. Inputting photons into the system in a Gaussian pulse as above, the photons are initially separated on average by less than the ground state separation, which initializes an oscillation in the two-photon wavefunction. The relative position of the two photons at the output depends on the ratio of the length over which the pulse travels in one oscillation period, L o , to the length of the atomic medium L s . An oscillation can then be observed at output by either adjusting the effective length of the system, for example by reducing the number of atoms by waiting for the atomic trap to decay, or by adjusting the strength G of the interaction potential, for example by adjusting the drive Rabi frequency Ω s . In the latter case, the oscillation length is proportional to 1/ √ ∆ M G. However, adjusting the strength G alone also changes the shape of the molecule, as the spatial extent of the ground state is approximately proportional to ∆ M /G. Reducing both ∆ M and G by the same factor ξ then results in the molecule keeping its shape while increasing L o by a factor of 1/ξ. In Fig. 5(d) we show a simulation of this situation, where we adjust G and ∆ M from the values used in Fig. 5(a)-(c) that initially lead to approximately one and a half oscillations over the length of the system. By reducing the values by a factor of three the exact simulation shows a full oscillation in the two photon output, i.e., L o changes from 3/2L s to 1/2L s , confirming the prediction of our simple model.

V. TECHNICAL REQUIREMENTS
Having described the main features of the two-photon molecule in the previous section, we now analyze more carefully the technical requirements in order to observe this physics. There are two main considerations, the first being the regime of validity of the effective equation (1) that describes polaritons as massive particles propagating in a potential, and the second, the extent to which dissipation affects the molecular state. Describing the polariton dynamics by Eq. (1) is equivalent to keeping only terms of up to quadratic order in an expansion of the photon dispersion relation about the EIT resonance. Neglecting higher order terms is a good approximation within the EIT window, which for |∆| > |Ω|, Γ , requires that the spread of frequencies of the relevant system dynamics remains within |Ω| 2 /|∆| of the EIT resonance. For a pulse that has length z 0 within the atomic medium (e.g., a Gaussian intensity distribution I(z) ∝ exp(−8z 2 /z 2 0 )) to be described by Eq. (1) it must have frequency width dω ∼ 4v g /z 0 < |Ω| 2 /|∆|, which can be recast as a requirement D p > 16|∆|/Γ for the amount of optical depth D p in the pulse length. This identifies a tradeoff, where we would like to increase the ratio |∆|/Γ to reduce the effect of dissipation in the effective mass term, however to do so requires an increase in the optical depth required in an experiment. Dissipation in our system comes from two main sources, one being spontaneous emission from the excited state |e that gives rise to the imaginary part of the effective mass term, and the other being the additional loss associated with introducing the photonic crystal mediated interactions. The former leads to loss as a pulse propagates over length L in the atomic medium, which in the limit of validity of the quadratic dispersion is approximately exp[−16L/(z 0 D p )]. This is the loss for linear propagation of a pulse and for two photons propagating in a molecule the loss is squared. To observe an oscillation of a photonic molecule we would then like to take a minimal value of the oscillation length L o to reduce loss. To see how L o is constrained, we consider the case where the interaction potential can be approximately described as a harmonic oscillator, with Gaussian ground state wavefunction exp[−2(r − r 0 ) 2 /(z 2 0 )] in relative coordinate space, in which case the photon molecule would oscillate with frequency ω M ≈ 32vg|∆| Γ1Dnz 2 0 . The oscillation length L o = 2πv g /ω M is reduced by increasing ω M , however this frequency is constrained to remain within the EIT window, which then results in the restriction L o > πz 0 /2. Taking the minimum value for L o and equating this with the system length L s , such that the photon molecule undergoes one full oscillation, then yields a propagation loss of exp[−16π/D p ].
As noted in Section II the loss related with the photonic crystal mediated interactions depends on the cooperativity of the atom-induced cavities in the photonic crystal. For a potential like that in Fig. 4(a), the loss rate is approximately βω M z 0 /(λC λ ) (taking r 0 ∼ z 0 ) where the proportionality constant β depends on the details of how the two potentials are combined (see Appendix C). Over one oscillation period of the photonic molecule we then expect the norm of the two photon wavefunction to be reduced by a factor of approximately exp[− 2πβ √ z0 √ λC λ ] as a result of this loss. Compared with the propagation loss, this loss increases with increasing pulse length z 0 , as larger constituent pulses require larger effective cavities to create the potential in the photonic crystal. The total loss then represents a balance between the two contributions and an optimal pulse length can be found z 0 = 2(4C λ Γ 2 /(βΓ 1D ) 2 ) 1/3 , at which point the total loss is exp(−6π(2β 2 Γ /(C λ Γ 1D )) 1/3 ) = exp(−48π/D p ). For parameters compatible with current experimental setups, Γ 1D = Γ , β ∼ 10 and C λ = 2 × 10 4 , we have z 0 ∼ 20a and the exp(−48π/D p ) ∼ 0.02. The total optical depth to observe a single oscillation with these parameters is then ∼ 140. In experimental systems to date only a few atoms have been trapped near the photonic crystals, however optical depths per atom of D a ∼ 2 have already been demonstrated [26], indicating that sufficient optical depth would be available with ∼ 70 atoms.
While we have focused on the loss of the photonic molecule above, signatures of the molecule and its oscillations can be observed in photon correlation measurements with greatly reduced requirements on the system parameters. In particular, for a weak coherent state input, the normalized second-order correlation function g (2) (τ ) will contain the salient features of the molecular dynamics even with significant two-photon losses, provided that the single-photon component decays in an even stronger manner. This strong single-photon loss can occur by exploiting the design of the correlated transparency window, and in particular enforcing that the transparency frequency for the two-photon molecule coincides with a region of the single-photon spectrum exhibiting strong absorption. This results in large loss of the single photon component, leading to the large values of g 2 (τ ) shown in Fig. 5(c), where our simulations include the loss mechanisms described above.
The behavior of g 2 (τ ) also depends on the nature of input and output. In particular, since the two-photon molecule is comprised of spatially separate photon pulses, this implies that there is a time ∼ r 0 /v g upon entering and exiting that only one of the photons is inside the medium. During this time, the two-photon component experiences loss as a single photon would. This additional loss to the two-photon wavepacket could be avoided if the pulses were first loaded into the medium and then the interactions switched when both photons had already entered the medium and a similar procedure applied at output.

VI. SUMMARY AND CONCLUSIONS
In summary we have demonstrated that the versatile platform of cold atoms coupled to photonic crystal waveguides leads to new possibilities in quantum optics. In particular, the ability to tune interactions between atoms over distances much greater than the wavelength of light allows the propagation of light through the atomic medium to be highly non-linear and non-local. In contrast with current experiments demonstrating nonlinearities due to Rydberg interactions between atoms [9][10][11][12][13][14][15], the interactions discussed here can be arranged to give level shifts of the order of the atomic linewidth over the entire atomic ensemble. This allows for correlated transparency of photons through the medium, where depending on the number of photons propagating in the system, the system is only transparent to particular frequencies and spatial configurations, for example a photon molecule. While the focus of our numerical studies here have mainly been in comparing single and two photon dynamics in the system, the intuition formed provides important insight into the solution of the multi-photon problem, which is expected to display rich many-body behavior [18,32].

VII. ACKNOWLEDGEMENTS
The authors thank E. Shahmoon In the main text we use the form D = 2N a Γ 1D /Γ to denote the resonant optical depth for N a two-level atoms. This expression can be derived using the transfer matrix formalism, where linear propagation in onedimension is well described by treating each atom as a refractive element with transmission and reflection coefficients t = 1 + r and r = −Γ 1D /(Γ 1D + Γ ) on resonance [40,44]. For Γ 1D Γ one expects that reflection is not important, and the transmitted intensity is given by the product of the individual atomic transmittance as ∼ |t| 2Na ∼ exp(−D). However when Γ 1D ∼ Γ a single atom can act as a strong reflector of light and the spacing of the atoms becomes crucial in how much light is transmitted through the sample [40,45]. In our simulations we choose the spacing between the atoms z a so that k 0 z a = 3π/2, in which case there is destructive interference between the reflections from each atom in the ensemble. The transmission is then reproduced by the expression exp(−D) for Γ 1D 0.2Γ . Under this condition the optical properties are the same for two samples with the same D = 2N a Γ 1D /Γ , and using this property we can may then reduce the number of atoms required to simulate a system with particular optical depth by increasing Γ 1D proportionally.
The linear properties in the EIT system can be even more robust, remaining the same for constant optical depth even when Γ 1D ∼ Γ . A single three-level atom has transmission and reflection coefficients t = 1 + r and r = −Γ 1D δ/[(Γ 1D + Γ − 2iδ)δ + 2iΩ 2 ] for ∆ = 0. We then find that after propagating through N a atoms with spacing k 0 z a = 3π/2 the probe field acquires a phase given by i exp[−N a Γ 1D δ/[(Γ − 2iδ)δ + 2iΩ 2 ], provided the two-photon detuning δ satisfies the condition |δ| 2|Ω| 2 /Γ 1D . As in the two level case above this phase fac-tor does not depend on Γ 1D and N a individually but on their product. Thus for probes with |δ| 2|Ω| 2 /Γ 1D we can again reduce the number of atoms and increase Γ 1D in a simulation to reduce computational requirements. From a practical point of view we can check this equivalence in our simulations by changing Γ 1D and N a appropriately and observing that the simulated behavior does not change.

Appendix B: Effective propagation of polaritons
When photons enter the EIT medium they are converted into polaritons, where the photonic excitation becomes a mixture of a photonic part and spin wave part with atomic population in the level |s . A single photon entering the system sees the EIT resonance when δ = 0 (see Fig. 1b). For more than one photon, the presence of the interaction between the |s levels shifts the EIT resonance. A shift in the EIT resonance will lead to changes in velocity and dispersion of each pulse, and with a spatially dependent potential as shown in Fig. 4a these effects will depend on the separation of polaritons in the media. However, it is not immediately clear from the Hamiltonian how and if the atom-atom interaction potential transforms into an interaction potential for polaritons. To understand more about the system we now derive an effective propagation equation for polaritons in the system.
To find the effective equation, we model the effective three-level media in the continuum limit, with linear atomic density n, coupled to the probe with slowly varying envelop E(t, z). In the frame rotating with the input frequency the system is described by the Hamiltonian [46] H = −n dz ∆ + i Γ 2 σ ee (z) + δσ ss (z)+ σ es (z)Ωe ikcz + cΓ 1D 2 σ eg (z)E(z)e ikpz + H.c.
In the slow light limit where v g = c cos 2 θ c and sin 2 θ ∼ 1, the bright state polariton has a pertubative effect on the propagation of the dark state polariton [47]. At lowest order in v g /c we get the effective propagation equation for the dark state We can make further approximations, neglecting the third spatial derivative term in the limit where the second spatial derivative terms dominate, that is when |∆ + iΓ /2| v g /z 0 for a polariton with pulse length z 0 . The ν(z) part of the second spatial derivative term can also be simplified in certain limits leading to an effective mass. Particularly, in the case of uniform interactions there is no spatial dependence and it may be taken outside of the derivative. In the photon molecule case discussed in the main text, the spatial dependence is approximately ω p z 2 /z 2 0 , which leads to terms that scale with ω p /D 2 p that can be neglected for large optical depth pulses. In this case we may replace ν(z) by ν(r 0 ).
The exact form of the effective mass then depends on the particular correlated EIT transparency we operate around, and is typically number dependent [32]. For example, when we have uniform interactions ν(z) = δ − 2U (N p −1) and when we are expanding around the transparency for N p photons, that is when δ = U (N p − 1), we have that the effective mass term m = − |Ω| 2 (2∆ M +iΓ )v 2 g depends on the renormalized detuning ∆ M = ∆+U (N p −1) (see for example Fig. 2(b)). Similarly in the photon molecule case we have ν(z) ≈ ν(r 0 ) = δ + 2V (r 0 ) and driving at the molecule transparency, i.e., δ ∼ −V (r 0 ), leads to the effective mass depending on ∆ M = ∆ − V (r 0 ). In both cases we have the renormalized detuning ∆ M = ∆+δ, and we arrive at the Hamiltonian in Eq. (1).
for G u = G we have that L /L u = α /α u . The strength of the potential is G = |Ωs| 2 ∆ 2 s g 2 u, ∆ u, ∼ 1.28Γ , which could be achieved for |Ω s | 2 /∆ 2 s ∼ 0.05, g u ∼ 2π × 2GHz, ∆ u ∼ 2π×31GHz and Γ ∼ 2π×5MHz. Furthermore the values of g u and L u are then consistent with a photonic crystal having g λ ∼ 2π × 10GHz (such as for the "alligator" structure described in Ref. [30]) and band curvature α ∼ 2 (such as for the structure described in Ref. [48]) for z a ∼ λ.