Interacting photon pulses in Rydberg medium

The understanding of dynamical evolutions of interacting photon pulses in Rydberg atomic ensemble is the prerequisite for realizing quantum devices with such system. We present an approach that efficiently simulates the dynamical processes, using a set of local functions we construct to reflect the profiles of narrowband pulses. For two counter-propagating photon pulses, our approach predicts the distinct phenomena from the widely concerned Rydberg blockade to the previously less noticed significant absorption in the anomalous dispersion regime, which can occur by respectively setting the pulse frequency to the appropriate values. Our numerical simulations also demonstrate how spatially extending photon pulses become deformed under realistic non-uniform interaction over their distributions.

The understanding of dynamical evolutions of interacting photon pulses in Rydberg atomic ensemble is the prerequisite for realizing quantum devices with such system. We present an approach that efficiently simulates the dynamical processes, using a set of local functions we construct to reflect the profiles of narrowband pulses. For two counter-propagating photon pulses, our approach predicts the distinct phenomena from the widely concerned Rydberg blockade to the previously less noticed significant absorption in the anomalous dispersion regime, which can occur by respectively setting the pulse frequency to the appropriate values. Our numerical simulations also demonstrate how spatially extending photon pulses become deformed under realistic non-uniform interaction over their distributions.
So far there had been no study about the real-time evolutions of photon pulses based on the full dynamics incorporating the general interaction between pulses as well as their realistic dissipation in the medium. A majority of the previous works about Rydberg EIT apply the steady-state propagating picture for continuous-wave (CW) light, according to which the time derivatives of the induced atomic excitations are assumed to be zero (see, e.g. the phenomenological superatom model in [7]). A process of finite sized pulses in the same medium can be much more complicated, since the atomic excitations induced by the pulses varying with time are surely timedependent too. Regarding the accompanying physical effects, the Rydberg blockade [38,39] used to be the theme of the previous researches. This well-known phenomenon takes place when atoms in the medium are under strong interaction to prevent their excitation to Rydberg levels. The processes of interacting photon pulses have also been approximated with a similar blockade model [28], i.e. the medium becomes a two-level system when photons meet inside a blockade radius and is under the EIT condition when they are away from each other. It is significantly meaningful to understand the photon-photon processes with a more realistic picture.
Here we present an approach based on the complete dynamical equations of the involved quantum fields to study the processes of interacting photon pulses in Rydberg medium. A main difference from the blockade potential model [28] and the quantum spin model [40], which describe the interaction mediated by Rydberg excitations as the potentials of either constant or discretely distributing values, is that we consider ever-changing interaction throughout pulse evolutions. In reality photon pulses evolve continuously in space and time without a clear-cut boundary like Rydberg blockade radius and, as we will show below, a variable interaction between pulses can lead to totally different phenomenon (significant pulse absorption) from the commonly concerned Rydberg blockade. More realistically, our approach goes beyond the assumption of uniform dissipation over pulses' spatial distributions [30] to capture the effects of inhomogeneous interaction magnitude.
The rest of the paper is organized as follows. In Sec. II we develop a theoretical approach using a type of local functions to the dynamical processes of interacting photon pulses in Rydberg EIT medium, starting from the exact dynamical equations of the involved quantum fields. The model used for illustration and the numerical simulation procedure are explained in details. The simulation results described in Sec. III mainly concern the dynamical evolutions of counter-propagating photon pulses. We illustrate the evolution processes of such single-photon pulse pairs, which manifest speedup propagation in the Rydberg blockade regime or significant absorption in the anomalous dispersion regime according to the sign of their detuning. The evolutions of polarization fields of the decaying atomic level and the change of pulse profiles under the realistic inhomogeneous interaction are also demonstrated with examples. Finally, Sec. IV contains the conclusions.

II. DYNAMICS OF INTERACTING PHOTONS
A. General two-photon process We start with two arbitrary weak light fieldsÊ l (x, t) (l = 1, 2) propagating in ensembles of the atomic level scheme in Fig. 1(a). Together with the pump beams or control fields with the Rabi frequency Ω c (t), they induce the atomic excitation distributions as the polarization fieldsP l (x, t) = √ Nσ l ge (x, t) and spinwave fieldŝ S l (x, t) = √ Nσ l gr (x, t). The flip operatorsσ µν = |µ ν| representing the transitions between the levels of atoms in an ensemble of high density N are treated as continuous fields. When the two light fields propagate in parallel along the z axis, the dynamical equations of the involved quantum fields read ( = 1) in a frame rotating at the central frequency ω p of the pulse fields l = 1 and 2. Here we consider the slowly varying fieldsÊ l (x, t) with their time derivatives much smaller than their multiplications by iω p , and they couple to the atoms with a constant defined as g = µ ge ω p / 0 (µ ge is the electric dipole matrix element and 0 the vacuum permittivity). The two-photon detuning δ = ∆ p +∆ c in Eq. (3) vanishes under the EIT condition ∆ p = −∆ c [see the definitions of these detunings in Fig. 1(a)] and, via a nonlocal potential ∆(x − x ), a spinwave fieldŜ l (x, t) experiences the interaction with the other oneŜ 3−l (x, t) as well as from itself. The quantum noise operatorsζ l (x, t),η l (x, t) are introduced to preserve the commutation relation for the quantum field operators,Ô l (x, t) =P l (x, t),Ŝ l (x, t) andÊ l (x, t), in the presence of the energy level decays at the rates γ, γ . Next we restrict the above process to what happens to two single-photon pulses, whose quantum state takes the form (4) before entering the medium, where the normalized functions f 1 (x 1 ), f 2 (x 2 ) are their snapshots at t = 0 in free space, and |0 is the vacuum state for the whole system plus reservoirs. The dynamical equations (1)-(3) can be obtained from the evolution operator U (t) = T exp{−i t 0 dτ H(τ )} of a total Hamiltonian H(t) including a stochastic part that accounts for the energy level decays (see Supplementary Material for details). Note that U (t) is not an ordinary unitary operator. Under the action U (t) the two-photon state in Eq. (4) will evolve to a general form as follows: Meanwhile, the quantum state of only one of the pulses will become (6) in the absence of the other, also exhibiting the possible loss with the converted noise components.

B. Dynamical equations for photon pulses
Any two-photon process can be described by the nine evolving two-particle functions EE(x, x , t), EP (x, x , t), · · · in Eq. (5). Though the Schrödinger equations governing their dynamical evolutions (see, e.g. the supplementary material of [19]) are linear compared with the coupled nonlinear Heisenberg-Langevin equations (1)-(3) of the field operators, it is not so straightforward to solve the dynamical equations from an initial profile EE(x, x , t = 0) = f 1 (x)f 2 (x ) in Eq. (4). A main difficulty in numerically solving this initial value problem of the evolved two-particle functions OO(x, x , t) = 00|Ô 1 (x, t)Ô 2 (x , t)|11 (Ô =Ê,P andŜ) is that the size of a pulse will undergo tremendous change in the process of becoming slowly propagating wavepacket in EIT medium; see, e.g. [36] for a discussion on the similar numerics for other functions. Many previous simulations have to start from the compressed pulses already inside EIT medium, and can hardly reflect their entrance process which is especially important to two pulses going together into the medium. So far, in different previous works, the two-particle functions have been calculated only with various simplifications, such as neglecting the photon losses [25,26], using the simplified dynamical equations from adiabatically eliminating the decaying level [28], working with their steady states (setting the time derivatives of the two-particle functions to be zero in the dynamical equations) [19], and adopting the analytical continuation from the results neglecting photon loss at high detuning ∆ p [41,42].
Different from all other works, we will apply the functions defined as whereÔ l (x, t) =Ê l (x, t),P l (x, t) andŜ l (x, t) (l = 1, 2), to deal with the dynamical problem. For one of the pulses the spinwave function of this type takes the exact form It is an inner product of the general two-body state Eq. (5) with the corresponding single-particle state Eq. (6) of the other pulse freely evolving in the absence of the concerned one (only three terms in (5) appear because the operatorŜ 1 (x) kills all others). The existence of the twoparticle functions in Eq. (5) and single-particle functions in (6) renders such functions well defined in any situation, and they truly reflect the profiles of evolved pulses except for their correlations. To find the dynamical equations for the defined functions, we multiply the vector |1, 1 to the right side of each term in Eqs. (1)-(3), and 0, 1| for l = 1 and 1, 0| for l = 2 to the left side of each term, resulting in the following set of equations: The derivations of the first two equations are straightforward, but the meaning of the third should be well explained. For the term on the third line in Eq. (3), the above procedure leads to where the operatorŜ 2 (x )Ŝ 1 (x) projects the component SS(x, x , t)|0, 0 out of the evolved state U (t)|1, 1 , while another operatorŜ 2 (x ) projects out S 0 2 (x , t)|0, 0 from U (t)|0, 1 , so that only the vacuum component in the inserted complete set of wavevectors will be left after taking the inner products with them. This is equivalent to adding a complex valued two-photon detuning as the potential given the system parameters that are capable of realizing slow light to have the first integral of the pure spinwave components dominating in Eq. (8). This time-dependent effective potential for two counter-propagating pulses obviously looks like the one shown in Fig. 1(b), with its magnitude (the absolute value) going continuously to a peak value when they are separated by the shortest distance. Meanwhile, performing the same procedure to the self-interaction term on the fourth line of Eq. (3) gives an exactly zero two-photon detuning, implying no selfinteraction for single photons.
In terms of our defined functions S l (x, t) in Eq. (8), the two-spinwave function 00|Ŝ 1 (x, t)Ŝ 2 (x , t)|11 can be approximated as SS(x, x , t) ≈ S l (x, t)S 0 3−l (x , t) (l = 1 or 2) if the pulses have sufficiently narrow bandwidths (see Supplementary Material). Substituting this approximate form into Eq. (13) gives This is the only major approximation we use, and the condition for the validity of the above separable form of SS(x, x , t) (see Supplementary Material) indicates that it works for arbitrary interaction between pulses as long as they are narrowband ones compatible with the EIT medium.
The independent evolutions of the functions S 0 (13) and (14) follow the equations which are found in a similar way to Eqs. (9)-(11). These functions evolve in the absence of the pulse interaction but are still under the dissipation due to a limited EIT width and other factors. The solution of the above equations gives the exact single-particle state in Eq. (6). The uniqueness of our approach is to simultaneously solve two sets of differential equations, (9)- (11) and (15)- (17), for finding the evolutions of the functions defined in Eq. (7). The advantage of this approach in numerical calculations will be discussed below. These local functions provide a substitute of the two-particle functions for studying the evolutions of interacting pulses, at the price of dispensing with their correlations such as the entanglement discussed in [26,30,43,44]. Moreover, instead of the simple point-point potential ∆(x − x ) appearing in the equations about the two-particle functions with more spatial variables, an effective potential V ef f l (x, t) determined with the distributions of pulses should be used, since there are less number of spatial derivatives to be integrated out in solving our equations.
With respect to the currently concerned single-photon pair states, the expectation values of the involved quantum field operatorsÔ l (x, t) =Ê l (x, t),P l (x, t) and S l (x, t) (l = 1, 2) always vanish. The functions in Eq. (7), however, give the amplitudes of the quantum fields of narrowband pulses through their average occupation numbers |O l (x, t)| 2 = 1, 1|Ô † lÔ l (x, t)|1, 1 (following Eq. (S17) in Supplement Material), the electromagnetic component among which is measurable in principle. For the narrowband interacting pulses with their entanglement neglected, their second-order correlation functions, for example, g (2) (τ ) = |EE(z = L, z = L − v g τ )| 2 (the unnormalized form of a function used in [19]), where L is the medium size and v g the pulse group velocity, can be approximated as g (2) These relations directly connect the functions defined in Eq. (7) with measurable quantities.

C. Model for illustration
The purpose of the current work is to apply the above general theory to photon pulses in Rydberg EIT media, where the nonlocal interaction potential is usually a vdW one ∆(x − x ) = C 6 /|x − x | 6 . Apart from what we have described above, there may exist some other factors relevant to the concerned processes in a realistic Rydberg atomic ensemble. One of them is the non-uniform density N (x) of the atoms, which is decided by how they are trapped. The extra decoherence from collision between atoms [34,45] and due to anisotropic interaction of the D levels of Rydberg atoms [46] can exist. Moreover, the transverse profiles of the pulses will change due to diffraction, which can be depicted with an additional kinetic term in the dynamical equations [26,43].
The most important question about a two-photon process in Rydberg medium is how the propagation of a pulse can be affected under the interaction with the other. Another essential point is how the pulse profiles should change if the interaction and dissipation over them are not uniform in reality. To answer these questions, we adopt a setup illustrated in Fig. 1(c). In this parallel waveguide setup two photon pulses either travel along the same direction (co-propagation) or respectively enter the opposite tips of two pencil-shaped ensembles (counterpropagation) and, given a time-dependent control field Ω c (t), they can also be stopped inside the atomic ensembles in which the evenly distributed atoms are assumed to be motionless. The prominent longitudinal extensions of the pulses make the setup suitable for illustrating the effects from the inhomogeneous pulse interaction between the different parts. Beyond the illustrative purpose, the technical advances toward the realization of the setup have been reported in [47,48]. The calculations with this model setup are also close to those for a process inside a single atomic ensemble, where the diffraction of the Gaussian beams can be neglected in certain domain (see Fig. S1(b) in Supplementary Material or the setup proposed in [29]).
Without loss of generality the profiles of the pulses at the entries to the ensembles (z = 0 or L) are supposed to be Ω p (ρ, where Ω M p is the maximum of the photons' Rabi frequency Ω p = gE l , and t p and τ p are the time scales indicating the peak arrival and pulse duration, respectively. We consider a single transverse mode J 0 (2ν 01 ρ/d), the Bessel function of order zero with its first zero point ν 01 , while more general transverse profiles can be used in the integral of Eq. (14) to find the effective potential V ef f l (x, t) over an ensemble separated from the other one by an adjustable distance. The field profiles on the ensemble axis ρ = 0 will be illustrated as the representation of those in the whole space. A dynamical process of the interacting three-dimensional pulses in the setup is thus reduced to a problem of finding the relevant functions O 1 (z, t) and O 2 (z, t) in 1 + 1 dimensional space and time.

D. Numerics in brief
Now we come to the practical numerical calculations with the model. Using Eqs. (10) and (11), one can expand the right-hand side of Eq. (9) as Under the condition g 2 N/|Ω c | 2 1 that is capable of realizing slow light, the time derivative on the left-hand side of Eq. (9) will be absorbed into the leading term of the above, reducing the equation to the one only with a spatial derivative. This rearranged form of Eq. (9) with only one spatial derivative is discretized for finding the spatial distribution of E l (z j , t i ) (over the lattice of z j for 1 ≤ j ≤ N s ) at a specified moment t i ≤ N t , given the distribution of the polarization field profile P l (z j , t i ) obtained with the discretized Eqs. (10) and (11).
We simply apply a fourth-order Runge-Kutta method in the iterative procedure toward the field profiles E l , P l and S l over the N s × N t space-time grid. The distribution of the constantly updated potential V ef f l (z i ) at a specific moment, which is used determine the evolution to the next moment, is found with another set of field profiles E 0 l , P 0 l and S 0 l from a similar iteration procedure with the discretized Eqs. (15)- (17). The group velocities v g,1 (t) and v g,2 (t) of the evolving pulses can be directly read from their simulated real-time trajectories, unless the wavepackets lose their distinct contours due to a significant group velocity dispersion that can also be well simulated (see Figs. 2(c) and 2(d) below).
At each temporal point t k on the side z 1 (the ensemble entry) of the space-time grid, we successively input the quantity Ω p,l (z 1 , t k ) = gE l (z 1 , t k ) from a given pulse's temporal profile, which lead to the distributions of Ω p,l (z i , t k ) (i > 1) and P l (z i , t k ), S l (z i , t k ) (i ≥ 1) over the further points z i through the iteration procedure. It is to solve the differential equations (9)-(11) and (15)-(17) as boundary value problems, and can clearly simulate the pulses' entry into the medium. In contrast it is not straightforward to deal with the evolving nonlocal functions in Eq. (5) as boundary value problem, and their numerical calculations as an initial value problem have to consider the pulse distributions outside the medium, which overwhelm the size of the EIT medium itself and thus make the simulation less efficient.

III. SIMULATION RESULTS
For simplicity only symmetric propagations of two identical pulses will be discussed in what follows. The generalization to the situations of pulses with different group velocities and different shapes is straightforward by using the different boundary conditions for each ensemble. The factors mentioned in the beginning paragraph of Sec. II(C) can also be included by the extensions of the numerical algorithm so as to apply the approach to more realistic situations.

A. Counter-propagating photon pulses
The effects of a gradually increasing interaction between pulses can be best seen from two counterpropagating photon pulses. Fig. 2 illustrates the dynamical evolutions of two photons passing by each other under their attractive interaction. As shown in Fig. 2(a), the red-detuned (∆ p > 0) photon pulses accelerate due to the mutual interaction. Instead, in Fig. 2(c), the pulses will be almost totally absorbed on the way approaching each other, if their detuning changes the sign. Such difference can be explained with the susceptibility χ(ω p ), defined as P (ω p ) = √ N χ(ω p )Ω p (ω p ), for the central frequency component under a constant interaction potential V 0 ; see the insets of Fig. 2. In the former situation the negative potential "pulls" the central frequency component initially under the EIT condition away from the absorption peak to the regime of effective two-level system, but it "pushes" the corresponding frequency component in the latter toward the peak of absorption in the anomalous dispersion regime. When they get closer, the red-detuned photon pulses entering the two-level regime will quickly escape from the medium due to much increased group velocity. However, the same medium becomes opaque to the blue-detuned (∆ p < 0) pulses because of the significantly enhanced absorption. It is evidenced by these different scenarios that, when two photon pulses approach each other in Rydberg EIT medium, there can exist richer phenomena than the wellknown Rydberg blockade that leads to a medium of effective two-level system.
The accuracy of the numerical simulations manifests with two features in the illustrated spinwave evolutions. In Fig. 2(b) about the Rydberg blockade scenario, the spinwave reappears around the exist of the medium, upon the restoration of the three-level system after the pulses separate. We particularly adopt a gradually switchingoff control field Ω c (t) which can stop the pulses in the medium. One sees a slight tendency of storing the spinwave in Fig. 2(d) (a series of parallel platforms of remnant spinwave after the control field is turned off), indicating that the medium is still a three-level one under the interaction between the negatively detuned pulses.
A potential application of the scenarios is to implement photon switch or photon transistor [33][34][35][36][37] with a slightly modified scheme of first storing one pulse in the medium. The stored pulse can easily block the blue-detuned ones coming into the medium and let the red-detuned ones go through. Such processes do not rely on the specific forms of interaction. For example, if the stored spinwave could be focused on only one point x 0 , the corresponding point-point potential ∆(x − x 0 ) can still give rise to the effects, since the approaching photon pulses nonetheless experience a gradually increasing interaction to modify their absorption and dispersion in the same ways. These effects can even be qualitatively captured by replacing the dynamical equation (11) with a corre-sponding Gross-Pitaevskii equation about the mean fields (see Supplementary Material). The interaction potential in the Gross-Pitaevskii equation considerably differs from the exact potential in Eq. (13) and the approximate potential in Eq. (14) by quantity, but also has a similar time-dependent pattern to the one shown in Fig. 1(b). Such flexibility with interaction makes the observation of the predicted phenomena more feasible.
The EIT width for the used photon pulses narrows down with their increased detuning ∆ p , impacting on their evolutions under mutual interaction. The scenarios in Fig. 2 happen when the interaction magnitude reaches the order of 10 −1 γ, and will take place under even lower interaction potentials as shown in Fig. 3 below about a higher detuning. To a pair of highly detuned pulses getting closer from where the mutual interaction is negligible, a perturbative interaction can easily alter their absorption and dispersion so that they are more likely to speed up or to be heavily damped as in Fig. 2. To fit into the narrow EIT width, the pulses with a high ∆ p should have sufficiently long duration τ p , and such widely spreading photon pulses in the medium induce much lower interaction potentials than the corresponding point-point vdW potentials obtained by shrinking them to single points. The evolutions of the pulses with considerable sizes will also become more complicated due to the inhomogeneity of interaction (see Sec. III(D) below). All these factors imply that forming the bound states of photons as recently proposed by employing the regime of high detuning [41,42] is experimentally difficult.

B. Further discussion on detuning signs
The different dynamical evolutions of the negatively and positively detuned pulses shown in Fig. 2 is one of our major predictions. A more interesting feature is that such difference can exist only under a gradually increasing rather than a suddenly increasing interaction potential, i.e. a potential like the solid curve instead of the dashed one in Fig. 1(b). To see the fact more clearly, we simplify the dynamical processes in Fig. 2

with a model of pulses propagating under constant external potentials.
Then the equation about the spinwaves [Eq. (11)] will reduce to the exact one with a constant V 0 . For a red-detuned pulse, the successively increased external potential magnitude turns its slow-light propagation under the EIT condition [ Fig. 3(a1)] into much faster propagation in the two-level regime [ Fig. 3(d1)]. As shown in Fig. 3(e1), the gradient of the dispersion curves at the detuning point ∆ p = 10γ lowers with the increased potential magnitude, indicating that the corresponding group velocity will go up to that of the effective two-level system. If the same series of potentials is applied to a blue-detuned pulse, its evolution can significantly differ. Fig. 3(b2) shows a complete absorption in contrast to the accelerated propagation in Fig.  3(b1). The dispersion curve for the central frequency of the pulse in Fig. 3(b2) has a negative gradient at ∆ p = −10γ [see Fig. 3(e2)], implying a heavy absorption accompanying the negative group velocity. Given the detuning |∆ p | = 10γ and the chosen EIT width in Fig. 3, rather low potentials (in the order of 10 −3 γ) are sufficient to see the difference.
The existing huge difference between the blue-and red-detuned pulses can not be predicted for two counterpropagation pulses, if their mutual interaction throughout the time is approximated by a blockade potential, the dashed one in Fig. 1(b). Interpreted with constant interaction potentials as in Fig. 3, an abrupt increase of the potential magnitude will directly turn the evolutions in Figs. 3(a1) and 3(a2) into the almost identical ones in 3(d1) and 3(d2). Under the highest interaction potential the dispersion curves of both blue and red detuning, which are illustrated in Fig. 3(e2) and Fig. 3(e1), respectively, stick together with that of the corresponding twolevel system, simply having Rydberg blockade. Therefore multi-photon CW beams, which create high and stable interaction within themselves, only exhibit the blockade behavior that can be analyzed in a steady-state framework [7]. In contrast we consider the completely dynamical processes of single-photon pulses, and the phenomena illustrated are very different. Like the spinwaves, the polarization fieldsP l (x, t) = √ Nσ l ge (x, t) are also evolving in the concerned processes. We present two numerical simulations of their evolutions in Fig. 4. The first one in Fig. 4(a) is about the ordinary EIT of one pulse going through the medium. The "gap" between two symmetric parts is the EIT window in which the dissipation from the induced polarization field is small. The polarization field inside the "gap" nonetheless changes with time (generally ∂P (x, t)/∂t = 0) as seen from the cross section views in Fig. 4(c). In the other example about two counter-propagating pulses as shown in Figs. 4(b) and 4(d), a peak value of the polarization field emerges at where the pulses are close to each other, leading to heavier damping due to their interaction. Here we consider resonant pulse with ∆ p = 0. For pulses with nonzero ∆ p = 0, their polarization field profiles P l (x, t) become asymmetric and vary with time more drastically, and it is also true to the pulses with rather narrow bandwidths like those used in Fig. 3.
The above results indicate that the approximation of adiabatically eliminating the degrees of freedom of the decaying level |e , i.e. setting ∂P (x, t)/∂t = 0 in Eq. (2), is not so suitable to pulses, though it works well for slow processes (compared with the time scale 1/γ) in atomic systems driven by CW light. The polarization field induced by a varying pulse is certainly time-dependent even with a considerable decay rate γ. A direct simulation using adiabatic elimination also shows the disappearance of the distinct evolutions for pulses with opposite-sign detunings (see Supplementary Material for an example). The complete dynamics involving the whole set of quantum fieldsÊ(x, t),P (x, t) andŜ(x, t) is therefore necessary to photon pulses.

D. Inhomogeneity of pulse interaction
Intuitively, atoms residing at different locations a pulse covers "feel" different long-range interaction from the other pulse due to their relative positions. The corresponding interaction potential V ef f l (x, t) in Eq. (13) or (14) is equivalent to a two-photon detuning that violates the EIT condition ∆ p + ∆ c = 0 at the location x and the moment t. Its non-uniformness leads to a space-time dependent dissipation of the pulses. In the concerned processes narrowband pulses are considered for reducing their losses in EIT medium and achieving good quantitative simulations with our approach described in Sec. 2B. Their large sizes in the medium make such inhomogeneity of interaction more obvious.
Here we illustrate the effects of inhomogeneous interaction with two examples in Fig. 5. The first group of plots in Figs. 5(a) and 5(b) is about two counter-propagating pulses. Since the pulses have considerable longitudinal extensions, the interaction between their front sides is much stronger than that between their back sides. We also apply a variable control field to stop the pulses. A large portion of their fronts can be absorbed under interaction, resulting in asymmetric shapes after they are stopped. The low decay rate of the Rydberg levels maintains the deformed spinwave profiles inside the medium after the control field is off. The other example in Figs. 5(c) and 5(d) illustrates the dynamical evolution of two co-propagating pulses until they are stopped together. A rather narrow pulse bandwidth corresponding to a large pulse size is used, so that the pulse dissipation is almost due to their interaction. Due to a longer interaction time for the co-propagating pulses, only a small portion of the initially induced spinwave will remain in the end.

IV. CONCLUSION
While promising possible applications in quantum information processing technology, interacting singlephoton pulses provide a clean channel to study the manybody physics of light in Rydberg medium, since they are without the self-interaction that makes their evolutions more complicated. However, the understanding of the completely dynamical processes is still rather challenging, as there exist the simultaneous evolutions of different types of quantum fields which determine the nonlocal interaction between and the dissipations of the pulses. To deal with the fully dynamical problem, we provide an approach based on the local functions defined in Eq. (7), which can well describe the pulse profiles though dispense with their correlations. Highly efficient numerical simulations of the complicated dynamical processes can be realized with the method, to capture the realistic photon losses in Rydberg EIT medium. It is also possible to extend the approach to the situations of more than two photons as recently discussed in [49,50].
An important application of the approach is to the processes of counter-propagating photon pulses. Our simulations give the complete dynamical pictures of how the pulses propagate under their mutual interaction in Rydberg EIT medium, showing that, in addition to the wellknown phenomenon of Rydberg blockade, a scenario of significant absorption in the anomalous dispersion regime can occur. This previously less noticed mechanism can manifest under perturbative interaction and adds to the methods of controlling photon transmission in Rydberg medium. One of its possible applications in quantum information processing technology is implementing a highly efficient photon switch. Furthermore, we have demonstrated how the profiles of pulses will change under realistic inhomogeneous interaction, as it could be relevant to a quantum memory storing photon pulses under their mutual interaction. These predictions about the dynamical pulse evolutions and photon-photon interaction could be valuable guide for the relevant experimental researches.

A. Dynamical Equations from Stochastic Hamiltonian
In an ensemble of the atoms with the level scheme shown in Fig. S1(a), an input light field (electromagnetic field)Ê(x, t) induces the polarization fieldP (x, t) = √ Nσ ge (x, t) and spinwave fieldŜ(x, t) = √ Nσ gr (x, t), in the presence of an control field Ω c (t). The flip operatorsσ µν = |µ ν| of the atomic excitations distributing over the ensemble with a high density N can be treated as continuous fields. The quantum fieldsÔ(x, t) = P (x, t),Ŝ(x, t) induced by weak light satisfy the commutation relation [Ô(x, t),Ô † (x , t)] = δ(x − x ) (the corresponding commutation relation for the fieldÊ(x, t) is up to a constant before the delta function), whereÔ(x) = 1 (2π) 3/2 dkĉ(k)e ik·x at t = 0 and [ĉ(k),ĉ † (k )] = δ(k−k ) for the wavevector mode operator.
A general process of two light fields that have entered a Rydberg atomic ensemble can be described by the following Hamiltonians. First, for two parallel propagating and slowly varying (the time derivatives of the fields are much smaller than their multiplications by iω p ) Gaussian beams with their diffraction negligible in the medium [see Fig. S1(b)], there is their kinetic Hamiltonian where "+" and "−" represent co-propagation and counter-propagation, respectively. Second, the coupling of the light fields with the atoms, which are of the level scheme in Fig. S1(a), is described by the following Hamiltonian in the rotating frame with respect to the central frequency ω p of the input pulses, where g = µ ge ω p / 0 (µ ge is the electric dipole matrix element and 0 the vacuum permittivity) is the atom-field coupling constant and the detuning ∆ p is defined in Fig. S1. A similar atom-field coupling Hamiltonian for a different level scheme is given in [51]. A narrowband pulse propagates with negligible absorption under the EIT condition ∆ p + ∆ c = 0. However, under the interaction between the induced spinwave fields, the EIT condition will be violated by shifting the levels |r of the relevant atoms. Generally the interaction Hamiltonian takes the form including both mutual and self interaction parts. The consequent dissipation from populating the levels |e that decay at the rate γ can be depicted by a stochastic Hamiltonian of the coupling between the system field operators and the quantum noise field operatorsζ l (x, t) andη l (x, t) of the environment, where the remnant decay of the levels |r is considered as well.
The Heisenberg-Langevin equations in Eqs. (1)-(3) of the main text can be derived with the total Hamiltonian H(t) = H p + H af + H int + H dis . The advantage of using the stochastic Hamiltonian in Eq. (S4) is that a concerned process can be simply depicted by an evolution operator as the time-ordered exponential U (t) = T exp{−i t 0 H(τ )dτ } of the total Hamiltonian. We write the noise operator increments as dB l (x, t) =ζ l (x, t)dt, which satisfy the Ito's rules generalized from the corresponding ones in [52]. An infinitesimal increment of the polarization fields, for example, will be found as where the above Ito's rules have been applied to the second order expansion of −iH dis (t)dt and the last term can be simply reduced to −γP (x, t)dt. Thus one will obtain the three dynamical equations as the Heisenberg-Langevin equations (1)-(3) in the main text, using the increments dÔ l (x, t) = U † (t + dt, t)Ô l (x, t)U (t + dt, t) − O l (x, t) of the system operatorsÔ l =Ê l ,P l andŜ l . These Heisenberg-Langevin equations are about an abstract dynamical process between two weak light fields. In the main text we restrict the dynamical processes to those of two single-photon pulses, and apply the waveguide setup in Fig. 1(c) of the main text to illustrate their dynamical properties.

B. Evolved Pulse Pair State and Approximate Two-spinwave Function
First, we take a look at how the initial photon pair state where we have used the evolution operator discussed in the last section. The factorization of the evolution operator U (t) is given by (2.189) in [53] or can be found in the appendices of [54,55] for its single mode version. Under the unitary action U I (t) the spinwave field operators are transformed to for l = 1, 2, so that the effective Hamiltonian U † I (t, 0){H(t) − H int }U I (t, 0) inside the time-ordered exponential of the above differs from the original form in the atom-field coupling part. When this timedependent effective Hamiltonian acts on the input, the output state will generally take the form considering the accumulated action I − iU † I (τ, 0){H(τ ) − H int }U I (τ, 0)dτ of the effective Hamiltonian at each moment, which converts the initial electromagnetic fields mostly to spinwave fields while displacing the pulses with the group velocities v g,1 and v g,2 . In this output we only show the dominant component of two spinwave fields in the currently concerned slow light regime. The succeeding operation U I (t) of the interaction Hamiltonian in Eq. (S7) will transform the spinwave field operators in the above as follows: which is via a procedure like U , 0 in the state vector. Two such transformed operators leads to an extra phase via the relation because e ±iφ l (x l ,t) |0, 0 = |0, 0 for l = 1, 2. Then one will have the general form of the evolved state The numerical calculation of this general two-body state is possible by using the effective Hamiltonian U † I (t, 0){H(t) − H int }U I (t, 0) including the part in (S9), but a closed form of the function s(x, x , t) is generally non-existing. Projecting out the two-spinwave component in U (t)|1, 1 with the operatorŜ 2 (x )Ŝ 1 (x) gives the formal form of the two-spinwave function (S15) On the other hand, the two-spinwave function can be written as (S16) For the narrowband pulses used in EIT medium, we take the approximation dk 2 |k 2 k 2 | ≈ |1 2 1| so that the above expression is reduced to This approximation based on the narrow bandwidth of the pulses loses the entanglement in the exact form (S15). By permuting the pulses, the above expression can also take the form S 0 1 (x , t)S 2 (x, t). These expressions are actually identical for pulses of identical profiles and symmetric propagations, but need not to take a symmetric form for different single-photon states |1 1 = |1 2 . Substituting the approximated form (S17) into Eq. (13) of the main text gives the potential in Eq. (14) there.
How good the approximation in Eq. (S17) is can be clearly seen from the definition of the field profile functions in Eq. (7) of the main text. According to the definition, a spinwave field function will be given as after substituting the approximate form (S17) into the above. This relation holds with the approximated normalization dx |S 0 2 (x , t)| 2 ≈ 1. The first condition for this approximate normalization is that the functions including S 0 l (x, t) (l = 1, 2), which evolve according to Eqs. (15)- (17) in the main text and independently of those of the functions defined in Eq. (7), should have low loss. It is true to these freely evolved functions under the EIT condition, given sufficiently narrow bandwidth of the used pulses. The second condition is automatically satisfied since the approximate equality on the third line in (S18) is highly close to an exact one with the system parameters in our concerned situations, as evidenced by the much higher spinwave magnitudes than those of the polarization field found with the same ensemble parameters; c.f. Fig. 2 and Fig. 4 in the main text. The first approximate equality in (S18), which leaves the dominant term only, is actually irrelevant to the evolved result such as Rydberg blockade or significant absorption in Fig. 2 of the main text. Realizing slow-light propagation when pulses first enter the medium, the chosen system parameters guarantee the validity of the equality with the fact |S 0 l (x, t)| |E 0 l (x, t)|, |P 0 l (x, t)|, while the functions SS(x, x , t), SE(x, x , t) and SP (x, x , t) in the same order of magnitudes (|SS(x, x , t)| |SE(x, x , t)|, |SP (x, x , t)| under the EIT condition for both pulses) will vanish together if the interaction could destroy the slow-light propagation.

C. Approximation with Mean Field
An intuitive approach to the concerned dynamical processes is to replace the quantum fields with their mean values as in the Gross-Pitaevskii equation or the Maxwell-Bloch equation. For a pair of pulses which are individually in any single-photon state |1 l (l = 1, 2), the mean values of the corresponding system field operators 1, 1|Ô l (x, t)|1, 1 are always vanishing, so it is necessary to use the two-particle functions in Eq. (5) or the field profiles defined in Eq. (7) of the main text. However, the mean values exist for coherent states and other photonic states whose average photon numbers are on the level of single photon. Here we study the dynamical processes in Fig. 2 of the main text with this mean value approximation by identifying the quantum field profiles O l (x, t) with the expectation values Ô l (x, t) of the quantum fields, while the self-interaction of the spinwave fields is neglected. Then we will only solve one set of equations ∂ t E l (x, t) + c∂ z E l (x, t) = ig √ N P l (x, t); (S19) ∂ t P l (x, t) = −(γ + i∆ p )P l (x, t) + iΩ * c (t)S l (x, t) + ig √ N E l (x, t); (S20) in the calculations of the field profiles. Fig. S2 shows the simulation of the processes based on this classical dynamics treatment. Because the effective potential dx ∆(x − x )|S l (x , t)| 2 in Eq. (S21) also has the pattern as the one in Fig. 1(b) of the main text, the similar effects to those in Fig. 2 of the main text can manifest as well. Such approximation with the mean fields evolved under interaction gives the lower bound of the interaction potential throughout the pulse evolutions, and the approximate potential in Eq. (14) of the main text provides an upper bound for the exact one.

D. Simulation with Adiabatic Elimination
A commonly used practice to simplify the dynamical equations of the similar processes is adiabatically eliminating the degrees of freedom for the decaying intermediate level |e or the polarization fieldsP (x, t). For the slow processes compared with the time scale 1/γ, this practice is applicable to driven-dissipation systems such as atoms in cavity. In the practical applications of the concerned processes, especially in quantum information processing, all used photonic states are pulses rather than CW light. We here check how the dynamical evolution of the pulses would be seen by the practice of adiabatic elimination. We let ∂ t P l (x, t) = 0 in Eq. (10) of the main text, to get the relation Substituting this relation into Eqs. (9) and (11) of the main text, one will obtain the reduced dynamical equations and .
The simulations of pulse evolution with these equations are shown in Fig. S3. One sees from these simulation results that the evolutions of the light fields with the opposite-sign photon detunings become almost identical. Here we consider the detunings of |∆ p | = 10γ. For a lower detuning the absorption becomes dominant, while the absorption can be reduced to very low level by higher detunings. A common feature of such reduced dynamics is the disappearance of the huge difference between opposite-sign photon detunings as shown in Fig.   2 of the main text. Moreover, the damping of a pulse will primarily depend on the magnitude of the detuning ∆ p and become insensitive to the pulse duration τ p , losing another property of the pulses propagating in EIT medium.