%\documentclass[aps,pra,twocolumn]{revtex4-2} 
\documentclass[aps,prl,twocolumn,superscriptaddress]{revtex4-2} 
%\documentclass[aps,prl,reprint]{revtex4-2}
\usepackage{mathtools}
\usepackage{graphicx}
\usepackage[colorlinks=true,linkcolor=blue,citecolor=blue]{hyperref}
\usepackage{enumerate}
\usepackage{siunitx}
\usepackage{comment}
\usepackage{setspace}
\usepackage{dsfont}
\usepackage{float}

%\usepackage[style=apa, backend=biber]{biblatex}
%\addbibresource{reference2.bib}

\begin{document}
\title{Narrowband four-photon states from spontaneous four-wave mixing}
\author{Yifan~Li}
\author{Justin Yu Xiang Peh}
\author{Chang Hoong Chow}
\author{Boon Long Ng}
\author{Vindhiya Prakash}
\affiliation{Centre for Quantum Technologies, National University of Singapore, 3 Science Drive 2, Singapore 117543}
\author{Christian~Kurtsiefer}
\affiliation{Centre for Quantum Technologies, National University of Singapore, 3 Science Drive 2, Singapore 117543}
\affiliation{Department of Physics, National University of Singapore, 2 Science Drive 3, Singapore 117551}
\email[]{christian.kurtsiefer@gmail.com}


\begin{abstract}
We observe time-correlated four photons within a correlation window of 20 ns  from spontaneous four-wave mixing via a double-$\Lambda$ scheme in a cold cloud of $^{87}$Rb atoms. In contrast to high-power pulsed pumping of nonlinear processes in crystals, our scheme generates correlated four-photon states by direct continuous-wave pumping at nominal powers. We verify the presence of genuinely correlated four-photon states over accidentals by higher-order intensity cross-correlation measurements and accidental subtraction. We obtain a time-correlated four-photon generation rate of $2.5(4) \times 10^6$ counts per second close to saturation. The photons produced are near-resonant with atomic transitions, and have a bandwidth in the order of MHz, making them readily compatible with quantum networking applications involving atoms. 
\end{abstract}

\maketitle

Multiphoton states, i.e. states with more than two photons entangled or correlated across single or multiple modes, are extremely useful resources for quantum sciences and technologies \cite{JWPan_Multiphoton}. At a foundational level, multiphoton Greenberger–Horne–Zeilinger (GHZ) states and W states have enabled powerful tests to disprove local realistic theories \cite{GHZOriginal,Pan2000quantum} and explore unique entanglement classes \cite{PhysRevLett.95.150404, Cirac2000}. 
Multiphoton states enable secure communication protocols \cite{Mark_SecretSharing, JWPan_Crypto} and also find application in quantum metrology \cite{Mitchell2004}. In the form of cluster states, they are essential for scalable and resource-efficient photonic quantum computing \cite{Raussendorf_OneWay, Browne_PRL}. States with four photons have also been used to encode decoherence-resistant quantum information \cite{Weinfurter_DecoherenceFree}. 
 
%\emph{(Part in Italics to be removed) Several types of entangled multiphoton states, such as three and four-photon GHZ states have been prepared by interfering and combining photon pairs from spontaneous parametric downconversion (SPDC) in independent crystals \cite{ BouwmeesterGHZ, Mikami_Multiphoton, Yongmin_Wstate, Zhong2018}. Another approach utilizes the onset of laser-like stimulation in SPDC, by retro-reflecting a pulsed pump and an initially generated pair to stimulate the generation of further pairs entangled with the reflected two-photon state \cite{Antia, Simon_EntanglementLaser}. However, such schemes require synchronized pump pulsing and complicated geometries for erasure of which-path information.  Without the need for such complicated schemes,   }.

%Such states have been prepared by interferometrically combining independent pairs from spontaneous parametric downconversion (SPDC) \cite{BouwmeesterGHZ, Mikami_Multiphoton, Yongmin_Wstate, Zhong2018}. However, a simpler alternative is the well-known method of directly producing multiphoton states by %As an %simpler alternative to such schemes, 
It is well known that some multiphoton states can be directly produced by 
strong pumping of non-linear processes like spontaneous parametric downconversion (SPDC), where the probability of producing more than one photon pair increases with the pump power \cite{PhysRevA.100.041802, dengGeneralizedTwomodeSqueezed1993}.  In directly pumping an SPDC process, the probability of producing entangled four photons is twice as high as producing two independent entangled pairs \cite{Weinfurter2001, Weinfurter2003}. Highly entangled W states have also been obtained from the higher-order component in a directly pumped SPDC process \cite{KieselWState}. 

Photons from SPDC typically have large bandwidths and correspondingly short coherence lengths, shorter than the length of the downconversion crystal itself.  Thus, in the above cases, pulsed pumps with large instantaneous powers and narrowband filters are used to isolate and analyze correlated multiphoton states from SPDC. This leads to losses. Direct filter-free analysis of the rich temporal structure of higher-order correlated photons from SPDC \cite{bettelliCommentCoherenceMeasures2010} has been challenging due to the jitter and response averaging of detectors \cite{razaviCharacterizingHeraldedSinglephoton2009}.  The large bandwidth of photons from SPDC also limits their use in quantum memory and repeater schemes that require efficient interfacing with material quantum systems. 

Here, we demonstrate a bright source for narrowband time-correlated photon quadruplets, matched to atomic transitions, based on direct continuous-wave (cw) pumping of spontaneous four-wave mixing (SFWM) in a cold atomic cloud. SFWM in atomic clouds is an excellent and bright \cite{Shiu2024, ChenPRR2022, Vuletic2006} alternative to SPDC for producing narrowband photon pairs with long coherence times \cite{duSubnaturalLinewidthBiphotons2008,DduNarrowbandBiphotonGeneration2008, kolchinElectromagneticallyinducedtransparencybasedPairedPhoton2007, srivathsanNarrowBandSource2013, balicGenerationPairedPhotons2005}. Photons from this process can be spectrally shaped to be narrower or wider than atomic transition linewidths, making them well-suited for quantum networking applications \cite{Kimble2008quantum}, such as memory, repeater \cite{sangouardQuantumRepeatersBased2011}, and entanglement distribution schemes involving atoms. Furthermore, their long coherence times, typically in the order of tens of nanoseconds, allow them to be well-resolved by off-the-shelf photon detection electronics.

We study correlated photon quadruplets from SFWM in a cold cloud of rubidium atoms using double Hanbury Brown and Twiss (HBT) type setups, one in each of the correlated modes. We introduce an efficient technique for identifying three-fold and four-fold coincidences of the photons generated from the nonlinear interaction. We analyze the temporal distribution of the detected triplet and quadruplet coincidences and observe that photon pairs bunch together in both measurements within a correlation window of 20 ns. The aggregate coincidences detected within this window is significantly larger than the sum of accidentals detected at longer delays, %The peak in three-fold coincidences is four times the maximum value at longer delays. 
indicating a strong contribution from quadruplets correlated in time, over accidental/uncorrelated four-photon states. % We estimate an instantaneous photon quadruplet generation rate of $2.3\times 10^{6}$ counts per second close to saturation of the nonlinear frequency conversion. 

In detail, our scheme is based on SFWM using a double-$\Lambda$ configuration of energy levels in a cold cloud of $^{87}$Rb atoms, similar to the systems reported in \cite{kolchinGenerationNarrowBandwidthPaired2006, kolchinElectromagneticallyinducedtransparencybasedPairedPhoton2007}.
The SFWM process is driven by a weak cw pump (of frequency $\omega_{p}$) detuned by $\Delta_p$ from $|5{S}_{1/2}, F=1\rangle  \rightarrow |5{P}_{3/2}, F'=2\rangle $ and a strong cw coupling laser (of frequency $\omega_{c}$) resonant to the $|5{S}_{1/2}, F=2\rangle  \rightarrow |5{P}_{1/2}, F'=2\rangle $ transition. Nonlinear interaction of the pump and coupling fields with the atomic medium generates correlated optical fields called Stokes and anti-Stokes by convention. The Stokes photons are generated at a frequency $\omega_s$ close to the $|5{P}_{3/2}, F'=2\rangle \rightarrow |5{S}_{1/2}, F=2\rangle$ transition and the anti-Stokes photons have a frequency $\omega_{a}$ resonant to the $|5{P}_{1/2}, F'=2\rangle \rightarrow |5{S}_{1/2}, F=1\rangle$ transition (refer Fig.\ref{fig:setup}(a)).  The pump and coupling fields are circularly polarized, orthogonal to each other, and are directed at an elongated magneto-optical trap (MOT) of cold $^{87}$Rb atoms, along the long axis in a counter-propagating configuration (Fig.\ref{fig:setup} (b)). The SFWM process is precluded by initializing atoms in the MOT to the $|5S_{1/2}, F=1\rangle$ hyperfine ground state via optical pumping. The MOT cooling beams are switched off during the SFWM measurement. The optical depth (OD) of the atomic cloud is $\sim 30$. The spatial modes for collecting the Stokes and anti-Stokes photons are focused on the atomic ensemble with a waist of $175\ \mu\text{m}$. The collection modes are at an angle of 1$^{\circ}$ to the pump and coupling fields, to reduce background scattering. Polarization filters and temperature-controlled etalon filters (bandwidth $\sim$ 100 MHz) are implemented in both Stokes and anti-Stokes collection arms to suppress unwanted photons. The photons collected in the Stokes and anti-Stokes arms are split using 50:50 fiber beamsplitters (BS) and detected using single photon detectors (D1-D4). A timestamp card with $2\ \text{ns}$ timing resolution records the photon arrival times in each of these four detectors. Second, third and fourth-order field correlations are analyzed using this data.
\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{Figures/setup_figure.pdf}
    \caption{\label{fig:setup}(a) Energy levels involved in the Double-$\Lambda$ spontaneous four-wave mixing in $^{87}$Rb. Solid blue and red arrows indicate cw pump and coupling fields respectively. Wiggly blue and red arrows indicate generated Stokes and anti-Stokes fields. Black dots indicate initialization of atoms in the $F=1$ hyperfine ground state. (b) Schematic of experimental setup. The pumping and coupling beams have a waist of $\sim 0.85\ \text{mm}$. The collection spatial mode is focused on the atomic ensemble with a waist of $175\ \mu\text{m}$. Detectors D1 and D2 detect the Stokes field, and D3 and D4 the anti-Stokes fields in a Hanbury-Brown and Twiss like setup. $\lambda/2$: half-wave plate, $\lambda/4$: quarter-wave plate, PBS: polarizing beamsplitter, D1-D4: single photon detectors. }
\end{figure}

A single frequency conversion process produces the following output state that can contain multiple Stokes and anti-Stokes photons \cite{Leonhardt, Lvovsky, sangouardQuantumRepeatersBased2011, deriedmattenTwoIndependentPhoton2004},

\begin{equation}\label{eq:two_mode_squeezing_state2}
    \begin{split}
    |\Psi\rangle & =  \frac{1}{\beta} \sum_{n=0}^{\infty} \left( \alpha \right)^{n} |n,n \rangle  \, \\
    %&= \frac{1}{\beta} \left( |0,0\rangle + \alpha |1,1\rangle + \alpha^2|2,2\rangle + ...\right).
    \end{split}
\end{equation}
 Here, $\beta\equiv \mathrm{cosh}\zeta$, $ \alpha \equiv \mathrm{tanh}\zeta $, $\zeta$ is the product of the strength of the nonlinear interaction and the duration of interaction, %($\hat{H}_{\mathrm{eff}}$) 
 and $|n,n\rangle$ indicates correlated Fock states with $n$ photons each in the Stokes and anti-Stokes modes. A complete expression for the interaction Hamiltonian and the nonlinear susceptibilities can be found in \cite{wenTransverseEffectsPairedphoton2006, DduNarrowbandBiphotonGeneration2008}. From Eq.~\ref{eq:two_mode_squeezing_state2}, it is evident that at small interaction strengths ($\zeta\ll 1$) the probability of generating states with four photons ($P_4$) relates to the probability of producing pairs ($P_2$) as $P_4=P^2_2$. In this case the four-photon state corresponds to two Stokes and two anti-Stokes photons that are correlated and entangled,  %scales perfectly quadratically with the amplitude for generating a single correlated Stokes-anti-Stokes pair, 
 generated within a single SFWM process \cite{deriedmattenTwoIndependentPhoton2004}. Furthermore, this is twice the probability of four-photon states present if the output contains a Poissonian distribution of photons \cite{deriedmattenTwoIndependentPhoton2004} (refer Supplementary Material). 
 
 %However, in the event of multiple independent frequency conversion processes occurring, four-photon states can contain uncorrelated and unentangled double pairs. In this case, the probability of creating $n$ pairs ($P'_n$) within a certain time window, is described by a Poisson distribution of mean $\mu$, $P'_n=e^{-\mu}\mu^n/n!$. For small $\mu$, $P'_4 \approx P^2_2/2$ \cite{deriedmattenTwoIndependentPhoton2004}. 
 In the following, we characterize the composition of correlated quadruplets versus uncorrelated double pairs from our source using higher-order intensity correlation measurements and verify the presence of twice as many correlated quadruplets as uncorrelated four-photons.  We also compare rates of pair production versus four-photon production for varying pump powers. 
 %From higher order intensity correlation measurements, we verify that the maximum quadruplets detected within a 20 ns correlation window is indeed Within a correlation window of 20 ns we observe a peak in three-fold coincidences that is twice the value at longer delays. This indicates that we observe quadruplets correlated in time within 20 ns, as opposed to a Poisson distribution of uncorrelated four photons detected outside the correlation time.


\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{Figures/cross_correlation_wo_fit.pdf}
    \caption{Normalized second-order correlation measurements. Main figure: Stokes-anti-Stokes cross-correlation as a histogram of coincidences for various delays $\tau$, normalized by the Stokes and anti-Stokes singles rates for a 2 ns bin size and an integration time of 150 s. Results averaged over 17 measurements. Oscillations are caused by the coupling field that drives the $|2\rangle \rightarrow |3\rangle$ transition at an effective Rabi frequency of $2\pi \times 55\, \text{MHz}\,$.  Insets: Unheralded autocorrelation measurements of Stokes photons (blue) with peak $g^{(2)}_{s,s}(0) = 2.07 \pm 0.02$  and anti-Stokes photons (red) with peak $g^{(2)}_{a,a}(0)=2.02 \pm 0.07$ (jointly labeled $g^{(2)}_{\text{auto}}(\tau)$). }
    \label{fig:crosscorrelation_w_auto}
\end{figure}


% $\Omega = 2\pi \times 55\ \text{MHz}$ and a decay rate of $\Gamma = 2\pi \times 10\ \text{MHz}$

We measure the second-order intensity correlations as a first step towards characterizing the statistical properties of the generated fields. The normalized second-order correlation between stationary fields $\hat{E}_i$ in mode $i$, detected at time $t_i$, and $\hat{E}_j$ detected at time $t_j = t_i+\tau_{ji}$ is \cite{glauberQuantumTheoryOptical1963}  

\begin{equation}\label{eq:normg2}
    g^{(2)}_{ji}(\tau_{ji})= \frac{\langle\hat{E}^{\dagger}_{i}(t_{i})\hat{E}^{\dagger}_{j}(t_{i}+\tau_{ji})\hat{E}_{j}(t_{i}+\tau_{ji})\hat{E}_{i}(t_{i})\rangle}{\langle \hat{E}^{\dagger}_{j}(t_{i}+\tau_{ji})\hat{E}_{j}(t_{i}+\tau_{ji})\rangle \langle \hat{E}^{\dagger}_{i}(t_{i})\hat{E}_{i}(t_{i})\rangle }\ , 
\end{equation}
where $ i,j \in\{s, a\}$ for the Stokes ($s$) and anti-Stokes ($a$) modes. 


The second-order autocorrelations $g^{(2)}_{{s,s}}(\tau)$, $g^{(2)}_{{a,a}}(\tau)$ and cross-correlation $g^{(2)}_{{s,a}}(\tau)$ were measured for pumping and coupling powers of about $800 \ \mu\text{W}$ and $ 10\ \text{mW}$ respectively and a pump detuning of $\Delta_{p} = 2\pi \times 40\ \text{MHz}$. Results are shown in Fig.\ref{fig:crosscorrelation_w_auto}. We infer a Stokes-anti-Stokes two photon correlation time of around $16$ ns from the time constant of an exponentially decaying fit to the $g^{(2)}_{{s,a}}(\tau)$ results. %The oscillations in the coincidences are due to the strong coupling beam which leads to an oscillation of population between energy levels $|2\rangle$ and $|3\rangle$ at an effective Rabi frequency of $2\pi \times 55\, \text{MHz}$\cite{kolchinElectromagneticallyinducedtransparencybasedPairedPhoton2007}. The pair generation rates inferred from these measurements are given in section \ref{sec:accidentcorrRates}.
The Stokes and anti-Stokes modes independently display thermal statistics as seen from their intensity autocorrelation at $\tau=0$ (inset in Fig.\ref{fig:crosscorrelation_w_auto}).% We measure $g^{(2)}_{s,s}(0) = 2.07 \pm 0.02$ for the Stokes mode and $g^{(2)}_{a,a}(0) = 2.02 \pm 0.07$ for the anti-Stokes mode. 


\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{Figures/triplet_count_normal.eps}
    \caption{ Normalized third-order correlation. (a) Normalized triple coincidences $g^{(3)}_{a,a,s}$ for various delays $\tau_{3s}$ and $\tau_{4s}$ between a detection in D3 and D4 respectively and a heralding Stokes photon in either of D1 or D2. Coincidences analyzed from data acquired over a measurement duration $T_m$ of 0.7 h, normalized by the accidental triplet rate $R_{s}R_{3}R_{4}\delta t^{2}T_{m}$, where the time bin $\delta t =2 \ \text{ns}$ and $R_i$ is the singles count in detector D$i$. The $g^{(3)}_{a,a,s}$ peak value of 18 indicates strongly correlated triplets.  (b) Comparison of the vertical ridge with $g^{(2)}_{a,s}$. Red dots: $g^{(3)}_{a,a,s}$ results averaged over $\tau_{4s}$ from 20 ns to 60 ns.  Solid line: normalized cross-correlation $g^{(2)}_{a,s}(\tau_{3s})$ between a Stokes detection in D1 or D2 and an anti-Stokes detection in D3. (c) Peak to ridge ratio. Blue dots: Trace at $\tau_{3s} = 10\,\mathrm{ns}$, normalized by average value at $\tau_{3s} = 10$ ns and $\tau_{4s  }= 20$ to $60$ ns, i.e., at long delays. The peak is close to 4 times the mean value in the ridge.}
    \label{fig:triplet_w_corr}
\end{figure}

We analyze the temporal distribution of coincidences involving more than two detections, to determine the ratio of correlated quadruplets to two independent pairs detected together by chance. %Both these scenarios lead to states containing multiples of Stokes-anti-Stokes pairs. 
Since there is no physical mechanism that generates states involving three photons, a measurement of triplet coincidences involving two Stokes and one anti-Stokes photons or two anti-Stokes and one Stokes photon provides similar information to a four-fold coincidence measurement of two anti-Stokes and two Stokes photons, while being faster to acquire and simpler to visualize.  



The normalized third-order correlation between the Stokes and anti-Stokes modes from two anti-Stokes detections at times $t_3$ and $t_4$ and a Stokes detection at time $t_s$ is 
\begin{equation} \label{eq:normg3}
    \begin{aligned}
        & g^{(3)}_{a,a,s}(t_3,t_4,t_s) =  \\ & \frac{ \langle \hat{E}^{\dagger}_{s}(t_{\text{s}})\hat{E}^{\dagger}_{a}(t_{\text{3}}) \hat{E}^{\dagger}_{a}(t_{4})\hat{E}_{a}(t_{4})\hat{E}_{a}(t_{\text{3}})\hat{E}_{s}(t_{s}) \rangle}{\langle \hat{E}^{\dagger}_{s}(t_{s})\hat{E}_{s}(t_{s}) \rangle \langle \hat{E}^{\dagger}_{a}(t_{3})\hat{E}_{a}(t_{3}) \rangle\langle \hat{E}^{\dagger}_{a}(t_{4})\hat{E}_{a}(t_{4}) \rangle}\ ,
    \end{aligned}
\end{equation}
where the numerator gives the triple-coincidence rate $G^{(3)}_{a,a,s}(t_3,t_4,t_s)$. This can be expressed in terms of second-order correlations as shown in Eq. 7 in %Eq.\ref{eq:TripletReduced} 
the Supplementary Material. 

Figure \ref{fig:triplet_w_corr} shows $g^{(3)}_{a,a,s}$ for triplets from an anti-Stokes detection each in D3 (at $t_3$) and D4 (at $t_4$) and a Stokes detection in either of D2 or D1 (at $t_s$), where the measurement was performed under the same conditions as the second order correlation measurements %\ref{section:second-order_correlation}. 
The results are represented in terms of relative delays $\tau_{3s}= t_3-t_s$ and  $\tau_{4s}= t_4-t_s$. The technique used to identify triplets from pair coincidences is described in the Supplementary Material %\textbf{\ref{app:algorithms}}.  

The features in Fig.\ref{fig:triplet_w_corr} can be understood intuitively or by analyzing Eq. 7 from the Supplementary Material %\ref{eq:TripletReduced} 
over various delays . Given a coherence time $\Delta t$ for the Stokes and anti-Stokes photons, when $\tau_{3s}, \tau_{4s}, \tau_{34} \gg \Delta t$, the triplet rate reduces to the background accidental rate which is normalized to 1 in Fig \ref{fig:triplet_w_corr}.
%In the measurement this translates to $R_3R_4R_s$ a product of the singles rates in D3, D4 and in the Stokes arm. 
%The diagonal in Fig \ref{fig:triplet_w_corr} is obtained for 
When $\tau_{3s}, \tau_{4s}, \gg \Delta t$ and $\tau_{34} \lesssim \Delta t$, the autocorrelation in the anti-Stokes mode dominates the result ($g^{(3)}_{a,a,s}(\tau_{3s},\tau_{4s},\tau_{34}) \rightarrow g^{(2)}_{a,a}(\tau_{34})$). In this case, the triplets are caused by the combination of an accidental click in the Stokes mode with a bunched thermal state in the anti-Stokes mode, forming the moderately bright diagonal in Fig.\ref{fig:triplet_w_corr}. 

When $\tau_{3s}, \tau_{34} \gg \Delta t$ but $\tau_{4s} \lesssim \Delta t$ (horizontal ridge) or when $\tau_{4s}, \tau_{34} \gg \Delta t\,$ but $ \tau_{3s} \lesssim \Delta t$ (vertical ridge) the strong cross-correlation between anti-Stokes (in D4 or D3 respectively) and Stokes photon pairs are the dominant contributions. Here, the triplets are formed by a combination of a correlated Stokes-anti-Stokes pair with an uncorrelated additional photon in the anti-Stokes mode. Thus, the maximum mean value in the horizontal and vertical ridges at long delays is equal to $g^{(2)}_{s, a}(0)$ as seen in Fig. \ref{fig:triplet_w_corr} (b). 

%$G^{(3)}_{a,a,s}(\tau_{3s},\tau_{4s},\tau_{34}) \rightarrow R^3(0) g^{(2)}_{a,s}(\tau_{4,s})$. Thus the horizontal ridge is proportional to the normalized intensity cross-correlation. Similarly the vertical ridge is $\propto R^3(0) g^{(2)}_{a,s}(\tau_{3,s})$ as the only contributing terms besides $R(0)$ come from the delay $\tau_{3s} \lesssim \Delta t$. This is illustrated in \ref{fig:triplet_w_corr} (b).  

In the region where $\tau_{3s}, \tau_{34}, \tau_{4s} \lesssim \Delta t$ the coincidences increase several-fold. Theoretically, the peak is expected to be 4 times the average in either the horizontal or vertical ridges (at delays longer than $\Delta t$) when the output contains highly correlated four-photon states (refer Supplementary Material). We see from Fig.\ref{fig:triplet_w_corr} (c) that in our measurement, the $g^{(3)}_{a,a,s}$ peak is about four times the mean along the vertical ridge (outside the central 20 ns window). Thus, we verify that the output of the SFWM process contains twice as many strongly correlated four-photons as uncorrelated double-pairs, that contribute to the high three-fold coincidences in the triplet measurement. 



%In the region where $\tau_{3s}, \tau_{34}, \tau_{4s} \lesssim \Delta t$ the coincidences increase several-fold. The theoretical maximum possible value of the normalized triplet rate is $g^{(3)}_{a,a,s}(0,0,0) = 4  g^{(2)}_{a,s}(0)-2 \approx 4  g^{(2)}_{a,s}(0)$ when $g^{(2)}_{a,s}(0)\gg 1 $, a condition true for highly non-classical pair sources (using $g^{(2)}_{a,a}(0)=2$ in Eq.7 from the Supplementary Material). Thus, theoretically, the peak is expected to be 4 times the average in either the horizontal or vertical ridges when the output contains highly correlated four-photon states. This is consistent with the discussion following Eq. \ref{eq:two_mode_squeezing_state2} since $P'_4$ includes the sum of average values in the ridges at long delays. 

%Due to the long coherence time of the Stokes and anti-Stokes photons, our triplet measurement is not limited by averaging effects due to detector resolution, which would have otherwise reduced the maximum of the triple-coincidence peak.  We see from Fig.\ref{fig:triplet_w_corr} (c) that in our measurement, the $g^{(3)}_{a,a,s}$ peak is about four times the mean along the vertical ridge (outside the central 20 ns window). Thus, we verify that the output of the SFWM process contains twice as many strongly correlated four-photons as uncorrelated double-pairs, that contribute to the high three-fold coincidences in the triplet measurement.

 %Since the horizontal and vertical ridges arise due to accidentals from Poissonian distributions of photons at longer delays, the coincidences from correlated four-photons will be twice the coincidences from accidentals following the discussion after Eq. \ref{eq:two_mode_squeezing_state2}.  The factor of 4 is equivalent to twice this contribution from the sum of both ridges. %When the output contains fewer correlated four-photons and more uncorrelated double-pairs leading to Poissonian distributions of photons, the triplet peak at zero delays would arise purely from the overlapping of the ridges and would be closer to 2 times the maximum in either of the horizontal or vertical ridges.

%In the region where $\tau_{3s}, \tau_{34}, \tau_{4s} \lesssim \Delta t$ the coincidences increase several-fold. The peak value of the triplet coincidence rate is $4 R^3(0) g^{(2)}_{a,s}(0)-2R^3(0)$. At a low accidental triplet rate $R^3(0)\rightarrow 0$, the peak triplet coincidence rate is 4 times the maximum triplet rate in the horizontal or vertical ridges $R^3(0) g^{(2)}_{a,s}(0)$. This is a signature that the SFWM output is dominated by correlated double pairs generated from the same SFWM event. In the presence of multiple uncorrelated SFWM events, the accidental triplet rate $R^3(0)$ would be non-negligible thereby the peak triplet rate would be $<4\times$ the rate in the horizontal and vertical ridges. The same result is arrived at in \cite{deriedmattenTwoIndependentPhoton2004}(refer eq. A16 and A17), where the extent of correlation between double pairs from pulsed SPDC is analyzed.   

 
%Thus we can better resolve the rich structures in the triple coincidences unlike previously reported measurements performed with broadband SPDC outputs \cite{razaviCharacterizingHeraldedSinglephoton2009}. 
%Thus we can directly verify the strong correlation between the double pairs that contribute to the triplet measurement, even with cw-pumping, from the four-fold increase in the peak measured triplet rate as compared to the maximum seen in the horizontal or vertical ridge outside the 20 ns Stokes-anti-Stokes correlation window. 

%	content...


%\subsection{Quadruple coincidences}\label{sec:FourCoincidence}


\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{Figures/four_photon_3d.eps}
    \caption{Quadruple-coincidence detection. Each slice shows unnormalized four-fold coincidences from a detection in each of D1 to D4 for a fixed delay $\tau_{12}$ with a 2 ns time bin, and a range of delays $\tau_{31}$ and $\tau_{41}$. Data acquired over a measurement duration of 0.7 h. The coincidences are peaked for $\tau_{12}=0\pm 2$ ns and $\tau_{31}$ and $\tau_{41}= 8 \pm 2$ ns.   }
    \label{fig:enter-label}
\end{figure}


We search for four-fold coincidences between detections of two photons in the anti-Stokes mode and two photons in the Stokes mode for further analysis of four-photon states produced from our SFWM source. 
%The quadruplet rate for four-fold coincidences from two Stokes detections at times $t_1$ and $t_2$ respectively and two anti-Stokes detections at time $t_3$ and $t_4$ respectively is described by the correlation function  
\begin{comment}
\begin{equation}\label{eq:quadrupletrate}
\begin{aligned}
 & G^{(4)}_{s,s,a,a}(t_1, t_2, t_3, t_4) = \\
 & \langle \hat{E}^{\dagger}_{a}(t_{\text{4}})
 \hat{E}^{\dagger}_{a}(t_{\text{3}}) \hat{E}^{\dagger}_{s}(t_{2})\hat{E}^{\dagger}_{s}(t_{1})\hat{E}_{s}(t_{1})\hat{E}_{s}(t_{2})\hat{E}_{a}(t_{\text{3}})\hat{E}_{a}(t_{\text{4}}) \rangle.\\
\end{aligned}
\end{equation}
\end{comment}
The normalized fourth-order cross-correlation is 

\begin{small}
\begin{equation}\label{eq:g4}
\begin{aligned}
 & {   g^{(4)}_{s,s,a,a}(t_1, t_2, t_3, t_4) =} \\
 & { \frac{\langle \hat{E}^{\dagger}_{a}(t_{\text{4}})
 \hat{E}^{\dagger}_{a}(t_{\text{3}}) \hat{E}^{\dagger}_{s}(t_{2})\hat{E}^{\dagger}_{s}(t_{1})\hat{E}_{s}(t_{1})\hat{E}_{s}(t_{2})\hat{E}_{a}(t_{\text{3}})\hat{E}_{a}(t_{\text{4}}) \rangle}{\langle \hat{E}^{\dagger}_{s}(t_{1})\hat{E}_{s}(t_{1}) \rangle \langle \hat{E}^{\dagger}_{s}(t_{2})\hat{E}_{s}(t_{2}) \rangle \langle \hat{E}^{\dagger}_{a}(t_{3})\hat{E}_{a}(t_{3}) \rangle\langle \hat{E}^{\dagger}_{a}(t_{4})\hat{E}_{a}(t_{4}) \rangle}},
\end{aligned}
\end{equation}
\end{small}where the numerator gives $G^{(4)}_{s,s,a,a}(t_1, t_2, t_3, t_4)$: the quadruplet rate for coincidences from two Stokes detections at times $t_1$ and $t_2$ respectively and two anti-Stokes detections at time $t_3$ and $t_4$ respectively. %This can be expressed in terms of the first-order auto and cross-correlations between the modes as shown in the Supplementary Material. 

We identify four-fold coincidences for detections at times $t_1$ to $t_4$ in detectors D1 to D4, under the same experimental conditions as used in previous measurements, with maximum delays up to 60 ns. We represent the data as sliced three dimensional histogram plots where each slice shows quadruplets for a fixed delay $\tau_{12}$ and various relative delays $\tau_{31}$ and $\tau_{41}$. We see the maximum density of quadruplets clustered around $\tau_{12}= 0 \pm 2$ ns and $\tau_{31}$ and $\tau_{41}= 810 \pm 2$ ns. Outside a 20 ns window centered at $(\tau_{12},\tau_{31},\tau_{41})=(0 \, \textrm{ns},10 \, \textrm{ns},10 \, \textrm{ns})$  the quadruplet count drops significantly indicating the presence of highly-correlated quadruplets within 20 ns. The horizontal and vertical ridges in slices as $\tau_{12}$ approaches $0$ arises from four-fold coincidences between accidentals and a correlated pair between D4-D1 or D3-D1 respectively. A relatively dull diagonal due to four-fold coincidences between accidentals and thermally bunched photons in D3-D4 can also be seen.

Due to the long coherence time of the Stokes and anti-Stokes photons, our triplet and quadruplet measurements are not limited by averaging effects due to detector resolution, which would have otherwise reduced the maximum of the triple and quadruple-coincidence peaks. %\textbf{VP:If anyone wants to say something more please suggest what can be added.}

%For easier visualisation, we represent the results as a projection on to the plane spanned by $\tau_{31}\equiv t_3 - t_1$ and $\tau_{42}\equiv t_4-t_2$ in Fig.\ref{fig:fourphoton_dis}(a) and by $\tau_{12}\equiv t_1 - t_2$ and $\tau_{34}\equiv t_3-t_4$ in Fig.\ref{fig:fourphoton_dis}(b). Fig.\ref{fig:fourphoton_dis}(a) shows bright horizontal and vertical ridges from four-fold coincidences between accidentals and a correlated pair between D4-D2 or D3-D1 respectively. Fig.\ref{fig:fourphoton_dis}(b) shows dull horizontal and vertical ridges from four-fold coincidences between accidentals and thermally bunched photons in D2-D1 or D3-D4 respectively. We see an increase in quadruplets by several-folds for relative delays within 20 ns. In both sub-figures, we see an increase in the quadruplet counts around zero delays that exceeds the combination of the horizontal and vertical ridges.  

%\subsection{Detection Rates}
\begin{figure}[h]
    \centering
    \includegraphics[width=0.95\linewidth]{Figures/rate_vs_pumping_power.eps}
    \caption{Top: Detection rates as a function of pump power. (a) Single count rates (right axis) of Stokes (blue circles) and anti-Stokes photons (red circles), and photon pair rate (green squares, left axis) as functions of pump power. (b) Total photon triplet rates (magenta triangles, left axis) and photon quadruplet rate (black diamonds, right axis) as functions of pump power. Bottom: Ratio of pairs, triplets and quadruplets to singles. The photon pair rate $R_{p}$ (green squares), photon triplet rate $R_{t}$ (magenta triangles), and photon quadruplet rate $R_{q}$ (black diamonds) from data in plots (a) and (b) represented in log-scale relative to the singles count rate $R_{\text{s}}$ in Stokes mode (c) and single count rate the anti-Stokes $R_{\text{a}}$ (d).  The change in single count rates is achieved by varying the pump power while keeping all other parameters constant. The detuning of the pump is $40 \ \text{MHz}$, while the coupling field is resonant with a fixed power of $10\ \text{mW}$. The atomic cloud has an $\text{OD}\simeq 30$.}
    \label{fig:rate_vs_pumping_power}
\end{figure}
We next examine the pump power dependency of multiphoton states. Fig.\ref{fig:rate_vs_pumping_power} (a), (b) show the total rates of Stokes/anti-Stokes singles ($R_{s/a}$), pairs ($R_{p}$), triplets ($R_{t}$) and quadruplets ($R_{q}$), aggregated over all detector combinations, defined within a coincidence window of $t_{c}=20\ \text{ns}$ and without subtraction of accidentals. $t_c$ of 20 ns is appropriate as the detection of correlated double-pairs is peaked within this window as seen from the correlation measurements

%From the data, we obtain singles rates in each detector ($R_i,\, i\in \{1,2,3,4 \}$) and singles in the Stokes (anti-Stokes) modes as $R_{s(a)}=R_{1(3)}+R_{2(4)}$.From the pair, triplet and quadruplet coincidence measurements we obtain the measured pair rate as $R_{p} = R_{13}+R_{14}+R_{23}+R_{24}$ ($R_{ij}$ is the rate for coincidences between  D$i$ and D$j$), triplet rate $R_{t}=R_{134}+R_{234}+R_{123}+R_{124}$ ($R_{ijk}$ is the rate for coincidences between D$i$, D$j$ and D$k$), and quadruplet rate ($R_{q}$) defined within a coincidence window of $t_{c}=20\ \text{ns}$, without subtraction of accidentals. $t_c$ of 20 ns is appropriate as the detection of correlated double-pairs is peaked within this window as seen from the correlation measurements. Fig.\ref{fig:rate_vs_pumping_power} (a) and (b) shows the singles, pairs, triplets and quadruplet rates as functions of pump power. The production of photon pairs and quadruplets saturates at pump powers over a mW. At low pump powers (approximately $<$ 200 $\mu$W) the pair rate scales linearly with the pump power while the triplet and quadruplet rate scales quadratically with the pump power.
%$R_p$ is obtained from the sum of pairs detected between a Stokes and an anti-Stokes photon from the following channel pairs:  D1-D3, D1-D4, D2-D3 and D2-D3. $R_t=R_{s34}+R_{as12}$ is from the total triplets detected in two possible combinations each of one Stokes with two anti-Stokes channels $\{$Ch1, D3, D4$\}$, $\{$Ch2, D3, D4$\}$  and one anti-Stokes with two Stokes channels $\{$Ch3, D1, D2$\}$, $\{$Ch4, D1, D2$\}$. $R_q$ is from the sum of quadruplet coincidences arising from a detection in each of the four channels within $t_c$. 


\begin{comment}
\begin{figure}
    \centering
    \includegraphics[width=0.9\linewidth]{Figures/rate_vs_single_count.pdf}
     \caption{ Ratio of pairs, triplets and quadruplets to singles. The photon pair rate $R_{p}$ (green dots), photon triplet rate $R_{t}$ (magenta dots), and photon quadruplet rate $R_{q}$ (black dots) from Fig. \ref{fig:rate_vs_pumping_power} represented in log-scale relative to the singles count rate $R_{\text{s}}$ in Stokes mode (a) and singles count rate the anti-Stokes $R_{\text{a}}$ (b).  Both axes are plotted on a logarithmic scale. The variation in the single count rates is achieved by varying the pump power while keeping all other parameters constant.}\label{fig:probability_vs_count_rate} 
\end{figure}
\end{comment}

A clear relationship between the pairs and quadruplets for increasing pump powers is better visualized in Fig.\ref{fig:rate_vs_pumping_power} (c) and (d), where the $R_p$, $R_t$ and $R_q$ are shown relative to $R_s$ and $R_a$, with axes in log scale. $R_{p}$ scales approximately linearly with $R_{s}$ and $R_{a}.$ %, with fitted slope parameters of $0.86 \pm 0.02$ and $0.95\pm 0.03$, respectively. 
%The slopes of $R_{t}$ and $R_{q}$ relative to $R_{s}$ are $1.87 \pm 0.02$ and $2.11 \pm 0.03$, respectively, while the slopes of $R_{t}$ and $R_{q}$ relative to $R_{a}$ are found to be $1.97 \pm 0.04$ and $2.34 \pm 0.03$, respectively. 
The slopes of $R_t$ and $R_q$ are both close to 2, which is expected from the fact that triplet and quadruplet photons originate from the same physical processes. Furthermore, these measurements indicate that the rate at which double pairs are detected scales close to quadratically with the rate of detecting photon-pairs. This confirms that the photons in the double pairs are produced from a higher-order process in frequency conversion. %However, for a more accurate relationship between the probability of generating states with four photons and states with pairs of photons, correct for accidentals and factor losses from inefficient optical paths, photon collection and detection in each detector, as described in the Supplementary Material. 
At a pump power of $800\,\mu$W, which is close to saturation, we detect pairs at the rate of $7.1(3)\times 10^{4} \mathrm{\, cps}$ and quadruplets at the rate of $21(3)\mathrm{\, cps}$. We perform accidental subtraction obtain a detection rate of $3(1)\,\mathrm{cps}$ for truly correlated quadruplets. Based on the detection rates and characterization of optical losses, we infer a quadruple generation rate of $ g_{q}=2.5 (4) \times 10^{6}\ \text{cps}$ and a pair generation rate of $g_{p}=1.3(3) \times 10^{7}\ \text{cps}$ at this power. Details on the procedure for accidental subtraction, channel losses and inferring generation rate can be found in the Supplementary Material. 

So far, narrowband multiphoton states have been demonstrated by spatially multiplexing two SFWM processes, spontaneous Raman events, and cascaded geometries \cite{dongExperimentalRealizationNarrowband2017, parkGenerationBrightFourphoton2022, Wu_2016, hubelDirectGenerationPhoton2010}. Our results show that direct pumping holds the potential to be a simpler alternative to producing multiphoton states. Furthermore, the correlated four photons generated here may potentially also be time-energy entangled \cite{Park_17, deriedmattenTwoIndependentPhoton2004}. To our knowledge, microscopic models of SFWM as a collective process involving individual atomic emitters, have dealt only with the generation of photon-pairs \cite{jiangQuantumLangevinTheory2023,raymondooiCorrelationPhotonPairs2007}. Our results pave the way for the extension of such models to better understand the microscopic origin of correlated multiphoton states, i.e., the physical origin of enhancement in the generation of a correlated second pair in the same mode as an initially generated pair. In the case of correlated four-photons generated from pulsed pumping of SPDC, such enhancement has been attributed to the bosonic nature of photons that leads to preferential bunching of otherwise indistinguishable photons into the same mode \cite{deriedmattenTwoIndependentPhoton2004, Weinfurter2001, Weinfurter2003}. While this may explain our observations, there may be other contributions as well. The long coherence length of the generated photon-pairs in our experiment exceeds the length of the nonlinear medium and may potentially lead to stimulation of other pairs into the same mode, a process similar to the observations in \cite{Antia}. We invite theoretical analyses of our experiment to explore the above processes and identify the underlying mechanism. 

\bibliographystyle{apsrev4-2}
\bibliography{reference2.bib}


%\typeout{} 

%\printbibliography



\end{document}