A Weakly-Interacting Many-Body System of Rydberg Polaritons Based on Electromagnetically Induced Transparency

Bongjune Kim, Ko-Tang Chen, Shih-Si Hsiao, Sheng-Yang Wang, Kai-Bo Li, Julius Ruseckas, Gediminas Juzeliūnas, Teodora Kirova, Marcis Auzinsh, Ying-Cheng Chen, Yong-Fan Chen, and Ite A. Yu1,7,∗ Department of Physics, National Tsing Hua University, Hsinchu 30013, Taiwan Institute of Theoretical Physics and Astronomy, Vilnius University, Saulėtekio 3, 10257 Vilnius, Lithuania Institute of Atomic Physics and Spectroscopy, University of Latvia, LV-1004 Riga, Latvia Laser Centre, University of Latvia, LV-1002, Riga, Latvia Institute of Atomic and Molecular Sciences, Academia Sinica, Taipei 10617, Taiwan Department of Physics, National Cheng Kung University, Tainan 70101, Taiwan Center for Quantum Technology, Hsinchu 30013, Taiwan

Here we proposed and experimentally demonstrated an innovative idea of utilizing a medium with a high optical depth (OD) and Rydberg atoms with a low principal quantum number to create a weakly-interacting manybody system based on the EIT effect. In this work, OD ≈ 80 and n = 32 were used. We studied the EIT polaritons also known as the dark-state polaritons [26][27][28], which, in the case of the Rydberg EIT, are quasi-particles representing superpositions of photons and Rydberg coherences. The Rydberg coherence is the coherence between the atomic ground and Rydberg states. Weakly-interacting Rydberg polaritons under an OD-enhanced interaction time can be employed in the study of manybody physics, such as the Bose-Einstein condensation (BEC) [29][30][31][32][33][34][35][36].
The system of the Rydberg EIT polaritons is an ensemble of bosonic quasiparticles [27,28]. The dispersion relation and momentum distribution of the Rydberg polaritons can be associated with their effective mass and the temperature [36]. The phase shift and attenuation of the output light induced by the DDI can be viewed as a consequence of the elastic and inelastic collisions among the polaritons. Hence, the change rates of the phase shift and attenuation are related to the elastic and inelastic collision rates, respectively. Finally, the propagation time of slow light or the diffusion time of stationary light in the EIT system is exactly the interaction time of the particles.
The high-OD medium provides for the Rydberg polaritons a long interaction time of a few µs, as demonstrated here and in our previous works [17,18]. On the other hand, the collision cross section of polaritons due to the DDI is around µm 2 in a low-n Rydberg state, leading to the elastic collision rate of a few MHz at the Rydbergpolariton density of 2×10 9 cm −3 used in this work. The low-n Rydberg state, together with a moderate polariton density, ensured a weakly-interacting regime in our experiment, where the mean distance between Rydbergstate atoms was much larger than the blockade radius. This made it possible to avoid the severe loss of Rydberg polaritons due to the blockade effect. In an ensemble with a weakly-interacting nature, a high OD can result in a sufficiently large product of the interaction time and collision rate for Rydberg polaritons, enabling them to exhibit many-body phenomena [36][37][38]. Up to now, the BECs of various kinds of polaritons have all been realized in cavity systems [32][33][34][35]. The present work deals with the cavity-free high-OD medium, in which the interaction time between the Rydberg polaritons is analogous to the storage time of polaritons in a cavity with a Q factor greater than 10 9 .
In Ref. [39], we developed a mean field theory to describe weakly-interacting Rydberg polaritons and predicted the DDI-induced phase shift and attenuation. In the present experiment we employed a laser-cooled atomic cloud with an OD of 80, and systematically studied the phase shift and attenuation of an output probe light for the EIT involving a Rydberg state of |32D 5/2 . Figure 1(a) shows the transitions driven by the probe and the coupling fields in the EIT system. The experimental data are in a good agreement with the theoretical predictions, revealing that the low-n Rydberg polaritons in the high-OD medium can form a weakly-interacting many-body system. Furthermore, we varied the DDI strength via the input photon flux and measured the transverse momentum distribution of the Rydberg polaritons. A larger DDI strength caused the width of the momentum distribution to become notably smaller, indicating the thermalization process was driven by elastic collisions [36,40,41]. The observed reduction of the momentum distribution width suggested that the combination of the µm-range DDI strength and the µs-long light-matter interaction time can make the BEC of the Rydberg polaritons feasible. The mean field model developed in Ref. [39] is summarized in this paragraph. Rydberg excitations are randomly distributed and can be considered approximately as particles of an ideal gas due to their weak interaction. Thus, the nearest-neighbor distribution (NND) [42] was utilized in our theory. Using the probability function of the NND and the atom-light coupling equations of an EIT system, we derived the following analytical formulas of the steady-state DDI-induced attenuation coefficient ∆β and phase shift ∆φ: where α is the OD of the system, Γ is the spontaneous decay rate of excited state |3 , C 6 is the van der Waals coefficient in SI units of Hz·m 6 , γ 0 is the intrinsic decoherence rate of the experiment, δ = ∆ p + ∆ c is the two-photon detuning, n atom is the atomic density, and ε is a phenomenological parameter. Utilizing ε, we related the average Rydberg-state population (ρ 22 ) to the input probe and coupling Rabi frequencies (Ω p and Ω c ) as ρ 22 = εΩ 2 p /Ω 2 c . In Eq. (2), the sign of + indicates C 6 < 0 and that of − indicates C 6 > 0. To reach the above formulas, we assumed Ω 2 p Ω 2 c and r 3 B r 3 a , where r B is the blockade radius and r a is the half mean distance between Rydberg polaritons.
The attenuation coefficient β and phase shift φ of the probe field are defined by the output-input (Ω p -Ω p ) relation of Ω p = Ω p exp(iφ − β/2). Thus, where β 0 and φ 0 are the attenuation coefficient and phase shift in the absence of DDI [39]. The mean field theory predicted that both β and φ depend on Ω 2 p linearly, and the dependence of Ω 2 p comes from the Rydberg-state population ρ 22 . In addition, the slope of β or φ versus Ω 2 p is asymmetric with respect to ∆ c = 0.
We carried out the experiment in cold 87 Rb atoms that had been produced by a magneto-optical trap (MOT). There were typically 5×10 8 trapped atoms with a temperature of about 350 µK in the MOT [43,44]. Before each measurement, we optically pumped all population to a single Zeeman state |F = 2, m F = 2 at the ground level |5S 1/2 [45]. Since the probe and coupling fields were σ + -polarized, each of the energy levels |1 , |2 and |3 shown in Fig. 1(a) was a single Zeeman state. The spontaneous decay rate, Γ, of |3 is 2π×6 MHz, and |1 ↔ |3 is a cycling transition. The life time, Γ −1 2 , of |2 is about 30 µs [46], and Γ 2 γ 0 in this work. Thus, the influence of Γ 2 was negligible. We estimated n atom ≈ 5×10 10 cm −3 in the experiment, and C 6 = −2π×260 MHz·µm 6 for |32D 5/2 , m J = 5/2 by using D ϕ = 0.343 [47]. Such C 6 together with n atom can make the DDI effect unobservable under a small Ω p or a low OD [48][49][50][51]. Figure 1(b) shows the experimental setup. The probe and coupling fields counter-propagated along the major axis of the atom cloud and completely overlapped with each other. The e −1 full width of the probe and coupling beams were 130 and 250 µm, respectively. Details of the experimental setup can be found in Sec. I of the Supplemental Material. The entire setup of red and brown optical paths shown in Fig. 1(b) formed the beat-note interferometer, which measured the phase shift of the output probe field. Details of the beat-note interferometer and the phase measurement can be found in Sec. II of the Supplemental Material and our previous work in Refs. [52,53]. The brown optical paths were present only in the phase shift measurement.
The parameters of Ω c , α (OD), and γ 0 were determined experimentally with the same method used in Ref. [43]. Details of the determination procedure can be found in Sec. III of the Supplemental Material. We set Ω c = 1.0Γ and were able to maintain γ 0 around 0.012(1)Γ throughout the experiment. This γ 0 included the effects of laser frequency fluctuation, Doppler shift, and other decoherence processes that appear in the Λ-type EIT system [43].
To verify the mean field theory of weakly-interacting Rydberg polaritons, we measured the attenuation coefficient β and phase shift φ as functions of Ω 2 p , as shown in Figs. 2(a) and 2(c). A nonzero two-photon detuning, δ, can significantly affect the phase shift and attenuation of the output probe light in the EIT medium. Thus, we utilized a beat-note interferometer to carefully determine the probe frequency for δ = 0. The determination procedure is illustrated in Sec. II of the Supplemental Material. The uncertainty of δ/2π was ±30 kHz, including the accuracy of the beat-note interferometer of ±10 kHz and the long-term drift of the two-photon frequency of ±20 kHz.
At δ = 0, we applied a square input probe pulse and measured the steady-state attenuation coefficient β and phase shift φ of the output probe pulse. The measure- ment procedure and representative data can be found in Sec. IV of the Supplemental Material. Figure 2(a) [and 2(c)] shows β (and φ) versus Ω 2 p at various ∆ c , where Ω 2 p corresponds to the peak or center intensity of the input Gaussian beam. At a given ∆ c , the data points approximately formed a straight line, as expected from the theory. We fitted the data with a linear function and plotted the slope of the best fit against ∆ c , as shown in Fig. 2(b) [and 2(d)]. It can be noticed that the y-axis interceptions of best fits were scattered in Fig. 2(c). Since a change of 30 kHz in δ/2π resulted in that of 0.4 rad in phase, the uncertainty of δ made the interceptions scatter around zero.
The strengths of the DDI effect on β and φ, i.e., the slope of β versus Ω 2 p and that of φ versus Ω 2 p (denoted as χ β and χ φ ), were asymmetric with respect to ∆ c = 0 as shown in Figs. 2(b) and 2(d). Such asymmetries are expected from the theory, as demonstrated by Eqs. (1) and (2). Considering the coupling detunings of −|∆ c | and |∆ c |, χ β (or χ φ ) of −|∆ c | was always larger (or smaller) than that of |∆ c |. The physical picture of the asymmetries is discussed in Ref. [39]. We fitted the data points in Figs. 2(b) and 2(d) with the functions of χ β and χ φ , where δ and γ 0 were set to 0 and 0.012Γ, and S DDI was the only fitting parameter. The S DDI of the best fit shown in Fig. 2(b) is in good agreement with that of Fig. 2(d). Using Eq. (3) and the values of S DDI , C 6 , n atom and Ω c , we obtained ε = 0.43±0.05. Considering that ε related the average Rydberg-state population to Ω 2 p and that Ω 2 p 160 mm 221 mm 57mm 632mm 229 mm 220mm FIG. 3: Sketch of key elements for images of probe beam on the EMCCD. Values of f indicate the focal lengths of lenses. PMF, CL b , L2 and L3 here are the same as those in Fig. 1 The probe beam coming out of CL b was collimated and had the e −1 full width of 0.92 mm. We adjusted the separation between L4 and L5 such that the beam was focused on nearly the center of the atomic cloud by L2, and had the width of 39 µm at the focal point.
corresponded to the center intensity of input Gaussian beam, the measured value of S DDI is also reasonable. Thus, the experimental data verified the mean field theory of weakly-interacting Rydberg polaritons.
To observe the thermalization process in this weaklyinteracting system, we varied the DDI strength and measured the transverse momentum distribution of the Rydberg polaritons. The moveable mirror M b shown in Fig. 1(b) was installed to direct the output probe beam to the EMCCD. Figure 3 shows the key elements for the EMCCD image of the probe beam. We reduced the waist of the input probe beam by adding the lens pair of L4 and L5, as shown in the figure, and made the input photons or initial Rydberg polaritons have a large transverse momentum distribution. The strengths of the DDI effect on β and φ (i.e., χ β and χ φ ) inferred the inelastic and elastic scattering rates of Rydberg polaritons in the medium, respectively. The elastic collisions produced the favored process of thermalization, while the inelastic collisions caused the disfavored process of decay. Hence, we chose ∆ c = 1Γ to have a suitable ratio of elastic to inelastic scattering rates for the study.
The elastic collision rate was estimated with the formula of phase shift per unit time, dφ/dt, being equal to the phase shift per collision, φ c , multiplied by the collision rate, R c . Using the propagation delay time and the data point shown in Fig. 2(c) of ∆ c = 1Γ and Ω p = 0.2Γ resulted in dφ/dt = 0.64 rad/µs. Note that the experimental condition of Ω p = 0.2Γ corresponded to the Rydberg-polariton density of 2×10 9 cm −3 . Considering the hard-sphere collision, φ c =ka where k is the root-mean-square value of the relative momentum between two colliding bodies, and a is the radius of a hard sphere. In view of the Rydberg polaritons, a can be treated as the blockade radius [54], which was about 2.1 µm in our case. The momentum distribution shown in Fig. 4(d) indicatedk = 0.075 µm −1 . Thus, φ c = 0.16 rad, inferring R c = 4.0 MHz. Under such an elastic collision rate, it was feasible to observe the thermalization effect in our experiment.
We measured the output probe beam size on the EM- CCD at ∆ c = 1Γ. To avoid the lensing effect [55], we set the two photon detuning, δ 0 , corresponding to the zero phase shift, i.e., φ = 0, in the measurement. Please note φ = 0 is the condition that the phase shift due to the ordinary EIT effect at δ 0 cancels out the phase shift due to the DDI effect. Details of the measurement of the probe beam size, i.e., the transverse momentum distribution of the Rydberg polaritons, and those of the study on the lensing effect can be found in Secs. V and VI of the Supplemental Material. The image of the probe beam, in the absence of the atoms, is shown in Fig. 4(a), while the two images in the presence of the atoms are shown in Figs. 4(b) and 4(c). It can be clearly observed that a larger value of Ω 2 p , i.e., a higher density of Rydberg polaritons because of n atom ρ 22 ∝ Ω 2 p , caused the beam size to become smaller. Since transverse momentums of the Rydberg polaritons were carried by the probe photons leaving the medium, the intensity profile of the output probe beam was able to be used to derive the transverse momentum distribution [32,36]. We fitted the intensity profiles of three images with a two-dimension Gaussian function, and utilized the results of the best fits to construct the momentum distributions shown in Figs. 4(d)-4(f). As the Rydberg-polariton density increased, the elastic collision rate also increased. Due to Ω c = 1.0Γ, the values of Ω p in the cases of Figs. 4(b) and 4(c) well satisfied the perturbation condition. Thus, the propagation times of the probe light, i.e., the interaction times of the Rydberg polaritons, were approximately the same in the two cases. Under the same interaction time, the higher collision rate due to the larger Rydberg-polariton density produced a smaller width of the transverse momentum distribution or a lower effective transverse temperature, which is the expected outcome of the thermalization process.
In conclusion, we proposed utilizing a high-OD EIT medium and a low-n Rydberg state to create a weaklyinteracting system of Rydberg polaritons, which can be a platform for the study of many-body physics, such as the polariton BEC. The experimental data of the attenuation coefficient and the phase shift influenced by the DDI are consistent with the predictions of the mean field theory developed in Ref. [39]. Furthermore, we measured the transverse momentum distribution of the Rydberg polaritons in the medium. A larger probe intensity at the input, i.e., a larger Rydberg-polariton density or a higher elastic collision rate, resulted in a smaller probe intensity profile at the output, i.e., a smaller momentum distribution width or a lower effective temperature of the Rydberg polaritons. The thermalization process was observed in this weakly-interacting many-body system. The DDI strength can be controlled and easily varied by the principal quantum number of the Rydberg state, the input probe intensity, and the coupling detuning and intensity in the EIT system, offering degrees of freedom in experiments. Our proposal is supported by the experimental demonstration and opens a new avenue in the study of quasi-particles in the field of many-body physics.

I. EXPERIMENTAL SETUP
In our experiment of the Rydberg-state electromagnetically induced transparency (EIT), the probe and coupling fields were generated by a homemade diode laser and a blue laser system (Toptica TA-SHG pro), respectively. The frequency stabilizations of the probe and coupling lasers are described as follows. We utilized the injection lock scheme to stabilize the frequency of the probe laser. In the scheme, the probe laser was used as the slave, and an external-cavity diode laser (ECDL) of Toptica DLC DL pro with the wavelength of about 780 nm was employed as the master. We utilized the Pound-Drever-Hall (PDH) scheme and the saturated absorption spectroscopy to lock the frequency of the ECDL with a heated vapor cell of Rb atoms. The blue laser had the wavelength of about 482 nm. We utilized the PDH scheme and the EIT spectroscopy to lock the frequency of the blue laser with another heated vapor cell. In this EIT spectroscopy, the ECDL and blue laser fields counter-propagated. We were able to lock the sum of the two laser frequencies to the EIT transition with the root-mean-square fluctuation of around 150 kHz [1].
As shown by the red path in Fig. 1(b) of the main text, the probe field came from the first-order beam of an acousto-optic modulator (AOM) and was coupled to a polarization maintained fiber (PMF). The PMF delivered the probe field to the cold atoms. Not shown in the figure, the probe field was twice diffracted from the AOM in the double-pass scheme. The AOM was used to shape the probe pulse, and its driving rf frequency and amplitude were precisely tuned or varied by a function generator. We used a rubidium atomic standard (SRS FS752) as the external clock of the function generator. The blue path in Fig. 1(b) of the main text represents the coupling field, which was switched on and off with another AOM not * Electronic address: yu@phys.nthu.edu.tw shown in the figure. Inside the atom cloud, the probe and coupling fields counter-propagated and completely overlapped to minimize the Doppler effect. The e −1 full widths of the probe and coupling beams were 130 and 250 µm, respectively.
Before each measurement, we first turned off the magnetic field of the magnetic-optical trap (MOT) and then performed the dark MOT for 2.5 ms to increase the optical depth (OD) of the system [2,3]. When the experimental condition of a low OD was needed, we were able to adjust the dark-MOT parameter to reduce the OD of the system by about 4 folds. The cigar-shaped atom cloud had the dimension of 1.8×1.8×6.0 mm 3 [4]. The probe and coupling fields propagated along the major axis of the atom cloud. After passing through the atoms, the probe field was detected by a photomultiplier tube (PMT). A digital oscilloscope (Agilent MSO6014A) acquired the signal from the PMT and produced raw data.

II. BEAT-NOTE INTERFEROMETER FOR THE PHASE MEASUREMENT
We employed the beat-note interferometer to measure the phase shift of the probe field induced by the atoms. The concept and illustration of the beat-note interferometer can be found in Refs. [5,6]. In the absence of the dipole-dipole interaction (DDI), the phase shift of the probe field also enabled us to precisely determine the probe frequency for the zero two-photon detuning, i.e., δ = 0, at a given coupling detuning, ∆ c . This is because the phase shift is equal to τ d δ according to the non-DDI EIT theory, where τ d is the propagation delay time. The data of attenuation coefficients and phase shifts under the DDI effect presented in Fig. 2 of the main text were taken at δ = 0.
As shown by the brown path in Fig. 1(b) of the main text, the zeroth-order beam of the AOM was unblocked when we performed the phase measurement. The firstorder beam of the AOM was the probe field. We combined the probe beam and the zeroth-order beam with a 50/50 beam splitter, BS a , in the figure to form the beat note, which had a fixed and stable frequency set by the AOM's driving frequency. The PMF can ensure both beams to completely overlap and have the same polarization. Coming out of the PMF, the beat-note signal was split into two by another 50/50 beam splitter, BS b . One, named the reference beat note, did not pass through the atoms and was detected by a photo detector (PD). The other, name the probe beat note, passed through the atoms and was detected by the PMT. Both beat notes were measured simultaneously. We compared the PMT and PD output signals to determine the phase difference between the probe and reference beat notes. Since the beat-note frequency (or wavelength) was around 133 MHz (or 2.3 m), the phase difference measured by the beat note interferometer was insensitive to the change or flucutation of the optical path [5]. The insensitivity was demonstrated by the result that the phase difference measured without the atoms was a stable constant. In addition, the frequency of the zero-order beam was far away from the transition frequency. Thus, the measured phase difference was mainly contributed from the phase shift of the probe field. We subtracted the phase difference measured with the atoms from that measured without the atoms to obtain the phase shift of the probe field induced by the atoms. We carefully set the two-photon detuning to zero at a given coupling detuning (∆ c ), before measuring the additional phase shift and attenuation induced by the DDI. To determine the probe frequency for δ = 0, the experimental condition was changed to α (OD) = 22∼26 and Ω p ≈ 0.05Γ such that the DDI effect was negligible, where Γ = 2π×6 MHz is the spontaneous decay rate of the excited state of the probe transition. During the beat-note measurement, the probe beam was a square pulse and the zeroth-order beam was a continuous wave. The zeroth-order beam was far detuned by about -22Γ and had the weak Rabi frequency of around 0.2Γ. Thus, the delay times of the probe field (as well as the ODs and decoherence rates of the experimental system) with and without the zeroth-order beam were nearly the same. Beat-note signals were acquired at a speed of 2×10 9 samples per second and stored in the segmented memory by the oscilloscope. We averaged beat-note data for 936 times and determined the phase difference between the probe and reference beat notes with the sinusoidal best fits of the data. For example, the high-frequency components of the lines in Fig. S1(a) and the circles in Fig. S1(b) demonstrate the beat-note data, while the red and black lines in Fig. S1(b) are the best fits. Once obtaining the probe frequency for δ = 0, we immediately switched back to the experimental condition of OD ≈ 81 and Ω p ≥ 0.1Γ, and took the data of the attenuation coefficient and phase shift under the DDI effect as shown by Fig. 2 of the main text.

III. DETERMINATION OF EXPERIMENTAL PARAMETERS
The experimental parameters of the coupling Rabi frequency, Ω c , the optical depth (OD), α, and the intrinsic decoherence rate, γ 0 , were determined in the same methods described in Ref. [1]. As the examples of the determination methods, Figs. 7(a)-7(c) of this reference demonstrated the data of a similar Rydberg-EIT experiment.
We determined Ω c by the frequency separation of two minima, i.e. the Autler-Towns splitting, in the probe transmission spectrum measured at the coupling detuning, ∆ c , of nearly zero and OD ≈ 1. An example of the spectrum is shown in Fig. 7(a) of Ref. [1]. Suppose the frequencies of the two minima are denoted as δ − and δ + . Then, Ω c was obtained by the relation of under the condition of ∆ c Ω c . The difference between the magnitudes of δ − and δ + is given by (S2) Using the above difference, we carefully minimized |∆ c | in the measurement such that the correction term ∆ 2 c /2Ω c in Eq. (S1) was about 1.3×10 −4 Γ. The condition of ∆ c = 0 was also determined in this step. To obtain the Autler-Towns splitting, we swept the probe frequency by varying the driving frequency of the AOM shown in Fig. 1(b) of the main text. The double-pass scheme of the AOM made the input probe power relatively stable during each run of the frequency sweeping. To prevent the spectrum from being distorted by the transient effect, a sufficiently slow sweeping rate of 240 kHz/µs was utilized.
After the value of Ω c was known, we determined the OD, i.e., α, in a designated study by measuring the propagation delay time of the probe pulse. An example of the delay time measurement is shown in Fig. 7(b) of Ref. [1]. A short Gaussian input pulse with e −1 full width of 0.66 µs was used in the measurement. To minimize the DDI effect, we employed a weak probe pulse with the peak Rabi frequency of 0.05Γ. According to the non-DDI EIT theory, the delay time, τ d is given by We used the measured value of τ d to obtain the OD of the experiment in the designated study. The intrinsic decoherence rate γ 0 was determined by measuring the probe transmission at both of the twophoton and one-photon resonances, i.e., not only δ = 0 but also ∆ c = 0. In the measurement, we set Ω c = 1.0Γ and α = 24. The input Gaussian pulse of the probe field had the e −1 full width of 7 µs and the peak Rabi frequency of 0.05Γ. The low OD and the weak probe pulse were used to made the DDI effect negligible. The long probe pulse made the measured transmission as the steady-state result. According to the non-DDI EIT theory, the peak transmission T max , i.e., the transmission at δ = 0, is given by The measured value of T max inferred the intrinsic decoherence rate of the experiment. We determined γ 0 before each study and its values maintained at (11∼13)×10 −3 Γ throughout this work.

IV. MEASUREMENTS OF TRANSMISSION AND PHASE SHIFT OF THE OUTPUT PROBE FIELD
In this section, we will describe the measurement procedure of the probe transmission and phase shift under the DDI effect. To make the DDI effect significant, a high OD of 81±3 and a moderate value of Ω c of (1.0±0.03)Γ were used and kept constant in all of the DDI studies presented in the paper. The input probe field was a long square pulse with 0.1Γ ≤ Ω p ≤ 0.2Γ. As compared with Ω c , such weak Ω p makes the probe field as the perturbation. Under the perturbation limit and without any DDI effect, the attenuation coefficient and phase shift of the output probe field do not change by varying Ω 2 p . Figures 2(a) and 2(c) of the main text show the attenuation coefficient and phase shift of the output probe field as functions of Ω 2 p . There are five different values of the coupling detuning, ∆ c , in each figure. We performed a set of consecutive measurements at each value of ∆ c . In each measurement set, we carried out the following steps: (i) The experimental parameters of Ω c , α, and γ 0 were verified by the methods described in Sec. III. (ii) We set and checked a designated value of ∆ c . (iii) The probe frequency for the two-photon resonance condition, i.e., δ = 0, was searched for and determined by the beat-note interferometer illustrated in Sec. II. (iv) We measured the probe transmissions or phase shifts at different values of Ω 2 p consecutively. (v) The steps (i), (ii), and (iii) were performed in the reverse order to ensure that the experimental condition did not change. Since it took about one hour to complete the steps (iii) and (iv), we confirmed the two-photon frequency had an uncertainty of merely ±20 kHz due to the long-term drift. Because of the above steps, it can be certain that the additional attenuation and phase shift were contributed mainly from the variation of Ω 2 p , i.e. from the DDI effect, but very little from the fluctuation or change of δ and γ 0 .
The representative data of the input and output probe transmissions are shown in Fig. S2. The data were averaged 512 times by the oscilloscope. We calculated the logarithm of the mean output transmission in the time interval between 5 and 6 µs to obtain the value of −β, where β is the attenuation coefficient. The attenuation coefficients at different values of Ω 2 p and ∆ c shown in Fig. 2(a) of the main text were determined in the same way. The representative data of the reference and probe beat notes are shown in Fig. S1(a) and illustrated in Sec. II. The data look like some high-frequency components added to the signals shown in Fig. S2. These highfrequency components are the beat notes. We zoomed in each time interval of 50 ns, and determined the phase difference between the probe and reference beat notes. In the absence of the atoms, the phase difference was a stable constant of 0.71 radians as shown by the blue circles in Fig. S1(c). In the presence of the atoms, the representative data of phase difference as a function of time are shown by the red circles in Fig. S1(c). The red line in the figure is the best fit of the red circles. We subtracted the steady-state value of the red line from 0.71 radians to obtain the phase shift of the probe field induced by the atoms. The phase shifts at different values of Ω 2 p and ∆ c shown in Fig. 2(c) of the main text were determined in the same way.

V. MEASUREMENT OF TRANSVERSE MOMENTUM DISTRIBUTION OF RYDBERG POLARITONS
To observe the evidence of the DDI-induced thermalization process in the Rydberg-polariton system, we studied the transverse momentum distribution of the Rydberg polaritons by taking images of the probe beam. After entering the medium of the atoms, the transverse momentums of the input probe photons were converted to those of the Rydberg polaritons. In the medium, the DDI can induce the elastic collisions among the Rydberg polaritons and reduce the width of their transverse momentum distribution. As the probe beam was exiting the medium, the transverse momentums of the Rydberg polaritons were converted back to those of the output photons. The images of the input and output probe beam profiles revealed the initial and final momentum distributions of the Rydberg polaritons. It can be seen as the consequence of cooling after the thermalization process that the output probe beam size becomes smaller than the input one, i.e., the width of the final transverse momentum distribution becomes narrower than that of the initial one. However, the lensing effect of the atoms can also change the probe beam size. The condition of the zero phase shift, i.e. φ = 0 can avoid the lensing effect, which will be explained in Sec. VI. For each of different values of Ω 2 p , we first searched for the two-photon detuning, δ 0 , corresponding to φ = 0, and then took the image of the output probe beam profile at δ 0 . To clearly observe the thermalization effect, we made the Rydberg polaritons have a large initial transverse momentum distribution. This was achieved by reducing the waist of the input probe beam by adding the lens pair of L4 and L5 shown in Fig. 3 of the main text. With (and without) the lens pair, the e −1 full width of the probe beam at the position of atom cloud center was 39 µm (and 130 µm). As shown in Fig. 1(b) of the main text, the moveable mirror M b directed the output probe beam to the electron-multiplying charge-coupled device (EM-CCD) camera instead of to the PMT. Let's call the experimental system without (and with) the installation of L4, L5, and M b as case I (and case II). The experimental parameters of Ω c and OD of the two cases under the same settings differed a little, i.e., less than 2%. We tuned the settings slightly to make the measured Ω c and OD of the two cases become the same. Then, the measured decoherence rates, γ 0 , had no observable difference. It took about a minute to switch from one case to another. We switched back and forth between the two cases in some measurements, and the experimental parameters of each case under the same settings were unchanged throughout the day.
Using nearly the same Ω c , OD, and γ 0 in the measurements of the two cases, we also confirmed that the slope, χ β , of the attenuation coefficient, β, versus Ω 2 p in the case I was consistent with that in the case II. Figure S3 shows β as a function of Ω 2 p which was measured in the case II at the zero two-photon detuning, i.e. δ = 0. In the measurement, the input probe field was a Gaussian pulse with the e −1 full width of 12 µs, since we took the EM-CCD image with the similar input pulse. Derived from the slope of the straight line in Fig. S3, the value of χ β was (16±1)/Γ 2 . According to the data in Fig. 2 of the main text, χ β under nearly the same experimental parameters in the case I was (17±4)/Γ 2 . Considering the uncertainties, the difference between the two χ β values is acceptable.
In the case II, the probe beam size was significantly re- We set the probe frequency for δ = 0 using the method in Sec. II. Black and red lines represent the linear best fits. We determined δ0 by the intersection between the best fit and the gray dashed line. As references, blue and green lines are the theoretical predictions calculated numerically [7].
duced. Even with the Rabi frequency of 0.2Γ, the probe power was still too weak to produce reasonable beat-note signals. Therefore, the beat-note interferometer was not able to provide a sufficient accuracy for the determination of the probe phase shift in the case II. Instead, we searched for the two-photon detuning, δ 0 corresponding to the zero phase shift, i.e, φ = 0, in the case I. Since the values of S DDI of the two cases were consistent, δ 0 found in the case I can be applied to the case II. The procedure of taking the image of the output probe beam profile is summarized in the following steps: (i) The experimental parameters were determined and set to the designed values. (ii) The probe frequency for δ = 0 was determined in the case I by the method described in Sec. II. (iii) At the above experimental parameters, we measured S DDI in the case II with the method similar to that shown in Fig. S3, and confirmed that the measured value was consistent with the value of S DDI in the case I. (iv) The two-photon detuning, δ 0 , corresponding to φ = 0 was determined in the case I. Representation data are shown in Fig. S4. We fitted the data with a straight line, and the best fit determined δ 0 . The uncertainty of δ 0 was ±2π×30 kHz. (v) At δ 0 , we took images of the output probe beam profile in the case II. We repeated the steps (iv) and (v) for different values of Ω p . (vi) The value of S DDI in the case II was measured again to verify it was unchanged during the steps (iv) and (v).

VI. STUDY ON THE LENSING EFFECT
In this section, we report the study on the lensing effect of the atoms. We experimentally demonstrated that condition of the zero phase shift, i.e., φ = 0, can avoid the lensing effect. Images of the output probe beam profiles were taken by the EMCCD for the study. Since the beam profile can also be affected by the thermalization process induced by the elastic collisions between the Rydberg polaritons, we performed the study with a low OD of 13 to minimize the DDI effect. Figure S5 shows the e −1 full width of the output probe beam profile as a function of the phase shift. The width was measured by the EMCCD image in the case II. The phase shift was converted from the two-photon detuning, δ, used in the measurement with the following formula: We verified the above formula and determined the probe frequency corresponding to δ with the beat-note interferometer in the case I. When taking images, we employed the Gaussian probe pulse of the e −1 full width of 12 µs.
To produce the images of a better signal-to-noise ratio, the peak Rabi frequency of the pulse was 0.15Γ.
The gray dashed line in Fig. S5 represents the beam width measured without the presence of the atoms. It can be clearly seen that the experimental data point of φ = 0 nearly locates on the gray dashed line, indicating that the condition of φ = 0 can avoid the lensing effect of the atoms. Such a result is also expected from the theory as shown by the black line in the figure. To calculate the black line, we utilized the ray-tracing method for Gaussian beams. In addition, we considered that the atoms behaved like a gradient index (GRIN) lens of which the refractive index, n, changes with the transverse distance, r, from the symmetric or central axis in the direction of light propagation. The ABCD matrix of a GRIN lens is given by    cos(gd) sin(gd)/g −g sin(gd) cos(gd) where d is the thickness of the lens and g is defined by n(r) = n 0 − g 2 r 2 /2 where n 0 the refractive index along the central axis [8]. Since the refractive index is proportional to the phase shift, φ, of the atoms, the following relation is used.
where 1/p is the proportionality. We substituted φ/p for g in the ABCD matrix to calculate the black line in Fig. S5. The black line is the best fit of the experimental data with the fitting parameters of d = 19.5 mm, p = 300 rad·mm 2 , and z 0 = 1.8 mm, where z 0 is the distance between the input beam waist and the center of the lens medium. Thus, the experimental data of the output probe beam width as a function of the phase shift can be well described by the effect of the GRIN lens made of the atom cloud. The DDI among the Rydberg polaritons produced a phase shift at δ = 0 and red-shifted the frequency of φ = 0. To avoid the lensing effect of the atoms as well as to observe the thermalization effect due to the DDI at the high OD, we searched for the two-photon detuning, δ 0 , corresponding to φ = 0 for each DDI strength or each value of Ω 2 p . The procedure of searching for δ 0 and then taking the image has been described in Sec. V. However, there was an uncertainty of ±2π×30 kHz in δ 0 due to the accuracy of the beat-note interferometer and the longterm drift of the two-photon frequency.
To check whether this uncertainty could be significant on the DDI effect, we measured the probe beam size not only at δ 0 but also at δ 0 ± 2π×50 kHz. Figure S6 shows the measured beam widths at these two-photon detunings. As Ω p = 0.1Γ, δ 0 = −2π×85 kHz. As Ω p = 0.2Γ, δ 0 = −2π×180 kHz. Both values of δ 0 were determined in the way described in Sec. V and shown by Fig. S4. In Fig. S6, the gray dashed line is the beam width measured without the presence of the atoms, and the black and red solid lines are the linear best fits of the experimental data. It is verified by the gentle slopes of the best fits that the uncertainty in δ 0 and the lensing effect played insignificant roles in the measurement of the probe beam width. Thus, the reduction of the probe beam width is the consequence of the transverse momentum distribution of the Rydberg polaritons being narrowed the thermalization process.