Photon propagation through dissipative Rydberg media at large input rates

We study the dissipative propagation of quantized light in interacting Rydberg media under the conditions of electromagnetically induced transparency. Rydberg blockade physics in optically dense atomic media leads to strong dissipative interactions between single photons. The regime of high incoming photon flux constitutes a challenging many-body dissipative problem. We experimentally study in detail the pulse shapes and the second-order correlation function of the outgoing field and compare our data with simulations based on two novel theoretical approaches well-suited to treat this many-photon limit. At low incoming flux, we report good agreement between both theories and the experiment. For higher input flux, the intensity of the outgoing light is lower than that obtained from theoretical predictions. We explain this discrepancy using a simple phenomenological model taking into account pollutants, which are nearly stationary Rydberg excitations coming from the reabsorption of scattered probe photons. At high incoming photon rates, the blockade physics results in unconventional shapes of measured correlation functions.

Generally, however, because of the many-body nature of the underlying open quantum system, the problem of strongly interacting photons is challenging to solve.Brute-force analytical or numerical approaches thus far remain limited to three or fewer photons [10,23,26].In recent years, progress has been made to develop effective theories for strongly interacting Rydberg polaritons in 1D.These theories are expected to be valid in various dispersive [20,24,27,29] and dissipative [30,31] regimes.Additionally, a promising numerical algorithm has emerged, based upon a matrix-product-state ansatz and the input-output formalism [32][33][34].Here, we show that the effective and numerical methods presented in Refs.[31] and [34], respectively, enable quantitative comparisons with rEIT experimental results, and together provide new insights into the microscopic workings of these experiments.
After the demonstration of dissipative Rydberg blockade at the single-photon level [8,9,35], a natural next step is the realization of a regular train of single photons, which could find many applications in quantum information and metrology [36][37][38].Here, we address this timely and exciting problem both theoretically and experimentally.
To be more specific, we consider photons propagating through a Rydberg medium, Fig. 1, in the regime in which a probe field E is on resonance with the |g − |e transition -the so-called dissipative regime [39,40].Van der Waals interactions between Rydberg levels lead to a blockade effect, where effectively only one atom may be excited to the Rydberg level |s within a blockade radius r b .The remaining atoms within the blockade radius then act as two-level atoms scattering incoming light.In the limit of large optical depth per blockade radius OD b = ODr b /L (where OD is the total optical depth for a medium of length L), only one photon per r b can enjoy EIT and propagate through the medium without loss, while other photons are scattered at the beginning of the medium [depicted by the solid wavy red arrow in Fig. 1 (where τ b = r b /v g is the blockade time and v g the EIT group velocity in the medium), a probe pulse shape with a well-defined beginning (sharp enough) can give rise to a train of single photons.The basic idea behind this train of photons is as follows.The first photon at the leading edge of the pulse forms a polariton in the beginning of the medium, r = 0, while a second photon can enter the medium only after the first polariton has propagated at least r b into the medium, r > r b .Hence, for higher R in , there is a high probability that one or more photons are scattered at the beginning of the medium leading to a projective measurement of the position of the polariton inside the medium, making this polariton shorter in time and hence wider in frequency.Due to the finite width of the EIT transparency window, these high-bandwidth polaritons [31] can decay in the medium [depicted by the dashed wavy red arrow in Fig. 1(b)], which puts additional constraints on OD b and OD required to observe an outgoing train of single photons [31].
In this work, for the first time, we experimentally demonstrate the time traces and correlation functions of the transmitted field in the regime of high incoming photon intensity and strong interactions.Up to now, Ryd-berg blockade physics in the dissipative regime resulted in the study of antibunching for photons separated by times smaller than the blockade time τ b , |t| < τ b .Here, we show experimentally and explain theoretically qualitatively new signatures of the blockade in the two-photon correlation function g (2) (τ ) as well as in the time traces R(t).In particular, Rydberg blockade leads to a local maximum in R(t) and g (2) (τ ) outside the blockade time τ b .This hump in output intensity [shown schematically in Fig. 1(d)] and correlations comes from the interplay of blockade physics, the finite width of the EIT transparency window, and the temporal shape of the input pulse.With this in mind, we extend the serialized hardsphere model introduced in Ref. [31] to include the temporal shape of the incoming photons as well as the decoherence of the Rydberg level.We show good agreement with output time traces predicted from exact numerics based on matrix product states (MPS) [34].We explore this regime experimentally and find qualitative signs of what the theories predict.Both the theoretical model and MPS numerics differ quantitatively from the experimentally observed time traces and correlations for high incoming photon rates.We believe that these deviations between theory and experiment are due to Rydberg pollutants, i.e., additional Rydberg excitations (in |s and other nearby Rydberg states) which are created by scattered probe photons.In order to capture the effect of pollutants, we describe a simple toy model for the dynamics of the pollutants in the system.These pollutants also prevent us, as well as other rEIT experiments, from seeing multiple subsequent humps in correlation functions [Fig.1(c-d)].In particular, pollutants prevent us from accessing higher rates for which humps would be more pronounced and therefore would lead to an output train of single photons.
The remainder of the paper is organized as follows.In Sec.II, we present two modeling approaches describing dissipative Rydberg EIT at large input rates.We first present a hard-sphere serialized model, then a model based on matrix product states, and finally compare their predictions.In Sec.III, we present experimental results, compare them with the theory, and discuss measurements suggesting that in order to explain observed data we need to include pollutants.In Sec.IV, we explain in detail the source and impact of the pollutants, as well as describe a numerically tractable toy model capturing the relevant physics.This leads to the quantitative agreement between the theory and the experiment.We summarize our work and give an outlook in Sec.V.

II. THEORY OF DISSIPATIVE RYDBERG EIT
The propagation of resonant light through a medium depends on the level structure of the atoms constituting the medium.In particular, a resonant two-level medium with levels |g and |e yields exponential attenuation of the transmission intensity by a factor T = exp(−OD).
Adding a third level |s and an appropriate control field makes the medium transparent, as interference suppresses population in |e , and a dark-state polariton with slow group velocity is generated.This effect is known as electromagnetically induced transparency (EIT) [41].
Let us now consider the propagation of photons through a dense medium of interacting three-level atoms under EIT conditions.Fig. 1(a) shows the level structure of the atoms with levels |g , |e , and |s .The control field has a Rabi frequency Ω [it takes time π/Ω to do a π pulse], and γ is the fullwidth of the level |e .The output intensity can be calculated using the following theory for dissipative dynamics developed in Refs.[30,31,42].A single photon incident under EIT conditions is converted into a Rydberg polariton (approximately a Rydberg spin wave) moving at a reduced group velocity v g .In the presence of strong Rydberg-Rydberg van der Waals interactions of the form C 6 /r 6 , this Rydberg polariton destroys EIT for any subsequent photon incident within a block- , where is the single-atom EIT linewidth [24].In the limit of large blockaded optical depth OD b = OD r b L , this leads to strong dissipation and absorption of all photons incident within a blockade time, τ b = r b /v g = OD b /(2γ EIT ), after the formation of a polariton.This is shown schematically in Fig. 1(b).Ideally, this would lead to the conversion of a continuous-wave input into a train of single photons separated in free space by the decompressed blockade radius, r b c/v g .However, the propagating polariton may decay because of the finite width of the EIT window, which washes out any spectral features sharper than 1/τ EIT = γ EIT / √ OD.Based on this intuition, the approximate output intensity [within the so-called hardsphere model] may be obtained through a serialized approach in which we first determine the output due to dissipative Rydberg-Rydberg interactions for perfect singlepolariton EIT conditions, and then frequency-filter the output with a filter of width 1/τ EIT [31].In contrast, the exact simulation using MPS does not rely on such an ansatz in treating the single-polariton EIT physics.We will now describe the hard-sphere and the MPS models in more detail.
A. Hard-sphere serialized model In Ref. [31], Zeuthen et al. develop a hard-sphere model to calculate outgoing photon rates and pulse shapes for incoming photon pulses that are longer than the medium.The basic assumptions in this model are as follows.Rydberg interactions are approximated by a hard-sphere potential of size r b .The medium is considered to be homogeneous with sharp boundaries, and polaritons only form in the beginning of the medium.Under these assumptions, it is possible to compute the output photon rate and the output time trace for a Poisson-distributed input at constant average input photon rate.Throughout the manuscript, by the time trace we mean the ensemble averaged time trace, i.e., the average over many experimental realizations of time traces.At perfect EIT, because of the hard-sphere dissipative interactions, the output rate for increasing incoming rate is saturated by one photon per blockade time.The finite EIT window can be accounted for by considering the effect of the scattered photons: Once the first polariton is formed at the beginning of the medium, the next photon arriving within a blockade time τ b of formation is scattered.This projects and localizes the first polariton wavefunction, with the time-width of the polariton being determined by the timing of the first scattering event.This means that higher input rates of photons will make the polariton wavefunction more localized in time.If the narrow polaritons do not fit in the EIT window [given by 1/τ EIT ], they may decay.This decay is governed by single-polariton physics, and we account for the EIT losses by using a Gaussian filter [31].This model was shown to be accurate in the limit of large OD b , where the predicted transmission rate was compared with exact numerical simulation of two-photon dynamics [31].We review the details of this approach in Appendix A.
The model discussed above assumes a constant input rate.Here, we extend this model to account for arbitrary input pulse shapes.We consider the input photons to be well-described by a coherent state, and the temporal shape is given by a real envelope h(t) satisfying dt h 2 (t) = 1.For the sake of brevity, we relegate technical details to the Appendix A. The Tukey pulse shape h(t), which is used in the experiment, is shown schematically in Fig. 4(b), where there is a ramp over the time t rise followed by a constant input rate.We first calculate the intensity G (1) (t) = E † (t)E(t) taking into account only blockade without EIT filtering.Then, we calculate the off-diagonal correlation function G (1) (t, t ) = E † (t)E(t ) and express it in terms of intensities G (1) (t).Finally, we convolve G (1) (t, t ) with a Gaussian filter function, which enables us to estimate the effect of a finite EIT window and leads to the intensity profiles shown in Fig. 3.
In the regime of t rise τ b , the output intensity predicted by the hard-sphere model is a train of single photons only in the limit of large input rates R in 1/τ b and large OD b .In terms of timescales, the condition on large OD b corresponds to a large blockade time, τ b τ EIT .Physically, this condition means that the photons in the train, which are necessarily each shorter than τ b , fit into the EIT transparency window, which has width 1/τ EIT .We define the ratio of these two time scales, ν = τ b /τ EIT = OD b /2 √ OD as a parameter quantifying whether it is possible to observe the photon train.In this scenario, G (1) (t) would exhibit pronounced oscillations as a function of time, as shown in Fig. 1(c), where we plot the predicted time trace for τ EIT = τ b /5, i.e. ν = 5.As long as ν > 1 and R in is appropriately chosen, the hard-sphere theory predicts oscillations in G (1) (t) with the separation of the peaks approximately given by τ b .However, in the experimentally relevant regime, we have ν ≈ 1.In this case, if we attempt to raise R in above 1/τ b to obtain the train, any oscillations in G (1) (t) are washed out due to strong filtering.
An interesting feature observed in the predicted time traces of the output intensity (Fig. 3) is the appearance of a hump at the start of the output time trace for larger input photon rates, in spite of the strong EIT filtering discussed above.This hump results from the interplay of two effects present for the parameters and pulse shapes relevant to the experiment: First, the incoming intensity |h(t)| 2 increases with time which naively would lead to the monotonous increase of the outgoing intensity.Second, the impact of EIT filtering is time-dependent because it depends on the input photon rate proportional to |h(t)| 2 ; and therefore, for greater h(t), each polariton is more localized due to the position-projecting scattering of photons at the beginning of the medium.In summary, the interplay of rising incoming intensity and stronger filtering at later times may (and for our parameters does) lead to a maximum in the outgoing intensity around the time when the amplitude of the output pulse settles to an approximately constant steady-state value (i.e., around approximately t rise + τ d , where τ d is the time delay of the transmitted pulse compared to the reference pulse).For a slower rise of h(t) (i.e., larger t rise ), the hump gets smaller.Note that the hump in the time trace indicates the existence, for a continuous-wave experiment, of an optimal input photon rate where the outgoing photon rate is maximum.One indeed sees a local maximum when plotting the outgoing steady state as a function of the input rate where the interplay between dissipative interactions and EIT gives rise to a hump [31].This is also consistent with experimental observations [see Fig. 2], however, the involved physics is more complex as we will discuss in Sec.III.

B. Detection of Rydberg pollutants
As a corollary, one might consider what happens at the end of the pulse.In this region, the incoming pulse rate R in (t) decays to zero in a time ∼ t fall .Using the same logic of weaker filtering for smaller intensities, one expects the presence of a hump at the end of the output pulse.However, as we discuss in Sec.III, the experimental measurements indicate the absence of any such hump at the end of the pulse.This leads us to conjecture the role of pollutants, which explains both the amplification of the hump in the beginning and the lack of a hump at the end of the output pulse.We discuss a simple model for the pollutants and its consequences in Sec.IV.

C. MPS method
In addition to the hard-sphere model described in the previous Section, we can also numerically obtain the output time traces using a novel time-evolution technique based on MPS introduced in Ref. [34].This method, presented in greater detail in Appendix B, relies on mapping the Maxwell-Bloch equations describing the original atomic ensemble to the propagation of a quantum field through a one-dimensional waveguide coupled to atoms.One key to this mapping is the use of a much smaller number of atoms [N 100] in the waveguide system (relative to the true number of atoms), while tuning the system parameters to ensure that macroscopic properties such as the optical depth and optical depth per blockade radius remain the same.Furthermore, all of the field properties are expressed in terms of the input field and correlation functions of the atoms alone via an input-output relation, while the dynamics of the atoms interacting with the field are encoded in a quantum spin model.As a final step, the dynamics are then solved using the MPS ansatz.The ansatz relies on the fact that, in many systems, the complete Hilbert space, which grows exponentially with atom number, is not necessary for a faithful representation of the physical states that occur, and, instead, a substantially restricted set of states, those formed from matrix products, is sufficient.This method has been extremely successful in studying condensed-matter many-body problems that would be intractable using direct diagonalization, and in Ref. [34] was applied to light propagation in atomic ensembles.
Here we extend the method in Ref. [34] to propagate the density matrix of the rEIT system in time, allowing for efficient numerical simulation of the highly dissipative system we study here.
The main benefit of the MPS method is that it allows a quantitative description beyond the hard-sphere model.Specifically, the nature of EIT in the Rydberg system is captured from first principles by using three-level atoms [as shown in Fig. 1(a)] directly in the simulation, rather than applying an approximate filter function to the photon wavepacket.Furthermore the full spatial form of the Rydberg interaction can be approximated to arbitrary precision by a sum of exponential interactions that are efficiently represented within the MPS method.Other details such as inhomogeneity in the atomic cloud, arbitrary time dependence of the input beam, and losses due to spontaneous emission and pure dephasing can also be implemented directly (see Appendix B).This allows us to check the results of the more intuitive hard-sphere model and to make qualitative comparisons with experimental results.Furthermore, we expect that this method will also be useful in other regimes of rEIT where effective models are not available.

D. Comparison between the MPS method and the
hard-sphere model Fig. 3 shows time traces from the MPS and the hardsphere models for a uniform atomic cloud with all other parameters as in the experiment.We fix OD and take L = 47.2 µm.We see good agreement between MPS and the hard-sphere model for small rates and/or initial times t < t rise .For higher rates, both methods agree qualitatively, with MPS giving a more pronounced hump.Note that, without the use of any fitting parameters, the absolute suppression of the incoming photon rate R out /R in in steady state is predicted by both theories to be at the 10% level.While in this sense the two theories agree well at the order-of-magnitude level, their predictions show appreciable relative deviations as seen in Fig. 3.This confirms that we can use the intuitive picture based on the hard-sphere model to explain qualitatively, but not quantitatively, MPS numerics and experimental data.
Note that, due to t rise τ b , the hump in output time traces is mostly due to the non-monotonic relationship between steady-state input and output intensities, postulated in Ref. [31] and discussed in section III A. The visibility of a train of photons depends on ν, which in our case is ≈ 1, making only the first hump in the train potentially visible.Furthermore, since t rise τ b , there is a large uncertainty in when the train actually begins, which further washes out the hump associated with the first photon in the train.As a result, the first photon in the train has only a minor contribution to the experimentally observed hump.respectively (here ph stands for photons).Note good agreement at initial times.Also note that the agreement is better for lower incoming rates than for higher ones.We do not see multiple humps because the blockade time τ b = 0.13 µs and filtering time τEIT = √ OD/γEIT = 0.11 µs are comparable, and because the rise time of the pulse trise = 0.8 µs is much greater than τ b .

III. EXPERIMENT
Next, we review the technical details of our experiment.We start a measurement by preparing 8 × 10 4 atoms of 87 Rb trapped in an optical dipole trap, producing a cigar-shaped atomic cloud at 4 µK with the density described by where σ R = 6.5 µm and σ = 23.6 µm characterize radial and longitudinal direction.All the atoms are optically pumped into the initial ground state |g = 5S 1/2 , F = 2, m F = 2 .We focus a weak 780 nm probe laser beam (Gaussian beam waist w 0,probe = 6.7 µm) into the cloud [Fig.4(a)], coupling the ground state |g and the intermediate state |e = 5P 3/2 , F = 3, m F = 3 .To establish EIT in the system, we add a strong 480 nm control laser beam (Gaussian beam waist w 0,control = 14 µm) coupling the intermediate state |e and the Rydberg state |s = 111S 1/2 , m J = 1/2 .The control Rabi frequency is measured to be Ω/2π = 10 MHz.From this, the Rydberg blockade radius is calculated to be r b = 18.7 µm [24].For these parameters, we observe a time delay τ d ≈ 0.31 µs of the weak probe pulses, from which we estimate the optical depth of our medium to be OD = 33.
The pulse sequence of a single experimental run is depicted in Fig. 4(b).To investigate the probe propagation at high photon rates, we send a Tukey-shaped probe pulse (2 µs uptime and 0.8 µs rise and fall times) with a varying amplitude into the medium, while the control light is on to maintain EIT conditions.The transmitted probe light is collected on a combination of four singlephoton counting modules (SPCMs).Our key experimental observations are the deformation of the probe pulse shapes transmitted through the cloud and the strong dependence of this deformation on the input probe photon  (c) and (d), respectively.In both cases, we observe the appearance of an initial hump in the transmitted pulses, the width of which is on the order of τ b .At a very low input photon rate of 0.6 ph/µs, this hump is completely absent [Fig.6(a)].At this rate, we only observe weak absorption caused by the Rydberg blockade and by the decoherence of the Rydberg level.We also observe the time delay of the transmitted pulse compared to the reference pulse.Besides the initial hump, we are interested in the steady-state transmission of the outgoing probe light.For this, we consider the orange-shaded regions indicated in Figs.4(c) and (d), where the transmission becomes approximately constant, as it does over a wide range of input photon rates we measure (the non-constant nature of this region in Fig. 4(d) is discussed in the next section).Fig. 2(a) shows the extracted steady-state transmission of the Tukey pulse as a function of the incoming rate.We find that after reaching a maximum for an input rate of R in ≈ 3 ph/µs, the output photon rate saturates to a constant value [depicted by the orange-shaded time window in Figs 4(c)  and (d)].Within this time window, we calculate from the experimental data the second-order correlation function g 2 (τ ) for the outgoing photons, as shown in Fig. 5.At low input photon rates, we find the previously observed anti-bunching at τ = 0 caused by the Rydberg-blockadeinduced nonlinearity of the medium [9].For higher input photon rates, the g (2) (τ ) correlation functions exhibit two striking features.First, the width of the anti-bunching dip shrinks, while at the same time we observe maximal bunching [g (2) (τ ) > g (2) (0)] of photons at separations τ approximately equal to the blockade time.Before we compare in Sec.III A our experimental data to the results of MPS numerics introduced in Sec.II, we briefly discuss the experimental observation of Rydberg pollutants in the optical medium and how they affect our experiment.
For the highest photon rates probed in our experiments, we observe that the outgoing probe photon rate, instead of reaching a steady-state value after the initial hump, continues to decrease on a timescale unrelated to the width of the hump [Fig.4(d)].We trace this effect back to the creation of stationary Rydberg excitations that are not accounted for in the theoretical models introduced in Sec.II.We quantify the number of these pollutant atoms and their effect on the probe photon transmission in two ways [Fig.4(b)].After each Tukey pulse, we probe the medium with a second Gaussian-shaped probe pulse (σ τ = 0.5 µs with constant peak amplitude of 2.4 ph/µs) to measure how the unwanted Rydberg excitations, created during the Tukey pulse and remaining in the cloud after the initial pulse has passed, reduce the transmission of this weak test pulse.Secondly, after the end of both pulses, we switch off the control light field and immediately ionize any remaining Rydberg atoms.The produced ions are collected on a microchannel plate detector (MCP).The intensity of the Gaussian test pulse is chosen so low and its length so short [compared with the lifetime of the Rydberg states on the order of ms] that the number of detected ions is unchanged by this test pulse.
The two observables characterizing the Rydberg pollutants which we extract from these additional measurements, namely the weak test pulse transmission and the number of detected ions after field ionization, are shown in Figs.2(b) and (c), respectively, as a function of the incoming probe photon rate.Specifically we find that, together with the growing number of detected ions, the test probe pulse transmission is reduced, meaning that the pollutant atoms affect the propagation of probe photons through the polluted medium.
It is important to note that, between the Tukey pulse and the field ionization, the control light is left on for multiple microseconds, which should depump stationary Rydberg excitations created during the probe pulse from the initial |s state.The fact that we still find a significant number of ions suggests that these Rydberg atoms have undergone a state change.Our field-ionization voltage is sufficiently high to ionize Rydberg states with n > 50, ensuring that we ionize atoms over a wide range of states near the original |s state.The claim that many of these atoms have transitioned to a state that interacts with |s only weakly (or have moved outside of the control and probe beams) is supported by the relatively weak suppression of test pulse transmission they cause.A stationary atom in state |s would block a significant part of the atomic cloud (OD b > 5), resulting in strong attenuation of probe photons.Finally, we notice that the the ion number grows with R in faster than linearly, which suggests that it is, at least partially, a two (or more)-body effect.
We discuss the possible origins of these pollutants in Sec.IV, where we also introduce an effective model to simulate their influence on pulse propagation.Fig. 2 suggests that this pollution effect becomes significant for large input photon rates R in >1 ph/µs.To further quantify when this pollution becomes important, in the follow-  ing section, we compare our experimental observations to the MPS theory developed in Sec.II.

A. Comparison of theory and experiment
To compare our experimental results quantitatively with theory, we use MPS simulations.The flexibility of the MPS model allows us to treat quantitatively crucial aspects of the experiment, such as the spatial dependence of the Rydberg interaction and the non-uniform cloud density n(z) along the probe beam direction z (see Appendix B for numerical details).Executing this model with the experimental parameters, we show in Fig. 6 the comparisons with the experimental results for time traces Experimental second-order correlation function (solid) compared with its MPS-simulated counterpart (dashed) for different input photon rates.Experimental parameters used in MPS simulations are as in all other figures except γss/2π = 120 kHz, because the g (2) data was taken in slightly different experimental conditions.The Rydberg interaction is modeled by three exponentials, as described in Appendix B 1. MPS density matrix simulations used time step 0.01/γ, N = 60, and bond dimension D = 140 for input rates 0.4 ph/µs = 0.052/τ b and 2.9 ph/µs = 0.38/τ b , and D = 180 for the input rate 7.7 ph/µs = 1/τ b .at various input rates, as well as the steady-state output rate as a function of input rate.
In the time trace shown in Fig. 6(a) and for low input rates in the steady state [Fig.6(d)], we see excellent agreement between the experiment and the MPS model.However, at higher input rates [Fig.6(b)-(c) and part of Fig. 6(d)], we see the presence of a much larger initial hump and lower steady state in the experimental output relative to our numerics.Furthermore, in Fig. 6(c), a second hump at the end of the output pulse is visible in the MPS simulation, but is absent in the experiment.This suggests that, for these higher rates, the pollution described above plays a role in determining both the size of the initial hump in the output pulse and the strength of the steady-state signal.
The pollution also plays a role in explaining the relation between g (2) measured in the experiment and the corresponding MPS simulations.In Fig. 7, we show this comparison between the theory prediction and the experimental observations for three different input photon rates.We see that the theory reproduces the qualitative feature of hump size increasing with input rate, however the humps are much larger in the experiment, suggesting once again that pollution is non-negligible at high input rates.

IV. POLLUTANTS
While the results of MPS simulations presented in the previous section qualitatively reproduce the experimen- tally observed effects both in the probe pulse shape and in the steady-state correlation functions, the lack of quantitative agreement suggests that the Rydberg pollutants we register in the experiment may have a strong effect on the probe pulse transmission even at low photon rates 4 ph/µs, which is lower than what Figs.2(b,c) may suggest.On the other hand, the relatively weak reduction of the test pulse transmission points towards the fact that the Rydberg pollutants have undergone a Rydberg state change and/or that there exists a process that removes them from the path of control and probe beams.As a possible explanation for the initial source of pollutants, we suggest radiation trapping [43] of scattered probe photons as an initial creation mechanism of Rydberg pollutants, followed by interaction-induced antiblockade and Rydberg-atom [44,45] collisions.
In this model, the pollutant creation proceeds as follows.Due to the finite extent of the cloud and the large waist of the control beam, photons scattered out of the probe mode do not necessarily escape the medium, but can instead be reabsorbed in state |s .Indeed, we estimate that our atomic cloud has a transverse optical depth of ∼ 13 at its center, and given that optical depth is a rough estimate of how many times a photon is scattered before leaving the medium, we expect that the lifetime of scattered photons could be enhanced by a factor of order 10.This radiation trapping leads to additional atoms in |s that are not part of polaritons propagating in the probe direction, but are however able to block the probe photon transmission.This effect in itself is not sufficient to explain the observation of ions, as even taking into account radiation trapping, such |s state excitations would still be expected to exit the system before the ionization pulse.Instead, through this process, atoms in |s with all possible angles between pairs of them are created.In this situation, both state-changing Rydberg collisions as well as direct anti-blockade excitation of other Rydberg levels can occur on the microsecond time scale of the experiments [44,45].Atoms in these additional states are not coupled to the control light and therefore are not de-pumped.Summarizing, the radiation trapping gives rise to both (a) the creation of the pollution atoms in the |s -state [which ultimately leave the medium] and (b) the creation of stationary pollutant Rydberg states (other than state |s ), which we observe as ions after field ionization [Fig.4(c)].

A. Effective pollutant model
Simulation of the full pollution process discussed above is prohibitively difficult.Specifically, the MPS model that we have used is only efficient in describing onedimensional propagation.Treatments of the full scattering problem in three-dimensions, so that radiation trapping is fully accounted for, are possible but currently only at the level of one or two total atomic excitations in the the system [46][47][48].Furthermore, taking into account the full family of Rydberg states and interactions would lead to an explosion of the computational Hilbert space.
Instead, we develop here a toy model that includes the basic features of the pollution process.We do so by modifying the existing MPS model to include an additional atomic 'pollutant' state |p .This state is populated by the decay from state |s at rate γ sp as shown in Fig. 8 and is assumed to induce the same Rydberg blockade as atoms in state |s (in future work, it may be interesting to consider extensions where states |p and |s have different blockade radii [49]).The population of state |p can then decay back to the ground state at rate γ pg .While state |p is not meant to represent any specific Rydberg state, we take it as a proxy for the pollution process.Atoms in state |p could represent atoms in state |s that are radiation trapped outside of the probe beam [but still inside the control beam] or atoms that have changed to new Rydberg states that still have a similar blockade radius.The decay γ pg then takes into account two phenomena: (a) a final escape from the control field of the radiatively-trapped |s excitations and (b) decay of the other Rydberg pollutant states [note that since the population of the ground state in the MPS model is essentially arbitrary, this decay can also represent decay to other long-lived Rydberg states that interact with |s only weakly].
In Fig. 9(a)-(b), we show time traces generated by the MPS model with and without pollution for two different input photon rates and compare these time traces to the experimental data.Choosing decay rates of γ sp /2π = γ pg /2π = 100 kHz, the modified MPS model provides much closer agreement with the experimental data than the original simple MPS model.Despite the simplicity of this toy model, we see that it can explain the much larger initial humps seen in the experimental time traces.Furthermore, in the g (2) correlation function shown in Fig. 9(c)-(d), the addition of pollutants to the theory also increases the size of the hump, however in this case for the parameters chosen the hump becomes larger than that seen in the experiment.
The success of this toy model leads us to conclude that pollutants do indeed play a major role in the observed output field, and may be the dominant determiner of the size of the humps we see both in the time traces and in g (2) .Meanwhile, this simple model neglects effects that are likely present in the system, such as the potential intensity dependence of γ sp and γ pg .The description of such effects requires a deeper understanding of which Rydberg processes take place and lead to pollution.Given the importance of these pollution effects at high intensity, we hope that this work will motivate further experimental and theoretical studies of this phenomenon and how it may be controlled and harnessed for applications.

V. CONCLUSIONS AND OUTLOOK
In this paper, we have discussed the physics of transmission of photons at high intensities through a Rydberg medium under the conditions of electromagnetically induced transparency.We have utilized a phenomenological model that produces reasonably good qualitative predictions for the time trace of the output intensities as well as for the steady state output rate.In addition, we utilized numerical MPS techniques to obtain a quantitative simulation of the system.The results of the two theoretical models qualitatively agree with each other.The discrepancy between these simulations and the observed experimental data points to the presence of pollutants.We extend the MPS model to include a simple treatment of pollutants consisting of an additional level.We tune this model to provide a better match to the experimental results.Our work motivates further investigation of high intensity rEIT.It highlights the importance of the role pollutants play in this strongly interacting many-body system, a role that requires additional theoretical and experimental studies and that may eventually be harnessed for applications.
We can now define a general density matrix describing this system of polaritons, where e R,R ({t R ; t R }) are the elements of the density matrix, and we have defined the shorthand for the time-ordered integral Dt Dt ≡ Note that the density matrix we consider here is the one that results after the entire pulse has entered the medium.
The expression for the density matrix obtained in the hard-sphere model is a generalization of that derived in Ref. [30].In Ref. [30], the authors develop a masterequation-type approach to calculate the outgoing photonic density matrix in the limit that the incident photons, after EIT compression, fit within the medium and that the full medium is blockaded.While the first photon forms a polariton, the subsequent photons get scattered and project the wavefunction of the first polariton.Given that the scattered photons are not detected (and hence traced out in our formalism), the density matrix of the outgoing single photon is no longer pure.In fact, the coherence in the single-polariton density matrix is related to the timing of the scattering.Ref. [31] generalized Ref. [30] to the case where the pulse size is larger than the blockade radius.We now provide an intuition for the output density matrix under the assumptions of the hardsphere model.The density matrix consists of coherences between many-polariton states.The polaritons must be at least one blockade time τ b apart from each other.Let I(t R ) denote the region R i=1 [t i , t i + τ b ) in which incoming photons are scattered by a polariton state.The quantum coherence (developed between the polariton states |t R and |t R ) arises from the fact that the projections of the Rydberg polariton wave function associated with a given set of scattering events does not fully determine the position of the polaritons.Put slightly differently, we can regard the wave function of the incoming light as a coherent superposition of arrival times for the photons; those coherences that are not destroyed by a given set of scattering events are then simply mapped into coherences of the Rydberg wave function.Since we assume that the scattered photons are not detected, they are traced out in our theory.Hence, we must perform an (incoherent) integral over all sets of scattering times.To this end, let us represent the temporal region in which scattering events can take place without destroying the coherence between states |t R and |t R as I(t R ) ∩ I(t R ).For m scatterings, the coherence factor becomes I(t R )∩I(t R ) h 2 (τ )dτ m in a straightforward generalization of Ref. [30].Utilizing the above insight, the elements of the density matrix for a Fock state input, |ψ in = |n in , can be expressed as where the factor nin! (nin−R)! is the number of ways in which one can pick an ordered set of R elements (polaritons) out of the n in incoming photons.Note that, for an input state with a definite photon number, it is not possible to have any coherence in the output between states with a different number of polaritons R = R since the environment knows the number of scattered photons and thus the number of remaining polaritons.Now we can generalize the result for an input coherent state |α = n α n √ n! e −|α| 2 /2 |n , where for simplicity we assume α ∈ R + .The general density matrix element is given by We can now utilize the expression for the general density matrix element to derive expressions for correlation functions such as G (1) (t, t ).
1. Expressions for G (1) (t, t) and G (1) (t, t ) for time-varying h(t) Here we present the expressions, used to plot Fig. 3 in the main text, for G (1) (t, t) and G (1) (t, t ) for time-varying h(t).By definition, we know that The above expression can be simplified further by carrying out the integrals over all t i > t.Using Eqs.(A4) and (A5) and assuming that the rise time of the pulse t rise fits at most R r [note that always R r ≥ 1] polaritons, and that the the presence of Rydberg levels and due to the dipolar interaction between them, may be added to the master equation without altering the above results.Solving the spin dynamics itself is challenging.Indeed, in the experiment, the atomic cloud contains tens of thousand of atoms, a number that cannot be feasibly modeled numerically.The dynamics may be found by evolving directly the density matrix in time or using the quantum jump formalism [53] for up to ∼ 20 atoms.To go beyond this, we first recognize that, in many propagation experiments, the number of atoms present is not of primary importance.Instead, the important quantity is the number of atoms multiplied by their coupling strength to the probe beam.The original large number of atoms in the ensemble can then be modeled by a much smaller chain of atoms where the overall optical depth of the system and optical depth per Rydberg blockade radius are conserved.Specifically, we can divide the atom cloud along the z axis into N slices of width a centered at positions z l .Then the collection of atoms within one slice |zj −z l |<a/2 Γ j /2σ j ge may be replaced by an effective atom Γ l /2σ l ge , whose coupling to the propagating field is the sum of the couplings of the original atoms Γ l = |zj −z l |<a/2 Γ j .In this way, the optical properties of the system are preserved as long as single-effectiveatom saturation effects are not present.To avoid saturation, the number of slices must be kept sufficiently large, which can be checked in numerical simulations by verifying that observables are invariant under changes in N provided that OD and OD b are maintained constant (see Appendix B 2). Grouping the atoms in this way also allows the non-uniform distribution of atoms in the cloud to be conveniently modeled by a non-uniform coupling to each atomic slice.
By reducing the full propagation dynamics to a model of tens of atoms, the dynamics can now be solved using matrix product states.There are then two possible treatments: (1) to represent the state of the system as a pure state MPS and propagate using the quantum jump formalism [54,55] or quantum state diffusion [56], or (2) to convert the density matrix of the system to an MPS and use the Liouvillian superoperator approach to find the time dynamics [57,58].In Ref. [34], light propagation was studied exclusively using the quantum jump method, and here we use that method to find the steadystate output shown in Fig. 6(d).
On the other hand, we find that using quantum jump trajectories to find the time traces in Fig. 6(a)-(c) is inefficient, as is quantum state diffusion, due to the large number of trajectories needed to reduce statistical noise.Instead, we propagate the density matrix using the master equation for the simple spin model above).To do so, we first map the density matrix ρ to a vector [57] by identifying local vector basis states, e.g., for two-level atoms {|g g| , |g e| , |e g| , |e e|} → {|gg), |ge), |eg), |ee)}.In this basis, the density operator can be rewritten as an , where the sum is over the basis states |β j ) for each atom j.MPS states are generalizations of product states, e.g., , where instead of having a complex coefficient c βj associated with the state of each atom we have a matrix A βj .This allows for entanglement to be introduced into the state in a controlled way by increasing the size of the matrices associated with each site.
Operators acting to the left or right of ρ, as required to represent the Liouvillian, can also be mapped into our new vector space using Kronecker products.Operator products such as O j ρI j , where O j and I j are the matrix representations in the original basis of operator O and the identity acting at site j, become O j ⊗I j |ρ).Here we have defined the Kronecker product A ⊗B = A ⊗ B T for notational convenience.
In this way, the Liouvillian for the spin-model described above becomes e ika(l−j) (I j ⊗σ j eg − σ j eg ⊗I j )σ l ge ⊗I l + σ j ge ⊗I j (I l ⊗σ l eg − σ l eg ⊗I l ) + e −ika(l−j) I j ⊗σ j eg (σ l ge ⊗I l − I l ⊗σ l ge ) + (σ j ge ⊗I j − I j ⊗σ j ge )I l ⊗σ l eg where t) e ikzj (σ j eg ⊗I j − I j ⊗σ j eg ) + e −ikzj (σ j ge ⊗I j − I j ⊗σ j ge ) .

(B6)
We can then express the entire Liouvillian as a matrix product operator with site matrices given by B7) for j = 2, . . ., N − 1, and The above operators can now be applied to the MPS representing the density operator to evolve the system in time.For example, a linear expansion of the master equation could be achieved by applying 1 + dtL to |ρ) for sufficiently small time step dt, where 1 + dtL is found from the matrix product operator above by simply multiplying all rates Γ j , γ, Γ j E in by dt and adding I j ⊗I j /N to each L j .Here, instead, to allow for larger time steps, we use a Runge-Kutta 4th order method [59].The steady state may also be found by minimizing (ρ|L † L|ρ) using traditional DMRG algorithms [60,61], however we find that time evolution in the quantum jump formalism yields the steady state more efficiently for the problem at hand.In this new vector space, we use the identity tr(A † ρ) = (A|ρ) to calculate expectation values, where |A) is the MPS vector mapping of the operator A.

Approximation of Rydberg power-law interactions by a series of exponentials
The MPO presented above describes the interaction and propagation of light in a 1D channel.We now show how this model is extended to include the case where the atoms also have a third level |s with Rydberg interactions of the form V (r) = C 6 /r 6 j<l σ j ss σ l ss .The classical driving of the control laser, H control = Ω j (σ es + σ se )/2, and decoherence of |s , L ss [ρ] = γ ss M j=1 (σ j ss ρσ j ss − σ j ss σ j ss ρ/2 − ρσ j ss σ j ss /2), are trivially included in the local part of the Liouvillian [Eq.(B6)].On the other hand, the interaction term is more complicated, as power-law decays have no known compact MPO representation.However, such decays can be approximated to arbitrary precision over finite distances by sums of exponentially decaying interactions of the form V (r) ≈ j η j λ r j [62][63][64].Here we use the technique described in Ref. [62] to find approximations of the Ryd-berg interaction over the range appropriate for the atomic cloud used in the experiment using 3-6 exponentials.In doing so, we also recognize that the large strength of the Rydberg interaction at short range can lead to numerical instabilities and instead fit the sum of exponentials to a fixed core strength at short range, as shown in Fig. 10(a).This is justified, as above a certain value the Rydberg interaction detunes the |s levels to the extent they no longer play a role in the physics.

FIG. 10. (a)
The Rydberg interaction is truncated at V = 20γ and is approximated using a sum of interactions with exponential form.The specific interaction being modeled has C6/2π = 1.8782 × 10 14 Hz µm −6 .Note that we only require a good approximation at the atomic positions, and in between the potential may take arbitrary values.Here we plot the approximate potential only at the positions of the atoms assuming a spin chain with interatomic distance a = 2.5 µm.In Fig. 10(b), we calculate the output intensity given a Tukey function input with incoming photon rate R in = 10.4 ph/µs, and for different approximations of the Rydberg interaction ranging from 3 to 6 exponentials.No great difference is seen between the curves produced with 3-6 exponentials.Furthermore, we have also checked that changing the core cutoff from 20γ (which is shown in the figure) to 30γ (not plotted) makes no difference to the dynamics.

Convergence with bond dimension and number of spins
The accuracy of the MPS methods depends on how well the quantum state of the system can be approximated by an MPS of constrained bond dimension D. Furthermore, we have modeled the system consisting of thousands of atoms by tens of atoms.To test that both of these approximations allow the nature of the light propagation to be faithfully represented, we test for convergence of the observable dynamics in both the bond dimension and in the number of atoms used in the simulations.
In Fig. 11, we show how increasing the bond dimension of the MPS used in the quantum jump simulations for finding the steady-state affects the observed output steady-state intensity.For input intensities below 10 ph/µs, the output is already well approximated by MPS with D = 10.For higher intensities, larger bond dimension is required: for the maximum input rate in the experiment of 71.8 ph/µs, convergence requires D > 50, Simulations using an MPS representation of the full density matrix require a larger bond dimension for convergence.However, despite this reduction in computational efficiency, this method may still be more efficient that the quantum jump approach as no summation over trajectories is required (where tens of thousands of trajectories are typically required for the convergence of a time trace).In Fig. 12, we show the convergence of the time traces for two different input rates.For an input rate of 4.2 ph/µs, a bond dimension of 100 is sufficient, while for 10.4 ph/µs a bond dimension of 180 is required.Finally, in Fig. 13, we show the convergence of the time trace for an input rate of 10.4 ph/µs with the number of atoms (slices) used in the simulations.For smaller number of atoms, e.g., N = 30, the time trace shows an overestimate of the output photon rate, however the qualitative behavior is still present.For higher number of atoms N = 60 the time traces converge.
FIG. 1.(a) Three-level rEIT scheme, where |g , |e , and |s are ground, excited, and long-lived Rydberg states, respectively.State |e decays spontaneously at rate γ, while γss describes the decoherence of |s .(b) Schematic representation of the theoretical model in position space at two time instances.A photon pulse, incident on the medium with velocity c, propagates as Rydberg polaritons with a group velocity vg inside the medium.However, due to Rydberg blockade, only one polariton per r b can propagate without loss.All other photons are scattered (represented by the solid wavy red arrow) at the beginning of the medium whenever a polariton is already inside the medium within r b .There are additional losses in the medium (represented by the dashed wavy red arrow) due to the finite width of the EIT transparency window.(c-d) Output pulse shapes as a function of time predicted by the theory (see Sec. II), for two different choices of incoming rates Rin, blockade time τ b , and EIT filtering time τEIT: (c) A time trace for τEIT = τ b /5 and Rin = 10/τ b , which gives rise to a train of photons [31].(d) A time trace for τEIT = τ b /2 and Rin = 3/τ b , which are closer to the parameters accessible in current experiments and in this work.Instead of well separated humps, the intensity exhibits oscillations with the peaks corresponding to the photon humps in (c).

FIG. 2 .
FIG. 2. As a function of the input photon rate, the figure shows (a) output photon rate in the steady-state region of the output pulse [see Fig. 4(c,d)], (b) transmission of the weak Gaussian test pulse through the medium, and (c) the number of detected ions.Lines connecting points are only a guide to the eye.

FIG. 3 .
FIG.3.Time traces of the output pulse for a uniform cloud with an input Tukey pulse of the same shape as the one used in the experiment.Comparison between MPS (dashed lines) and effective hard-sphere model (solid lines) without free parameters.The atomic cloud is taken as uniform in both models with L = 47.2 µm, chosen to be consistent with the length of the experimental cloud [described in Sec.III].Other parameters are OD = 33, γ/2π = 6.065MHz, Ω/2π = 10 MHz and C6/2π = 1.87820 × 10 14 Hz µm 6 for the n = 111 Rydberg state.Decoherence in the Rydberg level with full width γss/2π = 40 kHz, which we extract from the transmission at low incoming rates.MPS simulations used N = 60 and bond dimensions D = 80, 120, and 160 for rates 0.6 ph/µs = 0.078/τ b , 4.2 ph/µs = 0.55/τ b , and 10.5 ph/µs = 1.4/τ b , respectively (here ph stands for photons).Note good agreement at initial times.Also note that the agreement is better for lower incoming rates than for higher ones.We do not see multiple humps because the blockade time τ b = 0.13 µs and filtering time τEIT = √ OD/γEIT = 0.11 µs are comparable, and because the rise time of the pulse trise = 0.8 µs is much greater than τ b .

FIG. 6 .
FIG. 6.Comparison of the experimental photon output with the output simulated using MPS.(a)-(c) Time traces for various input photon rates: (a) 0.6 ph/µs = 0.078/τ b , (b) 4.2 ph/µs = 0.55/τ b , and (c) 10.5 ph/µs = 1.4/τ b .From (a), we see that the reference pulse is well described using a Tukey function.(d) Steady-state output as a function of the input rate.Experimental data and theoretical curves are shown with solid and dashed lines, respectively.Blue curves indicate input pulses while the orange ones depict the output.The input pulses indicated with ÷10 have been divided by a factor of 10 for easier viewing.The Rydberg interaction is modeled by a sum of five exponentials as described in Appendix B 1. In (a)-(c), MPS density matrix simulations use time step 0.01/γ, N = 60 effective atoms, and bond dimensions D equal to (a)-(b) 100, (c) 180.The steady-state results in (d) are from quantum jump MPS simulations with time step 0.01/γ, number of effective atoms N = 70, and bond dimensions dependent on the input rate as shown in Fig. 11 in Appendix B 2.
FIG. 7.Experimental second-order correlation function (solid) compared with its MPS-simulated counterpart (dashed) for different input photon rates.Experimental parameters used in MPS simulations are as in all other figures except γss/2π = 120 kHz, because the g(2) data was taken in slightly different experimental conditions.The Rydberg interaction is modeled by three exponentials, as described in Appendix B 1. MPS density matrix simulations used time step 0.01/γ, N = 60, and bond dimension D = 140 for input rates 0.4 ph/µs = 0.052/τ b and 2.9 ph/µs = 0.38/τ b , and D = 180 for the input rate 7.7 ph/µs = 1/τ b .

FIG. 8 .
FIG. 8. (a) Illustration of radiation trapping of scattered probe photons in the atomic cloud.Re-absorption of probe photons is possible within the larger control beam.The two three-level atoms that we zoom into schematically represent a process where the left atom emits a photon (red arrow), which is then absorbed by the right atom.(b) Level scheme for the effective model we introduce to incorporate the pollutant atoms in our numerics.
FIG. 10. (a)The Rydberg interaction is truncated at V = 20γ and is approximated using a sum of interactions with exponential form.The specific interaction being modeled has C6/2π = 1.8782 × 10 14 Hz µm −6 .Note that we only require a good approximation at the atomic positions, and in between the potential may take arbitrary values.Here we plot the approximate potential only at the positions of the atoms assuming a spin chain with interatomic distance a = 2.5 µm.(b) The number of exponentials used leads to only minor differences in the time trace shown here for an input photon rate of 10.4 ph/µs.Parameters used in MPS simulations: γ/2π = 6.065MHz, OD = 33, γss/2π = 40 kHz, N = 60, D = 100.The atomic cloud has a Gaussian distribution n(z) = exp[−z 2 /(2σ 2 )] with σ = 23.6 µm.

60 FIG. 11 .
FIG. 11.Convergence of the predicted steady-state output photon rate with bond dimension D of the MPS used in the quantum jump trajectories.Parameters used in MPS simulations: γ/2π = 6.065MHz, OD = 33, γss/2π = 40 kHz, N = 70.The Rydberg interaction is modeled by five exponentials as described in Appendix B 1. The atomic cloud has a Gaussian distribution n(z) = exp[−z 2 /(2σ 2 )] with σ = 23.6µm.The vertical bars denote the statistical errors s in the averaging of intensities over all trajectories.

FIG. 12 .
FIG. 12. Convergence of the time traces with bond dimension of the MPS for the input rate of 4.2 ph/µs (left) and 10.4 ph/µs (right).Higher input photon rates require larger bond dimension, where for an input rate of 4.2 ph/µs the time trace has already converged for D = 100, while at input rate of 10.4 ph/µs convergence is not seen until D = 180.Parameters used in MPS simulations: γ/2π = 6.065MHz, OD = 33, γss/2π = 40 kHz, N = 60.The Rydberg interaction is modeled by five exponentials as described in Appendix B 1. The atomic cloud has a Gaussian distribution n(z) = exp[−z 2 /(2σ 2 )] with σ = 23.6µm.