# Observation of Highly Correlated Ultrabright Biphotons Through Increased
Atomic Ensemble Density in Spontaneous Four-Wave Mixing

Jiun-Shiuan Shiu,1,2 Zi-Yu Liu,1,2 Chin-Yao Cheng,1,2 Yu-Chiao Huang,1,2 Ite
A. Yu,3,4 Ying-Cheng Chen,5 Chih-Sung Chuu,3,4 Che-Ming Li,2,6 Shiang-Yu
Wang,7 and Yong-Fan Chen1,2<EMAIL_ADDRESS>1Department of Physics,
National Cheng Kung University, Tainan 70101, Taiwan
2Center for Quantum Frontiers of Research $\&$ Technology, Tainan 70101,
Taiwan
3Center for Quantum Science and Technology, National Tsing Hua University,
Hsinchu 30013, Taiwan
4Department of Physics, National Tsing Hua University, Hsinchu 30013, Taiwan
5Institute of Atomic and Molecular Sciences, Academia Sinica, Taipei 10617,
Taiwan
6Department of Engineering Science, National Cheng Kung University, Tainan
70101, Taiwan
7Institute of Astronomy and Astrophysics, Academia Sinica, Taipei 10617,
Taiwan

(December 24, 2023)

###### Abstract

The pairing ratio, a crucial metric assessing a biphoton source’s ability to
generate correlated photon pairs, remains underexplored despite theoretical
predictions. This study presents experimental findings on the pairing ratio,
utilizing a double-$\Lambda$ spontaneous four-wave mixing biphoton source in
cold atoms. At an optical depth (OD) of 20, we achieved an ultrahigh biphoton
generation rate of up to $1.3\times 10^{7}$ per second, with a successful
pairing ratio of 61%. Increasing the OD to 120 significantly improved the
pairing ratio to 89%, while maintaining a consistent biphoton generation rate.
This achievement, marked by high generation rates and robust biphoton pairing,
holds great promise for advancing efficiency in quantum communication and
information processing. Additionally, in a scenario with a lower biphoton
generation rate of $5.0\times 10^{4}$ per second, we attained an impressive
signal-to-background ratio of 241 for the biphoton wavepacket, surpassing the
Cauchy-Schwarz criterion by approximately $1.5\times 10^{4}$ times.

Introduction.–Temporally correlated biphotons have recently garnered
considerable attention in the fields of optical quantum computing and quantum
communication, thanks to their exceptional nonclassical properties. Of
particular significance is their role as heralded single-photon sources, which
have found applications in diverse domains, including quantum cryptography [1,
2, 3, 4], quantum metrology [5, 6, 7, 8], and quantum imaging [9, 10, 11, 12].
Among the various biphoton sources, the spontaneous four-wave mixing (SFWM)
mechanism, distinguished by its operation near atomic resonance, stands out
for its ability to conveniently manipulate bandwidth and serves as a bridge
between different quantum devices, thus attracting significant interest.

The operation of SFWM near resonance accommodates diverse energy level
configurations, including the double-$\Lambda$ scheme [13, 14, 15] and
cascade-transition scheme [16, 17, 18]. This proximity to atomic resonance
enables the generation of bright biphotons with low optical power [19, 20], as
well as the production of narrowband biphotons [21, 22]. Especially for the
double-$\Lambda$ scheme, characterized by its intrinsic $\Lambda$-type
electromagnetically induced transparency (EIT) structure [23, 24, 25, 26, 27],
it not only significantly suppresses the generation of noise photons [28, 29,
30, 31] but also provides a wide bandwidth tuning capability [32, 33, 34].
This facilitates its direct application in conjunction with quantum devices
[35, 36, 37, 38] or reshaping of the biphoton waveforms [39, 40, 41].
Furthermore, the $\Lambda$ structure supports convenient implementation in
two- or three-level atomic systems [42, 43, 44, 45, 46].

Despite the numerous remarkable achievements of the double-$\Lambda$ SFWM
scheme, there is a frequently overlooked concern: the limited atomic density
hinders spontaneously emitted photons from achieving perfect coherence through
the four-wave mixing (FWM) process [47, 48]. Incoherent emissions not involved
in FWM processes can significantly diminish the biphoton pairing capacity,
referred to as the pairing ratio. Unfortunately, while the Heisenberg–Langevin
operator theory is predictive in understanding the concept of the pairing
ratio, there is currently a lack of relevant research and investigation in
experimental studies [49].

In this Letter, we present the results of a biphoton experiment conducted
using the double-$\Lambda$ configuration in a cold ensemble of 87Rb atoms via
SFWM. This configuration, chosen for its inherent EIT effect, allows for easy
control of the biphoton bandwidth. We achieved an exceptionally high biphoton
generation rate of approximately $1.3\times 10^{7}$ per second, surpassing
prior studies on the double-$\Lambda$ scheme. We thoroughly investigated the
pairing ratio under this ultrahigh generation rate, demonstrating its
enhancement through variations in optical depth (OD). Remarkably, we achieved
the highest pairing ratio of 0.89 at an OD of 120. Additionally, at a
relatively low biphoton generation rate of $5.0\times 10^{4}$ per second, the
signal-to-background ratio for the biphoton wavepacket reached 241, exceeding
the Cauchy-Schwarz criterion by approximately $1.5\times 10^{4}$ times.

Figure 1:  Diagram of the double-$\Lambda$ SFWM system and experimental setup.
M, mirror; L, lens; PBS, polarizing beam splitter; QWP, quarter-wave plate;
SMF, single-mode fiber; MMF, multimode fiber; EFS, etalon filter set; SPCM,
single-photon counting module. The inset shows the relevant energy levels of
the 87Rb atom.

Experimental setup.–We trapped cold 87Rb atoms using a standard magneto-
optical trap. After optically pumping them to the ground state
$|5S_{1/2},F=1\rangle$, as illustrated in Fig. 1, we irradiated the atomic
ensemble with a far-detuned driving field characterized by a Rabi frequency
$\Omega_{d}$ and a nearly resonant coupling field denoted by $\Omega_{c}$.
Synchronization between these fields was achieved through injection locking
with an external cavity diode laser (not shown). The driving field, with
detuning $\Delta_{d}$, operated on the $\sigma^{+}$-transition
$|5S_{1/2},F=1\rangle\rightarrow|5P_{3/2},F=2\rangle$, effectively suppressing
incoherent fluorescence from one-photon absorption. This allowed for primary
emission of Stokes photons via spontaneous Raman scattering [50].
Subsequently, the nearly resonant coupling field, with detuning $\Delta_{c}$,
acted on the $\sigma^{+}$-transition
$|5S_{1/2},F=2\rangle\rightarrow|5P_{3/2},F=2\rangle$, inducing the emission
of anti-Stokes photons. We used an elongated atomic ensemble to enhance the
FWM effect and boosting specific direction scattering probability.

In the SFWM experiment, biphotons were generated using a 10-$\mu$s driving
pulse in each 2.5-ms cycle. To prevent laser leakage, we adopted a backward
configuration where the driving and coupling beams counter-propagated, each
with $1/e^{2}$ full widths of 250 and 310 $\mu$m. The corresponding optical
powers of the 1$\Gamma$ driving and coupling fields were approximately 7.5 and
11.1 $\mu$W. Both fields drove $\sigma^{+}$-transitions, resulting in
generated photon pairs exhibiting $\sigma^{+}$ polarization, propagating in
opposite directions and passing through respective etalon filter sets (EFS).
The intersection angle between the Stokes (anti-Stokes) and driving (coupling)
beams was set at 1.7∘ for our experiment. Each EFS consisted of two etalons,
each with an extinction ratio of roughly 30 dB and an approximate bandwidth of
100 MHz, separated by an optical isolator. The total extinction ratios of the
Stokes and anti-Stokes channels were 114 dB and 124 dB, respectively.
Biphotons were detected using fiber-coupled single-photon counting modules
(SPCM-AQRH-13-FC). Upon detection, an 8-ns pulse was emitted from the SPCM
toward the time-of-flight multiscaler (TOF, MCS6A-4T8, not shown). In the
coincidence count experiment, we measured the time difference between Stokes
and anti-Stokes photons. When generated as a correlated pair, they arrived at
the SPCMs within the correlation time, contributing to a nonflat biphoton
wavepacket. The TOF generated a histogram of coincident counts based on these
data points, providing insight into the source of biphotons.

Theoretical model.–We used the Heisenberg–Langevin operator approach to
analyze the biphoton generation in double-$\Lambda$ SFWM [49]. More details
can be found in the Supplemental Material. The photon generation rate, given
by $R=\frac{c}{L}\langle\hat{a}^{\dagger}\hat{a}\rangle$ and derived from the
annihilation operator $\hat{a}$, leads to the following expressions:

$\displaystyle R_{s}=$
$\displaystyle\int\frac{d\omega}{2\pi}\left(|B|^{2}+\sum_{jk,j^{\prime}k^{\prime}}\int_{0}^{L}dz\,P_{jk}^{*}\mathcal{D}_{jk^{\dagger},j^{\prime}k^{\prime}}P_{j^{\prime}k^{\prime}}\right)$
$\displaystyle\equiv$ $\displaystyle\int d\omega\widetilde{R}_{s}(\omega),$
(1) $\displaystyle R_{as}=$
$\displaystyle\int\frac{d\omega}{2\pi}\left(|C|^{2}+\sum_{jk,j^{\prime}k^{\prime}}\int_{0}^{L}dz\,Q_{jk}\mathcal{D}_{jk,j^{\prime}k^{\prime\dagger}}Q_{j^{\prime}k^{\prime}}^{*}\right)$
$\displaystyle\equiv$ $\displaystyle\int d\omega\widetilde{R}_{as}(\omega),$
(2)

where $\mathcal{D}_{jk^{\dagger},j^{\prime}k^{\prime}}$ and
$\mathcal{D}_{jk,j^{\prime}k^{\prime\dagger}}$ are diffusion coefficients,
while $\widetilde{R}_{s(as)}$ represents the spectrum of Stokes (anti-Stokes)
photons. The total photon generation rates, $R_{s}$ and $R_{as}$, comprise two
components: correlated photons governed by the FWM process with coefficients
$B$ and $C$, and uncorrelated photons due to vacuum field fluctuations
represented by an integral term with diffusion coefficients. The pairing
ratio, $r_{p}$, denotes the ratio of the correlated photons to the total
generated photons. Under ideal conditions, $R_{s}$ and $R_{as}$ are nearly
equal. However, in experiments, $R_{as}$ is slightly smaller due to phase
mismatch and ground state decoherence. Therefore, the biphoton generation rate
$R_{\rm B}$ is contingent on $R_{as}$.

In the SFWM process within the atomic ensemble, a spontaneously emitted Stokes
photon from one atom may interact with nearby atoms, triggering stimulated
Raman scattering. Unlike spontaneous Raman scattering, the Stokes photon
generated by stimulated Raman scattering shares the same direction as the
incident Stokes photon, enhancing directionality. This collective enhancement
effect establishes paired correlation with the anti-Stokes photon through the
FWM process, reflected in coefficients $B$ and $C$. However, as SFWM relies on
vacuum field fluctuations, the generated photons exhibit isotropic
(uncorrelated) nature, represented by the integral term with diffusion
coefficients in Eqs. (1) and (2). While the paired correlations of these
biphotons can be established through the FWM process, it requires a
sufficiently high OD within the atomic ensemble.

In biphoton systems, the normalized Glauber second-order cross-correlation
function $g_{s\text{-}as}^{(2)}(\tau)$ is often used alongside the photon
generation rate. This function is a crucial parameter for evaluating the
temporal correlation between biphotons. The derived theoretical expression is
as follows:

$\displaystyle g_{s\text{-}as}^{(2)}(\tau)=$ $\displaystyle
1+\frac{1}{R_{s}R_{as}}\left|\int\frac{d\omega}{2\pi}e^{-i\omega\tau}\left(\vphantom{\frac{A}{B}}B^{*}D\right.\right.$
$\displaystyle+\sum_{jk,j^{\prime}k^{\prime}}\int_{0}^{L}dz\left.\left.\vphantom{\frac{A}{B}}P_{jk}^{*}\mathcal{D}_{jk^{\dagger},j^{\prime}k^{\prime}}Q_{j^{\prime}k^{\prime}}\right)\right|^{2}.$
(3)

The integral term on the right-hand side of Eq. (3) reveals the correlation of
biphotons. This term is equivalent to the wavepacket of the anti-Stokes single
photon, conditioned on the post-selection of a Stokes single photon. This
correlation provides valuable information for evaluating the biphoton source.
For instance, the peak signal-to-background ratio, denoted as $r_{\rm SB}$, is
defined as the maximum value of $[g_{s\text{-}as}^{(2)}(\tau)-1]$. It serves
as a standard metric for assessing the nonclassicality of a biphoton source.
In the case of a classical field, the Cauchy–Schwarz inequality universally
applies:
$\left[g_{s\text{-}as}^{(2)}(\tau)\right]^{2}\left[g_{s\text{-}s}^{(2)}(0)g_{as\text{-}as}^{(2)}(0)\right]^{-1}\leq
1.$ The normalized autocorrelation functions of the Stokes and anti-Stokes
fields can be derived as $g_{s\text{-}s}^{(2)}(\tau)=1+R_{s}^{-2}\left|\int
d\omega\widetilde{R}_{s}e^{-i\omega\tau}\right|^{2}$ and
$g_{as\text{-}as}^{(2)}(\tau)=1+R_{as}^{-2}\left|\int
d\omega\widetilde{R}_{as}e^{-i\omega\tau}\right|^{2}$. These equations
indicate that both the Stokes and anti-Stokes fields exhibit thermal states,
with $g_{s\text{-}s}^{(2)}(0)=g_{as\text{-}as}^{(2)}(0)=2$. Nonclassical
behavior is observed when $r_{\rm SB}>1$. Additional details and initial
proofs of the thermal field distributions for both the Stokes and anti-Stokes
fields can be found in the Supplemental Material.

Biphoton bandwidth.–Figure 2 presents the experimental biphoton coincidence
count rate $R_{\rm C}$ (or biphoton temporal wavepacket) for various coupling
field conditions. $R_{\rm C}$ is calculated as
$R_{s}R_{as}g_{s\text{-}as}^{(2)}(\tau)\Delta T+R_{\rm env}$ (refer to
Supplemental Material), where $R_{\rm env}$ accounts for environmental
background count rates, arising from laser leakage or SPCM dark counts. For
data processing, we used a time bin of $\Delta T=1/R_{s}$ to tally Stokes
photons, enabling the post-selection of a single Stokes photon. This ensures
that the background and correlated regions of the coincidence count rate
correspond to $R_{\rm B}+R_{\rm env}$ and $r_{p}$, respectively. The time bin
for detected anti-Stokes photons, $\Delta\tau=6.4$ ns, aligns with the time
interval between experimental data points in Fig. 2.

Figure 2:  Biphotons with controllable bandwidth. The red lines represent the
theoretical curves, while the black circles indicate the experimental data
points. The time bin for detecting the anti-Stokes photons is $\Delta\tau=6.4$
ns. Other parameters are ${\rm OD}=15$, $\Omega_{d}=1\Gamma$,
$\Delta_{d}=10\Gamma$, $\gamma_{21}=0.001\Gamma$, $\Delta kL=0.37\pi$, (a)
$\Omega_{c}=4\Gamma$, $\Delta_{c}=0\Gamma$, (b) $\Omega_{c}=1\Gamma$,
$\Delta_{c}=0\Gamma$, (c) $\Omega_{c}=1\Gamma$, $\Delta_{c}=1\Gamma$, (d)
$\Omega_{c}=1\Gamma$, $\Delta_{c}=3\Gamma$.

In Figs. 2(a) and 2(b), we set the coupling Rabi frequencies $\Omega_{c}$ to
4$\Gamma$ and 1$\Gamma$, respectively. With the OD fixed at 15, both cases
yielded a measured biphoton generation rate $R_{\rm B}$ of approximately
$3.4\times 10^{5}\text{s}^{-1}$. The delay time in Fig. 2(b) is significantly
longer than that in 2(a). This delay arises from two intrinsic properties of
the double-$\Lambda$ SFWM system: the damped Rabi oscillation with a period
denoted as $\tau_{\rm R}=2\pi/\sqrt{|\Omega_{c}|^{2}-\Gamma^{2}/4}$, and the
delay time attributed to the EIT effect denoted as $\tau_{\rm EIT}=\Gamma{\rm
OD}/|\Omega_{c}|^{2}$ [51]. Both characteristic times are influenced by the
coupling field. The damped Rabi oscillation periods in Figs. 2(a) and 2(b) are
calculated as 42 and 192 ns, respectively. As $\Omega_{c}$ decreases, the EIT
effect causes anti-Stokes photons to propagate slowly. The EIT delay times in
Figs. 2(a) and 2(b) are 25 and 398 ns, respectively. The overall delay time is
determined by the larger of $\tau_{\rm R}$ and $\tau_{\rm EIT}$, i.e.,
max$(\tau_{\rm R},\tau_{\rm EIT})$. Consequently, the behavior of the biphoton
wavepacket in Fig. 2(b), where $\tau_{\rm EIT}$ dominates, exhibits
characteristics reminiscent of slow light, with the slow light effect
noticeable in the trailing edge of the biphoton wavepacket. Conversely, in
Fig. 2(a), where $\tau_{\rm R}$ surpasses $\tau_{\rm EIT}$, subtle oscillatory
features are present within the biphoton wavepacket.

Figures 2(c) and 2(d) demonstrate how changes in coupling detuning
$\Delta_{c}$ affect the biphoton bandwidth. All experimental parameters were
consistent with those in Fig. 2(b), except for the $\Delta_{c}$. A shorter
delay time in Fig. 2(c) is observed due to the introduction of
$\Delta_{c}=1\Gamma$, which reduces the effective OD and shortens the EIT-
induced delay. Conversely, with $\Delta_{c}=3\Gamma$ in Fig. 2(d), the tail
lengthens again. This behavior is attributed to damped Rabi oscillations,
where a larger $\Delta_{c}$ weakens the interaction between the coupling field
and the atomic medium, requiring more time to convert spinwave excitations
into anti-Stokes photons. The values of $R_{\rm B}$ in Figs. 2(c) and 2(d) are
$3.4\times 10^{5}\text{s}^{-1}$ and $3.1\times 10^{5}\text{s}^{-1}$,
respectively. This demonstrates that by detuning the coupling field, we can
control the biphoton bandwidth without significantly reducing $R_{\rm B}$.
Further discussions can be found in the Supplemental Material.

Figure 3:  High-purity biphotons. The red lines represent the theoretical
curves, while the black circles indicate the experimental data points. The
time bin for detecting the anti-Stokes photons is $\Delta\tau=1.6$ ns. The
remaining parameters are set to ${\rm OD}=10$, $\Omega_{c}=4\Gamma$,
$\Delta_{d}=10\Gamma$, $\gamma_{21}=0.001\Gamma$, $\Delta kL=0.37\pi$, with
(a) $\Omega_{d}=0.5\Gamma$ and (b) $\Omega_{d}=3\Gamma$. (c) The peak signal-
to-background ratio $r_{\rm SB}$ versus $\Omega_{d}$. The black squares
represent the experimental data, and the black line is the curve fitted to
these experimental data points. (d) The biphoton generation rate $R_{\rm B}$
and pairing ratio $r_{p}$ as a function of $\Omega_{d}$. The experimental data
points for $R_{\rm B}$ and $r_{p}$ are represented by the unfilled blue and
solid magenta circles, respectively. The theoretical curves for $R_{\rm B}$
and $r_{p}$ are depicted by the blue and magenta lines, respectively.

High-purity biphotons.–Figure 3(a) illustrates high-purity biphoton generation
achieved with parameters $\Omega_{d}=0.5\Gamma$, $\Omega_{c}=4\Gamma$, and
${\rm OD}=10$. The theoretical $R_{\rm B}$ is calculated as $5.0\times
10^{4}\text{s}^{-1}$. In experiments, $R_{\rm B}$ was determined by
subtracting the total measured background count rate, $R_{\rm tot}=1.6\times
10^{5}\text{s}^{-1}$, from the environmental background count rate, $R_{\rm
env}=1.1\times 10^{5}\text{s}^{-1}$. This yielded an experimental $R_{\rm B}$
of approximately $5.0\times 10^{4}\text{s}^{-1}$, in close agreement with the
theoretical prediction. In Fig. 3(b), we theoretically calculated $R_{\rm B}$
to be $1.9\times 10^{6}\text{s}^{-1}$ at $\Omega_{d}=3\Gamma$. The
experimentally observed $R_{\rm B}$, obtained from measurements of $R_{\rm
tot}=2.0\times 10^{6}\text{s}^{-1}$ and $R_{\rm env}=1.2\times
10^{5}\text{s}^{-1}$, also closely matches theoretical prediction. As
$\Omega_{d}$ increases, both $R_{\rm B}$ and $R_{\rm tot}$ rise significantly.
However, this also leads to a notable decrease in $r_{\rm SB}$, as shown in
Fig. 3(c). At $\Omega_{d}=0.5\Gamma$, we observed an experimental $r_{\rm
SB}=241$, surpassing the Cauchy–Schwarz criterion by a factor of approximately
$1.5\times 10^{4}$. If the $R_{\rm env}$ in our experiment could be completely
eliminated, it would lead to a more pronounced violation of the Cauchy–Schwarz
criterion, exceeding the normal level by a factor of $5.9\times 10^{4}$.

Figure 3(d) shows the variation of $R_{\rm B}$ and $r_{p}$ with different
$\Omega_{d}$ values. At $\Omega_{d}=0.5\Gamma$ and $\Omega_{d}=3\Gamma$, the
corresponding $r_{p}$ values are 0.63 and 0.59, respectively. These
experimental $r_{p}$ values were determined by integrating the area under the
correlated biphoton wavepacket. In the SFWM process, atomic ensembles play a
crucial role in collectively enhancing the correlation between the Stokes and
anti-Stokes fields along the applied light direction. Therefore, with a fixed
OD, while increasing $\Omega_{d}$ can boost $R_{\rm B}$, the limited density
of atomic ensembles constrains their ability to produce correlated photon
pairs, leading to a slight decrease in $r_{p}$.

Figure 4:  Highly correlated ultrabright biphotons. The red lines represent
the theoretical curves, while the black circles indicate the experimental data
points. The time bin for detecting the anti-Stokes photons is $\Delta\tau=1.6$
ns. The remaining parameters are set to $\Omega_{d}=3\Gamma$,
$\gamma_{21}=0.001\Gamma$, $\Delta kL=0.37\pi$, with (a) ${\rm OD}=20$,
$\Omega_{c}=4\Gamma$, $\Delta_{d}=5\Gamma$, and (b) ${\rm OD}=120$,
$\Omega_{c}=8.8\Gamma$, $\Delta_{d}=14.9\Gamma$. (c) The peak signal-to-
background ratio $r_{\rm SB}$ versus OD. The black squares represent the
experimental data, and the black line is the curve fitted to these
experimental data points. (d) The biphoton generation rate $R_{\rm B}$ and
pairing ratio $r_{p}$ as a function of OD. The experimental data points for
$R_{\rm B}$ and $r_{p}$ are represented by the unfilled blue and solid magenta
circles, respectively. The theoretical curves for $R_{\rm B}$ and $r_{p}$ are
depicted by the blue and magenta lines, respectively.

Highly correlated biphotons.–Figure 4(a) showcases the generation of
ultrabright biphotons using specific parameters: $\Omega_{d}=3\Gamma$,
$\Omega_{c}=4\Gamma$, $\Delta_{d}=5\Gamma$, and ${\rm OD}=20$, resulting in a
remarkable theoretical $R_{\rm B}$ of $1.3\times 10^{7}\text{s}^{-1}$. This
exceeds rates reported in the literature for the double-$\Lambda$ SFWM scheme.
Under these conditions, the experimental total background count rate was
$R_{\rm tot}=1.3\times 10^{7}\text{s}^{-1}$, with environmental background at
$R_{\rm env}=2.3\times 10^{5}\text{s}^{-1}$, accounting for only 1.8$\%$ of
the total. Thus, in this high $R_{\rm B}$ scenario, the primary source of
background count arises from the high photon generation rate rather than
environmental factors. Furthermore, in this scenario, the measured $r_{\rm
SB}$ was 2.4, surpassing the Cauchy-Schwarz criterion by a factor of 2.9,
while $r_{p}$ was only 0.61. Although increasing the coupling power can
enhance $r_{\rm SB}$, as demonstrated in Fig. 2, it does not lead to
corresponding improvements in $r_{p}$. To enhance both $r_{p}$ and $r_{\rm
SB}$, we further increased the OD. In addition to ${\rm OD}=20$, we measured
the biphoton wavepacket at ${\rm OD}=40$, $60$, $80$, $100$, and $120$. We
fine-tuned $\Omega_{c}$ to maintain a consistent biphoton bandwidth, while
keeping $\Omega_{d}=3\Gamma$ constant and adjusting $\Delta_{d}$ to maintain a
theoretucal $R_{\rm B}$ of $1.3\times 10^{7}\text{s}^{-1}$. Specific
parameters can be found in the Supplemental Material.

In Fig. 4(b), we present the scenario with ${\rm OD}=120$,
$\Omega_{c}=8.8\Gamma$, and $\Delta_{d}=14.9\Gamma$. Here, the measured
$r_{\rm SB}$ at 4.2 exceeds the Cauchy–Schwarz criterion by 6.8 times. An
evident positive correlation emerges between increased OD and enhanced
$r_{p}$, resulting in a higher $r_{\rm SB}$ due to augmented coincidence
receptions [Fig. 4(c)]. This enhancement stems from the increased accumulation
of biphoton correlations along a specific direction at higher OD values.
Photons generated at higher OD levels are more likely to encounter subsequent
atoms, amplifying the collective enhancement through the FWM process. Figure
4(d) illustrates the relationships between $R_{\rm B}$ and $r_{p}$ with OD. At
${\rm OD}=120$, we observed the highest $r_{p}$ of 0.89, indicating a
significant improvement in correlated photon pair generation. The experimental
$R_{\rm B}=1.3\times 10^{7}\text{s}^{-1}$ signifies the successful generation
of approximately $1.2\times 10^{7}$ pairs of correlated photons per second.
Additionally, the Fourier transform of $\left(R_{\rm C}-R_{\rm tot}\right)$
reveals a biphoton bandwidth of approximately 24 MHz, resulting in a spectral
brightness of the biphoton source at $5.4\times
10^{5}\text{s}^{-1}\text{MHz}^{-1}$, surpassing the highest achieved by sub-
megahertz biphoton sources [52]. These results highlight the crucial role of
high OD in SFWM-based biphoton sources, allowing for higher values of $R_{\rm
B}$ and $r_{p}$. This enables the generation of a large quantity of high-
quality correlated photon pairs for use in various quantum systems.

Conclusion.–Our investigation into the biphoton pairing ratio, utilizing the
double-$\Lambda$ SFWM in cold 87Rb atoms, revealed a marginal decrease with
higher biphoton generation rates. Nonetheless, this trend can be effectively
addressed by elevating the atomic ensemble density. The highest pairing ratio
observed was 0.89 at an OD of 120, accompanied by an ultrabright biphoton
generation rate of up to $1.3\times 10^{7}\text{s}^{-1}$, surpassing
previously reported rates achieved via the double-$\Lambda$ SFWM scheme.
Furthermore, our experiment demonstrated the highest signal-to-background
ratio of the biphoton wavepacket at 241, achieved at a low biphoton generation
rate of $5.0\times 10^{4}\text{s}^{-1}$. This outstanding performance exceeded
the Cauchy–Schwarz criterion by approximately $1.5\times 10^{4}$ times. These
results underscore the capability of the double-$\Lambda$ SFWM scheme in
advancing biphoton sources for future quantum technologies.

We thank Meng-Jung Lin, Chih-Min Yang, I-Chia Huang, and Ting-Ho Wu for their
contributions to the initial setup of the experimental system. This work was
supported by the National Science and Technology Council of Taiwan under Grant
Nos. of 112-2112-M-006-034, 111-2639-M-007-001-ASP, and 111-2119-M-007-007.

## References

  * [1] C. Liu, S. Zhang, L. Zhao, P. Chen, C.-H. Fung, H. F. Chau, M. M. T. Loy, and S. Du, Differential-phase-shift quantum key distribution using heralded narrow-band single photons, Opt. Express 21, 9505 (2013).
  * [2] C. H. Bennett and G. Brassard, Quantum cryptography: public key distribution and coin tossing, Theor. Comput. Sci. 560, 7 (2014).
  * [3] F. Xu, X. Ma, Q. Zhang, H.-K. Lo, and J.-W. Pan, Secure quantum key distribution with realistic devices, Rev. Mod. Phys. 92, 025002 (2020).
  * [4] E. Fitzke, L. Bialowons, T. Dolejsky, M. Tippmann, O. Nikiforov, T. Walther, F. Wissel, and M. Gunkel, Scalable network for simultaneous pairwise quantum key distribution via entanglement-based time-bin coding, PRX Quantum 3, 020341 (2022).
  * [5] P. R. Tapster, S. F. Seward, and J. G. Rarity, Sub-shot-noise measurement of modulated absorption using parametric down-conversion, Phys. Rev. A 44, 3266 (1991).
  * [6] S. Slussarenko, M. M. Weston, H. M. Chrzanowski, L. K. Shalm, V. B. Verma, S. W. Nam, and G. J. Pryde, Unconditional violation of the shot-noise limit in photonic quantum metrology, Nat. Photonics 11, 700 (2017).
  * [7] Luca Pezzè, Augusto Smerzi, Markus K. Oberthaler, Roman Schmied, and Philipp Treutlein, Quantum metrology with nonclassical states of atomic ensembles, Rev. Mod. Phys. 90, 035005 (2018).
  * [8] Y. Chen, L. Hong, and L. Chen, Quantum interferometric metrology with entangled photons, Front. Phys. 10, 892519 (2022).
  * [9] T. B. Pittman, Y. H. Shih, D. V. Strekalov, and A. V. Sergienko, Optical imaging by means of two-photon quantum entanglement, Phys. Rev. A 52, R3429 (1995).
  * [10] Y. Shih, Quantum imaging, IEEE J. Sel. Top. Quantum Electron. 13, 1016 (2007).
  * [11] D.-S. Ding, Z.-Y. Zhou, B.-S. Shi, and G.-C. Guo, Single-photon-level quantum image memory based on cold atomic ensembles, Nat. Commun. 4, 2527 (2013).
  * [12] P.-A. Moreau, E. Toninelli, T. Gregory, and M. J. Padgett, Imaging with Quantum States of Light, Nat. Rev. Phys. 1, 367 (2019).
  * [13] V. Balić, D. A. Braje, P. Kolchin, G. Y. Yin, and S. E. Harris, Generation of paired photons with controllable waveforms, Phys. Rev. Lett. 94, 183601 (2005).
  * [14] C. Shu, X. Guo, P. Chen, M. M. T. Loy, and S. Du, Narrowband biphotons with polarization-frequency-coupled entanglement, Phys. Rev. A 91, 043820 (2015).
  * [15] C. Shu, P. Chen, T. K. A. Chow, L. Zhu, Y. Xiao, M. M. T. Loy, and S. Du, Subnatural-linewidth biphotons from a Doppler-broadened hot atomic vapour cell, Nat. Commun. 7, 12783 (2016).
  * [16] T. Chanelière, D. N. Matsukevich, S. D. Jenkins, T. A. B. Kennedy, M. S. Chapman, and A. Kuzmich, Quantum telecommunication based on atomic cascade transitions, Phys. Rev. Lett. 96, 093604 (2006).
  * [17] B. Srivathsan, G. K. Gulati, B. Chng, G. Maslennikov, D. Matsukevich, and C. Kurtsiefer, Narrow band source of transform-limited photon pairs via four-wave mixing in a cold atomic ensemble, Phys. Rev. Lett. 111, 123602 (2013).
  * [18] J. Park, H. Kim, and H. S. Moon, Polarization-entangled photons from a warm atomic ensemble using a Sagnac interferometer, Phys. Rev. Lett. 122, 143601 (2019).
  * [19] D.-S. Ding, Z.-Y. Zhou, B.-S. Shi, X.-B. Zou, and G.-C. Guo, Generation of non-classical correlated photon pairs via a ladder-type atomic configuration: theory and experiment, Opt. Express 20, 11433 (2012).
  * [20] A. Bruns, C.-Y. Hsu, S. Stryzhenko, E. Giese, L. P. Yatsenko, I. A. Yu, T. Halfmann, and T. Peters, Ultrabright and narrowband intra-fiber biphoton source at ultralow pump power, Quantum Sci. Technol. 8, 015002 (2023).
  * [21] L. Zhao, X. Guo, C. Liu, Y. Sun, M. M. T. Loy, and S. Du, Photon pairs with coherence time exceeding 1 $\mu$s, Optica 1, 84 (2014).
  * [22] Y.-S. Wang, K.-B. Li, C.-F. Chang, T.-W. Lin, J.-Q. Li, S.-S. Hsiao, J.-M. Chen, Y.-H. Lai, Y.-C. Chen, Y.-F. Chen, C.-S. Chuu, and I. A. Yu, Temporally ultralong biphotons with a linewidth of 50 kHz, APL Photonics 7, 126102 (2022).
  * [23] S. E. Harris, J. E. Field, and A. Imamoğlu, Nonlinear optical processes using electromagnetically induced transparency, Phys. Rev. Lett. 64, 1107 (1990).
  * [24] K.-J. Boller, A. Imamoğlu, and S. E. Harris, Observation of electromagnetically induced transparency, Phys. Rev. Lett. 66, 2593 (1991).
  * [25] L. V. Hau, S. E. Harris, Z. Dutton, and C. H. Behroozi, Light speed reduction to 17 metres per second in an ultracold atomic gas, Nature 397, 594 (1999).
  * [26] M. D. Lukin, Colloquium: Trapping and manipulating photon states in atomic ensembles, Rev. Mod. Phys. 75, 457 (2003).
  * [27] M. Fleischhauer, A. Imamoglu, and J. P. Marangos, Electromagnetically induced transparency: optics in coherent media, Rev. Mod. Phys. 77, 633 (2005).
  * [28] M. Fleischhauer and M. D. Lukin, Dark-state polaritons in electromagnetically induced transparency, Phys. Rev. Lett. 84, 5094 (2000).
  * [29] A. Peng, M. Johnsson, W. P. Bowen, P. K. Lam, H.-A. Bachor, and J. J. Hope, Squeezing and entanglement delay using slow light, Phys. Rev. A 71, 033809 (2005).
  * [30] C.-Y. Cheng, J.-J. Lee, Z.-Y. Liu, J.-S. Shiu, and Y.-F. Chen, Quantum frequency conversion based on resonant four-wave mixing. Phys. Rev. A 103, 023711 (2021).
  * [31] H. Hsu, C.-Y. Cheng, J.-S. Shiu, L.-C. Chen, and Y.-F. Chen, Quantum fidelity of electromagnetically induced transparency: the full quantum theory, Opt. Express 30, 2097 (2022).
  * [32] S. Du, P. Kolchin, C. Belthangady, G. Y. Yin, and S. E. Harris, Subnatural linewidth biphotons with controllable temporal length, Phys. Rev. Lett. 100, 183603 (2008).
  * [33] C.-Y. Hsu, Y.-S. Wang, J.-M. Chen, F.-C. Huang, Y.-T. Ke, E. K. Huang, W. Hung, K.-L. Chao, S.-S. Hsiao, Y.-H. Chen, C.-S. Chuu, Y.-C. Chen, Y.-F. Chen, and I. A. Yu, Generation of sub-MHz and spectrally-bright biphotons from hot atomic vapors with a phase mismatch-free scheme, Opt. Express 29, 4632 (2021).
  * [34] S.-S. Hsiao, W.-K. Huang, Y.-M. Lin, J.-M. Chen, C.-Y. Hsu, and I. A. Yu, Temporal profile of biphotons generated from a hot atomic vapor and spectrum of electromagnetically induced transparency, Phys. Rev. A 106, 023709 (2022).
  * [35] S. Zhou, S. Zhang, C. Liu, J. F. Chen, J. Wen, M. M. T. Loy, G. K. L. Wong, and S. Du, Optimal storage and retrieval of single-photon waveforms, Opt. Express 20, 24124 (2012).
  * [36] Y. Wang, J. Li, S. Zhang, K. Su, Y. Zhou, K. Liao, S. Du, H. Yan, and S.-L. Zhu, Efficient quantum memory for single-photon polarization qubits, Nat. Photonics 13, 346 (2019).
  * [37] C.-Y. Cheng, Z.-Y. Liu, P.-S. Hu, T.-N. Wang, C.-Y. Chien, J.-K. Lin, J.-Y. Juo, J.-S. Shiu, I. A. Yu, Y.-C. Chen, and Y.-F. Chen, Efficient frequency conversion based on resonant four-wave mixing, Opt. Lett. 46, 681 (2021).
  * [38] Z.-Y. Liu, J.-S. Shiu, C.-Y. Cheng, and Y.-F. Chen, Controlling frequency-domain Hong-Ou-Mandel interference via electromagnetically induced transparency, Phys. Rev. A 108, 013702 (2023).
  * [39] P. Kolchin, C. Belthangady, S. Du, G.Y. Yin, and S. E. Harris, Electro-optic modulation of single photons, Phys. Rev. Lett. 101, 103601 (2008).
  * [40] J. F. Chen, S. Zhang, H. Yan, M. M. T. Loy, G.K. L. Wong, and S. Du, Shaping biphoton temporal waveforms with modulated classical fields, Phys. Rev. Lett. 104, 183604 (2010).
  * [41] L. Zhao, X. Guo, Y. Sun, Y. Su, M. M. T. Loy, and S. Du, Shaping the Biphoton Temporal Waveform with Spatial Light Modulation, Phys. Rev. Lett. 115, 193601 (2015).
  * [42] P. Kolchin, S. Du, C. Belthangady, G. Y. Yin, and S. E. Harris, Generation of narrow-bandwidth paired photons: use of a single driving laser, Phys. Rev. Lett. 97, 113602 (2006).
  * [43] J. Wen, S. Du, and M. H. Rubin, Biphoton generation in a two-level atomic ensemble, Phys. Rev. A 75, 033809 (2007).
  * [44] S. Du, J. Wen, M. H. Rubin, and G. Y. Yin, Four-wave mixing and biphoton generation in a two-level system, Phys. Rev. Lett. 98, 053601 (2007).
  * [45] M. O. Araújo, L. S. Marinho, and D. Felinto, Observation of nonclassical correlations in biphotons generated from an ensemble of pure two-level atoms, Phys. Rev. Lett. 128, 083601 (2022).
  * [46] L. S. Marinho, M. O. Araújo, W. Martins, and D. Felinto, Enhancing nonclassical correlations for light scattered by an ensemble of cold two-level atoms, Opt. Lett. 48, 3323 (2023).
  * [47] C.-K. Chiu, Y.-H. Chen, Y.-C. Chen, I. A. Yu, Y.-C. Chen, and Y.-F. Chen, Low-light-level four-wave mixing by quantum interference, Phys. Rev. A 89, 023839 (2014).
  * [48] J.-Y. Juo, J.-K. Lin, C.-Y. Cheng, Z.-Y. Liu, I. A. Yu, and Y.-F. Chen, Demonstration of spatial-light-modulation-based four-wave mixing in cold atoms, Phys. Rev. A 97, 053815 (2018).
  * [49] P. Kolchin, Electromagnetically-induced-transparency-based paired photon generation, Phys. Rev. A 75, 033814 (2007).
  * [50] M. G. Raymer and J. Mostowski, M. G. Raymer and J. Mostowski, Stimulated Raman scattering: unified treatment of spontaneous initiation and spatial propagation, Phys. Rev. A 24, 1980 (1981).
  * [51] S. Du, J. Wen, and M. H. Rubin, Narrowband biphoton generation near atomic resonance, J. Opt. Soc. Am. B 25, C98 (2008).
  * [52] J.-M. Chen, C.-Y. Hsu, W.-K. Huang, S.-S. Hsiao, F.-C. Huang, Y.-H. Chen, C.-S. Chuu, Y.-C. Chen, Y.-F. Chen, and I. A. Yu, Room-temperature biphoton source with a spectral brightness near the ultimate limit, Phys. Rev. Res. 4, 023132 (2022).

Supplemental Material for Observation of Highly Correlated Ultrabright
Biphotons Through Increased Atomic Ensemble Density in Spontaneous Four-Wave
Mixing

Jiun-Shiuan Shiu,1,2 Zi-Yu Liu,1,2 Chin-Yao Cheng,1,2 Yu-Chiao Huang,1,2 Ite
A. Yu,3,4 Ying-Cheng

Chen,5 Chih-Sung Chuu,3,4 Che-Ming Li,2,6 Shiang-Yu Wang,7 and Yong-Fan
Chen1,2

1Department of Physics, National Cheng Kung University, Tainan 70101, Taiwan
2Center for Quantum Frontiers of Research $\&$ Technology, Tainan 70101,
Taiwan
3Center for Quantum Science and Technology, National Tsing Hua University,
Hsinchu 30013, Taiwan
4Department of Physics, National Tsing Hua University, Hsinchu 30013, Taiwan
5Institute of Atomic and Molecular Sciences, Academia Sinica, Taipei 10617,
Taiwan
6Department of Engineering Science, National Cheng Kung University, Tainan
70101, Taiwan
7Institute of Astronomy and Astrophysics, Academia Sinica, Taipei 10617,
Taiwan

This supplemental material provides a comprehensive explanation of the
theoretical model, offering detailed insights into experimental parameters and
data analysis. In Section I, we introduce the theoretical model, covering the
derivation of the field operator evolution in the double-$\Lambda$ spontaneous
four-wave mixing atomic ensemble, biphoton generation rate, and corresponding
correlation functions. Section II presents a proof that Stokes and anti-Stokes
photons adhere to a thermal state distribution. Finally, Section III
elaborates on experimental measurements, including the coincidence count rate
and biphoton collection efficiency. Additionally, it discusses the biphoton
temporal profile and the experimental parameters presented in Figure 4 of the
main text.

##  I Theoretical Model

###  I.1 Derivation of Field Operators

In this research, we used the Heisenberg–Langevin operator approach to
characterize the biphotons generated with the spontaneous four-wave mixing
(SFWM) process [1]. The slowly varying collective atomic operator
$\hat{\sigma}_{jk}(z,t)=\sum_{l=1}^{N_{z}}\hat{\sigma}^{l}_{jk}(z,t)/N_{z}$
was used to model the intricate atomic dynamics during the interaction, where
$N_{z}$ is the number of atoms within the infinitesimal one-dimensional (1D)
spatial interval $dz$. Here, $\hat{\sigma}^{l}_{jk}$ denotes the adiabatic
atomic operator of the $l$th atom that has been transformed by the slowly
varying amplitude. The dynamics of slowly varying collective atomic operators
are governed by the Heisenberg–Langevin equations (HLEs):
$\frac{\partial}{\partial
t}\hat{\sigma}_{jk}=\frac{i}{\hbar}[\hat{H},\hat{\sigma}_{jk}]+\hat{R}_{jk}+\hat{F}_{jk}$,
where $\hat{R}_{jk}$ and $\hat{F}_{jk}$ are the adiabatic relaxation term and
the Langevin noise, respectively. The interaction Hamiltonian $\hat{H}$ of the
double-$\Lambda$ SFWM atomic ensemble is

$\displaystyle\hat{H}=$ $\displaystyle-\frac{\hbar
N}{2L}\int_{0}^{L}dz\left(\Omega_{d}\hat{\sigma}_{31}+\Omega_{c}\hat{\sigma}_{42}+2g_{s}\hat{a}_{s}\hat{\sigma}_{32}\right.\left.+2g_{as}\hat{a}_{as}\hat{\sigma}_{41}e^{-i\Delta
kz}+\Delta_{d}\hat{\sigma}_{33}+\Delta_{c}\hat{\sigma}_{44}+{\rm
H.c.}\right),$ (S1)

where $N$ and $L$ are the total number of atoms and the atomic medium length,
respectively. The applied high-intensity coupling and driving fields,
characterized by the semiclassical Rabi frequency $\Omega_{c(d)}$, induce the
generation of Stokes and anti-Stokes photons, which are treated as quantum
fields. Here, $g_{s}=d_{32}\sqrt{\bar{\omega}_{s}/2\hbar\epsilon_{0}V}$ is the
coupling constant between the Stokes field and atoms, where $d_{32}$ indicates
the dipole moment of the $|2\rangle\leftrightarrow|3\rangle$ transition,
$\bar{\omega}_{s}$ is the center angular frequency of Stokes field,
$\epsilon_{0}$ is the vacuum permittivity, and $V$ is the interaction volume.
The coupling constant
$g_{as}=d_{41}\sqrt{\bar{\omega}_{as}/2\hbar\epsilon_{0}V}$ between the anti-
Stokes field and atoms is described in the same manner. The operator
$\hat{a}_{s(as)}$ is the slowly varying annihilation operator of the Stokes
(anti-Stokes) field. The phase shift $\Delta
kz\approx(\vec{k}_{d}-\vec{k}_{s}+\vec{k}_{c}-\vec{k}_{as})\cdot\vec{z}$
represents the phase-mismatching value in the 1D approximation. All relevant
HLEs for the double-$\Lambda$ SFWM system can be derived as follows:

$\displaystyle\frac{\partial\hat{\sigma}_{11}}{\partial t}=$ $\displaystyle
i\left[\frac{\Omega_{d}^{*}}{2}\hat{\sigma}_{13}-\frac{\Omega_{d}}{2}\hat{\sigma}_{31}+g_{as}^{*}\hat{a}_{as}^{\dagger}\hat{\sigma}_{14}e^{i\Delta
kz}-g_{as}\hat{a}_{as}\hat{\sigma}_{41}e^{-i\Delta
kz}\right]+\Gamma_{31}\hat{\sigma}_{33}+\Gamma_{41}\hat{\sigma}_{44}+\hat{F}_{11},$
(S2) $\displaystyle\frac{\partial\hat{\sigma}_{22}}{\partial t}=$
$\displaystyle
i\left[\frac{\Omega_{c}^{*}}{2}\hat{\sigma}_{24}-\frac{\Omega_{c}}{2}\hat{\sigma}_{42}+g_{s}^{*}\hat{a}_{s}^{\dagger}\hat{\sigma}_{23}-g_{s}\hat{a}_{s}\hat{\sigma}_{32}\right]+\Gamma_{32}\hat{\sigma}_{33}+\Gamma_{42}\hat{\sigma}_{44}+\hat{F}_{22},$
(S3) $\displaystyle\frac{\partial\hat{\sigma}_{33}}{\partial t}=$
$\displaystyle
i\left[\frac{\Omega_{d}}{2}\hat{\sigma}_{31}-\frac{\Omega_{d}^{*}}{2}\hat{\sigma}_{13}+g_{s}\hat{a}_{s}\hat{\sigma}_{32}-g_{s}^{*}\hat{a}_{s}^{\dagger}\hat{\sigma}_{23}\right]-(\Gamma_{31}+\Gamma_{32})\hat{\sigma}_{33}+\hat{F}_{33},$
(S4) $\displaystyle\frac{\partial\hat{\sigma}_{44}}{\partial t}=$
$\displaystyle
i\left[\frac{\Omega_{c}}{2}\hat{\sigma}_{42}-\frac{\Omega_{c}^{*}}{2}\hat{\sigma}_{24}+g_{as}\hat{a}_{as}\hat{\sigma}_{41}e^{-i\Delta
kz}-g_{as}^{*}\hat{a}_{as}^{\dagger}\hat{\sigma}_{14}e^{i\Delta
kz}\right]-(\Gamma_{41}+\Gamma_{42})\hat{\sigma}_{44}+\hat{F}_{44},$ (S5)
$\displaystyle\frac{\partial\hat{\sigma}_{13}}{\partial t}=$ $\displaystyle
i\left[\frac{\Omega_{d}}{2}(\hat{\sigma}_{11}-\hat{\sigma}_{33})-g_{as}\hat{a}_{as}\hat{\sigma}_{43}e^{-i\Delta
kz}+g_{s}\hat{a}_{s}\hat{\sigma}_{12}\right]-\frac{\gamma_{31}-2i\Delta_{d}}{2}\hat{\sigma}_{13}+\hat{F}_{13},$
(S6) $\displaystyle\frac{\partial\hat{\sigma}_{24}}{\partial t}=$
$\displaystyle
i\left[\frac{\Omega_{c}}{2}(\hat{\sigma}_{22}-\hat{\sigma}_{44})-g_{s}\hat{a}_{s}\hat{\sigma}_{34}+g_{as}\hat{a}_{as}\hat{\sigma}_{21}e^{-i\Delta
kz}\right]-\frac{\gamma_{42}-2i\Delta_{c}}{2}\hat{\sigma}_{24}+\hat{F}_{24},$
(S7) $\displaystyle\frac{\partial\hat{\sigma}_{21}}{\partial t}=$
$\displaystyle
i\left[\frac{\Omega_{d}^{*}}{2}\hat{\sigma}_{23}-\frac{\Omega_{c}}{2}\hat{\sigma}_{41}-g_{s}\hat{a}_{s}\hat{\sigma}_{31}+g_{as}^{*}\hat{a}_{as}^{\dagger}\hat{\sigma}_{24}e^{i\Delta
kz}\right]-\frac{\gamma_{21}}{2}\hat{\sigma}_{21}+\hat{F}_{21},$ (S8)
$\displaystyle\frac{\partial\hat{\sigma}_{23}}{\partial t}=$ $\displaystyle
i\left[g_{s}\hat{a}_{s}(\hat{\sigma}_{22}-\hat{\sigma}_{33})-\frac{\Omega_{c}}{2}\hat{\sigma}_{43}+\frac{\Omega_{d}}{2}\hat{\sigma}_{21}\right]-\frac{\gamma_{32}-2i\Delta_{d}}{2}\hat{\sigma}_{23}+\hat{F}_{23},$
(S9) $\displaystyle\frac{\partial\hat{\sigma}_{41}}{\partial t}=$
$\displaystyle i\left[-g_{as}^{*}\hat{a}_{as}^{\dagger}e^{i\Delta
kz}(\hat{\sigma}_{11}-\hat{\sigma}_{44})+\frac{\Omega_{d}^{*}}{2}\hat{\sigma}_{43}-\frac{\Omega_{c}^{*}}{2}\hat{\sigma}_{21}\right]-\frac{\gamma_{41}+2i\Delta_{c}}{2}\hat{\sigma}_{41}+\hat{F}_{41},$
(S10) $\displaystyle\frac{\partial\hat{\sigma}_{43}}{\partial t}=$
$\displaystyle
i\left[\frac{\Omega_{d}}{2}\hat{\sigma}_{41}-\frac{\Omega_{c}^{*}}{2}\hat{\sigma}_{23}+g_{s}\hat{a}_{s}\hat{\sigma}_{42}-g_{as}^{*}\hat{a}_{as}^{\dagger}\hat{\sigma}_{13}e^{i\Delta
kz}\right]-\frac{\gamma_{43}-2i\Delta_{d}+2i\Delta_{c}}{2}\hat{\sigma}_{43}+\hat{F}_{43}.$
(S11)

Here $\Gamma_{jk}$ is the spontaneous decay rate from the excited state
$|j\rangle$ to the ground state $|k\rangle$. $\gamma_{jk}$ is the decoherence
rate corresponding to the states $|j\rangle$ and $|k\rangle$, which is due to
the dephasing from the spontaneous emission. Of these, although both the
states $|1\rangle$ and $|2\rangle$ are ground states, the term $\gamma_{21}$
is introduced because of some effects such as atomic collisions and magnetic
inhomogeneity. In our system,
$\Gamma_{31}=\Gamma_{32}=\Gamma_{41}=\Gamma_{42}=\frac{\Gamma}{2}$, where
$\Gamma$ is the spontaneous decay rate of the rubidium-87 D2 line, and
therefore we have $\gamma_{31}=\gamma_{32}=\gamma_{41}=\gamma_{42}=\Gamma$ and
$\gamma_{43}=2\Gamma$.

To simplify the solution of the HLEs, we employed perturbation theory. By
neglecting the influence of the quantum field operator $\hat{a}_{s(as)}$ in
the HLEs, we can derive the zeroth-order steady-state collective atomic
operators in the form of
$\hat{\sigma}_{jk}^{(0)}=\langle\hat{\sigma}_{jk}^{(0)}\rangle+\sum_{mn}\epsilon_{jk}^{mn}\hat{F}_{mn}$,
with expectation values:
$\langle\hat{\sigma}_{11}^{(0)}\rangle=\frac{|\Omega_{c}|^{2}(\Gamma^{2}+4\Delta_{d}^{2})+|\Omega_{d}|^{2}|\Omega_{c}|^{2}}{M}$,
$\langle\hat{\sigma}_{22}^{(0)}\rangle=\frac{|\Omega_{d}|^{2}(\Gamma^{2}+4\Delta_{c}^{2})+|\Omega_{d}|^{2}|\Omega_{c}|^{2}}{M}$,
$\langle\hat{\sigma}_{33}^{(0)}\rangle=\langle\hat{\sigma}_{44}^{(0)}\rangle=\frac{|\Omega_{d}|^{2}|\Omega_{c}|^{2}}{M}$,
$\langle\hat{\sigma}_{13}^{(0)}\rangle=\frac{i(\Gamma+2i\Delta_{d})|\Omega_{c}|^{2}\Omega_{d}}{M}$,
and
$\langle\hat{\sigma}_{24}^{(0)}\rangle=\frac{i(\Gamma+2i\Delta_{c})|\Omega_{d}|^{2}\Omega_{c}}{M}$.
Here,
$M=|\Omega_{d}|^{2}(\Gamma^{2}+4\Delta_{c}^{2})+|\Omega_{c}|^{2}(\Gamma^{2}+4\Delta_{d}^{2})+4|\Omega_{d}|^{2}|\Omega_{c}|^{2}$.
Despite our assumption that both input fields operate in a single mode
described by steady-state zeroth-order HLEs, the solutions of the zeroth-order
collective atomic operators still exhibit time-varying terms due to the
presence of the vacuum field. These time-varying terms lead to a dynamic
equilibrium, which is reflected in the linear combination of $\hat{F}_{mn}$.
By substituting the derived zeroth-order adiabatic collective atomic operators
back into the HLEs, we can obtain the first-order HLEs as follows:

$\displaystyle\frac{\partial\hat{\sigma}_{21}^{(1)}}{\partial
t}=i\left[\frac{\Omega_{d}^{*}}{2}\hat{\sigma}_{23}^{(1)}-\frac{\Omega_{c}}{2}\hat{\sigma}_{41}^{(1)}-g_{s}\hat{a}_{s}\langle\hat{\sigma}_{31}^{(0)}\rangle+g_{as}^{*}\hat{a}_{as}^{\dagger}e^{i\Delta
kz}\langle\hat{\sigma}_{24}^{(0)}\rangle\right]-\frac{\gamma_{21}}{2}\hat{\sigma}_{21}^{(1)}+\hat{F}_{21},$
(S12) $\displaystyle\frac{\partial\hat{\sigma}_{23}^{(1)}}{\partial
t}=i\left[g_{s}\hat{a}_{s}\left(\langle\hat{\sigma}_{22}^{(0)}\rangle-\langle\hat{\sigma}_{33}^{(0)}\rangle\right)-\frac{\Omega_{c}}{2}\hat{\sigma}_{43}^{(1)}+\frac{\Omega_{d}}{2}\hat{\sigma}_{21}^{(1)}\right]-\frac{\Gamma-2i\Delta_{d}}{2}\hat{\sigma}_{23}^{(1)}+\hat{F}_{23},$
(S13) $\displaystyle\frac{\partial\hat{\sigma}_{41}^{(1)}}{\partial
t}=i\left[-g_{as}^{*}\hat{a}_{as}^{\dagger}e^{i\Delta
kz}\left(\langle\hat{\sigma}_{11}^{(0)}\rangle-\langle\hat{\sigma}_{44}^{(0)}\rangle\right)+\frac{\Omega_{d}^{*}}{2}\hat{\sigma}_{43}^{(1)}-\frac{\Omega_{c}^{*}}{2}\hat{\sigma}_{21}^{(1)}\right]-\frac{\Gamma+2i\Delta_{c}}{2}\hat{\sigma}_{41}^{(1)}+\hat{F}_{41},$
(S14) $\displaystyle\frac{\partial\hat{\sigma}_{43}^{(1)}}{\partial
t}=i\left[\frac{\Omega_{d}}{2}\hat{\sigma}_{41}^{(1)}-\frac{\Omega_{c}^{*}}{2}\hat{\sigma}_{23}^{(1)}+g_{s}\hat{a}_{s}\langle\hat{\sigma}_{42}^{(0)}\rangle-
g_{as}^{*}\hat{a}_{as}^{\dagger}e^{i\Delta
kz}\langle\hat{\sigma}_{13}^{(0)}\rangle\right]-(\Gamma-i\Delta_{d}+i\Delta_{c})\hat{\sigma}_{43}^{(1)}+\hat{F}_{43}.$
(S15)

When substituting the zeroth-order results into the first-order HLEs, we
omitted the product of the field ladder operators and the Langevin noises
$\hat{F}_{mn}$ due to their negligible effects. To investigate biphoton
spectra, we applied the Fourier transform on Eqs. (S12)–(S15). We chose the
angular frequency shift, represented as $\omega$, relative to the central
angular frequency of the Stokes field, $\bar{\omega}_{s}$, as the frequency-
domain variable in the Fourier transform. In other words,
$\widetilde{a}_{s}(z,\omega)=\frac{1}{2\pi}\int dt\hat{a}_{s}(z,t)e^{i\omega
t}\equiv\mathcal{F}\\{\hat{a}_{s}(z,t)\\}$. To maintain the factor of the
integrand as $e^{i\omega t}$, the Fourier transform of the creation operator
should be adjusted to
$\widetilde{a}_{s}^{\dagger}(z,-\omega)=\frac{1}{2\pi}\int
dt\hat{a}_{s}^{\dagger}(z,t)e^{i\omega t}$. Similarly, we have
$\widetilde{a}_{as}(z,\omega)=\mathcal{F}\\{\hat{a}_{as}(z,t)\\}$,
$\widetilde{a}_{as}^{\dagger}(z,-\omega)=\mathcal{F}\\{\hat{a}_{as}^{\dagger}(z,t)\\}$,
$\widetilde{\sigma}_{jk}^{(1)}(z,\omega)=\mathcal{F}\\{\hat{\sigma}_{jk}^{(1)}(z,t)\\}$,
and $\widetilde{F}_{jk}(z,\omega)=\mathcal{F}\\{\hat{F}_{jk}(z,t)\\}$.
Therefore, the solutions for the frequency-domain collective atomic operators
$\widetilde{\sigma}_{23}$ and $\widetilde{\sigma}_{41}$ can be expressed as
$\widetilde{\sigma}_{23(41)}=\epsilon_{23(41)}g_{s}\widetilde{a}_{s}+\eta_{23(41)}g_{as}^{*}\widetilde{a}_{as}^{\dagger}e^{i\Delta
kz}+\sum_{jk}\xi_{jk}^{s(as)}\widetilde{F}_{jk}$, where
$\\{jk\\}=\\{21,23,41,43\\}$. We combine these results with the frequency-
domain Maxwell–Schrödinger equations used to describe our experimental
backward double-$\Lambda$ SFWM system as follows:

$\displaystyle\left(-\frac{i\omega}{c}+\frac{\partial}{\partial
z}\right)\widetilde{a}_{s}=\frac{iNg_{s}^{*}}{c}\widetilde{\sigma}_{23}^{(1)},$
(S16) $\displaystyle\left(-\frac{i\omega}{c}-\frac{\partial}{\partial
z}\right)\widetilde{a}_{as}^{\dagger}=-\frac{iNg_{as}}{c}\widetilde{\sigma}_{41}^{(1)}e^{-i\Delta
kz}.$ (S17)

We can simplify the expressions by introducing a change of variables, defining
$\widetilde{a}_{as,\Delta
k}^{\dagger}\equiv\widetilde{a}_{as}^{\dagger}e^{i\Delta kz}$ to remove the
exponential phase variation. This allows us to reorganize Eqs. (S16) and (S17)
into the following matrix form:

$\displaystyle\frac{\partial}{\partial z}\begin{pmatrix}\widetilde{a}_{s}\\\
\widetilde{a}_{as,\Delta
k}^{\dagger}\end{pmatrix}=\begin{pmatrix}g_{R}&\kappa_{s}\\\
\kappa_{as}&\Gamma_{as}\end{pmatrix}\begin{pmatrix}\widetilde{a}_{s}\\\
\widetilde{a}_{as,\Delta
k}^{\dagger}\end{pmatrix}+\frac{iN}{c}\sum_{jk}\begin{pmatrix}g_{s}^{*}\xi_{jk}^{s}\\\
g_{as}\xi_{jk}^{as}\end{pmatrix}\widetilde{F}_{jk}.$ (S18)

The Raman gain coefficient, denoted as $g_{R}$, is expressed as:
$g_{R}=\frac{iN|g_{s}|^{2}}{c}\epsilon_{23}+\frac{i\omega}{c}$. The first
component, $\frac{iN|g_{s}|^{2}}{c}\epsilon_{23}$, characterizes the Raman
process, while the second component, $\frac{i\omega}{c}$, represents the phase
term arising from free evolution. The electromagnetically induced transparency
(EIT) profile coefficient, denoted as $\Gamma_{as}$, takes the form:
$\Gamma_{as}=\frac{iN|g_{s}|^{2}}{c}\eta_{41}+i\Delta k-\frac{i\omega}{c}$. In
this expression, the first component, $\frac{iN|g_{s}|^{2}}{c}\eta_{41}$,
describes the EIT characteristics, and the second term, $i\Delta k$,
represents the phase-mismatching effect. Furthermore, we have the Stokes
(anti-Stokes) coupling coefficients $\kappa_{s}$ and $\kappa_{as}$, defined
as: $\kappa_{s}=\frac{iNg_{s}^{*}g_{as}^{*}}{c}\eta_{23}$ and
$\kappa_{as}=\frac{iNg_{s}g_{as}}{c}\epsilon_{41}$. The coefficient
$\kappa_{s}$ ($\kappa_{as}$) quantifies how the anti-Stokes (Stokes) field
affects the Stokes (anti-Stokes) field.

The solution to Eq. (S18) can be obtained by directly integrating from $z=0$
to $z=L$, as illustrated below:

$\displaystyle\begin{pmatrix}\widetilde{a}_{s}(L)\\\ \widetilde{a}_{as,\Delta
k}^{\dagger}(L)\end{pmatrix}=$ $\displaystyle{\rm
exp}\left[\begin{pmatrix}g_{R}&\kappa_{s}\\\
\kappa_{as}&\Gamma_{as}\end{pmatrix}L\right]\begin{pmatrix}\widetilde{a}_{s}(0)\\\
\widetilde{a}_{as,\Delta
k}^{\dagger}(0)\end{pmatrix}+\frac{iN}{c}\sum_{jk}\int_{0}^{L}dz\,{\rm
exp}\left[\begin{pmatrix}g_{R}&\kappa_{s}\\\
\kappa_{as}&\Gamma_{as}\end{pmatrix}(L-z)\right]\begin{pmatrix}g_{s}^{*}\xi_{jk}^{s}\\\
g_{as}\xi_{jk}^{as}\end{pmatrix}\widetilde{F}_{jk}$ $\displaystyle\equiv$
$\displaystyle\begin{pmatrix}A^{\prime}&B^{\prime}\\\
C^{\prime}&D^{\prime}\end{pmatrix}\begin{pmatrix}\widetilde{a}_{s}(0)\\\
\widetilde{a}_{as,\Delta
k}^{\dagger}(0)\end{pmatrix}+\sqrt{\frac{N}{c}}\sum_{jk}\int_{0}^{L}dz\,\begin{pmatrix}P^{\prime}_{jk}\\\
Q^{\prime}_{jk}\end{pmatrix}\widetilde{F}_{jk}.$ (S19)

Since the Stokes and anti-Stokes fields are co-propagating outputs in the
backward double-$\Lambda$ SFWM system, we now reorganize Eq. (S19) to describe
this biphoton output as follows:

$\displaystyle\begin{pmatrix}\widetilde{a}_{s}(L)\\\ \widetilde{a}_{as,\Delta
k}^{\dagger}(0)\end{pmatrix}=$
$\displaystyle\begin{pmatrix}A^{\prime}-\frac{B^{\prime}C^{\prime}}{D^{\prime}}&\frac{B^{\prime}}{D^{\prime}}\\\
-\frac{C^{\prime}}{D^{\prime}}&\frac{1}{D^{\prime}}\end{pmatrix}\begin{pmatrix}\widetilde{a}_{s}(0)\\\
\widetilde{a}_{as,\Delta
k}^{\dagger}(L)\end{pmatrix}+\sqrt{\frac{N}{c}}\sum_{jk}\int_{0}^{L}dz\,\begin{pmatrix}P^{\prime}_{jk}-\frac{B^{\prime}}{D^{\prime}}Q^{\prime}_{jk}\\\
-\frac{1}{D^{\prime}}Q^{\prime}_{jk}\end{pmatrix}\widetilde{F}_{jk}$
$\displaystyle\equiv$ $\displaystyle\begin{pmatrix}A&B\\\
C&D\end{pmatrix}\begin{pmatrix}\widetilde{a}_{s}(0)\\\
\widetilde{a}_{as,\Delta
k}^{\dagger}(L)\end{pmatrix}+\sqrt{\frac{N}{c}}\sum_{jk}\int_{0}^{L}dz\,\begin{pmatrix}P_{jk}\\\
Q_{jk}\end{pmatrix}\widetilde{F}_{jk}.$ (S20)

$A$ and $D$ represent coefficients for mode preservation, elucidating the
influence of the Stokes and anti-Stokes inputs on their respective outputs. On
the other hand, $B$ and $C$ represent coefficients for mode conversion,
describing the energy conversion between the Stokes and anti-Stokes fields,
which is the process of photon mutual conversion in the double-$\Lambda$ four-
wave mixing (FWM) [2, 3].

###  I.2 Biphoton Generation Rate and Correlation Functions

The biphoton generation rate is related to the photon flux passing through the
interaction cross-section area with the atoms. The photon flux can be obtained
using the Poynting vector, denoted as
$\vec{S}=\frac{1}{\mu_{0}}\vec{E}\times\vec{B}^{*}$, which represents the
energy flux per unit area of an electromagnetic wave. By dividing the surface
integral of the Poynting vector by the energy of a single photon, and taking
into account the quantized electric field in the form of
$\vec{E}=\hat{e}\sqrt{\frac{\hbar\bar{\omega}}{2\epsilon_{0}V}}\hat{a}e^{i(\vec{k}\cdot\vec{r}-\omega
t)}+{\rm H.c.}$, we can determine the quantum photon generation rate as
follows:

$\displaystyle R=\frac{1}{\hbar\bar{\omega}}\int\vec{S}\cdot
d\vec{\sigma}=\frac{c}{L}\langle\hat{a}^{\dagger}\hat{a}\rangle+\frac{c}{2L}.$
(S21)

The parameter $\hat{e}$ is the unit vector of the polarized direction. The
term $\frac{c}{2L}$ on the right-hand side of Eq. (S21) arises from the zero-
point energy of the vacuum field and does not contribute to the actual count
of generated photons. Thus, the photon generation rate can be expressed as
$R=\frac{c}{L}\langle\hat{a}^{\dagger}\hat{a}\rangle$. By employing the
frequency-domain commutation relation
$[\widetilde{a}(\omega^{\prime}),\widetilde{a}^{\dagger}(-\omega)]=\frac{L}{2\pi
c}\delta(\omega+\omega^{\prime})$ and Einstein relation
$\langle\widetilde{F}_{jk}^{\dagger}(z,-\omega)\widetilde{F}_{j^{\prime}k^{\prime}}(z^{\prime},\omega^{\prime})\rangle=\frac{L}{2\pi
N}\mathcal{D}_{jk^{\dagger},j^{\prime}k^{\prime}}\delta(z-z^{\prime})\delta(\omega+\omega^{\prime})$,
we can derive the photon generation rates for the Stokes ($R_{s}$) and anti-
Stokes ($R_{as}$) fields, as given by Eqs. (1) and (2) presented in the main
text. In our SFWM system, it is worth noting that $g_{p}=g_{s}\equiv g$. To
conform to commonly used experimental parameters, we perform the substitution
$\frac{g^{2}N}{c}=\frac{\rm{OD}\Gamma}{4L}$, where OD denotes the optical
depth of the atomic medium. In addition to the photon generation rate, we can
also derive the normalized Glauber second-order cross-correlation function,
denoted as $g_{s\text{-}as}^{(2)}(\tau)$, as depicted in Eq. (3) in the main
text. Similarly, we can obtain the normalized autocorrelation functions
$g_{s\text{-}s}^{(2)}(\tau)$ and $g_{as\text{-}as}^{(2)}(\tau)$ for the
generated Stokes and anti-Stokes photons.

###  I.3 Diffusion Coefficients and Einstein Relation

The HLE for the $l$th atom is given by $\frac{\partial}{\partial
t}\hat{\sigma}_{jk}^{l}=\frac{i}{\hbar}[\hat{H}^{l},\hat{\sigma}_{jk}^{l}]+\hat{R}_{jk}^{l}+\hat{F}_{jk}^{l}$.
According to the Markovian approximation, the correlation between Langevin
noises follows a Dirac delta function as
$\langle\hat{F}_{jk}^{l}(t)\hat{F}_{j^{\prime}k^{\prime}}^{l}(t^{\prime})\rangle=\mathcal{D}_{jk,j^{\prime}k^{\prime}}\delta(t-t^{\prime})$,
where $\mathcal{D}_{jk,j^{\prime}k^{\prime}}$ represents the diffusion
coefficient satisfying the Einstein relation [4]. In this case, the Einstein
relation can be derived from the following equation:

$\displaystyle\frac{\partial}{\partial
t}\langle\hat{\sigma}_{jk}^{l}\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}\rangle=$
$\displaystyle\left\langle\left(\frac{i}{\hbar}[\hat{H}^{l},\hat{\sigma}_{jk}^{l}]+\hat{R}_{jk}^{l}+\hat{F}_{jk}^{l}\right)\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}\right\rangle+\left\langle\hat{\sigma}_{jk}^{l}\left(\frac{i}{\hbar}[\hat{H}^{l},\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}]+\hat{R}_{j^{\prime}k^{\prime}}^{l}+\hat{F}_{j^{\prime}k^{\prime}}^{l}\right)\right\rangle$
$\displaystyle=$
$\displaystyle\frac{i}{\hbar}\langle[\hat{H}^{l},\hat{\sigma}_{jk}^{l}\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}]\rangle+\langle\hat{R}_{jk}^{l}\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}\rangle+\langle\hat{\sigma}_{jk}^{l}\hat{R}_{j^{\prime}k^{\prime}}^{l}\rangle+\langle\hat{F}_{jk}^{l}\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}\rangle+\langle\hat{\sigma}_{jk}^{l}\hat{F}_{j^{\prime}k^{\prime}}^{l}\rangle,$
(S22)

where the term of
$\langle\hat{F}_{jk}^{l}\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}\rangle$ can be
expressed as

$\displaystyle\langle\hat{F}_{jk}^{l}(t)\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}(t)\rangle=$
$\displaystyle\langle\hat{F}_{jk}^{l}(t)\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}(t-\Delta
t)\rangle+\int_{t-\Delta
t}^{t}dt^{\prime}\left\langle\hat{F}_{jk}^{l}(t)\frac{d\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}(t^{\prime})}{dt^{\prime}}\right\rangle$
$\displaystyle=$ $\displaystyle\frac{i}{\hbar}\int_{t-\Delta
t}^{t}dt^{\prime}\langle\hat{F}_{jk}^{l}(t)[\hat{H}^{l}(t^{\prime}),\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}(t^{\prime})]\rangle+\int_{t-\Delta
t}^{t}dt^{\prime}\left\langle\hat{F}_{jk}^{l}(t)\hat{R}_{j^{\prime}k^{\prime}}^{l}(t^{\prime})\right\rangle+\int_{t-\Delta
t}^{t}dt^{\prime}\left\langle\hat{F}_{jk}^{l}(t)\hat{F}_{j^{\prime}k^{\prime}}^{l}(t^{\prime})\right\rangle.$
(S23)

The term
$\langle\hat{F}_{jk}^{l}(t)\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}(t-\Delta
t)\rangle$ yields zero as the Langevin noises at time $t$ are incapable of
influencing the atom at time $t-\Delta t$, owing to the nonanticipatory
property. Additionally, given that the terms
$[\hat{H}^{l}(t^{\prime}),\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}(t^{\prime})]$
and $\hat{R}_{j^{\prime}k^{\prime}}^{l}(t^{\prime})$ evolve at a slower rate
compared to the Langevin noises, the integrals of the first two terms in Eq.
(S23) both approach to zero. Consequently, Eq. (S23) simplifies to
$\langle\hat{F}_{jk}^{l}\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}\rangle=\mathcal{D}_{jk,j^{\prime}k^{\prime}}/2$.
A similar derivation leads to the relation
$\langle\hat{\sigma}_{jk}^{l}\hat{F}_{j^{\prime}k^{\prime}}^{l}\rangle=\mathcal{D}_{jk,j^{\prime}k^{\prime}}/2$.
Employing these relations, Eq. (S22) can be further expressed as follows:

$\displaystyle\mathcal{D}_{jk,j^{\prime}k^{\prime}}=$
$\displaystyle\frac{\partial}{\partial
t}\langle\hat{\sigma}_{jk}^{l}\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}\rangle-\frac{i}{\hbar}\langle[\hat{H}^{l},\hat{\sigma}_{jk}^{l}\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}]\rangle-\langle\hat{R}_{jk}^{l}\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}\rangle-\langle\hat{\sigma}_{jk}^{l}\hat{R}_{j^{\prime}k^{\prime}}^{l}\rangle$
$\displaystyle=$
$\displaystyle\delta_{kj^{\prime}}\langle\hat{R}_{jk^{\prime}}^{l}\rangle-\langle\hat{R}_{jk}^{l}\hat{\sigma}_{j^{\prime}k^{\prime}}^{l}\rangle-\langle\hat{\sigma}_{jk}^{l}\hat{R}_{j^{\prime}k^{\prime}}^{l}\rangle,$
(S24)

and this is commonly referred to as the Einstein relation. After establishing
these relationships, we can extend our findings from the single-atom scenario
to collective situations. In the context of an atomic ensemble system, we
define the collective Langevin noise as
$\hat{F}_{jk}=\sum_{l=1}^{N_{z}}\hat{F}_{jk}^{l}/N_{z}$, where
$N_{z}=N{dz}/L$. With this definition in place, we proceed to formulate

$\displaystyle\langle\hat{F}_{jk}(z,t)\hat{F}_{j^{\prime}k^{\prime}}(z^{\prime},t^{\prime})\rangle=\frac{1}{N_{z}^{2}}\sum_{l,l^{\prime}}\langle\hat{F}_{jk}^{l}\hat{F}_{j^{\prime}k^{\prime}}^{l^{\prime}}\rangle\delta_{ll^{\prime}}\delta_{zz^{\prime}}=\frac{1}{N_{z}^{2}}\sum_{l}\langle\hat{F}_{jk}^{l}\hat{F}_{j^{\prime}k^{\prime}}^{l}\rangle\delta_{zz^{\prime}}=\frac{\delta_{zz^{\prime}}}{N_{z}}\mathcal{D}_{jk,j^{\prime}k^{\prime}}\delta(t-t^{\prime}).$
(S25)

By using the relation $\delta_{zz^{\prime}}=\delta(z-z^{\prime}){dz}$, we can
obtain the Einstein relations in both the time and frequency domains for the
atomic ensemble system as follows:

$\displaystyle\langle\hat{F}_{jk}(z,t)\hat{F}_{j^{\prime}k^{\prime}}(z^{\prime},t^{\prime})\rangle=\frac{L}{N}\mathcal{D}_{jk,j^{\prime}k^{\prime}}\delta(z-z^{\prime})\delta(t-t^{\prime}).$
(S26)
$\displaystyle\langle\widetilde{F}_{jk}(z,\omega)\widetilde{F}_{j^{\prime}k^{\prime}}(z^{\prime},\omega^{\prime})\rangle=\frac{L}{2\pi
N}\mathcal{D}_{jk,j^{\prime}k^{\prime}}\delta(z-z^{\prime})\delta(\omega+\omega^{\prime}).$
(S27)

###  I.4 Proof of Thermal States for Stokes and Anti-Stokes Fields

We employ the reduced density operator approach to derive the density matrices
of the Stokes and anti-Stokes fields generated from the double-$\Lambda$ SFWM
system. This method transforms the obtained field operator into the form of a
density matrix [5, 6, 7]. Consider an open system with a Stokes field
$\rho^{s}(t)$, anti-Stokes field $\rho^{as}(t)$, and vacuum reservoir
$\rho^{R}(t)$. Initially, these fields are assumed to be uncorrelated, meaning
$\rho(0)=\rho^{s}(0)\otimes\rho^{as}(0)\otimes\rho^{R}(0)$. The dynamics of
this open system can be described by $\rho(t)=U(t)\rho(0)U^{\dagger}(t)$,
where $U(t)$ is the time evolution operator of the total density matrix. The
Stokes field $\rho^{s}(t)$ can be obtained by tracing over the degree of
freedom of $\rho^{as}(t)$ and $\rho^{R}(t)$ from the total density matrix,
i.e. $\rho^{s}(t)={\rm Tr}_{as,R}\\{\rho(t)\\}$. The density matrix element of
the generated Stokes field in the Fock-state basis can be given by

$\displaystyle\rho^{s}_{mn}(t)$ $\displaystyle=\langle m|{\rm
Tr}_{as,R}\left\\{U(t)\rho(0)U^{\dagger}(t)\right\\}|n\rangle$
$\displaystyle={\rm Tr}_{s}\left\\{\vphantom{\frac{A}{B}}|n_{s}\rangle\langle
m_{s}|\otimes{\rm
Tr}_{as,R}\left\\{U(t)\rho(0)U^{\dagger}(t)\right\\}\right\\}$
$\displaystyle={\rm
Tr}\left\\{\vphantom{\frac{A}{B}}\left(\vphantom{A^{B^{B}}_{C_{C}}}|n_{s}\rangle\langle
m_{s}|\otimes I_{as}\otimes I_{R}\right)U(t)\rho(0)U^{\dagger}(t)\right\\}$
$\displaystyle={\rm
Tr}\left\\{\vphantom{\frac{A}{B}}U^{\dagger}(t)\left(\vphantom{A^{B^{B}}_{C_{C}}}|n_{s}\rangle\langle
m_{s}|\otimes I_{as}\otimes I_{R}\right)U(t)\rho(0)\right\\}\equiv{\rm
Tr}\left\\{\vphantom{\frac{A}{B}}\hat{\rho}^{s}_{mn}(t)\rho(0)\right\\}.$
(S28)

Consider the vacuum input $\rho(0)=|0_{s}\rangle\langle
0_{s}|\otimes|0_{as}\rangle\langle 0_{as}|\otimes\rho^{R}(0)$ in the
double-$\Lambda$ SFWM system. By rewriting the vacuum density matrix with the
annihilation and creation operators $|0\rangle\langle
0|=\sum_{l=0}^{\infty}\frac{(-1)^{l}}{l!}[\hat{a}^{\dagger}(0)]^{l}[\hat{a}(0)]^{l}$
[8], we can derive the density matrix element of the Stokes field as follows:

$\displaystyle\rho^{s}_{mn}(t)=\frac{1}{\sqrt{m!n!}}\sum_{l=0}^{\infty}\frac{(-1)^{l}}{l!}\langle
0_{s}|[\hat{a}_{s}^{\dagger}(t)]^{n+l}[\hat{a}_{s}(t)]^{m+l}|0_{s}\rangle.$
(S29)

Using the inverse Fourier transform and combining Eq. (S20), the Stokes field
operator in the time domain can be obtained as follows:

$\displaystyle\hat{a}_{s}(L,t)=$ $\displaystyle\int d\omega e^{-i\omega
t}\left[A\widetilde{a}_{s}(0,\omega)+B\widetilde{a}_{as}^{\dagger}(L,-\omega)+\widetilde{F}_{s}(\omega)\right],$
(S30)

where $\widetilde{F}_{s}(\omega)$ is a noise-correlated term that is a linear
combination of $\widetilde{F}_{21}$, $\widetilde{F}_{23}$,
$\widetilde{F}_{41}$, and $\widetilde{F}_{43}$. After substituting Eq. (S30)
into Eq. (S29) and considering the initial vacuum input $\rho(0)$, we find
that only the condition where $m=n$ contributes to a non-zero value of the
density matrix element of the Stokes field. With these results, we can derive
the density matrix element of the Stokes field as follows:

$\displaystyle\rho^{s}_{nn}=\left\langle\sum_{l=0}^{\infty}\frac{(-1)^{l}}{l!n!}\left\\{\int
d\omega e^{-i\omega
t}\left[B^{*}\widetilde{a}_{as}(L,\omega)+\widetilde{F}_{s}^{\dagger}(-\omega)\right]\right\\}^{l+n}\left\\{\int
d\omega e^{-i\omega
t}\left[B\widetilde{a}_{as}^{\dagger}(L,-\omega)+\widetilde{F}_{s}(\omega)\right]\right\\}^{l+n}\right\rangle.$
(S31)

To further simplify Eq. (S31), we define $\hat{\zeta}\equiv\int d\omega
e^{-i\omega t}B^{*}\widetilde{a}_{as}(L,\omega)$ and $\hat{\eta}\equiv\int
d\omega e^{-i\omega t}\widetilde{F}_{s}(\omega)$, resulting in the following
expression:

$\displaystyle\rho^{s}_{nn}$
$\displaystyle=\left\langle\sum_{l=0}^{\infty}\frac{(-1)^{l}}{l!n!}\left(\hat{\zeta}+\hat{\eta}^{\dagger}\right)^{l+n}\left(\hat{\zeta}^{\dagger}+\hat{\eta}\right)^{l+n}\right\rangle$
$\displaystyle=\left\langle\sum_{l=0}^{\infty}\frac{(-1)^{l}}{l!n!}\sum_{\alpha=0}^{l+n}C^{l+n}_{\alpha}\hat{\zeta}^{\alpha}(\hat{\eta}^{\dagger})^{l+n-\alpha}\sum_{\beta=0}^{l+n}C^{l+n}_{\beta}(\hat{\zeta}^{\dagger})^{\beta}\hat{\eta}^{l+n-\beta}\right\rangle$
$\displaystyle=\sum_{l=0}^{\infty}\frac{(-1)^{l}}{l!n!}\sum_{\alpha=0}^{l+n}\sum_{\beta=0}^{l+n}C^{s+n}_{\alpha}C^{s+n}_{\beta}\left\langle\hat{\zeta}^{\alpha}(\hat{\zeta}^{\dagger})^{\beta}\right\rangle\left\langle(\hat{\eta}^{\dagger})^{l+n-\alpha}\hat{\eta}^{l+n-\beta}\right\rangle,$
(S32)

where we apply the binomial coefficient $C^{a}_{b}=\frac{a!}{b!(a-b)!}$, and
only when $\beta$ equals $\alpha$, there is a non-zero contribution, leading
to the following expression:

$\displaystyle\rho^{s}_{nn}=\sum_{l=0}^{\infty}\frac{(-1)^{l}}{l!n!}\sum_{\alpha=0}^{l+n}(C^{l+n}_{\alpha})^{2}\left\langle\hat{\zeta}^{\alpha}(\hat{\zeta}^{\dagger})^{\alpha}\right\rangle\left\langle(\hat{\eta}^{\dagger})^{l+n-\alpha}\hat{\eta}^{l+n-\alpha}\right\rangle.$
(S33)

The operator product term
$\left\langle\hat{\zeta}^{\alpha}(\hat{\zeta}^{\dagger})^{\alpha}\right\rangle$
can be simplified using the generalized Wick’s theorem. Along with the
commutation relation
$[\widetilde{a}_{as}(L,\omega),\widetilde{a}_{as}^{\dagger}(L,-\omega^{\prime})]=\frac{L}{2\pi
c}\delta(\omega+\omega^{\prime})$, this term can be expressed as:

$\displaystyle\left\langle\hat{\zeta}^{\alpha}(\hat{\zeta}^{\dagger})^{\alpha}\right\rangle=$
$\displaystyle\alpha!\left\langle\hat{\zeta}\hat{\zeta}^{\dagger}\right\rangle^{\alpha}=\alpha!\left[\frac{L}{c}\int\frac{d\omega}{2\pi}|B|^{2}\right]^{\alpha}.$
(S34)

Similarly, the operator product term
$\left\langle(\hat{\eta}^{\dagger})^{l+n-\alpha}\hat{\eta}^{l+n-\alpha}\right\rangle$
can be simplified as follows:

$\displaystyle\left\langle(\hat{\eta}^{\dagger})^{l+n-\alpha}\hat{\eta}^{l+n-\alpha}\right\rangle=(l+n-\alpha)!\left[\int\frac{d\omega}{2\pi}\langle\widetilde{F}_{s}^{\dagger}(-\omega)\widetilde{F}_{s}(\omega)\rangle\right]^{l+n-\alpha}.$
(S35)

Thus, Eq. (S32) yields

$\displaystyle\rho^{s}_{nn}$
$\displaystyle=\sum_{l=0}^{\infty}\frac{(-1)^{l}}{l!n!}\sum_{\alpha=0}^{l+n}(C^{l+n}_{\alpha})^{2}\alpha!\left[\frac{L}{c}\int\frac{d\omega}{2\pi}|B|^{2}\right]^{\alpha}(l+n-\alpha)!\left[\int\frac{d\omega}{2\pi}\langle\widetilde{F}_{s}^{\dagger}(-\omega)\widetilde{F}_{s}(\omega)\rangle\right]^{l+n-\alpha}$
$\displaystyle=\sum_{l=0}^{\infty}\frac{(-1)^{l}(l+n)!}{l!n!}\sum_{\alpha=0}^{l+n}C^{l+n}_{\alpha}\left[\frac{L}{c}\int\frac{d\omega}{2\pi}|B|^{2}\right]^{\alpha}\left[\int\frac{d\omega}{2\pi}\langle\widetilde{F}_{s}^{\dagger}(-\omega)\widetilde{F}_{s}(\omega)\rangle\right]^{l+n-\alpha}$
$\displaystyle=\sum_{l=0}^{\infty}\frac{(-1)^{l}(l+n)!}{l!n!}\left[\frac{L}{c}\int
d\omega\widetilde{R}_{s}(\omega)\right]^{l+n}=\frac{\left(\frac{L}{c}R_{s}\right)^{n}}{\left(1+\frac{L}{c}R_{s}\right)^{n+1}},$
(S36)

where the term
$\frac{L}{c}R_{s}=\langle\hat{a}^{\dagger}_{s}(L,t)\hat{a}_{s}(L,t)\rangle$
represents the average photon number of the Stokes photons $\bar{n}_{s}$.
Hence, the density matrix of the Stokes field becomes

$\displaystyle\rho^{s}=\sum_{n=0}^{\infty}\frac{\bar{n}_{s}^{n}}{\left(1+\bar{n}_{s}\right)^{n+1}}|n\rangle\langle
n|,$ (S37)

which shows the characteristic of thermal field distribution. Similarly, we
can consider the anti-Stokes field operator in the time domain to have the
following form:

$\displaystyle\hat{a}^{\dagger}_{as}(0,t)=$ $\displaystyle\int d\omega
e^{-i\omega
t}\left[C\widetilde{a}_{s}(0,\omega)+D\widetilde{a}_{as}^{\dagger}(L,-\omega)+\widetilde{F}_{as}(\omega)\right].$
(S38)

This identical structure enables us to deduce its density matrix, denoted as
$\rho^{as}$, utilizing the previously mentioned method. Consequently, we
acquire the density matrix for the anti-Stokes field as outlined below:

$\displaystyle\rho^{as}=\sum_{n=0}^{\infty}\frac{\bar{n}_{as}^{n}}{\left(1+\bar{n}_{as}\right)^{n+1}}|n\rangle\langle
n|,$ (S39)

also exhibiting the characteristics of a thermal distribution.

##  II Experimental Details

###  II.1 Coincidence Count Rate

The coincidence count rate characterizes the detection rate in the anti-Stokes
channel for a specific number of received Stokes photons. This rate is closely
related to the normalized Glauber second-order cross-correlation function,
defined as $R_{\rm C}(\tau)=R_{s}R_{as}g_{s\text{-}as}^{(2)}(\tau)\Delta T$,
where $\Delta T$ represents the time bin for detecting the Stokes photons. To
ensure a fixed value of unity for the post-selected Stokes photon, we chose a
time bin of $\Delta T=1/R_{s}$. This choice leads us to the following
expressions:

$\displaystyle R_{\rm C}(\tau)=$ $\displaystyle
R_{as}+\frac{1}{R_{s}}\left|\int\frac{d\omega}{2\pi}e^{-i\omega\tau}\left(B^{*}D+\sum_{jk,j^{\prime}k^{\prime}}\int_{0}^{L}dz\,P_{jk}^{*}\mathcal{D}_{jk^{\dagger},j^{\prime}k^{\prime}}Q_{j^{\prime}k^{\prime}}\right)\right|^{2}.$
(S40)

By selecting a total of $2^{18}$ receptions in the Stokes channel, each with a
purity $P_{s}$, we can effectively convert the coincidence count rate $R_{\rm
C}$ into coincidence counts $N_{\rm C}$ using the following relations:

$\displaystyle N_{\rm C}=$ $\displaystyle 2^{18}P_{s}\times R_{\rm
C}\eta_{as}\Delta\tau+2^{18}P_{s}\times R_{\rm
noise}^{as}\Delta\tau+2^{18}(1-P_{s})\times(R_{\rm
noise}^{as}+R_{as}\eta_{as})\Delta\tau.$ (S41)

These relations enable us to quantitatively determine the number of
coincidence counts based on the coincidence count rate from the reception
characteristics $P_{s}$ in the Stokes channel, the detection efficiency
$\eta_{as}$ in the anti-Stokes channel, and the relevant noise rate $R_{\rm
noise}^{as}$. The first term in Eq. (S41) represents the detection of
biphotons, corresponding to the registration of both Stokes and anti-Stokes
photons. This term includes the contribution of the background from the
biphoton generation rate. The second and last terms are environmental
background, which are respectively caused by the pure and impure detections in
the Stokes channel. If Stokes photons are detected, the impure registers in
anti-Stokes channel results in an undesired background; this is represented by
the second term. By contrast, an uncorrelated background is necessarily the
case if the Stokes channel receives noised photons; this is indicated by the
last term. Finally, the experimental biphoton coincidence count rate can be
expressed as $N_{\rm C}/(2^{18}P_{s}\eta_{as}\Delta\tau)\equiv R_{\rm
C}+R_{\rm env}$, where $R_{\rm env}$ represents the count rate of the
environmental background.

###  II.2 Collection Efficiency

For a collection time $T_{1}$, the total number of received Stokes photons
during $T_{1}$ is given by $R_{s}\eta_{s}T_{1}$, where $\eta_{s}$ denotes the
collection efficiency of the Stokes photons. Due to the anti-Stokes collection
efficiency $\eta_{as}$ and the pairing ratio $r_{p}$, only
$R_{s}\eta_{s}T_{1}\times\eta_{as}r_{p}$ Stokes photons are successfully
paired with corresponding anti-Stokes photons. On the other hand, when Stokes
and anti-Stokes photons from different pairs are detected, their mutual delay
times become stochastic, leading to the appearance of a uniform background
signal. Within the collection time $T_{1}$, for each individual detected
Stokes photon (total of $R_{s}\eta_{s}T_{1}$ photons), the corresponding anti-
Stokes photon detection rate is $R_{as}\eta_{as}$. Consequently, the average
background count within each time spacing $\Delta\tau$ can be expressed as
$R_{s}R_{as}\eta_{s}\eta_{as}T_{1}\Delta\tau$. We refer this $\Delta\tau$ to
the time bin of the detected anti-Stokes photons.

To ensure sufficient counts at low generation rates, we adopted an alternative
approach in our measurements, avoiding the fixed collection time $T_{1}$.
Instead, we maintained a constant count of $2^{18}$ in the Stokes channel to
ensure an ample number of detections. As a result, the correlated coincidence
counts and uncorrelated background counts were adjusted to
$2^{18}P_{s}\eta_{as}r_{p}$ and $2^{18}P_{s}R_{as}\eta_{as}\Delta\tau$,
respectively. Here, $P_{s}$ represents the purity of Stokes detection and is
calculated as $P_{s}=\frac{R_{s}\eta_{s}}{R_{s}\eta_{s}+R_{\rm noise}^{s}}$,
where $R_{\rm noise}^{s}$ denotes the count rate of the Stokes channel
originating from laser leakage, environmental photons, and dark counts from
the SPCM. This parameter was determined by conducting single-channel detection
on the Stokes channel over a time interval $T_{2}$, resulting in a count of
$R_{s}\eta_{s}T_{2}+R_{\rm noise}^{s}T_{2}$. Furthermore, by theoretically
calculating the Stokes generation rate of this process, we obtained the Stokes
collection efficiency $\eta_{s}$. Similarly, the collection efficiency
$\eta_{as}$ for the anti-Stokes channel was determined using the same
methodology, yielding a count of $R_{as}\eta_{as}T_{2}+R_{\rm
noise}^{as}T_{2}$ within the time interval $T_{2}$, where $R_{\rm noise}^{as}$
represents the noise count rate of the anti-Stokes channel. By incorporating
these experimentally derived collection efficiencies $\eta_{s}$ and
$\eta_{as}$ into the correlated coincidence counts and uncorrelated background
counts, respectively, we were able to obtain the experimental pairing ratio
and biphoton generation rate.

The collection efficiencies obtained from our experiments using guiding beams
differ from those acquired by measuring the Stokes and anti-Stokes channels in
conjunction with the theoretical biphoton generation rate. Specifically, the
collection efficiencies obtained from the guiding beams for the Stokes and
anti-Stokes fields were 9$\%$ and 6$\%$, respectively, while the corresponding
values measured from the biphoton experiments were approximately 2$\%$ and
1$\%$, respectively. For convenience, we refer to the method of determining
collection efficiency using the guiding beam as the ”guiding method.” On the
other hand, the approach that relies on the measured photon counts from
biphoton experiments and the theoretical generation rate is termed the
”experimental method.” The disparity between these two approaches may arise
from differences in the propagation of the guiding beam and the biphoton
field, potentially leading to variations in the coupling efficiency with
etalons and optical fibers. In this research, we consider the experimental
method to be more suitable. This choice is reinforced by the alignment between
our experimental results and the theoretical predictions for the signal-to-
background ratio and biphoton wavepacket characteristics. Conversely, using
the collection efficiencies of 9$\%$ and 6$\%$ obtained from the guiding
method would lead to deviations from the theoretically predicted generation
rates. This, in turn, would result in differences between the theoretical
signal-to-background ratios and biphoton wavepacket characteristics and their
experimental counterparts.

###  II.3 Pairing Ratio and Heralding Efficiency

After determining the collection efficiencies $\eta_{s}$ and $\eta_{as}$, we
divided Eq. (S41) by a factor of $2^{18}P_{s}\eta_{as}\Delta\tau$ to transform
the histogram of coincidence counts $N_{\rm C}$ into the experimental
coincidence count rate $R_{\rm C}+R_{\rm env}$. This transformation aims to
correct for the impact of collection efficiencies, thereby exploring the
intrinsic characteristics of biphotons. Specifically, under this
transformation, the photon number of the Stokes field was post-selected to be
one. Therefore, the area of the region temporally correlated with the post-
selected Stokes photon exactly represents the experimental pairing ratio. This
also indicates that, under the same biphoton generation rate, an increase in
the pairing ratio enlarges the area of the temporally correlated region while
maintaining the same background coincidence count rate. In other words, with a
fixed biphoton wavepacket bandwidth, an increase in the pairing ratio enhances
the signal-to-background ratio, thereby improving the purity of the heralded
single photon.

In contrast to the pairing ratio, the heralding efficiency addresses the
practical application of collection efficiency [9]. It focuses on measuring
anti-Stokes photons and determining the proportion heralded by Stokes photons,
essentially representing the collection efficiency $\eta_{s}$ of the Stokes
channel [10]. For instance, assuming a biphoton generation rate of $R_{\rm B}$
and a pairing ratio of one ($r_{p}=1$), the detection rates on both sides,
$R_{1}=R_{\rm B}\eta_{s}$ and $R_{2}=R_{\rm B}\eta_{as}$, can be obtained by
separately measuring each channel. Through coincidence detection experiments,
the temporally correlated biphoton detection rate $R_{3}=R_{\rm
B}\eta_{s}\eta_{as}$ can be obtained. These temporally correlated photons
occupy a fraction of the detected photons in the anti-Stokes channel with a
proportion of $\eta_{s}$. Specifically, when detecting the anti-Stokes
channel, only the rate of $R_{3}$ is temporally correlated, representing a
single-photon Fock state, while the remaining rate of $R_{2}-R_{3}$
corresponds to a thermal state.

###  II.4 Biphoton Temporal Profile

We introduce an alternative method to characterize the generated biphoton by
directly calculating the biphoton wave function [11]. The SFWM is a process of
third-order susceptibility $\chi^{(3)}$:
$H_{I}=\epsilon_{0}A\int_{-L/2}^{L/2}dz\chi^{(3)}E_{d}^{(+)}E_{c}^{(+)}\hat{E}_{s}^{(-)}\hat{E}_{as}^{(-)}+{\rm
H.c.}$, where $E_{d}^{(+)}=E_{d}e^{i(k_{d}z-\omega_{d}t)}$ and
$E_{c}^{(+)}=E_{c}e^{-i(k_{c}z+\omega_{c}t)}$ are positive-frequency driving
and coupling electric fields, respectively. The generated Stokes field is
treated using the method of quantized field with the form
$\hat{E}_{s}^{(-)}=\sqrt{\frac{\hbar\bar{\omega}_{s}}{2\epsilon_{0}V}}e^{-i(k_{s}z-\bar{\omega}_{s}t)}\int
d\omega_{1}\widetilde{a}_{s}^{\dagger}(-\omega_{1})e^{-i\left[\frac{k_{s}\chi_{s}^{*}(-\omega_{1})}{2}z+\omega_{1}t\right]}$,
whereas the anti-Stokes field is described with the form
$\hat{E}_{as}^{(-)}=\sqrt{\frac{\hbar\bar{\omega}_{as}}{2\epsilon_{0}V}}e^{i(k_{as}z+\bar{\omega}_{as}t)}\int
d\omega_{2}\widetilde{a}_{as}^{\dagger}(-\omega_{2})e^{i\left[\frac{k_{as}\chi_{as}^{*}(-\omega_{2})}{2}z-\omega_{2}t\right]}$.
In our calculations, we have supposed
$k(\omega)=k_{0}\sqrt{1+\chi(\omega)}\approx k_{0}+\frac{\chi(\omega)}{2}$.
Upon substituting these electric fields into the interaction Hamiltonian, it
can be obtained that

$\displaystyle H_{I}=$
$\displaystyle\frac{\hbar\sqrt{\bar{\omega}_{s}\bar{\omega}_{as}}}{2L}E_{d}E_{c}\int
d\omega_{1}\int d\omega_{2}\int_{-L/2}^{L/2}dz\,e^{i\Delta
kz}e^{-\frac{i}{2}\left[k_{s}\chi_{s}^{*}(-\omega_{1})-k_{as}\chi_{as}^{*}(-\omega_{2})\right]z}e^{-i(\omega_{1}+\omega_{2})t}\chi^{(3)}(\omega_{1},\omega_{2})\widetilde{a}_{s}^{\dagger}(-\omega_{1})\widetilde{a}_{as}^{\dagger}(-\omega_{2})$
$\displaystyle+{\rm H.c.}$ $\displaystyle=$
$\displaystyle\frac{\hbar\sqrt{\bar{\omega}_{s}\bar{\omega}_{as}}}{2L}E_{d}E_{c}\int
d\omega_{1}\int
d\omega_{2}\int_{-L/2}^{L/2}dz\,e^{i\kappa(\omega_{1},\omega_{2})z}e^{-i(\omega_{1}+\omega_{2})t}\chi^{(3)}(\omega_{1},\omega_{2})\widetilde{a}_{s}^{\dagger}(-\omega_{1})\widetilde{a}_{as}^{\dagger}(-\omega_{2})+{\rm
H.c.}$ $\displaystyle=$
$\displaystyle\frac{\hbar\sqrt{\bar{\omega}_{s}\bar{\omega}_{as}}}{2}E_{d}E_{c}\int
d\omega_{1}\int d\omega_{2}\,{\rm
sinc}\left[\frac{\kappa(\omega_{1},\omega_{2})L}{2}\right]e^{-i(\omega_{1}+\omega_{2})t}\chi^{(3)}(\omega_{1},\omega_{2})\widetilde{a}_{s}^{\dagger}(-\omega_{1})\widetilde{a}_{as}^{\dagger}(-\omega_{2})+{\rm
H.c.},$ (S42)

where we define $\kappa(\omega_{1},\omega_{2})=\Delta
k-\frac{k_{s}\chi_{s}^{*}(-\omega_{1})}{2}+\frac{k_{as}\chi_{as}^{*}(-\omega_{2})}{2}$.
The biphoton wave function $|\psi\rangle$ is the solution of the Schrödinger
equation: $\frac{\partial}{\partial
t}|\psi\rangle=\frac{1}{i\hbar}H_{I}|\psi\rangle$, which yields the results of

$\displaystyle|\psi(t)\rangle=$
$\displaystyle|0\rangle+\frac{1}{i\hbar}\int_{-\infty}^{t}dt^{\prime}H_{I}(t^{\prime})|\psi(t^{\prime})\rangle$
$\displaystyle=$
$\displaystyle|0\rangle+\frac{1}{i\hbar}\int_{-\infty}^{t}dt^{\prime}H_{I}(t^{\prime})|0\rangle+\frac{1}{(i\hbar)^{2}}\int_{-\infty}^{t}dt^{\prime}\int_{-\infty}^{t^{\prime}}dt^{\prime\prime}H_{I}(t^{\prime})H_{I}(t^{\prime\prime})|\psi(t^{\prime\prime})\rangle.$
(S43)

The first term in Eq. (S43) can be safely neglected, as the vacuum field is
inherently undetectable. Moreover, we can disregard the integrand
$H_{I}(t^{\prime})H_{I}(t^{\prime\prime})$ in the last terms, given its
minimal impact attributed to $\chi^{(3)}$. Consequently, the expression for
the steady-state biphoton wave function $|\psi(t\rightarrow\infty)\rangle$ can
be stated as follows:

$\displaystyle|\psi(t\rightarrow\infty)\rangle=$
$\displaystyle-\frac{iE_{d}E_{c}\sqrt{\bar{\omega}_{s}\bar{\omega}_{as}}}{2}\int_{-\infty}^{\infty}dt^{\prime}\int
d\omega_{1}\int d\omega_{2}\,{\rm
sinc}\left[\frac{\kappa(\omega_{1},\omega_{2})L}{2}\right]e^{-i(\omega_{1}+\omega_{2})t}\chi^{(3)}(\omega_{1},\omega_{2})\widetilde{a}_{s}^{\dagger}(-\omega_{1})\widetilde{a}_{as}^{\dagger}(-\omega_{2})|0\rangle$
$\displaystyle=$ $\displaystyle-i\pi
E_{d}E_{c}\sqrt{\bar{\omega}_{s}\bar{\omega}_{as}}\int d\omega_{1}\,{\rm
sinc}\left[\frac{\kappa(\omega_{1},-\omega_{1})L}{2}\right]\chi^{(3)}(\omega_{1},-\omega_{1})\widetilde{a}_{s}^{\dagger}(-\omega_{1})\widetilde{a}_{as}^{\dagger}(\omega_{1})|0\rangle.$
(S44)

It is worth noting that in the calculation of this biphoton wave function,
$\widetilde{a}_{s}^{\dagger}$ and $\widetilde{a}_{as}^{\dagger}$ appear
together. This implies that this method can only calculate correlated photons
participating in FWM and, as a result, cannot directly provide information
about the pairing ratio $r_{p}$.

By utilizing the steady-state biphoton wave function, we can calculate the
Glauber second-order correlation function as
$G^{(2)}_{s\text{-}as}(\tau)=\langle\psi(t\rightarrow\infty)|\hat{a}_{s}^{\dagger}(t)\hat{a}_{as}^{\dagger}(t+\tau)\hat{a}_{as}(t+\tau)\hat{a}_{s}(t)|\psi(t\rightarrow\infty)\rangle$.
Since the biphoton wave function is associated with
$\widetilde{a}_{s}^{\dagger}\widetilde{a}_{as}^{\dagger}|0\rangle$, we can
insert $|0\rangle\langle 0|$ between $\hat{a}_{as}^{\dagger}(t+\tau)$ and
$\hat{a}_{as}(t+\tau)$, where $|0\rangle$ denotes the vacuum state. This
yields $G^{(2)}_{s\text{-}as}(\tau)=|\langle
0|\hat{a}_{as}(t+\tau)\hat{a}_{s}(t)|\psi(t\rightarrow\infty)\rangle|^{2}$.
The term $\langle
0|\hat{a}_{as}(t+\tau)\hat{a}_{s}(t)|\psi(t\rightarrow\infty)\rangle$ can be
expressed as follows:

$\displaystyle\langle
0|\hat{a}_{as}(t+\tau)\hat{a}_{s}(t)|\psi(t\rightarrow\infty)\rangle$
$\displaystyle=-i\pi E_{d}E_{c}\sqrt{\bar{\omega}_{s}\bar{\omega}_{as}}\int
d\omega_{1}\,{\rm
sinc}\left[\frac{\kappa(\omega_{1},-\omega_{1})L}{2}\right]\chi^{(3)}(\omega_{1},-\omega_{1})\langle
0|\hat{a}_{as}(t+\tau)\hat{a}_{s}(t)\widetilde{a}_{s}^{\dagger}(-\omega_{1})\widetilde{a}_{as}^{\dagger}(\omega_{1})|0\rangle$
$\displaystyle=-i\pi E_{d}E_{c}\sqrt{\bar{\omega}_{s}\bar{\omega}_{as}}\int
d\omega_{1}\,{\rm
sinc}\left[\frac{\kappa(\omega_{1},-\omega_{1})L}{2}\right]\chi^{(3)}(\omega_{1},-\omega_{1})e^{i(k_{s}\frac{L}{2}-\bar{\omega}_{s}t)}e^{i(k_{as}\frac{L}{2}-\bar{\omega}_{as}t-\bar{\omega}_{as}\tau)}$
$\displaystyle\quad\times\int
d\omega_{3}e^{i\left[\frac{k_{s}\chi_{s}(\omega_{3})}{2}\frac{L}{2}-\omega_{3}t\right]}\int
d\omega_{4}e^{i\left[\frac{k_{as}\chi_{as}(\omega_{4})}{2}\frac{L}{2}-\omega_{4}t-\omega_{4}\tau\right]}\langle
0|\widetilde{a}_{as}(\omega_{4})\widetilde{a}_{s}(\omega_{3})\widetilde{a}_{s}^{\dagger}(-\omega_{1})\widetilde{a}_{as}^{\dagger}(\omega_{1})|0\rangle$
$\displaystyle=-i\pi
E_{d}E_{c}\sqrt{\bar{\omega}_{s}\bar{\omega}_{as}}\frac{L^{2}}{4\pi^{2}c^{2}}e^{i(k_{s}\frac{L}{2}-\bar{\omega}_{s}t)}e^{i(k_{as}\frac{L}{2}-\bar{\omega}_{as}t-\bar{\omega}_{as}\tau)}$
$\displaystyle\quad\times\int d\omega_{1}\,{\rm
sinc}\left[\frac{\kappa(\omega_{1},-\omega_{1})L}{2}\right]\chi^{(3)}(\omega_{1},-\omega_{1})e^{\frac{i}{2}\left[k_{s}\chi_{s}(-\omega_{1})+k_{as}\chi_{as}(\omega_{1})\right]\frac{L}{2}}e^{-i\omega_{1}\tau}$
$\displaystyle=\frac{-iE_{d}E_{c}\sqrt{\bar{\omega}_{s}\bar{\omega}_{as}}}{2}\frac{L^{2}}{c^{2}}e^{i\phi(t)}\int\frac{d\omega_{1}}{2\pi}e^{-i\omega_{1}\tau}{\rm
sinc}\left[\frac{\kappa(\omega_{1},-\omega_{1})L}{2}\right]\chi^{(3)}(\omega_{1},-\omega_{1})e^{\frac{i}{2}\left[k_{s}\chi_{s}(-\omega_{1})+k_{as}\chi_{as}(\omega_{1})\right]\frac{L}{2}},$
(S45)

where
$\phi(t)=\frac{k_{s}+k_{as}}{2}L-\bar{\omega}_{as}\tau-(\bar{\omega}_{s}+\bar{\omega}_{as})t$.
By employing the semiclassical model and assuming a large driving detuning
$\Delta_{d}$, we derive the forms of susceptibilities as follows:

$\displaystyle\chi_{s}(\omega)=$
$\displaystyle\frac{n|d_{32}|^{2}}{\epsilon_{0}\hbar}\frac{|\Omega_{d}|^{2}}{\Delta_{d}^{2}}\frac{\omega-i\Gamma/2}{|\Omega_{c}|^{2}-4(\omega-i\gamma_{21}/2)(\omega-i\Gamma/2)},$
(S46) $\displaystyle\chi_{as}(\omega)=$
$\displaystyle\frac{n|d_{41}|^{2}}{\epsilon_{0}\hbar}\frac{4(\omega-\Delta_{c}+i\gamma_{21}/2)}{|\Omega_{c}|^{2}-4(\omega-\Delta_{c}+i\gamma_{21}/2)(\omega+i\Gamma/2)},$
(S47) $\displaystyle\chi^{(3)}(\omega)=$
$\displaystyle\frac{nd_{14}d_{23}d_{31}d_{42}}{\epsilon_{0}\hbar^{3}}\frac{1}{\Delta_{d}}\frac{4}{|\Omega_{c}|^{2}-4(\omega-\Delta_{c}+i\gamma_{21}/2)(\omega+i\Gamma/2)}.$
(S48)

On substituting the $\chi^{(3)}$ into $G_{s\text{-}as}^{(2)}(\tau)$, we have

$\displaystyle
G_{s\text{-}as}^{(2)}(\tau)=A_{0}\left|\int\frac{d\omega}{2\pi}e^{-i\omega\tau}{\rm
sinc}\left[\frac{\kappa(\omega,-\omega)L}{2}\right]A_{1}(\omega)e^{\frac{i}{2}\left[k_{s}\chi_{s}(-\omega)+k_{as}\chi_{as}(\omega)\right]\frac{L}{2}}\right|^{2},$
(S49)

where
$A_{0}=\frac{L^{4}n^{2}|d_{14}|^{2}|d_{23}|^{2}\bar{\omega}_{s}\bar{\omega}_{as}}{4c^{4}\epsilon_{0}^{2}\hbar^{2}}=\frac{L^{2}}{c^{2}}\frac{\Gamma^{2}{\rm
OD}^{2}}{16}$ and
$A_{1}(\omega)=\frac{\Omega_{d}}{\Delta_{d}}\frac{\Omega_{c}}{|\Omega_{c}|^{2}-4(\omega-\Delta_{c}+i\gamma_{21}/2)(\omega+i\Gamma/2)}$.
In the case of a large $\Delta_{d}$, i.e., $|\chi_{s}|\ll|\chi_{as}|$, we can
further simplify $G_{s\text{-}as}^{(2)}(\tau)$ as follows:

$\displaystyle
G_{s\text{-}as}^{(2)}(\tau)=\frac{L^{2}}{c^{2}}\left|\frac{\Gamma{\rm
OD}}{4}\int\frac{d\omega}{2\pi}e^{-i\omega\tau}{\rm sinc}\left[\frac{\Delta
kL}{2}+\frac{k_{as}L\chi_{as}^{*}(-\omega)}{4}\right]A_{1}(\omega)e^{i\frac{k_{as}L\chi_{as}(\omega)}{4}}\right|^{2}.$
(S50)

The coincidence count rate can be determined using the expression $R_{\rm
C}(\tau,\Delta T)=\frac{c^{2}}{L^{2}}G_{s\text{-}as}^{(2)}(\tau)\Delta T$.
Here, we use a time bin of $\Delta T=1$ second to calculate the correlated
biphoton generation rate by integrating $R_{\rm C}(\tau)$, resulting in
$R_{\rm B}r_{p}$ (i.e., $\int d\tau R_{\rm C}(\tau,1\,{\rm s})=R_{\rm
B}r_{p}$). However, the formula in Eq. (S50) does not account for accidental
receptions, which can result from biphotons originating from different pairs
or those that remain uncorrelated.

After obtaining Eq. (S50), we proceed to explore the regime of damped Rabi
oscillations [11]. This regime tends to dominate in cases of low OD or when
the coupling field strength $\Omega_{c}$ is large. In this scenario, the term
of ${\rm sinc}\left[\frac{\Delta
kL}{2}+\frac{k_{as}L\chi_{as}^{*}(-\omega)}{4}\right]e^{i\frac{k_{as}L\chi_{as}(\omega)}{4}}$
in the integrand of Eq. (S50) can be approximated as ${\rm sinc}(\Delta
kL/2)$. To further simplify the expression, we assume $\gamma_{21}=0\Gamma$,
leading to the following form:

$\displaystyle G_{s\text{-}as}^{(2)}(\tau)=$
$\displaystyle\frac{L^{2}}{c^{2}}\left|\frac{\Gamma\Omega_{d}\Omega_{c}{\rm
OD}}{4\Delta_{d}}{\rm sinc}\left(\frac{\Delta
kL}{2}\right)\right|^{2}\left|\int\frac{d\omega}{2\pi}e^{-i\omega\tau}\frac{1}{|\Omega_{c}|^{2}-4(\omega-\Delta_{c})(\omega+i\Gamma/2)}\right|^{2}$
$\displaystyle=$
$\displaystyle\frac{L^{2}}{c^{2}}\left|\frac{\Gamma\Omega_{d}\Omega_{c}{\rm
OD}}{16\Delta_{d}}{\rm sinc}\left(\frac{\Delta
kL}{2}\right)\right|^{2}\left|\int\frac{d\omega}{2\pi}e^{-i\omega\tau}\frac{1}{(\omega+i\gamma_{e}/4-\Omega_{e}/2)(\omega+i\gamma_{e}/4+\Omega_{e}/2)}\right|^{2}$
$\displaystyle=$
$\displaystyle\frac{L^{2}}{c^{2}}\left|\frac{\Gamma\Omega_{d}\Omega_{c}{\rm
OD}}{8\Omega_{e}\Delta_{d}}{\rm sinc}\left(\frac{\Delta
kL}{2}\right)\right|^{2}e^{-\frac{\Gamma}{2}\tau}\left|{\rm
sin}\left({\frac{\Omega_{e}}{2}\tau}\right)\right|^{2}.$ (S51)

Here, we introduce the parameters $\gamma_{e}=\Gamma+i\Delta_{c}$ and
$\Omega_{e}=\sqrt{|\Omega_{c}|^{2}+\Delta_{c}^{2}-\Gamma^{2}/4+i\Gamma\Delta_{c}}$.
Subsequently, we further expand $\Omega_{e}$ as $\alpha+i\beta$, where both
$\alpha$ and $\beta$ are real. As a result, the term $\left|{\rm
sin}(\Omega_{e}\tau/2)\right|^{2}$ can be expressed as $[{\rm
sin}^{2}(\alpha\tau/2)+{\rm sinh}^{2}(\beta\tau/2)]$, where ${\rm
sinh}x=\frac{e^{x}-e^{-x}}{2}$. Because the term of ${\rm
sin}^{2}(\alpha\tau/2)$ represents the oscillating behavior, the decay
processes can be described by
$e^{-\frac{\Gamma}{2}\tau}\frac{(e^{\beta\tau/2}-e^{-\beta\tau/2})^{2}}{4}$.
Furthermore, given that $e^{-\beta\tau/2}$ decays more severely, the overall
decay behaviors depend on the term $e^{-(\frac{\Gamma}{2}-\beta)\tau}$, which
yields the $1/e$ decay time of $1/(\frac{\Gamma}{2}-\beta)$. Keeping other
parameters constant, increasing $\Delta_{c}$ extends the tail of the biphoton
wavepacket. This is because the heightened $\Delta_{c}$ weakens the coupling
between the coupling field and the atoms, making it more challenging for the
coupling field to efficiently convert the pre-established spinwave excitation
$\widetilde{\sigma}_{21}$ into anti-Stokes photons. Under the condition of
$\Delta_{c}=0$, the $1/e$ decay time and oscillating period are $2/\Gamma$ and
$2\pi/\sqrt{|\Omega_{c}|^{2}-\Gamma^{2}/4}$, respectively.

Figure S1:  EIT slow light effect. Parameters: ${\rm OD}=15$,
$\Omega_{c}=1\Gamma$, $\gamma_{21}=0\Gamma$. (a) EIT delay time $\tau_{\rm
EIT}$ as a function of detuning $\omega$. $\tau_{\rm EIT}$ is given by Eq.
(S53). The solid black, red, and blue curves represent cases with
$\Delta_{c}=0\Gamma$, $\Delta_{c}=1\Gamma$, and $\Delta_{c}=3\Gamma$,
respectively. The magenta dashed line indicates a delay time of 398 ns
calculated using $\tau_{\rm EIT}=\Gamma\,{\rm OD}/|\Omega_{c}|^{2}$. (b)
Temporal profiles passing through the EIT medium. The black dashed curve
represents the input Gaussian pulse with a $1/e^{2}$-full width of 400 ns. The
solid black, red, and blue curves correspond to cases with
$\Delta_{c}=0\Gamma$, $\Delta_{c}=1\Gamma$, and $\Delta_{c}=3\Gamma$,
respectively.

The group delay time in the steady-state condition can be obtained from
$\tau_{\rm EIT}=\frac{L}{{\rm
Re}[d\omega/dk]}-\frac{L}{c}=L\frac{d}{d\omega}{\rm
Re}[\frac{\omega}{c}\sqrt{1+\chi}]-\frac{L}{c}=\frac{L}{c}{\rm
Re}[\sqrt{1+\chi}]+\frac{L}{c}\omega\frac{d}{d\omega}{\rm
Re}[\sqrt{1+\chi}]-\frac{L}{c}$. The real part of $\sqrt{1+\chi}$ can be dealt
with the following method:

$\displaystyle{\rm Re}[\sqrt{1+\chi}]=$ $\displaystyle{\rm
Re}\left\\{\sqrt{1+{\rm Re}[\chi]+i{\rm Im}[\chi]}\right\\}$ $\displaystyle=$
$\displaystyle{\rm Re}\left\\{\sqrt{\sqrt{(1+{\rm Re}[\chi])^{2}+({\rm
Im}[\chi])^{2}}\left[\frac{1+{\rm Re}[\chi]}{\sqrt{(1+{\rm
Re}[\chi])^{2}+({\rm Im}[\chi])^{2}}}+i\frac{{\rm Im}[\chi]}{\sqrt{(1+{\rm
Re}[\chi])^{2}+({\rm Im}[\chi])^{2}}}\right]}\right\\}$ $\displaystyle\equiv$
$\displaystyle{\rm Re}\left\\{\left[(1+{\rm Re}[\chi])^{2}+({\rm
Im}[\chi])^{2}\right]^{\frac{1}{4}}\sqrt{e^{i\theta}}\right\\}$
$\displaystyle=$ $\displaystyle\left[(1+{\rm Re}[\chi])^{2}+({\rm
Im}[\chi])^{2}\right]^{\frac{1}{4}}{\rm
cos}(\frac{\theta}{2})=\sqrt{\frac{\sqrt{(1+{\rm Re}[\chi])^{2}+({\rm
Im}[\chi])^{2}}}{2}+\frac{1+{\rm Re}[\chi]}{2}}.$ (S52)

Here, we aim to analyze the propagation characteristics of the anti-Stokes
field. In this context, $\omega$ represents the angular frequency denoted as
$\omega_{as}$. Consequently, the group delay time can be expressed as:

$\displaystyle\tau_{\rm EIT}(\omega)=\frac{L}{c}{\rm
Re}[\sqrt{1+\chi_{as}(\omega)}]+k_{as}L\frac{d}{d\omega}{\rm
Re}[\sqrt{1+\chi_{as}(\omega)}]-\frac{L}{c}.$ (S53)

We consider the specific scenario where $\Delta_{c}=0\Gamma$ and
$\omega=0\Gamma$, allowing us to derive the EIT group delay time $\tau_{\rm
EIT}=\frac{\Gamma{\rm OD}}{|\Omega_{c}|^{2}}$ [12]. It is important to note
that under two-photon resonance conditions ($\omega=\Delta_{c}$), $\tau_{\rm
EIT}$ remains constant regardless of variations in $\Delta_{c}$. Therefore,
relying solely on Eq. (S53), we cannot fully elucidate the influence of
$\Delta_{c}$ on the delay time.

Figure S1(a) shows $\tau_{\rm EIT}$ corresponding to different values of
$\Delta_{c}$ and $\omega$. It is evident from the figure that, regardless of
the value of $\Delta_{c}$, the EIT group delay time is equal to
$\frac{\Gamma\,{\rm OD}}{|\Omega_{c}|^{2}}$ when the condition
$\omega=\Delta_{c}$ is met, as mentioned above. Note that for larger
$\Delta_{c}$ values near $\omega=\Delta_{c}$, there is a significant variation
in $\tau_{\rm EIT}$, which can lead to significant waveform distortion.
Conversely, in the vicinity of $\omega=\Delta_{c}=0\Gamma$, there is a
relatively flat region of $\tau_{\rm EIT}$, allowing waveforms to propagate
through the medium without severely distortion. Figure S1(b) illustrates the
transmission of a Gaussian pulse through the EIT medium. The Gaussian pulse
has a $1/e^{2}$ full width of 400 ns. Although Eq. (S53) states that
$\tau_{\rm EIT}$ remains constant at the two-photon resonance condition, the
pulse behaviors after passing through the EIT medium are noticeably different
for various $\Delta_{c}$ values. It is observed that the tail of the EIT pulse
decreases as $\Delta_{c}$ increases, especially in the case of
$\Delta_{c}=3\Gamma$. This observation further supports the claim made about
the significant variation of $\tau_{\rm EIT}$ in Fig. S1(a).

By analyzing the damped Rabi oscillation and the EIT group delay, we can gain
insights into the behaviors depicted in Figs. 2(c) and 2(d) of the main text.
In Fig. 2(c) of the main text, we introduce $\Delta_{c}=1\Gamma$, which leads
to a reduction in $\tau_{\rm EIT}$, as indicated by the red curve in Fig.
S1(b). Increasing $\Delta_{c}$ to $3\Gamma$, as demonstrated in Fig. 2(d) of
the main text, further decreases $\tau_{\rm EIT}$ due to the reduced effective
OD. However, introducing $\Delta_{c}=3\Gamma$ simultaneously extends the $1/e$
decay time $1/(\frac{\Gamma}{2}-\beta)$ of the damped Rabi oscillation. As a
result, the elongated tail in Fig. 2(d) of the main text is attributed to the
regime of damped Rabi oscillation.

###  II.5 Experimental Parameters of Figure 4

The driving Rabi frequency $\Omega_{d}$ in Fig. 4(c) and 4(d) of the main text
is 3$\Gamma$. The other parameters are listed in Table S1.

OD |  $\Omega_{c}$ $(\Gamma)$ |  $\Delta_{d}$ $(\Gamma)$
---|---|---
20 | 4.0 | 5.0
40 | 5.0 | 8.1
60 | 6.2 | 10.2
80 | 7.2 | 12.0
100 | 8.0 | 13.5
120 | 8.8 | 14.9
Table S1: Parameters in Fig. 4(c) and 4(d) of the main text.

## References

  * [1] P. Kolchin, Electromagnetically-induced-transparency-based paired photon generation, Phys. Rev. A 75, 033814 (2007).
  * [2] C.-Y. Cheng, J.-J. Lee, Z.-Y. Liu, J.-S. Shiu, and Y.-F. Chen, Quantum frequency conversion based on resonant four-wave mixing, Phys. Rev. A 103, 023711 (2021).
  * [3] Z.-Y. Liu, J.-S. Shiu, C.-Y. Cheng, and Y.-F. Chen, Controlling frequency-domain Hong-Ou-Mandel interference via electromagnetically induced transparency, Phys. Rev. A 108, 013702 (2023).
  * [4] J. C. Garrison and R. Y. Chiao, Quantum Optics (Oxford University Press, New York, 2008).
  * [5] M. O. Scully and M. S. Zubairy, Quantum Optics (Cambridge University Press, Cambridge, England, 1997).
  * [6] F. Kheirandish, Quantum dynamics of a driven damped harmonic oscillator in Heisenberg picture: exact results and possible generalizations, Eur. Phys. J. Plus 135, 243 (2020).
  * [7] H. Hsu, C.-Y. Cheng, J.-S. Shiu, L.-C. Chen, and Y.-F. Chen, Quantum fidelity of electromagnetically induced transparency: the full quantum theory, Opt. Express 30, 2097 (2022).
  * [8] W. H. Louisell, Quantum statistical properties of radiation (John Wiley and Sons, Inc., New York, 1973).
  * [9] R. H. Hadfield, Single-photon detectors for optical quantum information applications, Nat. photonics 3, 696 (2009).
  * [10] S. Ramelow, A. Mech, M. Giustina, S. Gröblacher, W. Wieczorek, J. Beyer, A. Lita, B. Calkins, T. Gerrits, S. W. Nam, A. Zeilinger, and R. Ursin, Highly efficient heralding of entangled single photons, Opt. Express 21, 6707 (2013).
  * [11] S. Du, J. Wen, and M. H. Rubin, Narrowband biphoton generation near atomic resonance, J. Opt. Soc. Am. B 25, C98 (2008).
  * [12] M. Fleischhauer, A. Imamoglu, and J. P. Marangos, Electromagnetically induced transparency: Optics in coherent media, Rev. Mod. Phys. 77, 633 (2005).