#  SUPPLEMENTAL MATERIAL
Experimental Demonstration of Stationary Dark-State Polaritons Dressed by
Dipole-Dipole Interaction

Bongjune Kim,1 Ko-Tang Chen,1 Kuei-You Chen,1 Yu-Shan Chiu,1 Chia-Yu Hsu,1 Yi-
Hsin Chen,2,3 and Ite A. Yu1,3 1Department of Physics, National Tsing Hua
University, Hsinchu 30013, Taiwan
2Department of Physics, National Sun Yat-sen University, Kaohsiung 80424,
Taiwan
3Center for Quantum Science and Technology, National Tsing Hua University,
Hsinchu 30013, Taiwan

##  I Experimental Setup

In the system for the creation of the stationary dark-state polaritons (DSP)
SLP , there are two ground states $|1\rangle$
($|5S_{1/2},F$=$1,m_{F}$=$1\rangle$) and $|2\rangle$
($|5S_{1/2},F$=$2,m_{F}$=$1\rangle$) and an excited state $|3\rangle$
($|5P_{3/2},F$=$2,m_{F}$=$2\rangle$) to form the $\Lambda$-type configuration
of the electromagnetically induced transparency (EIT) effect. Since laser-
cooled 87Rb atoms were employed in the experiment, all the energy levels are
referred to those of 87Rb. The coupling (or probe) fields propagated in the
forward and backward directions, and their Rabi frequencies are denoted as
$\Omega_{c+}$ and $\Omega_{c-}$ (or $\Omega_{p+}$ and $\Omega_{p-}$),
respectively. Both coupling fields drove the transition from $|2\rangle$ to
$|3\rangle$, and both probe fields drove that from $|1\rangle$ to $|3\rangle$.
All the fields had the $\sigma+$ polarization. The peak Rabi frequency of the
input Gaussian probe pulse (denoted as $\Omega_{p0}$) ranged from
$0.07$$\Gamma$ to $0.40$$\Gamma$, where $\Gamma$ ($=2$$\pi$$\times$6 MHz) is
the spontaneous decay rate of $|3\rangle$. Throughout this work,
$\Omega_{c+}=$ 0.54$\Gamma$ during the input and output stages of the probe
pulse, and $\Omega_{c+}=\Omega_{c-}=$ 0.44$\Gamma$ during the stationary DSP.
The determination method of $\Omega_{c+}$ and $\Omega_{c-}$ will be described
in  Sec. IV.

The lasers that generated the coupling and probe fields were injection-locked
to a master laser. Because of the injection lock, the contribution of the
laser frequency fluctuations to the decoherence rate in the $\Lambda$-type EIT
system is negligible OurPRA2019 . In the forward direction, we maintained the
frequencies of $\Omega_{p+}$ and $\Omega_{c+}$ around the resonance
frequencies of their transitions ($\Delta_{+}\approx 0$), and kept their
frequency difference, i.e., their two-photon frequency, to the two-photon
resonance. In the backward direction, we also kept the two-photon frequency of
$\Omega_{p-}$ and $\Omega_{c-}$ to the two-photon resonance, but red-detuned
their frequencies by about 1.0$\Gamma$ with respect to their transition
frequencies ($\Delta_{-}=$ $-1.0$$\Gamma$). A sufficiently large difference
between the one-photon detuning in the forward direction, i.e., $\Delta_{+}$,
and that in the backward direction, i.e., $\Delta_{-}$, is necessary for the
formation of stationary DSPs in cold media SLP .

In the system of the two-photon transition (TPT) that induced the dipole-
dipole interaction (DDI), there are a ground state $|2\rangle$ which is the
same state $|2\rangle$ in the $\Lambda$-type EIT system, an excited state
$|4\rangle$ ($|5P_{3/2},F$=$3,m_{F}$=$2\rangle$), and a Rydberg state
$|5\rangle$ ($|32D_{5/2},m_{J}$=$3/2\,\&\,5/2\rangle$). The spontaneous decay
rate of $|4\rangle$ is equal to $\Gamma$, and that of $|5\rangle$ is
2$\pi$$\times$7.9 kHz or 1.3$\times$$10^{-3}$$\Gamma$. The two laser fields
drove the transition from $|2\rangle$ to $|4\rangle$ and that from $|4\rangle$
to $|5\rangle$, and they had the Rabi frequencies $\Omega_{a}$ and
$\Omega_{b}$, respectively. Throughout this work, $\Omega_{a}=$ 2.0$\Gamma$
and $\Omega_{b}=$ 1.6$\Gamma$ during the TPT. The determination method of
$\Omega_{a}$ and $\Omega_{b}$ will be described in  Sec. IV.

The frequencies of the TPT fields of $\Omega_{a}$ and $\Omega_{b}$ were
stabilized with the Pound-Drever-Hall scheme. The root-mean-square value of
the two-photon frequency fluctuation is about 150 kHz. We utilized acousto-
optic modulators (AOMs) and electro-optic modulators (EOMs) to tune
frequencies and to modulate amplitudes of the laser fields. The rf or
microwave frequency sources of the AOMs and EOMs were locked to a rubidium
frequency standard (SRS FS725). Details of the stabilization and tuning of the
laser field frequencies can be found in Ref. OurPRA2019 . We fixed the sum of
the frequencies, i.e., two-photon frequency, of $\Omega_{a}$ and $\Omega_{b}$
to the two-photon resonance, while detuning their frequencies by 5.0$\Gamma$
with respect to their transition frequencies ($\Delta_{a}=-\Delta_{b}=$
+5.0$\Gamma$).

In the $\Lambda$-type EIT system, the $e^{-1}$ full width of the input probe
beam was 130 $\mu$m, and those of the forward and backward coupling beams were
3.1 mm. The forward and backward coupling fields completely overlapped and
covered the entire atom cloud.

In the TPT system, the field of $\Omega_{a}$ propagated in the forward
direction and that of $\Omega_{b}$ in the backward direction. This counter-
propagation configuration minimized the Doppler effect in the TPT. The
$e^{-1}$ full width of the $\Omega_{a}$ beam was 4.9 mm, and that of the
$\Omega_{b}$ beam was 250 $\mu$m. The beams of $\Omega_{a}$ and $\Omega_{b}$
fully covered the region of stationary DSPs in the atom cloud.

The output probe field was coupled to a single-mode optical fiber to minimize
the background signals from the other light fields, and detected by a
photomultiplier tube (PMT). A digital oscilloscope (Agilent MSO6014A) acquired
the signal from the PMT. All the timings of $\Omega_{c+}^{2}$,
$\Omega_{c-}^{2}$, $\Omega_{a}^{2}$, and $\Omega_{b}^{2}$ were monitored by
fast photo detectors. In the measurement of the phase shift of the output
probe pulse, we employed the beat-note interferometer beatnote . The
measurement method of the phase shift was the same as that described in Ref.
OurCommunPhys2021 .

In the experiment, the optical depth (OD) of the system was maintained around
36. Using a known value of $\Omega_{c+}$ in the $\Lambda$-type EIT system, we
measured the propagation delay time, $\tau_{d}$, of the slow-light probe pulse
to determine the OD, $\alpha$, where $\tau_{d}=\alpha\Gamma/\Omega_{c+}^{2}$
OurPRA2019 . The ground-state decoherence rate, $\gamma_{\Lambda}$, affects
the coherence involving the DSP, and its value was about
9$\times$$10^{-4}$$\Gamma$ or 2$\pi$$\times$5.4 kHz in this work. We measured
the output transmission of the probe field, $T_{\rm SL}$, to determine the
value of $\gamma_{\Lambda}$, where $T_{\rm
SL}=\exp(-2\alpha\gamma_{\Lambda}\Gamma/\Omega_{c+}^{2})$ OurPRA2019 .
Examples of data for the determinations of the coupling Rabi frequency,
optical density, and decoherence rate in a $\Lambda$-type EIT system can be
found in Fig. 3 of Ref. OurPRA2019 .

##  II Theoretical Model

To make theoretical predictions, we numerically solved the optical Bloch
equations (OBEs) of the density-matrix elements of the atoms,
$\rho_{ij}(t,z)$, and Maxwell-Schrödinger equations (MSEs) of the
forward/backward probe pulses (with the Rabi frequencies of $\Omega_{p+}(t,z)$
and $\Omega_{p-}(t,z)$), where $t$ and $z$ represent the temporal and spatial
variables. There are five energy levels in the calculation consisting of two
ground states $|1\rangle$ and $|2\rangle$, two excited states $|3\rangle$ and
$|4\rangle$, and a Rydberg state $|5\rangle$. The system of the $\Lambda$-type
electromagnetically induced transparency (EIT) is formed by $|1\rangle$,
$|2\rangle$ and $|3\rangle$. It is driven by the two probe pulses and the
forward/backward coupling fields ($\Omega_{c+}(t)$ and $\Omega_{c-}(t)$) to
generate the stationary dark-state polariton (DSP) SLP . The coherences of
$\rho_{31+}$ and $\rho_{31-}$ correspond to the forward/backward probe
transitions, and those of $\rho_{32+}$ and $\rho_{32-}$ to the
forward/backward coupling transitions. In the system of the two-photon
transition (TPT), the two fields of $\Omega_{a}(t)$ and $\Omega_{b}(t)$ drove
the transition of $|2\rangle\rightarrow$ $|4\rangle\rightarrow$ $|5\rangle$
with a large one-photon detuning. The coherences of $\rho_{42}$ and
$\rho_{52}$ correspond to the transitions of $\Omega_{a}$ and $\Omega_{b}$,
respectively. The TPT excites and deexcites the population between $|2\rangle$
and $|5\rangle$ or, equivalently, converts the coherence of $\rho_{21}$ to
that of $\rho_{51}$ and vice versa.

The OBEs describe how the optical fields influence the density-matrix
operator, and they are given by

$\displaystyle\frac{\partial\rho_{21}}{\partial
t}=\frac{i}{2}(\Omega^{*}_{c+}\rho_{31+}+\Omega^{*}_{c-}\rho_{31-}-\Omega_{p+}\rho^{*}_{32+}-\Omega_{p-}\rho^{*}_{32-}+\Omega^{*}_{a}\rho_{41})+i\delta_{\Lambda}\rho_{21}-\gamma_{\Lambda}\rho_{21},$
(S1) $\displaystyle\frac{\partial\rho_{31\pm}}{\partial
t}=\frac{i}{2}[\Omega_{p\pm}(\rho_{11}-\rho_{33})+\Omega_{c\pm}\rho_{21}]+i\Delta_{p\pm}\rho_{31\pm}-\frac{\Gamma_{3}}{2}\rho_{31\pm},$
(S2) $\displaystyle\frac{\partial\rho_{41}}{\partial
t}=\frac{i}{2}(\Omega_{a}\rho_{21}+\Omega^{*}_{b}\rho_{51}-\Omega_{p+}\rho_{43+}-\Omega_{p-}\rho_{43-})+i(\Delta_{a}+\delta_{\Lambda})\rho_{41}-\frac{\Gamma_{4}}{2}\rho_{41},$
(S3) $\displaystyle\frac{\partial\rho_{51}}{\partial
t}=\frac{i}{2}(\Omega_{b}\rho_{41}-\Omega_{p+}\rho_{53+}-\Omega_{p-}\rho_{53-})+i(\delta_{R}+\Delta_{\rm
DDI}+\delta_{\Lambda})\rho_{51}-\left(\frac{\Gamma_{5}}{2}+\gamma_{R}+\gamma_{\rm
DDI}+\gamma_{\Lambda}\right)\rho_{51},$ (S4)
$\displaystyle\frac{\partial\rho_{32\pm}}{\partial
t}=\frac{i}{2}[\Omega_{c\pm}(\rho_{22}-\rho_{33})+\Omega_{p\pm}\rho^{*}_{21}-\Omega_{a}\rho^{*}_{43\pm}]+i\Delta_{c\pm}\rho_{32\pm}-\frac{\Gamma_{3}}{2}\rho_{32\pm},$
(S5) $\displaystyle\frac{\partial\rho_{42}}{\partial
t}=\frac{i}{2}[\Omega_{a}(\rho_{22}-\rho_{44})+\Omega^{*}_{b}\rho_{52}-\Omega_{c+}\rho_{43+}-\Omega_{c-}\rho_{43-}]+i\Delta_{a}\rho_{42}-\frac{\Gamma_{4}}{2}\rho_{42},$
(S6) $\displaystyle\frac{\partial\rho_{52}}{\partial
t}=\frac{i}{2}(\Omega_{b}\rho_{42}-\Omega_{a}\rho_{54}-\Omega_{c+}\rho_{53+}-\Omega_{c-}\rho_{53-})+i(\delta_{R}+\Delta_{\rm
DDI})\rho_{52}-\left(\frac{\Gamma_{5}}{2}+\gamma_{R}+\gamma_{\rm
DDI}\right)\rho_{52},$ (S7) $\displaystyle\frac{\partial\rho_{43\pm}}{\partial
t}=\frac{i}{2}(\Omega_{a}\rho^{*}_{32\pm}+\Omega^{*}_{b}\rho_{53\pm}-\Omega^{*}_{p\pm}\rho_{41}-\Omega^{*}_{c\pm}\rho_{42})+i(\Delta_{a}-\Delta_{c\pm})\rho_{43\pm}-\frac{\Gamma_{3}+\Gamma_{4}}{2}\rho_{43\pm},$
(S8) $\displaystyle\frac{\partial\rho_{53\pm}}{\partial
t}=\frac{i}{2}(\Omega_{b}\rho_{43\pm}-\Omega^{*}_{c\pm}\rho_{52}-\Omega^{*}_{p\pm}\rho_{51})+i(\delta_{R}+\Delta_{\rm
DDI}-\Delta_{c\pm})\rho_{53\pm}-\frac{\Gamma_{3}+\Gamma_{5}}{2}\rho_{53\pm},$
(S9) $\displaystyle\frac{\partial\rho_{54}}{\partial
t}=\frac{i}{2}[\Omega_{b}(\rho_{44}-\rho_{55})-\Omega^{*}_{a}\rho_{52}]+i(\Delta_{\rm
DDI}+\Delta_{b})\rho_{54}-\frac{\Gamma_{4}+\Gamma_{5}}{2}\rho_{54},$ (S10)
$\displaystyle\frac{\partial\rho_{11}}{\partial
t}=\frac{i}{2}(\Omega^{*}_{p+}\rho_{31+}-\Omega_{p+}\rho^{*}_{31+}+\Omega^{*}_{p-}\rho_{31-}-\Omega_{p-}\rho^{*}_{31-})+\frac{\Gamma_{3}}{2}\rho_{33},$
(S11) $\displaystyle\frac{\partial\rho_{22}}{\partial
t}=\frac{i}{2}(\Omega^{*}_{c+}\rho_{32+}-\Omega_{c+}\rho^{*}_{32+}+\Omega^{*}_{c-}\rho_{32-}-\Omega_{c-}\rho^{*}_{32-}+\Omega_{a}^{*}\rho_{42}-\Omega_{a}\rho^{*}_{42})+\frac{\Gamma}{2}\rho_{33}+\Gamma_{4}\rho_{44},$
(S12) $\displaystyle\frac{\partial\rho_{44}}{\partial
t}=\frac{i}{2}(\Omega_{a}\rho^{*}_{42}-\Omega_{a}^{*}\rho_{42}+\Omega^{*}_{b}\rho_{54}-\Omega_{b}\rho^{*}_{54})-\Gamma_{4}\rho_{44},$
(S13) $\displaystyle\frac{\partial\rho_{55}}{\partial
t}=\frac{i}{2}(\Omega_{b}\rho^{*}_{54}-\Omega^{*}_{b}\rho_{54})-\Gamma_{5}\rho_{55},$
(S14)

where $\Delta_{c\pm}$ and $\Delta_{p\pm}$ are the one-photon detunings of
$\Omega_{c\pm}$ and $\Omega_{p\pm}$, respectively, $\delta_{\Lambda}$
($=\Delta_{p+}-\Delta_{c+}=\Delta_{p-}-\Delta_{c-}$) is the two-photon
detuning of the $\Lambda$-type EIT transitions, $\gamma_{\Lambda}$ is the
ground-state decoherence rate, $\Delta_{a}$ and $\Delta_{b}$ are the one-
photon detunings of $\Omega_{a}$ and $\Omega_{b}$, respectively, $\delta_{R}$
($=\Delta_{a}+\Delta_{b}$) is the two-photon detuning of the TPT, $\gamma_{R}$
is the Rydberg-state decoherence rate, $\Delta_{\rm DDI}$ and $\gamma_{\rm
DDI}$ are the DDI-induced frequency shift and decoherence rate, and
$\Gamma_{3}$ ($=\Gamma$ or 2$\pi$$\times$6 MHz), $\Gamma_{4}$ ($=\Gamma$), and
$\Gamma_{5}$ (= 1.3$\times$$10^{-3}$$\Gamma$) represent the spontaneous decay
rates of $|3\rangle$, $|4\rangle$, and $|5\rangle$, respectively.

Since $\gamma_{\Lambda}$, $\gamma_{R}$, $\gamma_{\rm DDI}$ $\ll$ $\Gamma$ in
the experiment, these three decoherence rates do not appear in the time-
derivative equations of $\rho_{53\pm}$ and $\rho_{54}$ in the calculation.
Furthermore, $\Gamma_{5}\ll\Gamma$, $\Gamma_{5}\ll\gamma_{R}$, and
$\Gamma_{5}^{-1}\gg$ the TPT duration in the experiment. Thus, we set
$\Gamma_{5}=0$ in the calculation, which did not change the outcome. The
forward probe and coupling fields had
$\Delta_{c+}=\Delta_{p+}\equiv\Delta_{+}\approx 0$, and the backward probe and
coupling fields had $\Delta_{c-}=\Delta_{p-}\equiv\Delta_{-}=-1.0$$\Gamma$.
Such an arrangement of the one-photon detunings suppressed possible high-order
coherences, which could be established by the counter-propagating fields of
$\Omega_{p+}$ and $\Omega_{c-}$ ($\Omega_{p-}$ and $\Omega_{c+}$) SLP .
Consequently, we do not need to consider the high-order coherences in the
calculation.

In the present work, the experimental condition satisfied the weak-interaction
regime, i.e., $(r_{B}/r_{a})^{3}\ll 1$ where $r_{B}$ is the blockade radius
and $r_{a}$ is the half mean distance between Rydberg atoms. Hence,
$\Delta_{\rm DDI}$ and $\gamma_{\rm DDI}$ in the OBEs are given by OurOEMFT ;
OurCommunPhys2021 ; OurPRR2022

$\displaystyle\Delta_{\rm DDI}=A\rho_{55},$ (S15) $\displaystyle\gamma_{\rm
DDI}=A\rho_{55},$ (S16)

where $A$ represents the frequency shift (also the decoherence rate) per unit
population of $|5\rangle$, and its value depends on the atomic density in the
system and the DDI strength. The appearances of $\Delta_{\rm DDI}$ and
$\gamma_{\rm DDI}$ make the OBEs nonlinear. The OBEs involving the DDI effect
are very similar to those in our previous works of Refs. OurCommunPhys2021 ;
OurPRR2022 , which were derived based on the mean field model studied in Ref.
OurOEMFT . In those works, the theoretical predictions are in agreement with
the experimental data.

The MSEs describe how the forward and backward probe fields are affected by
the atomic coherences, and they are given by

$\displaystyle\frac{1}{c}\frac{\partial}{\partial
t}\Omega_{p+}+\frac{\partial}{\partial
z}\Omega_{p+}=i\frac{\alpha\Gamma}{2L}\rho_{31+},$ (S17)
$\displaystyle\frac{1}{c}\frac{\partial}{\partial
t}\Omega_{p-}-\frac{\partial}{\partial
z}\Omega_{p-}+i\Delta_{k}\Omega_{p-}=i\frac{\alpha\Gamma}{2L}\rho_{31-},$
(S18)

where $\alpha$ and $L$ are the optical depth and length of the atom cloud, and
$\Delta_{k}$
$\equiv(\vec{k}_{p-}-\vec{k}_{c-}-\vec{k}_{p+}+\vec{k}_{c+})$$\cdot$$\hat{z}$
and $\vec{k}_{p\pm}$ and $\vec{k}_{c\pm}$ are the wave vectors of
$\Omega_{p\pm}$ and $\Omega_{c\pm}$, respectively. The value of $L\Delta_{k}$
is the degree of phase mismatch in the system. Although $\Delta_{k}$ only
appears in the MSE of $\Omega_{p-}$, it reduces the energies of $\Omega_{p+}$
and $\Omega_{p-}$ equally due to the coupling between the forward and backward
fields, i.e., the four-wave mixing process Ldk2004 ; SLP2012 ; FBS2021 .

##  III Illustration of Figures 2(a)and 2(b) of the main text

We utilize  Figs. 2(a) and 2(b) of the main text to illustrate the procedure
and method for forming the stationary DSPs and driving the two-photon
transition (TPT) from the ground state to the Rydberg state. In the figures,
the black and red circles are the experimental data of the input and output
probe signals as functions of time, the red lines are the theoretical
predictions, the black line is the Gaussian function used as the input probe
pulse in the theoretical calculation, and the blue, green, magenta, and cyan
lines represent the timing sequences of the forward coupling (the Rabi
frequency denoted as $\Omega_{c+}$), backward coupling ($\Omega_{c-}$), and
two TPT fields ($\Omega_{a}$ and $\Omega_{b}$), respectively. The theoretical
predictions were calculated numerically with the equations presented in Sec.
II, and they are consistent with the experimental data well.

In  Fig. 2(a) of the main text, we formed the stationary DSP and did not drive
the TPT. The duration before $t\approx$ 2 $\mu$s is the Input period specified
on the figure. The probe pulse with the peak Rabi frequency $\Omega_{p0}$ was
sent into the atom cloud in the forward direction under the presence of only
$\Omega_{c+}$ ($=0.54$$\Gamma$). As the majority of the probe pulse had
entered the atom cloud and propagated slowly, we quickly switched off
$\Omega_{c+}$ as indicated by the falling edge of the blue line at $t\approx$
2 $\mu$s. At this moment, the part of the pulse centering around its peak
inside the cloud was stored. This is the end of the Input period. The black
circles after $t\approx$ 2 $\mu$s and the red circles before $t\approx$ 2
$\mu$s were the parts of the pulse outside the cloud, and they were not
stored.

After the Input period shown on  Fig. 2(a) of the main text, there was a fixed
storage time. Then, the period of Stationary DSP started. We simultaneously
switched on $\Omega_{c+}$ and $\Omega_{c-}$ (both are equal to 0.44$\Gamma$),
as indicated by the green and blue pulses, to form the stationary DSP for a
duration of about 0.7 $\mu$s. Within the stationary-DSP duration, the leakage
of the probe pulse in the forward direction is shown by the red circles
underneath the blue square pulse. There was also the leakage in the backward
direction, but we did not measure it. After keeping $\Omega_{c+}$ and
$\Omega_{c-}$ on for about 0.7 $\mu$s, we switched them off. This is the end
of the Stationary DSP period. The remaining probe energy was stored in the
atoms for another fixed time. Next, the Output period started. We switched on
$\Omega_{c+}$ (= 0.54$\Gamma$) again and retrieved the probe pulse. The red
circles after the rising edge of the blue line at $t\approx$ 6.5 $\mu$s
represent the retrieved probe pulse. Other details of forming the stationary
DSP can be found in Ref. SLP .

In  Fig. 2(b) of the main text, we formed the stationary DSP and also drove
the TPT simultaneously. Since the Input period here was completely the same as
that of  Fig. 2(a) of the main text, it is not drawn again. The TPT fields
$\Omega_{a}$ (= 2.0$\Gamma$) and $\Omega_{b}$ (= 1.6$\Gamma$) were applied to
induce the DDI during the Stationary DSP period. The field $\Omega_{a}$ is
indicated by the magenta pulse. Without the presence of $\Omega_{a}$, the
field of $\Omega_{b}$ does not affect the effects of EIT, slow light, and
storage of light at all. Hence, for simplicity it was constantly present.

Figure 2(b) of the main text shows the representative data of the stationary
DSP dressed by the DDI, and the magenta square pulse corresponds to the TPT
duration. The TPT fields produced the Rabi oscillation between $|2\rangle$ and
$|5\rangle$. Consequently, in the figure the red circles of the probe leakage
within the duration exhibit the oscillation. The node (anti-node) of the
oscillation corresponds to all the population in $|5\rangle$ (in $|2\rangle$).
We set the TPT frequency to the two-photon resonance by minimizing the anti-
node. Due to the DDI-induced attenuation, the retrieved probe energy in  Fig.
2(b) of the main text is obviously less than that in  Fig. 2(a) of the main
text.

##  IV Determination of the Rabi Frequencies

In this section, we illustrate the methods that determine the Rabi frequencies
of the forward and backward coupling fields $\Omega_{c+}$ and $\Omega_{c-}$,
and that of the two-photon-transition (TPT) field $\Omega_{a}$. We also
checked whether the Rabi frequency of the TPT field $\Omega_{b}$ determined by
the data in  Fig. 2(b) of the main text is reasonable.

The Autler-Townes splitting (ATS) in the EIT transmission spectrum was used to
determine the Rabi frequencies of $\Omega_{c+}$ and $\Omega_{c-}$. Figure S1
shows the representative EIT spectra. We followed the determination method
described in Ref. OurPRA2019 . The frequency difference between two
transmission minima, $\delta_{+}$ and $\delta_{-}$, gave the Rabi frequency.
As mentioned in the reference, the optical depth (OD) of the atom cloud was
reduced to a value between 1 and 2 in the measurement of the spectrum. At a
low OD, the minima can be clearly identified, and any asymmetry in the
spectrum due to a nonzero one-photon detuning of the coupling field can be
noticed. In the EIT system,

$\delta_{+}-\delta_{-}\approx\Omega+\Delta^{2}/2\Omega,$ (S19)

where $\Omega$ is the coupling Rabi frequency and $\Delta$ is the one-photon
detuning. As long as $\Delta^{2}\ll\Omega^{2}$, the angular-frequency
separation between $\delta_{+}$ and $\delta_{-}$ is nearly the same as
$\Omega$.

Figure S1: (a) [or (b)] Probe transmission spectrum as a function of the two-
photon detuning, $\delta$, under the constant presence of the forward (or
backward) coupling field of $\Omega_{c+}$ (or $\Omega_{c-}$) in the
$\Lambda$-type EIT system, and (c) that under the constant presence the TPT
field of $\Omega_{b}$ in the Rydberg-EIT system. In the measurement of each
spectrum, the one-photon detuning was nearly zero. Black circles are the
experimental data. Red curves are the best fits to determine the frequencies
of minima. The frequency separation between the two minima gives the Autler-
Townes splitting of (a) 0.44$\Gamma$, (b) 0.44$\Gamma$, or (c) 2.0$\Gamma$.

In the $\Lambda$-type EIT system, we employed the forward probe field
($\Omega_{p+}$) to measure the values of the forward and backward coupling
Rabi frequencies, $\Omega_{c+}$ and $\Omega_{c-}$, as shown in Figs. S1(a) and
S1(b), respectively. The sweeping of the probe frequency deviated the probe
propagation direction very little. Since the propagation direction of
$\Omega_{c-}$ was opposite to that of $\Omega_{p+}$, the Doppler effect made
the decoherence rate in Fig. S1(b) much larger than that in Fig. S1(a).
Consequently, the EIT peak transmission in Fig. S1(b) is significantly
reduced, and that in Fig. S1(a) is nearly 100%. Nevertheless, the frequency
separation between the two transmission minima is not affected by the
decoherence rate.

Figure S2: Main plot: In the system of all-optical switching (AOS),
transmission of the retrieved probe energy versus the pulse duration of the
field of $\Omega_{a}$. Circles are the experimental data taken at
$\Delta_{a}=$ 5.0$\Gamma$. Black line is the best fit with the fitting
function given by Eq. (S20), which determines $\Omega_{a}=$ 2.0$\Gamma$.
Inset: Black and red lines are the signals of the probe pulses without and
with the atoms, respectively. The black line is scaled down by a factor of
0.5. Blue and green dashed lines represent the time sequences of the coupling
field and the field of $\Omega_{a}$. Filled area in red represents the
retrieved probe energy after the AOS effect.

To determine the value of $\Omega_{a}$ of the TPT field, we utilized the all-
optical switching (AOS) effect in the four-level $N$-type system AOS ; XPM ,
which consists of states $|1\rangle$, $|2\rangle$, $|3\rangle$, and
$|4\rangle$ in this work. In the AOS, the forward probe ($\Omega_{p+}$) and
coupling ($\Omega_{c+}$) fields form the $\Lambda$-type EIT configuration, and
the presence of the field of $\Omega_{a}$ causes an additional attenuation of
the output probe field. As shown in the inset of Fig. S2, we first stored the
probe pulse in the atoms, then applied the pulse of $\Omega_{a}$ during the
storage time, and finally retrieved the probe pulse. The one-photon detuning,
$\Delta_{a}$, of the $\Omega_{a}$ pulse was kept to $+5.0$$\Gamma$ in the
measurement. Except that the fields of $\Omega_{c-}$ and $\Omega_{b}$ were
absent, all the other experimental arrangements were the same as those in the
measurement with the TPT. In Fig. S2, the circles and black line are the
experimental data and best fit of the retrieved probe energy transmission,
$T_{\rm AOS}$, as a function of the $\Omega_{a}$ pulse duration, $\tau_{a}$.
The fitting function is given by

$T_{\rm
AOS}=T_{0}\exp\left(-\frac{\Omega_{a}^{2}\Gamma}{\Gamma^{2}+4\Delta_{a}^{2}}\tau_{a}\right),$
(S20)

where $T_{0}$ is the transmission without the presence of $\Omega_{a}$. The
best fit determined the value of $\Omega_{a}$ used in the TPT.

We measured the ATS induced by the TPT field $\Omega_{b}$ in a Rydberg-EIT
system OurPRA2019 . In the system, the coupling field $\Omega_{b}$ and a weak
probe field drove the transition of $|5P_{3/2},F=3,m_{F}=3\rangle$
$\rightarrow$ $|32D_{5/2},m_{J}=5/2\rangle$ and that of
$|5S_{1/2},F=2,m_{F}=2\rangle$ $\rightarrow$ $|5P_{3/2},F=3,m_{F}$=$3\rangle$,
respectively, and the two fields counter-propagated. The input probe Rabi
frequency was set to about 0.05$\Gamma$ such that the DDI effect was
negligible. Figure S1(c) shows the Rydberg-EIT spectrum taken at a nearly zero
one-photon detuning. The ATS in the spectrum is 2.0$\Gamma$. However, the
transition of $\Omega_{b}$ in Fig. S1(c) was not the same as that of
$\Omega_{b}$ in the TPT, which was from $|5P_{3/2},F=3,m_{F}=2\rangle$ to
$|32D_{5/2},m_{J}=3/2~{}\&~{}5/2\rangle$. According to the ATS and the
Clebsch-Gordan coefficients, the Rabi frequency of the transition to
$|32D_{5/2},m_{J}=5/2\rangle$ and that to $|32D_{5/2},m_{J}=3/2\rangle$ should
be 1.4$\Gamma$ and 1.1$\Gamma$, respectively. Their root-mean-square value is
1.8$\Gamma$. It is reasonable that the TPT field $\Omega_{b}$ had the Rabi
frequency of 1.6$\Gamma$ as determined by  Fig. 2(b) of the main text.

##  V Degree of Phase Mismatch During the Stationary Dark-State Polariton

The stationary dark-state polariton (DSP) was formed by the balance between
the two four-wave mixing (FWM) processes of
$\Omega_{p+}-\Omega_{c+}+\Omega_{c-}\rightarrow\Omega_{p-}$ and
$\Omega_{p-}-\Omega_{c-}+\Omega_{c+}\rightarrow\Omega_{p+}$ SLP2012 , where
$\Omega_{p\pm}$ and $\Omega_{c\pm}$ denote the forward/backward probe and
coupling fields, respectively. The former process consists of absorbing a
forward probe photon, emitting a forward coupling photon, and absorbing a
backward coupling photon to finally generate a backward probe photon. The
latter process is just the reversed former process to generate a forward probe
photon from a backward probe photon. In the FWM, the phase mismatch causes an
additional attenuation of the stationary DSP, and its measure, i.e., the
degree of the phase mismatch, is given by $L\Delta_{k}$, where $L$ is the
medium length and $\Delta_{k}$ represents the sum of the wave vectors of the
four fields projected on the propagation direction. Please see Eqs. (S17) and
(S18) for the role of $L\Delta_{k}$ in the Maxwell-Schrödinger equations. The
definition and details of $\Delta_{k}$ are given right after the equations.

To determine the value of $L\Delta_{k}$ in the stationary-DSP experiment, we
measured the retrieved probe energy as a function of the stationary-DSP
duration. The experimental condition or parameters of the measurement was the
same as that of  Fig. 2(a) of the main text. As shown in the inset of Fig. S3,
we first stored the probe pulse in the atoms, then applied the pulse of the
forward and backward coupling fields to generate the stationary DSP for a
given duration, $\tau_{\rm DSP}$, and finally retrieved the probe pulse after
the execution of $\tau_{\rm DSP}$. The time length from storing to retrieving
the probe pulse was fixed to about 10 $\mu$s to accommodate the variation of
the stationary-DSP duration. The main plot of Fig. S3 shows the retrieved
probe energy transmission versus $\tau_{\rm DSP}$ at three different values of
the peak Rabi frequency, $\Omega_{p0}$, of the input probe pulse. In the
figure, the circles are the experimental data. The retrieved probe energy
transmissions of $\Omega_{p0}=$ 0.07$\Gamma$ and 0.24$\Gamma$ are about the
same, but is higher than that of $\Omega_{p0}=$ 0.40$\Gamma$. As
$\Omega_{p0}\leq$ 0.24$\Gamma$, the probe pulse is weak enough to be treated
as the perturbation, and the result does not depend on $\Omega_{p0}$. On the
other hand, as $\Omega_{p0}$ becomes larger, the probe pulse is no longer the
perturbation and its intensity influences the outcome.

Figure S3: Main plot: Transmission of the retrieved probe energy versus the
stationary-DSP duration. Black, red, and blue circles are the experimental
data taken at $\Omega_{p0}=$ 0.07$\Gamma$, 0.24$\Gamma$, and 0.40$\Gamma$,
respectively, where $\Omega_{p0}$ is the peak Rabi frequency of the input
probe pulse. Gray ($\Omega_{p0}=$ 0.07$\Gamma$ and 0.24$\Gamma$) and blue
($\Omega_{p0}=$ 0.40$\Gamma$) areas are the theoretical predictions calculated
with $L\Delta_{k}$ = 0.90 rad. The experimental parameters were $\alpha$ =
36$\pm$1, $\Omega_{c+}=$ 0.54$\Gamma$ before the storage and after the
retrieval, $\Omega_{c+}=\Omega_{c-}=$ 0.44$\Gamma$ during the stationary DSP,
and $\gamma_{\Lambda}$ = (9$\pm$1)$\times$$10^{-4}$$\Gamma$. Thicknesses of
the gray and blue areas represent the uncertainties or fluctuations in the OD
and ground-state decoherence rate. Inset: Representative data of the probe
signals without (black circles) and with (red circles) the atom cloud. Black
(red) line is the best fit (theoretical prediction) of the black (red)
circles. The black circles and line are scaled down by a factor of 0.5. Blue
and green lines represent the timing sequence of $\Omega_{c+}$ and
$\Omega_{c-}$, the two lines completely overlap during the stationary DSP (4
$\mu$s $<t<$ 9 $\mu$s here). The red circles underneath the green-line pulse
reveal that during the stationary DSP the probe light leaked out of the atoms
in the forward direction.

In Fig. S3, the colored areas represent the theoretical predictions. We used
the experimentally-determined values of the optical density ($\alpha$), the
forward and backward coupling Rabi frequencies ($\Omega_{c+}$ and
$\Omega_{c-}$), and the ground-state decoherence rate ($\gamma_{\Lambda}$) in
the theoretical calculation. The determination methods of $\alpha$ and
$\gamma_{\Lambda}$ can be found in the last paragraph of  Sec. I, and that of
$\Omega_{c+}$ and $\Omega_{c-}$ in the second and third paragraphs of  Sec.
IV. Due to the fluctuations in $\alpha$ and $\gamma_{\Lambda}$, the prediction
of each $\Omega_{p0}$ is displayed as an area. A single value of $L\Delta_{k}$
of 0.9 rad was used in all the predictions such that the consistency between
the data and predictions is satisfactory as shown in Fig. S3. This nonzero
$L\Delta_{k}$ can be caused by a little misalignment among the fields in the
FWM processes.

##  VI Estimation of the Coefficient A of Dipole-Dipole Interaction

We estimate the strength of dipole-dipole interaction (DDI) or the coefficient
$A$ with a simplified model in this section, where $A$ is defined as the DDI-
induced decoherence rate, $\gamma_{\rm DDI}$, and frequency shift,
$\Delta_{\rm DDI}$, per unit Rydberg-state population, $\rho_{55}$, as shown
by Eqs. (S15) and (S16). The coefficient $A$ was used in the equations of the
coherences $\rho_{51}$, $\rho_{52}$, $\rho_{53}$, and $\rho_{54}$ relating to
the Rydberg state $|5\rangle$, i.e., Eqs. (S4), (S7), (S9), and (S10),
respectively. In the experiment, we utilized a series of measurements of the
output probe energy transmission against the input probe intensities or
$\Omega_{p0}^{2}$ to determine $A=$ 0.60$\Gamma$ as shown by  Figs. 2(a) and
2(b) of the main text.

The formation of stationary dark-state polaritons (DSP) is not relevant to the
DDI strength. Consequently, state $|3\rangle$ and the probe and coupling
fields are not present in the model here. To focus on the DDI effect, we only
consider states $|1\rangle$, $|2\rangle$, and $|5\rangle$. An effective field
with the Rabi frequency of $\Omega_{\rm eff}$ drives the one-photon transition
(OPT) between $|2\rangle$ and $|5\rangle$. Due to a large one-photon detuning,
$\Delta_{a}$, in the experiment, the intermediate state $|4\rangle$ of the
two-photon transition (TPT) can be eliminated. Thus, the TPT of
$|2\rangle\leftrightarrow|4\rangle\leftrightarrow|5\rangle$ with the driving
fields of $\Omega_{a}$ and $\Omega_{b}$ in the experiment is equivalent to the
OPT of $|2\rangle\leftrightarrow|5\rangle$ here, and

$\Omega_{\rm eff}\approx\frac{\Omega_{a}\Omega_{b}}{2\Delta_{a}}.$ (S21)

Considering $\Omega_{a}=$ 2.0$\Gamma$, $\Omega_{b}=$ 1.6$\Gamma$, and
$\Delta_{a}=$ 5.0$\Gamma$ of the TPT in the experiment, $\Omega_{\rm eff}=$
0.32$\Gamma$ in the OPT.

The OPT drives the Rabi oscillation between $|2\rangle$ and $|5\rangle$, and
the equations of motion are given by

$\displaystyle\frac{\partial\rho_{55}}{\partial t}=\frac{i}{2}(\Omega_{\rm
eff}\rho^{*}_{52}-\Omega^{*}_{\rm eff}\rho_{52}),$ (S22)
$\displaystyle\frac{\partial\rho_{52}}{\partial t}=\frac{i}{2}\Omega_{\rm
eff}(\rho_{22}-\rho_{55})+iA\rho_{55}\rho_{52}-A\rho_{55}\rho_{52},$ (S23)
$\displaystyle\rho_{0}=\rho_{22}+\rho_{55},$ (S24)

where $A\rho_{55}$ is substituted for $\gamma_{\rm DDI}$ and $\Delta_{\rm
DDI}$ and $\rho_{0}$ is the initial total population in the system. In the
above equations, we neglect the spontaneous decay rate of $|5\rangle$ (i.e.,
$\Gamma_{5}$) and all the other decoherence rates (i.e., $\gamma_{\Lambda}$
and $\gamma_{R}$), which appear in Eqs. (S1)-(S14).

The OPT also converts the coherences $\rho_{21}$ and $\rho_{51}$ back and
forth, and the equations of the conversion are given by

$\displaystyle\frac{\partial\rho_{21}}{\partial t}=\frac{i}{2}\Omega_{\rm
eff}\rho_{51},$ (S25) $\displaystyle\frac{\partial\rho_{51}}{\partial
t}=\frac{i}{2}\Omega_{\rm eff}\rho_{21}+i\omega\rho_{51}-\gamma\rho_{51},$
(S26)

where $\omega$ and $\gamma$ correspond to the detuning and decoherence rate of
the OPT. Since the output probe energy transmission is proportional to
$|\rho_{21}|^{2}$, we are interested in the solution of $\rho_{21}$ given by

$\displaystyle{\rm Re}[\rho_{21}(\omega,\gamma)]$ $\displaystyle\approx$
$\displaystyle\sqrt{\rho_{0}}e^{-\gamma t/2}\left[\cos\left(\frac{\omega
t}{2}\right)\cos\left(\frac{\Omega^{\prime}t}{2}\right)+\frac{\omega}{\Omega^{\prime}}\sin\left(\frac{\omega
t}{2}\right)\sin\left(\frac{\Omega^{\prime}t}{2}\right)+\frac{\gamma}{\Omega^{\prime}}\cos\left(\frac{\omega
t}{2}\right)\sin\left(\frac{\Omega^{\prime}t}{2}\right)\right],$ (S27)
$\displaystyle{\rm Im}[\rho_{21}(\omega,\gamma)]$ $\displaystyle\approx$
$\displaystyle\sqrt{\rho_{0}}e^{-\gamma t/2}\left[\sin\left(\frac{\omega
t}{2}\right)\cos\left(\frac{\Omega^{\prime}t}{2}\right)-\frac{\omega}{\Omega^{\prime}}\cos\left(\frac{\omega
t}{2}\right)\sin\left(\frac{\Omega^{\prime}t}{2}\right)+\frac{\gamma}{\Omega^{\prime}}\sin\left(\frac{\omega
t}{2}\right)\sin\left(\frac{\Omega^{\prime}t}{2}\right)\right],$ (S28)
$\displaystyle|\rho_{21}(\omega,\gamma)|^{2}$ $\displaystyle\approx$
$\displaystyle\rho_{0}e^{-\gamma
t}\left\\{\frac{1}{2}[1+\cos(\Omega^{\prime}t)]+\frac{\gamma}{\Omega^{\prime}}\sin(\Omega^{\prime}t)+\frac{\gamma^{2}}{2\Omega^{\prime
2}}[1-\cos(\Omega^{\prime}t)]+\frac{\omega^{2}}{2\Omega^{\prime
2}}[1-\cos(\Omega^{\prime}t)]\right\\},$ (S29) $\displaystyle\Omega^{\prime}$
$\displaystyle\equiv$ $\displaystyle\sqrt{\Omega_{\rm eff}^{2}+\omega^{2}}.$
(S30)

To estimate the value of $A$, we utilized the mean-field approach based on the
probability function of nearest-neighbor distribution (NND) OurOEMFT ;
NNDistribution . Considering Rydberg atoms are randomly distributed, the
probability of finding the DDI-induced frequency shift within $\omega$ and
$\omega+d\omega$ is given by

$\displaystyle P(\omega)\\!\\!\\!$
$\displaystyle=\frac{1}{\omega_{a}}\frac{\left[1+\sqrt{1+4(\omega/\omega_{a})}\right]^{2}}{4(\omega/\omega_{a})^{2}\sqrt{1+4(\omega/\omega_{a})}}$
(S31)
$\displaystyle\times\exp{\left[-\frac{1+\sqrt{1+4(\omega/\omega_{a})}}{2(\omega/\omega_{a})}\right]},$

where

$\omega_{a}\equiv|C_{6}|/r_{a}^{6},$ (S32)

$C_{6}$ is the van der Waals coefficient, $r_{a}$ is the half-mean distance
between the Rydberg atoms given by

$(4\pi/3)r_{a}^{3}\equiv(n_{\rm atom}\rho_{55})^{-1},$ (S33)

and $n_{\rm atom}$ is the atomic density. The derivation of $P(\omega)$ using
the NND can be found in Ref. OurOEMFT . We performed the following integrals
to simulate the DDI effect:

$\displaystyle\bar{\rho}_{21,r}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\omega P(\omega)~{}{\rm
Re}[\rho_{21}(\omega,\gamma=0)],$ (S34) $\displaystyle\bar{\rho}_{21,i}(t)$
$\displaystyle=$ $\displaystyle\int_{0}^{\infty}d\omega P(\omega)~{}{\rm
Im}[\rho_{21}(\omega,\gamma=0)],$ (S35) $\displaystyle|\bar{\rho}_{21}|^{2}$
$\displaystyle\equiv$
$\displaystyle[\bar{\rho}_{21,r}(t)]^{2}+[\bar{\rho}_{21,i}(t)]^{2}.$ (S36)

Because of the distribution function $P(\omega)$, $|\bar{\rho}_{21}|^{2}$ is a
function of $\rho_{55}$ as indicated by Eqs. (S31)-(S33).

Figure S4: Estimation of the DDI strength or coefficient $A$. (a) The solution
of $\rho_{55}(t)$ by numerically solving Eqs. (S22)-(S24). (b) The ensemble
average of the microscopic DDI effect, i.e., $|\bar{\rho}_{21}|^{2}$
calculated from Eqs. (S34)-(S36). Black circles, red circles, blue squares,
and green diamonds represent the results from the first to the fourth
iterations. Lines represent the corresponding best fits with the fitting
function of $|\rho_{21}(\omega=A\rho_{55},\gamma=A\rho_{55})|^{2}$ given by
Eq. (S29), in which $A$ is the fitting parameter.

The process of the estimation of $A$ is illustrated as the followings: (i) We
set $A=0$ and numerically solved Eqs. (S22)-(S24) to obtain $\rho_{55}(t)$ as
shown by the black circles in Fig. S4(a). (ii) With $\rho_{55}(t)$, we used
Eqs. (S34)-(S36) to evaluate the temporal data points of
$|\bar{\rho}_{21}|^{2}$ represented by the black circles in Fig. S4(b). (iii)
The data points in Fig. S4(b) are fitted with the function of
$|\rho_{21}(\omega=A\rho_{55},\gamma=A\rho_{55})|^{2}$ shown by Eq. (S29),
where $A$ is the only fitting parameter. The fitting means to relate the
ensemble average of the microscopic DDI effect to the macroscopic quantity
$A$, i.e.,

$|\bar{\rho}_{21}|^{2}=|\rho_{21}(\omega=A\rho_{55},\gamma=A\rho_{55})|^{2}.$
(S37)

The black line in Fig. S4(b) is the best fit, which determined the value of
$A$ in the first iteration, denoted as $A^{(1)}$. (iv) In the second
iteration, we started with the value of $A^{(1)}$ in Eqs. (S22)-(S24) and
solved the equations to obtain $\rho_{55}^{(2)}(t)$ as shown by the red
circles in Fig. S4(a). Steps (ii) and (iii) were repeated, and the red circles
and line in Fig. S4(b) are the results, giving $A^{(2)}$. (v) The process of
(i) $\rightarrow$ (ii) $\rightarrow$ (iii) was performed several times until
the results converged as demonstrated by the blue squares/line and green
diamonds/line in Figs. S4(a) and S4(b).

In the estimation, we set $C_{6}=$ –2$\pi$$\times$130 MHz$\cdot$$\mu$m6,
$n_{\rm atom}=0.05$ $\mu$m-3, and $\Omega_{\rm eff}=$ 0.32$\Gamma$ according
to the experimental condition. Considering the population of $|2\rangle$
achievable in this work, the value of $A$ estimated by the method in this
section depends on $\rho_{0}$ very little. As the estimation gives $A$ of
0.54$\Gamma$, the experimentally determined $A$ of 0.60$\Gamma$ sounds
reasonable.

##  VII The Input Probe Pulse with a Large Rabi Frequency

We employed the input probe pulses with five different peak Rabi frequencies,
$\Omega_{p0}$, of 0.07$\Gamma$, 0.2$\Gamma$, 0.24$\Gamma$, 0.3$\Gamma$,
0.4$\Gamma$ in the experiment. As $\Omega_{p0}\leq$ 0.24$\Gamma$, the
experimental data of slow light were well consistent with the theoretical
predictions. As $\Omega_{p0}\geq$ 0.3$\Gamma$, the experimental data of slow
light deviated from the theoretical predictions. A representative example of
the discrepancy taken at $\Omega_{p0}$ of 0.4$\Gamma$ is shown in Fig. S5(a).
In the theoretical calculations, we used the parameters of the coupling Rabi
frequency ($\Omega_{c+}$), optical depth or OD ($\alpha$), and ground-state
decoherence rate ($\gamma_{\Lambda}$), which were determined experimentally by
the methods described in  Sec. IV and in the last paragraph of  Sec. I. The
determinations were carried out at the experimental condition that the typical
electromagnetically induced transparency (EIT) criterion, i.e.,
$\Omega_{p0}^{2}\ll\Omega_{c+}^{2}$ and $2\gamma\Gamma\ll\Omega_{c+}^{2}$, is
satisfied.

Figure S5: (a) Slow light under the constant presence of the coupling field.
(b, c) The formation of the stationary DSP without and with the TPT or
equivalently DDI, respectively. $\Omega_{p0}=$ 0.4$\Gamma$, and all the other
experimental parameters are the same as those shown in the caption of  Figs.
2(a) and 2(b) of the main text. Black and red circles are the experimental
data of the input and output probe pulses. Black lines represent the Gaussian
best fits and red lines are the theoretical predictions. Blue, green, magenta,
and cyan lines show the timing sequences of the coupling fields of
$\Omega_{c+}$ and $\Omega_{c-}$ as well as the TPT fields of $\Omega_{a}$ and
$\Omega_{b}$. The part of the probe pulse between the two vertical gray dashed
lines in (a) was stored in the atoms at $t\approx$ 6.0 $\mu$s (when
$\Omega_{c+}$ was switched off) in (b) and (c).

The comparison between the experimental data and theoretical predictions in
Fig. S5(a) leads to the following observations: (i) The experimental output
pulse peak is smaller than the theoretical one, (ii) the experimental output
pulse width is broader than the theoretical one, and (iii) the experimental
propagation delay time is approximately the same as the theoretical one. These
observations can explain the results of the output probe energies without and
with the two-photon transition (TPT).

Figure S5(b) shows the input and output probe pulses without the presence of
the TPT. The timing sequences of the forward ($\Omega_{c+}$) and backward
($\Omega_{c-}$) coupling fields are nearly the same as those shown in  Fig.
2(a) of the main text. As we switched off $\Omega_{c+}$ at $t\approx$ 2.0
$\mu$s, the stored probe energy of the experimental data was less than that of
the theoretical prediction due to the observation (i). This first factor could
make the experimental probe energy less than the predicted one at the output
stage, i.e., $t\gtrsim$ 6.6 $\mu$s. Nevertheless, the stored probe pulse shape
of the experimental data had a smoother spatial variation than that of the
theoretical prediction due to the observation (ii). A smoother pulse shape
indicates a narrower frequency bandwidth, which makes the energy loss less
during the formation of stationary dark-state polaritons (DSPs) SLP ; SLPCJP ,
i.e., the duration of the green plus and blue pulse in Fig. S5(b). A narrower
frequency bandwidth also makes the energy loss less during the pulse
propagation at the output stage OurPRL2013 . This second factor could make the
experimental output probe energy more than the predicted. Finally, the medium
of EIT plus stationary DSP can be viewed as a frequency plus spatial bandpass
filter. This third factor can make the output probe pulse shape is
approximately independent of the input one. Putting the three factors
together, we can expect that the experimental and theoretical output probe
pulses are consistent as demonstrated by Figs. S5(b) and S5(c), and also by
the rightmost blue- and red-circle data points and their corresponding
predictions in  Fig. 2(c) of the main text.

When the probe Rabi frequency is comparable to the coupling Rabi frequency as
well as the spontaneous decay rate in the EIT system, the EIT effect is no
longer independent of the probe Rabi frequency. We think a three-dimensional
calculation of the optical Bloch equations and the Maxwell-Schrödinger
equation is necessary to make a precise theoretical prediction of the slow
light. Nevertheless, the three-dimensional calculation is beyond the scope of
this work. The discrepancy between the experimental data and theoretical
prediction at a large input probe Rabi frequency can be explained, and it does
not affect the conclusion of this work at all.

##  VIII The DDI-Induced Phase Shift under the Rabi Oscillation

When we applied the two-photon transition (TPT) in the experiment, it not only
produced the Rabi oscillation between the ground state $|2\rangle$ and the
Rydberg state $|5\rangle$ but also induced the dipole-dipole interaction
(DDI). The Rabi oscillation alone resulted in a phase shift of the ground-
state coherence $\rho_{21}$, and the DDI added an extra phase shift to
$\rho_{21}$. Note that the phase shift of $\rho_{21}$ was transferred to that
of the output probe pulse after the retrieval XPM .

Figure S6: Simulations of the phase shift of the ground-state coherence
$\rho_{21}$, or equivalently the phase shift of the output probe pulse, as a
function of time. Black and red lines are the results without and with the DDI
of $A=$ 0.60$\Gamma$, respectively. In the simulations, $\Omega_{p0}=$
0.24$\Gamma$, $\Delta_{-}=0$, $\Omega_{a}=\Omega_{b}=$ 2.0$\Gamma$, and the
remaining parameters are the same as those in the caption of  Figs. 2(a) and
2(b) of the main text. All decay and decoherence mechanisms are switched off
in (a) and intact in (b). Dashed lines represent the baselines of the
oscillations. In each subfigure, the difference between the two dashed lines
reveals the DDI-induced phase shift.

Figure S6(a) shows the theoretical predictions, elucidating the phase shift of
$\rho_{21}$ induced by the Rabi oscillation and also that induced by the DDI.
To enhance the DDI-induced phase shift, the decay rates $\gamma_{\Lambda}$,
$\gamma_{R}$, $\gamma_{\rm DDI}$, and $\Gamma_{4}$ as well as the degree of
phase mismatch $L\Delta_{k}$ were set to zero in the calculation of the
predictions. We varied the value of $A$ to switch off or on the DDI-induced
phase shift, i.e., DDI-induced frequency shift or $\Delta_{\rm
DDI}=A\rho_{55}$. In Fig. S6(a), the black solid line represents the result of
$A=0$, i.e., the DDI is switched off, and it reveals the oscillation of the
phase shift, which is the consequence of the Rabi oscillation. The black
dashed is the baseline of the oscillation. At $A=$ 0.60$\Gamma$, the red solid
line represents the result of the phase shift produced by the DDI together
with the Rabi oscillation. The red dashed is the baseline of the oscillation.
Thus, the phase difference between the red and black dashed lines is solely
caused by the DDI. As the interaction time gets longer, the DDI-induced phase
shift becomes larger as expected.

We next consider a more realistic case by restoring all the decay terms in the
calculation. The values of the decay rates $\gamma_{\Lambda}$, $\gamma_{R}$,
$\Gamma_{4}$, and $L\Delta_{k}$ were set to those specified in the caption of
Figs. 2(a) and 2(b) of the main text, and the DDI-induced attenuation was
brought back, i.e., $\gamma_{\rm DDI}=A\rho_{55}$. In Fig. S6(b), the black
and red solid lines are the results without and with the DDI of $A=$
0.60$\Gamma$. The difference between the black and red dashed lines
corresponds to the DDI-induced phase shift. Due to the decay terms, the amount
of the DDI-induced phase shift is reduced.

In the simulations for the results shown in Figs. S6(a) and S6(b), we made the
two TPT fields have the same Rabi frequency, i.e., $\Omega_{a}=\Omega_{b}=$
2.0$\Gamma$, and set the backward coupling field frequency to the resonance of
its transition, i.e., $\Delta_{-}=0$. The remaining calculation parameters are
the same as the experimental parameters shown in the caption of  Figs. 2(a)
and 2(b) of the main text. In the experiment, $\Omega_{a}$ (= 2.0$\Gamma$)
$\neq\Omega_{b}$ (= 1.6$\Gamma$) and $\Delta_{-}=-1$$\Gamma$. Both
$\Omega_{a}\neq\Omega_{b}$ and $\Delta_{-}=-1$$\Gamma$ can induce AC Stark
shifts, which produces a phase shift to $\rho_{21}$. Nevertheless, the phase
shifts produced by the causes of the Rabi oscillation,
$\Omega_{a}\neq\Omega_{b}$, and $\Delta_{-}=-1$$\Gamma$ depend on the input
probe Rabi frequency a little. In  Fig. 3(b) of the main text, the data points
are the results that the measured phase shifts at different values of the
input probe Rabi frequency, $\Omega_{p0}$, were subtracted by that at the
smallest value of $\Omega_{p0}$. Therefore, the phase shifts revealed by these
data points are mainly contributed from the DDI.

##  IX Estimations of the BEC transition temperature, stationary-DSP
temperature, and elastic collision rate

In Ref. Fleischhauer2008 , Fleischhauer et al. proposed to realize the Bose-
Einstein condensation (BEC) with stationary dark-state polaritons (DSP). We
employ the formulas given by Ref. Fleischhauer2008  to estimate the BEC
transition temperature and the stationary-DSP temperature, and use the hard-
sphere model to evaluate the elastic collision rate in this section. The
estimations are based on the experimental condition in the present work. The
BEC transition temperature, $T_{c}$, is given by

$T_{c}=\frac{2\pi\hbar^{2}}{[\zeta(3/2)]^{2/3}k_{B}}\frac{n_{\rm
DSP}^{2/3}}{(m_{\perp}^{2}m_{\parallel})^{1/3}}$ (S38)

with

$\displaystyle m_{\perp}$ $\displaystyle=$ $\displaystyle\frac{\hbar
k_{p}\Gamma(\alpha/L)}{\Omega_{c}^{2}},$ (S39) $\displaystyle m_{\parallel}$
$\displaystyle=$
$\displaystyle\frac{\hbar\Gamma^{2}(\alpha/L)^{2}}{8\Omega_{c}^{2}|\Delta|},$
(S40)

where $\zeta(3/2)$ is the Riemann zeta function evaluated at 3/2, $\hbar$ and
$k_{B}$ are the Planck and Boltzmann constants, $n_{\rm DSP}$ is the density
of DSPs, $m_{\perp}$ and $m_{\parallel}$ represent the effective masses in the
transverse and longitudinal directions, $k_{p}$ is wave vector of the probe
field, $\alpha$ and $L$ are the optical depth and length of the medium,
$\Omega_{c}$ and $\Delta$ are the total (forward plus backward) coupling Rabi
frequency and one-photon detuning for the creation of stationary DSPs, and
$\Gamma$ is the spontaneous decay rate of the excited state $|3\rangle$. In
the experiment, $\alpha=$ 36, $L=$ 6 mm, $\Omega_{c}=$
$\sqrt{2}$$\times$0.44$\Gamma$, and $|\Delta|$ = 1.0$\Gamma$. We obtain
$m_{\perp}/m_{\rm Rb}=2.4\times 10^{-6}$ and $m_{\parallel}/m_{\rm
Rb}=2.3\times 10^{-10}$, where $m_{\rm Rb}$ is the mass of the 87Rb atom.
Based on the numerical simulation for the data shown in  Fig. 2(c) of the main
text, the peak value of $n_{\rm DSP}(=n_{\rm atom}\rho_{22})$ is
3.7$\times$109 cm-3 of $\Omega_{p0}=$ 0.4$\Gamma$, where $n_{\rm atom}=$
4.1$\times 10^{10}$ cm-3 is the atom density and $\rho_{22}=$ 0.09 is the
maximum population in $|2\rangle$ near the end of the DDI during the
stationary DSP. Thus, we estimate $T_{c}=$ 4.0 mK.

We utilize the following formula to estimate the effective temperature,
$T_{p}$, of the stationary DSP Fleischhauer2008

$T_{p}=\frac{\hbar^{2}}{k_{B}}\left(\frac{k_{\perp}^{2}}{2m_{\perp}}+\frac{k_{\parallel}^{2}}{4m_{\parallel}}\right),$
(S41)

where $k_{\perp}$ and $k_{\parallel}$ represent the root-mean-square values of
the DSP momentums in the transverse and longitudinal directions. According to
the $e^{-1}$ full width of the input probe beam of 0.13 mm, we obtained
$k_{\perp}=$ 31 mm-1. Based on the numerical simulation for the data shown in
Fig. 2(c) of the main text, the full width of the stationary-DSP distribution
was about 6 mm, giving $k_{\parallel}=$ 0.67 mm-1. Hence, $T_{p}=$ 3.8 $\mu$K,
which has already been well below $T_{c}$.

Following the idea that the phase shift is a consequence of the elastic
collision Fleischhauer2008 , we employ the formula of $\Delta\phi/\Delta
t=R_{c}\phi_{c}$ to estimate the elastic collision rate, $R_{c}$, where
$\Delta\phi/\Delta t$ is the average phase shift per DDI interaction time,
i.e., the TPT pulse duration, and $\phi_{c}$ is the phase shift per collision.
Under the hard-sphere model, $\phi_{c}=\bar{k}r_{B}$, where $\hbar\bar{k}$ is
the root-mean-square value of the relative momentum between the colliding
spheres, and $r_{B}$ is the radius of the spheres. We assume $r_{B}=$ 2.0
$\mu$m, which is the blockade radius of $32D_{5/2}$ ($C_{6}=2\pi$$\times$ 130
MHz$\cdot$$\mu$m-6) under the effective TPT Rabi frequency of 0.32$\Gamma$.
According to the data in  Fig. 3(b) of the main text, $\Delta\phi/\Delta t=$
1.0 rad/$\mu$s at $\Omega_{p0}=$ 0.4$\Gamma$. Thus, $R_{c}$ is 33 $\mu$s-1 and
the thermal equilibrium of stationary DSPs is feasible OurCommunPhys2021 .

## References

  * (1) Y.-W. Lin, W.-T. Liao, T. Peters, H.-C. Chou, J.-S. Wang, H.-W. Cho, P.-C. Kuan, and I. A. Yu, “Stationary Light Pulses in Cold Atomic Media and without Bragg Gratings,” Phys. Rev. Lett. 102, 213601 (2009).
  * (2) B. Kim, K.-T. Chen, C.-Y. Hsu, S.-S. Hsiao, Y.-C. Tseng, C.-Y. Lee, S.-L. Liang, Y.-H. Lai, J. Ruseckas, G. Juzeliūnas, and I. A. Yu, “Effect of laser-frequency fluctuation on the decay rate of Rydberg coherence,” Phys. Rev. A 100, 013815 (2019).
  * (3) Y.-F. Chen, Y.-C. Liu, Z.-H. Tsai, S.-H. Wang, and I. A. Yu, “Beat-note interferometer for direct phase measurement of photonic information,” Phys. Rev. A 72, 033812 (2005).
  * (4) B. Kim, K.-T. Chen, S.-S. Hsiao, S.-Y. Wang, K.-B. Li, J. Ruseckas, G. Juzeliūnas, T. Kirova, M. Auzinsh, Y.-C. Chen, Y.-F. Chen, and I. A. Yu, “A weakly-interacting many-body system of Rydberg polaritons based on electromagnetically induced transparency,” Commun. Phys. 4, 101 (2021).
  * (5) K.-T. Chen, B. Kim, C.-C. Su, S.-S. Hsiao, S.-J. Huang, W.-T. Liao, and I. A. Yu, “Increasing the decoherence rate of Rydberg polaritons due to accumulating dark Rydberg atoms,” Phys. Rev. Research 4, 023024 (2022).
  * (6) S.-S. Hsiao, K.-T. Chen, and I. A. Yu, “Mean field theory of weakly-interacting Rydberg polaritons in the EIT system based on the nearest-neighbor distribution,” Opt. Express 28, 28414 (2020).
  * (7) D. A. Braje, V. Balić, S. Goda, G. Y. Yin, and S. E. Harris, “Frequency Mixing Using Electromagnetically Induced Transparency in Cold Atoms,” Phys. Rev. Lett. 93, 183601 (2004).
  * (8) Y.-H. Chen, M.-J. Lee, W. Hung, Y.-C. Chen, Y.-F. Chen, and I. A. Yu, “Demonstration of the Interaction between Two Stopped Light Pulses,” Phys. Rev. Lett. 108, 173603 (2012).
  * (9) K.-F. Chang, T.-P. Wang, C.-Y. Chen, Y.-H. Chen, Y.-S. Wang, Y.-F. Chen, Y.-C. Chen, and I. A. Yu, “Low-loss high-fidelity frequency beam splitter with tunable split ratio based on electromagnetically induced transparency,” Phys. Rev. Research 3, 013096 (2021).
  * (10) C.-Y. Wang, Y.-F. Chen, S.-C. Lin, W.-H. Lin, P.-C. Kuan, and I. A. Yu, “Low-light-level all-optical switching,” Opt. Lett. 31, 2350 (2006).
  * (11) Y.-F. Chen, C.-Y. Wang, S.-H. Wang, and I. A. Yu, “Low-Light-Level Cross-phase-Modulation Based on Stored Light Pulses,” Phys. Rev. Lett. 96, 043603 (2006).
  * (12) S. Chandrasekhar, “Stochastic problems in physics and astronomy,” Rev. Mod. Phys. 15, 1 (1943).
  * (13) W.-T. Liao, T. Peters, E.-C. Shen, and I. A. Yu, “Propagation, Broadening, and Energy Decay of Quasi-Stationary Light Pulses in Thermal Atoms,” Chin. J. Phys. 47, 043603 (2009).
  * (14) Y.-H. Chen, M.-J. Lee, I.-C. Wang, S. Du, Y.-F. Chen, Y.-C. Chen, and I. A. Yu, “Coherent Optical Memory with High Storage Efficiency and Large Fractional Delay,” Phys. Rev. Lett. 110, 083601 (2013).
  * (15) M. Fleischhauer, J. Otterbach, and R. G. Unanyan, “Bose-Einstein Condensation of Stationary-Light Polaritons,” Phys. Rev. Lett. 101, 163601 (2008).