%\documentclass[aps,pra,twocolumn]{revtex4-2} 
\documentclass[aps,pra,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}

\begin{document}
\title{Direct generation of correlated four photons 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 report the direct observation of multiphoton states generated from a bright narrowband correlated photon source based on double-$\Lambda$ four-wave mixing in a cold atomic ensemble. We investigate the temporal correlations among twin modes using second-order, joint third-order, and joint fourth-order correlation functions. The instantaneous generation rate of photon quadruplets is estimated at $2.7(1.2)\times 10^{5}$ counts per second, along with an instantaneous pair generation rate of around $3.4(7) \times 10^{6}$ counts per second. These findings provide valuable insights into the effects of the double-pair emission in the four-wave mixing. Furthermore, the high generation rates also indicate that double-pair emission from four-wave mixing can serve as a promising source for narrowband multiphoton states. \textbf{VP: Abstract to be worked on at the very end}
\end{abstract}


\maketitle

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}. Here, the third-order nonlinear interaction of a weak continuous-wave (cw) \textit{pump} field and a strong cw \textit{coupling} field with the atomic medium generates narrowband correlated Stokes (frequency $\omega_s)$ and anti-Stokes (frequency $\omega_{as}$) fields. This process produces a two-mode squeezed state of the correlated Stokes and anti-Stokes fields described by \cite{Leonhardt, Lvovsky, sangouardQuantumRepeatersBased2011}

\begin{equation}\label{eq:two_mode_squeezing_state2}
\begin{align}
    |\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{align}
\end{equation}
 where $\beta\equiv \mathrm{cosh}\zeta$, $ \alpha \equiv \mathrm{tanh}\zeta $, $\zeta$ is the integral of the effective interaction Hamiltonian 
 %($\hat{H}_{\mathrm{eff}}$) 
 over the duration of the interaction and is proportional to the pump power, $n$ is the number of photons and $|n,n\rangle$ indicates 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 the probability amplitude for generating double-pairs of Stokes and anti-Stokes photons ($p_4$) scales quadratically with the amplitude for generating a single correlated Stokes-anti-Stokes pair ($p_2$), within a single SFWM event (i.e $p_4=(p_2)^2 = \alpha^2 $). In this case the 4 photons in the double pair are correlated in time and frequency. 
 
 However, when $N$ multiple independent SFWM events occur the overall output state can be described as a tensor product of the output of individual SFWM events giving $|\psi\rangle ^{\otimes N}$. In this case, the overall amplitude for detecting double pairs would have contributions from: a) $N$ independent processes that each degenerate a correlated double pair and b)  $N(N-1)/2$ combinations where a single pair is generated each in two separate events \cite{deriedmattenTwoIndependentPhoton2004}. It is important to note that the double pairs that contribute to case (b) are not correlated. \cite{deriedmattenTwoIndependentPhoton2004} analyses the distinction between generating a correlated double pair vs two independent double pairs from a spontaneous parametric process and shows that the limit of a large $N$, $p_4$ and $p_2$ become Poissonian distributions and relate with each other as $p_4 = p_2^2/2$. 

 Here we demonstrate a source of highly correlated narrow-band photon quadruplets or double pairs where the four-photons are generated simultaneously from a single SFWM event. We demonstrate this by ...

 

.. section to be expanded and improved later

\section{Experimental Setup}
\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{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 initialisation 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}$. Channels 1 and 2 analyse the Stokes field, and Channels 3 and 4 the anti-Stokes fields in a Handbury-Brown-Twiss like setup. $\lambda/2$: Half wave-plate, $\lambda/4$: quarter wave-plate, PBS: polarising beam-splitter, Ch: channel, APD: avalanche photo-diode. }
\end{figure}
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 field  (of frequency $\omega_{c}$) resonant to the $|5{S}_{1/2}, F=2\rangle \rightarrow |5{P}_{1/2}, F=2\rangle$ transition. 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_{as}$ resonant to the $|5{P}_{1/2}, F=2\rangle \rightarrow |5{S}_{1/2}, F=1\rangle$ transition (refer \ref{fig:setup}(a)).  The fields are circularly polarised, 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 (\ref{fig:setup} (b)). The atoms in the MOT are initialized in the $|5S_{1/2}, F=1\rangle$ hyperfine ground state and the MOT trapping beams are switched off during the SFWM measurement. The optical depth (OD) of the atomic cloud is $\sim 30$. The generated correlated photon-pairs are collected from a waist of 175 $\mu$m, at an angle of 1$^{\circ}$, to reduce background scattering from the cw fields. Polarisation filters and temperature controlled etalon filters (bandwidth $\sim$ 100 MHz) are implemented in both the 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 beam-splitters (BS) and detected using avalanche photo-diodes (APDs). A timestamp card with $2\ \text{ns}$ time resolution records the photon arrival times in each of these four channels. Second, third and fourth order field correlations are analysed using this information.


\section{Results}
\subsection{Second-Order Correlation}\label{section:second-order_correlation}


\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{cross_correlation_wo_fit.pdf}
    \caption{Normalized second-order correlation measurements. (a) The Stokes-anti-Stokes cross-correlation as a histogram of coincidences for various delays $\tau$, normalised 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 of periodicity 18 ns are caused by the Rabi frequency of the coupling field. Insets: Unheralded auto-correlation measurements of Stokes photons $g^{(2)}_{s,s}(\tau)$ (blue) and anti-Stokes photons $g^{(2)}_{as,as}(\tau)$ (red) (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 correlations and statistical properties of the generated fields. The normalised 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, as\}$ for the Stokes($s$) and anti-Stokes ($as$) mode. 


The second-order auto correlations $g^{(2)}_{{s,s}}(\tau)$, $g^{(2)}_{{as,as}}(\tau)$ and cross-correlation $g^{(2)}_{{s,as}}(\tau)$ were measured for pump and coupling powers of about $800 \ \mu\text{W}$ and $ 10\ \text{mW}$ respectively and a pump detuning of $\Delta_{p} = 40\ \text{MHz}$. From the $g^{(2)}_{{s\text{-}as}}(\tau)$ results shown in Fig.\ref{fig:crosscorrelation_w_auto} we infer a correlation time of around $16$ ns between the Stokes and anti-Stokes photons. The oscillations in the coincidences are due to the strong coupling beam which leads to population transfer between $|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:generationrate}.

The Stokes and anti-Stokes modes independently display thermal statistics as seen from their intensity auto-correlation 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)}_{as,as}(0) = 2.02 \pm 0.07$ for the anti-Stokes mode. 


\subsection{ Triple coincidences } \label{sec:TripleCoincidence}
\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{triplet_count.pdf}
    \caption{ Triple-coincidence measurement. (a) Triple coincidences for various delays $\tau_{3s}$ and $\tau_{4s}$ between a detection in Ch3 and Ch4 respectively and a heralding Stokes photon, over a measurement duration $T_m$ of 0.7 h. (b) Comparison of the vertical ridge with the normalized second-order cross-correlation. Red dots: Mean triplet count, averaged over $\tau_{4\text{-}s}$ from 20 ns to 60 ns and normalized by the accidental triplet counts. The accidental triplet here is $R_{s}R_{3}R_{4}\delta t^{2}T_{m}$, where $\delta t =2 \ \text{ns}$ is the time bin and $R_i$ is the singles count in channel $i$. Solid line: normalized cross-correlation function $g^{(2)}_{\text{as-s}}(\tau_{3\text{-}s})$ between Stokes and Ch3. (c) Ratio of triple coincidence peak value with the average value along the vertical ridge away from the peak. Blue dots: Triple coincidences at $\tau_{3\text{-}s} = 8$ ns divided by mean triple coincidences at $\tau_{3\text{-}s} = 8$ ns and $\tau_{4\text{-}s  }= 20-60$ ns. The peak is approximately 4 times the value in the ridge.}
    \label{fig:triplet_w_corr}
\end{figure}

We analyse the temporal distribution of coincidences involving more than two detections, to distinguish correlated double pairs from two independent pairs detected together by chance. Both these scenarios lead to states containing multiples of Stokes-anti-Stokes pairs and there is no physical mechanism that generates states involving three photons. Thus a measurement of triplet coincidences involving two Stokes and one anti-Stokes photons or two anti-Stokes and one Stokes photon, provides the same information as a four-fold coincidence measurement of two anti-Stokes and two Stokes photons, while being simpler to acquire and visualise.  


%Analysis of triple coincidences between two photon detections in the anti-Stokes (Stokes) mode and a single detection in the Stokes (anti-Stokes) mode, is a direct measure of the rate of generating more than one Stokes-anti-Stokes pair from the SFWM process within the coincidence window under consideration. 


%The triple-coincidence rate for finding anti-Stokes photons at times $t_3$ and $t_4$ and a Stokes photon at time $t_s$ is 


The normalised 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)}(t_3,t_4,t_s) =  \\ & \frac{ \langle \hat{E}^{\dagger}_{s}(t_{\text{s}})\hat{E}^{\dagger}_{as}(t_{\text{3}}) \hat{E}^{\dagger}_{as}(t_{4})\hat{E}_{as}(t_{4})\hat{E}_{as}(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}_{as}(t_{3})\hat{E}_{as}(t_{3}) \rangle\langle \hat{E}^{\dagger}_{as}(t_{4})\hat{E}_{as}(t_{4}) \rangle}\ 
    \end{aligned}
\end{equation}

where the numerator gives the triple-coincidence rate $P_T(t_3,t_4,t_s)$. This can be expressed in terms of the second-order correlations as shown in the Appendix in Eq. \ref{eq:TripletReduced}. 

\begin{comment}
    
Similar to SPDC, the joint Stokes-anti-Stokes output state of SFWM is a zero-mean Gaussian state. Hence the Gaussian moment factoring theorem (Isserli's theorem) can be applied \cite{shapiroSemiclassicalQuantumBehavior1994} to express the triple coincidences in terms of the only non-zero second-order moments which are the temporal autocorrelation $R(\tau)$ 
\begin{equation}
    R(\tau) \propto \langle \hat{a}^\dagger_k(t+\tau) \hat{a}_k(t)\rangle \, k \in\{s, as\} ,
\end{equation}
and the phase sensitive cross-correlation
\begin{equation}
    C(\tau) \propto \langle \hat{a}^\dagger_s(t+\tau) \hat{a}_i(t)\rangle,
\end{equation}
to give 
\end{comment}
 

A measurement of the triple coincidences for an anti-Stokes detection each in channels 3 (at $t_3$) and 4 (at $t_4$) and a Stokes detection in either of the channels 2 or 1 (at $t_s$), for a duration of 0.7 hours is shown in  Fig.\ref{fig:triplet_w_corr}, where the measurement was performed for the same conditions as in section \ref{section:second-order_correlation}. Here we have represented the results in terms of relative delays $\tau_{3s}= t_3-t_s, \tau_{4s}= t_4-t_s$ and $\tau_{34}=t_3-t_4$. The technique used to identify triplets from double coincidences is described in Appendix \ref{app:algorithms}.  

The features in Fig. \ref{fig:triplet_w_corr} are understood intuitively or from analyzing Eq.\ref{eq:TripletReduced} over various delays. For a coherence time $\Delta t$ for the Stokes and anti-Stokes photons, when $\tau_{3s}, \tau_{4s}, \tau_{34} \gg \Delta t$, $g^{(3)}(\tau_{3s},\tau_{4s},\tau_{34})$ reduces to the background accidental rate $R^3(0)$ normalised to 1.
%In the measurement this translates to $R_3R_4R_s$ a product of the singles rates in ch3, ch4 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$ but $\tau_{34} \lesssim \Delta t$, the autocorrelation in the anti-Stokes mode dominates the result  ($g^{(3)}(\tau_{3s},\tau_{4s},\tau_{34}) \rightarrow g^{(2)}_{as,as}(\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-Stoked mode, forming the moderately bright diagonal in Fig. \ref{fig:triplet_w_corr}. 


For $\tau_{3s}, \tau_{34} \gg \Delta t$ but $\tau_{4s} \lesssim \Delta t$ (horizontal ridge) or when $\tau_{4s}, \tau_{34} \gg \Delta t\, \& \, \tau_{3s} \lesssim \Delta t$ (vertical ridge) the cross-correlation between anti-Stokes (in Ch4 or Ch3 ) 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 photon from a separate SFWM event in the other anti-Stokes channel. Thus the maximum mean value in the horizontal and vertical ridges is equal to $g^{(2)}_{s, as}(0)$ as seen in Fig. \ref{fig:triplet_w_corr} (b). 

%$P_T(\tau_{3s},\tau_{4s},\tau_{34}) \rightarrow R^3(0) g^{(2)}_{as,s}(\tau_{4,s})$. Thus the horizontal ridge is proportional to the normalised intensity cross-correlation. Similarly the vertical ridge is $\propto R^3(0) g^{(2)}_{as,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. From Eq.\ref{eq:TripletReduced}, using $g^{(2)}_{as,as}(0)=2$, the theoretical peak value of the normalised triplet rate is $g^{(3)}(0,0,0) = 4  g^{(2)}_{as,s}(0)-2 \approx 4  g^{(2)}_{as,s}(0)$ when $g^{(2)}_{as,s}(0)\gg 1 $ which is true for highly non-classical pair sources. Thus, theoretically the peak is 4x the maximum in either of the horizontal or vertical ridges when the output contains highly correlated four-photon states. When the output contains fewer correlated four-photons and more uncorrelated double pairs from multiple SFWM events, the triplet peak at zero delays would arise from the overlapping of the ridges and would be closer to 2x 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)}_{as,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)}_{as,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 analysed.   

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 triplet peak is four times the mean along the vertical ridge (outside the central 20ns window). Thus, we verify that the output of the SFWM process contains strongly correlated double pairs, that contribute to the high three-fold coincidences in the triplet measurement. 
%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. 


\subsection{Quadruple coincidences}\label{sec:FourCoincidence}
\begin{figure}
\includegraphics[width=\linewidth]{four_photon_coind_dis.pdf}
 \caption{ Quadruple-coincidence measurement. (a) Photon quadruplet events from a detection recorded in each of the four channels within a $60 \ \text{ns}$ time window, shown as a distribution of delays $\tau_{\text{3-1}}$ (between Ch3 and Ch1) and $\tau_{\text{4-1}}$ (between Ch4 and Ch2), over a measurement period of $0.7$ hour. Red square indicates the peak within a $20$ ns coincidence window. (b) Same data plotted against the Stokes delay $\tau_{\text{1-2}}$ (between Ch1 and Ch2) and anti-Stokes  $\tau_{3\text{-}4}$ (between Ch3 and Ch4).}
    \label{fig:fourphoton_dis}
\end{figure}


We look for four-fold coincidences between detections of two photons in the anti-Stokes mode and two photons in the Stokes mode to obtain the double pair generation rate and as an additional check to verify the direct generation of correlated double-pairs 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{equation}\label{eq:quadrupletrate}
\begin{aligned}
 & Q(t_1, t_2, t_3, t_4) = \\
 & \langle \hat{E}^{\dagger}_{as}(t_{\text{4}})
 \hat{E}^{\dagger}_{as}(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}_{as}(t_{\text{3}})\hat{E}_{as}(t_{\text{4}}) \rangle.
\end{aligned}
\end{equation}
$Q$ can be expressed in terms of the first-order auto and cross correlations between the modes as shown in Appendix \ref{app:coincidencerates}. 

We measure four-fold coincidences for detections at times $t_1$ to $t_4$ in Channels Ch$_1$ to Ch$_4$ respectively, under the same conditions as in section \ref{section:second-order_correlation}, within a window of 100 ns. 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 Ch4-Ch2 or Ch3-Ch1 respectively. Fig.\ref{fig:fourphoton_dis}(b) shows dull horizontal and vertical ridges from four-fold coincidences between accidentals and thermally bunched photons in Ch2-Ch1 or Ch3-Ch4 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}
\textcolor{blue}{From the data we identify singles rates in each channel ($R_i,\, i\in \{1,2,3,4 \}$) and singles in the Stokes (anti-Stokes) modes as $R_{s(as)}=R_{1(3)}+R_{2(4)}$.}
From the pair, triplet and quadruplet coincidence measurements we calculate the measured pair rate ($R_{P}$), triplet rate ($R_{T}$), and quadruplet rate ($R_{Q}$) defined within a $t_{c}=20\ \text{ns}$ coincidence window, without subtraction of accidentals. \textcolor{red}{$t_c$ of 20 ns is appropriate as the detection of correlated-double pairs is peaked within this window as seen from sections \ref{sec:TripleCoincidence} and \ref{sec:FourCoincidence}}. \textcolor{blue}{$R_P$ is obtained from the sum of pairs detected between a Stokes and an anti-Stokes photon from the following channel pairs:  Ch1-Ch3, Ch1-Ch4, Ch2-Ch3 and Ch2-Ch3. $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, Ch3, Ch4$\}$, $\{$Ch2, Ch3, Ch4$\}$  and one anti-Stokes with two Stokes channels $\{$Ch3, Ch1, Ch2$\}$, $\{$Ch4, Ch1, Ch2$\}$. $R_Q$ is from the sum of quadruplet coincidences arising from a detection in each of the four channels within $t_c$.} Fig.\ref{fig:rate_vs_pumping_power} shows the singles, pairs, triplets and quadruplet rates as functions of pump power. The photon pair and quadruplet production saturate at pump powers over a mW. At low pump powers (approximately $<$ 200 $\mu$W) the pair rate scales linearly with the pump power and the triplet and quadruplet rate scales quadratically with the pump power.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{rate_vs_pumping_power.pdf}
    \caption{ Detection rates a function of pump power (a) Single count rates (left axis) in Stokes (blue dots) and anti-Stokes channels (red dots) and correlated photon pair rate (green dots, right axis) as functions of pumping field power. (b) Correlated photon triplet rates across three detector channels (magenta dots, left axis) and correlated photon quadruplet rate (grey dots, right axis) as functions of pumping field power. The detuning of the pumping field is $40 \ \text{MHz}$, while the coupling field is resonant with a fixed power of $10\ \text{mW}$. The atomic cloud has $\text{OD}\simeq 30$.}
    \label{fig:rate_vs_pumping_power}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{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{as}}$ (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}

\textcolor{red}{This is better visualised in Fig. \ref{fig:probability_vs_count_rate}}, where the pair, triplet and quadruplet rates are shown relative to singles rates in the Stokes and anti-Stokes modes, with axes in log scale. $R_{P}$ scales approximately linearly with $R_{s}$ and $R_{as}$, 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 s}$ are found to be $1.97 \pm 0.04$ and $2.34 \pm 0.03$, respectively. The slopes of $R_T$ and $R_Q$ match within the error bars as expected from the fact that triplet and quadruplet photons originate from \textcolor{red}{the same physical processes}. Furthermore, these measurements indicate that that the rate at which double pairs are detected scales close to quadratically with the rate of detecting correlated photon-pairs, which shows that \textcolor{red}{the photons in the double pairs are produced from a higher-order process within a frequency conversion event.} However, to identify the exact relationship between the probability of generating states with four photons and states with pairs of photons, we need to correct for accidentals and factor losses from inefficient optical paths, photon collection and detection in each channel. 

\subsection{Accidental corrected detection rates}\label{sec:accidentcorrRates}

At fixed pump and coupling powers of $800\, \mu$W and 10 mW respectively, and $\Delta_p=40$ MHz, the mean measured $R_s= (1.04 \pm 0.07)\times 10^{6} \text{\ cps}$ and $R_{as}= (1.10 \pm 0.06) \times 10^{6} \text{\ cps}$. To determine the accidental corrected values for the detected pair, triplet, and quadruplet rates, we performed accidental correction as described in the Appendix \ref{app:accidental}. This gives a corrected correlated pair detection rate $R_{cp}= (4.8 \pm 0.3) \times 10^{4} \text{\ cps}$. 

\textcolor{red}{As we can see in Appendix \ref{app:accidental}, the correction of accidentals for triplets depends on which mode is used as the herald. Thus, we report channel-specific accidental-corrected triplet-rates.} The corrected rate of detected photon triplets consisting of one Stokes photon and two anti-Stokes photons in Ch3 and Ch4 is $R_{cs34}=(330 \pm 80)\text{\ cps}$, while the rate for triplets consisting of two Stokes photons in Ch1 and Ch2 and one anti-Stokes photon is $ R_{cas12}=(299 \pm 84)\text{\ cps}$. \textcolor{blue}{The total accidental-corrected triplet-rate is $R_{ct}=R_{cas12}+R_{cs34}.$} Finally, the accidental subtracted correlated photon quadruplets is found to be %R_{cq}=(3.1 \pm 1.1) \times 10^{4}\text{\ cph}$ 
$R_{cq}=(8.6 \pm 3.0) \text{\ cps}$. 

\subsection{Channel losses}
\textcolor{blue}{We carefully characterise the losses in each channel to estimate the rate of correlated pairs and double-pairs directly generated from the SFWM process. }
The total transmission and detection probability in each channel $k \in\{1,2,3,4\}$, which includes transmission of the collection and filtering setup, splitting efficiency of the fiber based 50:50 beamsplitter and quantum efficiency of the detector in the respective channel is denoted by $\eta_{k}$. \textcolor{red}{The optical losses in each channel is determined by measuring the transmission of a laser beam (at the target wavelength) from outside the vacuum chamber to just before the detector in each channel. We measure transmissions of  $\approx 11.5\%$ each for Ch1 and Ch2 and $ 12.5\%$ each for Ch3 and Ch4, which include the contributions from the filter Etalon and fiber beam-splitter.} Factoring in the quantum efficiencies of the APDs in each channel \textcolor{red}{(about 60-70$\%$ per APD)} the total efficiencies $\eta_{k}$ for each channel are given as 0.078, 0.083, 0.080, and 0.067 for channels 1,2,3 and 4 respectively. These values provide an estimate of the upper bound for effective efficiencies, as they do not account for absorption in the atomic ensemble or spatial mode mismatch between the photons and the collection optics. %In the absence of such additional losses, the measured triplet rate between three channels should directly relate to the quadruplet rate by the factor $\eta_k$ associated with the fourth channel. 
From the discrepancy between expected and measured ratios between the triplet and quadruplet rates, we infer additional losses \textcolor{blue}{that are} $38 \%$ for the channels in the Stokes arm and $25 \%$ for the anti-Stokes channels, which we attribute to \textcolor{red}{a combination of the above-mentioned factors.}

\subsection{Generation rates and probabilities for pairs vs double pairs}\label{sec:generationrate}

\textcolor{blue}{
We wish to estimate the probabilities $p_2 =|\alpha_2|^2$ and $p_4=|\alpha_4|^2$ of generating pairs and double-pairs from the accidental-corrected coincidence-rates reported in section \ref{sec:accidentcorrRates} and the efficiencies $\eta_k$, assuming the following general form for the state at the output of the $^{87}$Rb cloud 
\begin{equation}
    \alpha_0|0,0\rangle+\alpha_2|1,1\rangle+\alpha_4|2,2\rangle.
\end{equation}
We also wish to estimate rates $R_{g,2}$ and $R_{g,4}$ for generating pairs and double-pairs respectively, from the SFWM process. }

\textcolor{red}{We note that the generation of a pair or a double pair is a random event in time and consequently their detections are also random events in time in the absence of information about the generation.} Thus, the probability ($P$) of either generating (or detecting) a certain number of photons within $t_c$ is related to the respective rate ($R$) of generating (or detecting) that many photons by $P = Rt_{c}e^{-Rt_{c}}$. \textcolor{blue}{Here $P$ is $P_p$ or $P_t$ or $P_q$ and $R$ is $R_{cp}$ or $ R_{ct}$ or $ R_{cq}$ for the detected probabilities and rates while $P$ is $ p_2$ or $p_4 $ and $R $ is $ R_{gp}$ or $R_{gq}$ for the probabilities and rates associated with generation of pairs or double-pairs.} 

%Using the accidental-corrected detection rates mentioned above, we wish to estimate a) the rate of generating pairs and quadruplets from the SFWM process and b) check if the probability of generating double-pairs scales quadratically with the probability of generating pairs. The latter is a strong indication that the quadruplets from the SFWM process are correlated and are not caused by the merging of pairs from multiple SFWM. For this, we ident

%The detection of a pair, a triplet, or a quadruplet is a random event in time, and the probability of these processes follows a Poisson distribution in time. The probability of detecting genuine pairs ($P_2$) triplets ($P_3$) and quadruplets ($P_4$) within the coherence time $t_c$ is related to the respective accidental corrected detection rates for pairs ($R_{cp}$), triplets ($R_{ct}$) and quadruplets ($R_{cq}$) through $P = Rt_{c}e^{-Rt_{c}}$. 



\textcolor{blue}{Using the detection rates from \ref{sec:accidentcorrRates}, we estimate the probabilities for detecting coincidences between genuine pairs, triplets and quadruplets $\{P_p,P_t, P_q\}$.} Following the procedure in Appendix \ref{app:imperfectdetection}, we then relate these with the probabilities of generating pairs and quadruplets $\{p_2, p_4\}$ by factoring in imperfect detectors and losses in each channel. We then once again use the Poissonian probability to rate relation to arrive at  $R_{gp}$ and $R_{gq}$, the generation rates for pairs and double-pairs.  



%Under the assumption that the pair and quadruplet generation processes are each random processes in time following Poissonian distributions, the probability for generating a correlated pair ($p_2$) or quadruplet ($p_4$) relates to the generation rate of pairs $R_{g,2}$ and quadruplets $R_{g,4}$ through $p_j=R_{g,j}t_{c}e^{-R_{g,j}t_{c}}$ for $j\in{2,4}$.
Based on this, we calculate the double-pair probability $p_{4}=0.0055(24)$, corresponding to $ R_{g,4}=2.7 (1.2) \times 10^{5}\ \text{cps}$. We obtain a pair generation probability  $p_2=0.0681(140)$, from $R_{g,2}=3.4(7) \times 10^{6}\ \text{cps}$. \textbf{VP: I think if we get $p_2$ and $p_4$ values for 2 more pump powers thenwe can say something about the factor $\kappa$ in $p_4=\kappa p_2^2$. We can say if it is closer to 1 or further from 1. I think with just 1 data point (also large error bars) it is hard to make the statement that was here previously which was ``The quadratic relationship holds". }
% The quadratic relation holds between photon pair probability and four-photon stat probability. 

\section{Conclusion}
Main points to make in conclusion \textbf{please add you thoughts in points here}
\begin{itemize}
    \item 
\end{itemize}


\appendix

\section{Double, Triple and Quadruple coincidences}\label{app:coincidencerates}

Here we express the double, triple and four-fold coincidences between modes in terms of the phase sensitive first-order cross-correlation,
\begin{equation}
C(\tau_{ij}) = \langle \hat{E}_i(t+\tau_{ij}) \hat{E}_{j}(t)\rangle
\end{equation}
and the first-order autocorrelation
\begin{equation}
    R(\tau_{ij}) = \langle \hat{E}^\dagger_i(t+\tau_{ij}) \hat{E}_j(t)\rangle
\end{equation}
where $\{i,j\}\in \{s, as\}$ to represent the Stokes or anti-stokes modes or $\{i,j\}\in\{1,2,3,4\}$ to represent one of the four detection channels. 
Note $ R(0)$ is the pair generation rate. % and $R_{s,as}(\tau)=C_{s,s}(\tau)=C_{as,as}(\tau)=0$ 


It is well known that when the state under consideration is the output of a parametric photon-pair production such as SPDC or SFWM, the Gaussian moment factoring theorem (Isserli's theorem) can be applied to the intensity correlations in Eq.\ref{eq:normg2} to give the following expression for the normalised intensity cross-correlation \cite{shapiroSemiclassicalQuantumBehavior1994}
\begin{equation}\label{eq:normcross}
    g^{(2)}_{s,as}(\tau)=1+\frac{|C(\tau_{s,as})|^2}{|R(0)|^2},
\end{equation}

and the normalised intensity auto-correlation as
\begin{equation}\label{eq:normauto}
    g^{(2)}_{i,i}(\tau)=1+\frac{|R(\tau_{i,i})|^2}{|R(0)|^2}
\end{equation}

%For the output of SPDC and SFWM the triple coincide rate depends on $C(\tau)$ and $R(\tau)$ as follows.
 
 The triple-coincidence rate $P_T$ for two anti-Stokes detections and a single Stokes detection, 
 \begin{equation}
     P_T(t_3,t_4,t_s)=\langle \hat{E}^{\dagger}_{s}(t_{\text{s}})\hat{E}^{\dagger}_{as}(t_{\text{3}}) \hat{E}^{\dagger}_{as}(t_{4})\hat{E}_{as}(t_{4})\hat{E}_{as}(t_{\text{3}})\hat{E}_{s}(t_{s}) \rangle
 \end{equation}
 can be similarly expanded and expressed in terms of relative delays to give \cite{razaviCharacterizingHeraldedSinglephoton2009,bettelliCommentCoherenceMeasures2010}
\begin{equation}\label{eq:tripletrate}
    \begin{aligned}
        P_T(\tau_{3s},\tau_{4s},\tau_{34})  = \, &   R(0)[R(0)^2 + |R(\tau_{34})|^{2}] \\ 
        & + R(0)[|C(\tau_{3s})|^{2}+|C(\tau_{4s})|^{2}] \\
        & + 2\text{Re}\{C(\tau_{3s})C^{\ast}(\tau_{4s})R(\tau_{34})\}. \\
        =\, & R(0)^3 g^{(2)}_{as, as}(\tau_{34})+ R(0)^3 g^{(2)}_{as, s}(\tau_{3s})\\ 
        &+R(0)^3g^{(2)}_{as, s}(\tau_{4s}) - 2R(0)^3  \\
        &  +2 R(0)^3 \sqrt{g^{(2)}_{as, s}(\tau_{3s})-1} \sqrt{{g^{(2)}_{as, s}(\tau_{4s})-1}} \\
        & \times \sqrt{{g^{(2)}_{as, as}(\tau_{34})-1}}.
    \end{aligned}
\end{equation}
This equation shows that accidental triplet events are caused by singles detected in each of the three channels, a correlated pair between the Stokes mode and Ch3 (Ch4) with an accidental in Ch4 (Ch3) and thermally bunched photons in the anti-Stokes (causing coincidences in Ch3 and Ch4) with an accidental in the Stokes mode (Ch1 or Ch2). 
The normalised third-order correlation can be expressed in terms of the normalised intensity auto and cross correlations using  \ref{eq:normauto} and \ref{eq:normcross} as
\begin{equation}
    \begin{aligned}
    g^{(3)}(\tau_{3s},\tau_{4s},\tau_{34})=&\frac{P_T(\tau_{3s},\tau_{4s},\tau_{34})}{R(0)^3}\\
        =\, & g^{(2)}_{as, as}(\tau_{34})+ g^{(2)}_{as, s}(\tau_{3s})\\ 
        &+g^{(2)}_{as, s}(\tau_{4s}) - 2  \\
        &  +2 \sqrt{g^{(2)}_{as, s}(\tau_{3s})-1} \sqrt{{g^{(2)}_{as, s}(\tau_{4s})-1}} \\
        & \times \sqrt{{g^{(2)}_{as, as}(\tau_{34})-1}}.
    \end{aligned} \label{eq:TripletReduced}
\end{equation}

Similarly the quadruplet rate for detecting two Stokes photons at times $t_1,$ and $t_2$ and two anti-Stokes photons at times $t_3$ and $t_4$ is given in Eq.\ref{eq:quadrupletrate}. The moment factoring theorem is used to obtain the following expression in terms of relative delays. 

\begin{equation}\label{eq:quadrupletrate}
    \begin{aligned}
        P_Q(\tau_{34},& \tau_{24}, \tau_{14},\tau_{23} ,\tau_{21},\tau_{13})=R(0)^4 \\
        &+R(0)^2[ |R(\tau_{43})|^2+|R(\tau_{21})|^2]\\
        &+R(0)^2[ |C(\tau_{23})|^2+|C(\tau_{24})|^2+|C(\tau_{13})|^2+|C(\tau_{14})|^2]\\
        &+2R(0)\mathrm{Re}\{R(\tau_{34})C^*(\tau_{24})C(\tau_{23})\}\\
        &+2R(0)\mathrm{Re}\{R(\tau_{34})C^*(\tau_{14})C(\tau_{13})\}\\
        &+2R(0)\mathrm{Re}\{R(\tau_{21})C^*(\tau_{14})C(\tau_{24})\}\\
        &+2R(0)\mathrm{Re}\{R(\tau_{21})C^*(\tau_{13})C(\tau_{23})\}\\
        &+ |R(\tau_{43})|^2|R(\tau_{21})|^2 \\ 
        &+|C(\tau_{13})|^2|C(\tau_{24})|^2+|C(\tau_{14})|^2|C(\tau_{23})|^2\\&+2\mathrm{Re}\{C^*(\tau_{24})C^*(\tau_{13})C(\tau_{14})C(\tau_{23})\}]\\
        &+2 \mathrm{Re}\{R(\tau_{34})R(\tau_{12})C^*(\tau_{24})C(\tau_{13})\}\\
        &+2 \mathrm{Re}\{R(\tau_{34})R(\tau_{21})C(\tau_{23})C^*(\tau_{14})\}\\
    \end{aligned}
\end{equation}

Out of the 17 terms that sum up to give $Q$ in Eq \ref{eq:quadrupletrate}, all terms other than the last three are contributions due to accidentals. Terms 2-7, are due to two accidentals combined with either a correlated Stokes-anti-Stokes pair or bunched photons in one of the Stokes or anti-Stokes modes. Terms 8-11 are caused by an accidental combined with a correlated triplet in three channels. Term 12 is from bunching in both the Stokes and anti-Stokes modes. Terms 13 and 14 are correlated pairs from separate SFWM events contributing to four-fold coincidences. Terms 15, 16 and 17 are due to correlated double-pairs from the same SFWM event. 

\section{Event searching algorithms}\label{app:algorithms}
The search for multiple coincidences over timestamp data of four channels is a computationally resource intensive task. 
We employ the following strategy to simplify the search.
We first identify pair coincidences at various relative delays for each pair from the following possible pairs of Stokes-anti-Stokes channels  1-3, 1-4, 2-3 and 2-4. %We note the respective photon arrival timestamp data.  
The triplet coincidences are then identified based on the pair detections that share a photon arrival timestamp. Quadruplet events are identified from triplet detections that share a common timestamp with a pair detection.

For instance, to identify triplet events involving Ch1, Ch3, and Ch4, we first locate pairs between Ch1 and Ch3 detected at timestamps $t_{1}$ and $t_{3}$ respectively. Similarly, we locate pairs between Ch1 and Ch4 detected at timestamps $t_{1}'$ and $t_{4}$. Pair events where $t_{1}=t_{1}'$ are taken to form the triplet event $(t_{1}, t_{3}, t_{4})$. This information can be used to identify quadruplets between the four channels. For this we compare pair events between between Ch2 (at $t_2$) and Ch4 (at time $t_4^*$)  with the previously identified triplet event with times $t_{1}, t_{3}, t_{4}$. When $t_{4}=t_{4}^*$, the events are combined to form a quadruplet detection with the timestamp $t_1, t_2, t_3,t_4$. Thus, we can efficiently search for pair, triplet, and quadruplet events from timestamp data and plot the temporal distribution of these coincidences.


\section{Accidental Correction Procedure}\label{app:accidental}
Subtraction of accidental coincidences is performed for $R_P, R_T$ and $R_Q$ which are the pair, triplet and quadruplet rates within a 20 ns coincidence window ($t_c$). 
The corrected correlated photon pair rate $R_{cp}= R_P- R_{ap}$ where the accidental pair rate $R_{ap} = \sum_{i,j} R_{i}R_{j}t_{c}$ for $\{i,j\} \in \{\{1,3\},\{1,4\},\{2,3\}, \{2,4\}\}$. 

Accidental correction for triplets is more involved. We see from \ref{app:coincidencerates} that only the last term in Eq.\ref{eq:TripletReduced} is a consequence of double pairs and all other terms (corresponding to the diagonal, horizontal, vertical ridges and general background in Fig \ref{fig:triplet_w_corr}) contribute to accidentals and have to be corrected for. Let us consider triplet detections in  Ch1, Ch3 and Ch4. The measured equivalent of the overall background $R^3(0)$ is $R_{1}R_{3}R_{4}t_{c}^{2}$. The combined accidentals from the background, cross and auto-correlations is then $R_{a134} =\int_{0}^{t_{c}}\int_{0}^{t_{c}}R_{1}R_{3}R_{4}g^{(2)}_{as,as}(\tau_{34})d\tau_{31}d\tau_{41}+ \int_{0}^{t_{c}}\int_{0}^{t_{c}} R_{1}R_{3}R_{4}[g^{(2)}(\tau_{31}) + g^{(2)}(\tau_{41})] d\tau_{31} d\tau_{41}- 2R_{1}R_{3}R_{4}t_{c}^{2}$. This procedure is repeated for triplets involving channel combinations $\{2,3,4\}, \{3,1,2\}$ and $\{4,1,2\}$. The accidental-corrected triplet-rate within $t_c$ for triplets between Stokes-Ch3-Ch4 is $R_{cs34}=R_{s34}-(R_{a134}+R_{a234})$, where $R_{s34}=R_{134}+R_{234}$ is the sum of detected triplet rates between $\{1,3,4\}$ and $\{2,3,4\}$.  The accidental-corrected triplet rate within $t_c$ for triplets between anti-Stokes-Ch1 and Ch2 is $R_{cas12}=R_{as12}-(R_{a312}+R_{a412})$, where $R_{as12}=R_{312}+R_{412}$ is the sum of detected triplet rates between $\{3,1,2\}$ and $\{4,1,2\}$ .  

A similar procedure is used to estimate the accidentals in the quadruplets. Terms 1- 14 in Eq.\ref{eq:quadrupletrate} are accidentals. Their contributions are evaluated using the relations $R(\tau)= R(0)\int_{0}^{t_c}\sqrt{g^{(2)}_{\textrm{auto}}(\tau)-1}\,d \tau$ and $C(\tau)= R(0)\int_{0}^{t_c}\sqrt{g^{(2)}_{s,as}(\tau)-1}\,d \tau$

%equivalent of the intensity cross correlation contribution $R^3(0)g^{(2)}_{as, s}(\tau_{3s})+R^3(0)g^{(2)}_{as, s}(\tau_{4s})$ is $\int_{0}^{t_{c}}\int_{0}^{t_{c}} R_{1}R_{3}R_{4}[g^{(2)}(\tau_{31}) +g^{(2)}(\tau_{41})] d\tau_{31} d\tau_{41}$ and equivalent of the autocorrelation contribution $R^3(0)g^{(2)}_{as, as}(\tau_{34})$ is $\int_{0}^{t_{c}}\int_{0}^{t_{c}}R_{1}R_{3}R_{4}g^{(2)}_{as,as}(\tau_{34})d\tau_{31}d\tau_{41}$. All these contributions are eliminated as in Eq.\ref{eq:TripletReduced} to obtain the corrected correlated triplet rate involving channels 1,3 and 4. This procedure is repeated for triplets involving channel combinations $\{2,3,4\}, \{1,2,3\}$ and $\{1,2,4\}$. These are all summed to obtain $R_{at}$- the total accidental corrected triplet rate within $t_c$.  


%Quadruplet events are identified by combining triplet events in Ch1, Ch3 and Ch4 with pairs between Ch2 and Ch4 as described in \ref{app:algorithms}. To simplify our estimation of quadruplet accidentals, we treat the triplets in channels 1,3 and 4 as singles in a combined channel 134. It is straightforward to see that the background accidental quadruplet rate is then $R_{T}R_{4}t_{c}$. Similar to the triplet case, Another contribution to the accidentals comes from The other comes from cross-correlations between the Stokes-anti-Stokes pairs in the accidental triplets in Ch1, Ch3 and Ch4 and th



% Quadruplet events refer to coincidences that involve both triplet events (Ch1-Ch2-Ch3) and pairs (Ch1-Ch4). To estimate the accidental quadruplet rate $R_{aq}$, we employ an algorithm-dependent method to simplify the analysis of correlations across four channels. One contribution to the accidental quadruplet rate comes from the triplet events $R_{t}$ combined with the single photon rate $R_{4}$, expressed as $R_{t}R_{4}t_{c}$. The other arises from correlated photon pairs intersecting with uncorrelated photon triplets, given by $\int_{0}^{t_{c}}R_{at}R_{4}(g^{(2)}_{\text{as-s}}(\tau_{\text{4-1}})-1)d\tau_{\text{4-1}}$. These two contributions serve effective evaluation for the accidental quadruplet rate.  



\section{}\label{app:imperfectdetection}

In this section we derive a relationship between the probabilities of detecting coincidences between genuine pairs ($P_p$), triplets ($P_t$) and quadruplets ($P_q$) with the probabilities of generating pairs ($p_2$) and double pairs ($p_4$), by accounting for the transmission and detection efficiencies in each channel $\eta_k\, k\in\{1,2,3,4\}$. 

The detection operator of a photon detector that detectes a photon in the Stokes mode with a non-unit efficiency $\eta$ and which does not resolve the photon number is modeled by \cite{Sekatski_2012}
\begin{equation}
    \hat{D}(\eta_{k}) = \mathds{1} -(1-\eta_{k})^{\hat{a}^{\dagger}\hat{a}} . 
\end{equation}
where $\hat{a}$ represents the annihilation operator for a photon in the Stokes mode and $\eta_{k}$ is the total detection efficiency including transmission efficiency in channel $k$. A detection event in the anti-stokes mode is expressed similarly with $\hat{b}$ indicating annihilation of a photon in the anti-stokes mode. Note that we have neglected the contribution from darkcounts in the detector as this and incoherent scattering photons from the atomic cloud contribute solely to the flat background of accidental counts in the coincidence measurements. The probability of recording a detection event with this detector is $\text{Tr}[\hat{D}(\eta_{k})\hat{\rho}]$ where $\hat{\rho}$ is the density operator representing the initial photon state. 

%To characterize the generation rates of photon pairs and photon quadruplets before they traverse lossy collection channels, it is essential to account for imperfections in collection and detection processes. Dark counts from the detectors and incoherent scattering photons from the atomic cloud contribute solely to the flat background of accidental counts in the coincidence measurements between the twin modes. Consequently, the measured correlated photon pair rates remain unaffected by dark counts and incoherent noise photons. Referring to the theoretical treatment of imperfections in single-photon detectors \cite{Sekatski_2012}, we can model the detection of a single photon in channel $k$ using the following measurement operator:
%\begin{equation}
 %   \hat{D}(\eta_{k}) = \mathds{1} -\eta_{k})^{\hat{a}^{\dagger}\hat{a}} . 
%\end{equation}
%where $\hat{a}$ represents the annihilation operator for either the Stokes or anti-Stokes mode, $\eta_{k}$ denotes the total transmission and detection probability in channel $k$, which incorporate the transmission in the collection and filtering setup, $50\%$ splitting efficiency of the beamsplitter and quantum efficiency of the detector in this channel. Using this model, the probability of getting a click in this detector is expressed as $\text{Tr}[\hat{D}(\eta_{k})\hat{\rho}]$ where $\hat{\rho}$ is the density operator representing the initial photon state. 

---editing in progress

The detection of a correlated photon pair in is described by a detection  is modeled as 
\begin{equation}
\begin{aligned}
& \hat{D}_{\text{s}}(\eta_{k})\hat{D}_{\text{as}}(\eta_{l}) = \hat{D}_{\text{s}}(\eta_{k}) \otimes \hat{D}_{\text{as}}(\eta_{l}) \\
& = (\mathds{1} -(1-\eta_{k})^{\hat{a}^{\dagger}\hat{a}} )\otimes (\mathds{1} -(1-\eta_{l})^{\hat{b}^{\dagger}\hat{b}} )\ ,
\end{aligned}
\end{equation}
where $\hat{a}$ and $\hat{b}$ are the annihilation operators of Stokes and anti-Stokes mode, respectively. 
The normal-ordered form of $f(\hat{a}^{\dagger}\hat{a})$ is expressed as:
\begin{equation}
\begin{aligned}
f(\hat{a}^{\dagger}\hat{a}) &= 1 - :(1-\eta)^{\hat{a}^{\dagger}\hat{a}}: \\
&= 1 - :e^{\hat{a}^{\dagger}\hat{a} \ln{(1-\eta)}}: \
&\approx 1 - :e^{-\eta \hat{a}^{\dagger}\hat{a}}:
\end{aligned}
\end{equation}
Thus, the measurement operator for a pair in two channels of the same mode is given by:
\begin{equation}
\begin{aligned}
& \hat{D}_{\text{s}}(\eta_{1})\hat{D}_{\text{s}}(\eta_{2}) = \\ & \mathds{1} - (1-\eta_{1})^{\hat{a}^{\dagger}\hat{a}}- (1-\eta_{2})^{\hat{a}^{\dagger}\hat{a}} + (1-\eta_{1}-\eta_{2})^{\hat{a}^{\dagger}\hat{a}} \\
& \hat{D}_{\text{as}}(\eta_{3})\hat{D}_{\text{as}}(\eta_{4}) = \\  &\mathds{1} - (1-\eta_{3})^{\hat{a}^{\dagger}\hat{a}}- (1-\eta_{4})^{\hat{a}^{\dagger}\hat{a}} + (1-\eta_{3}-\eta_{4})^{\hat{a}^{\dagger}\hat{a}} \\
\end{aligned}
\end{equation}
For the photon triplets consisting of two photons in the Stokes mode and one photon in the anti-Stokes mode, the detection operator is given by:
\begin{equation}
\begin{aligned}
& \hat{D}_{\text{s}}(\eta_{1})\hat{D}_{\text{s}}(\eta_{2})\hat{D}_{\text{as}}(\eta_{l}) = \hat{D}_{\text{s}}(\eta_{1})\hat{D}_{\text{s}}(\eta_{2})\otimes \hat{D}_{\text{as}}(\eta_{l}) \\
& = \mathds{1} - (1-\eta_{1})^{\hat{a}^{\dagger}\hat{a}}- (1-\eta_{2})^{\hat{a}^{\dagger}\hat{a}} + (1-\eta_{1}-\eta_{2})^{\hat{a}^{\dagger}\hat{a}} \\
& \otimes  (\mathds{1} -(1-\eta_{l})^{\hat{b}^{\dagger}\hat{b}}
\end{aligned}
\end{equation}
On the other hand, for a photon triplet consisting of two photons in the anti-Stokes mode and one photon in the Stokes mode, the detection operator is given by:
\begin{equation}
\begin{aligned}
& \hat{D}_{\text{s}}(\eta_{k})\hat{D}_{\text{as}}(\eta_{3})\hat{D}_{\text{as}}(\eta_{4}) =  \hat{D}_{\text{s}}(\eta_{k})\hat{D}_{\text{as}}(\eta_{3})\hat{D}_{\text{as}}(\eta_{4}) \\
& =  (\mathds{1} -(1-\eta_{k})^{\hat{a}^{\dagger}\hat{a}} \\
& \otimes \mathds{1} - (1-\eta_{3})^{\hat{b}^{\dagger}\hat{b}}- (1-\eta_{4})^{\hat{b}^{\dagger}\hat{b}} + (1-\eta_{3}-\eta_{4})^{\hat{b}^{\dagger}\hat{b}} \\
\end{aligned}
\end{equation}
The detection of a photon quadruplet across four channels is modeled as follows:
\begin{equation}
\begin{aligned}
&\hat{D}_{\text{s}}(\eta_{1})\hat{D}_{\text{s}}(\eta_{2})\hat{D}_{\text{as}}(\eta_{3})\hat{D}_{\text{as}}(\eta_{4})\\  & =\hat{D}_{\text{s}}(\eta_{1})\hat{D}_{\text{s}}(\eta_{2})\otimes\hat{D}_{\text{as}} (\eta_{3})\hat{D}_{\text{as}}(\eta_{4})  \\
& =  \mathds{1} - (1-\eta_{1})^{\hat{a}^{\dagger}\hat{a}}- (1-\eta_{2})^{\hat{a}^{\dagger}\hat{a}} + (1-\eta_{1}-\eta_{2})^{\hat{a}^{\dagger}\hat{a}}  \\
& \otimes \mathds{1} - (1-\eta_{3})^{\hat{b}^{\dagger}\hat{b}}- (1-\eta_{4})^{\hat{b}^{\dagger}\hat{b}} + (1-\eta_{3}-\eta_{4})^{\hat{b}^{\dagger}\hat{b}}
\end{aligned}
\end{equation}
When applying the general form of photon state
\begin{equation}
    \sqrt{a_{0}}|0,0\rangle + \sqrt{a_{1}}|1,1\rangle + \sqrt{a_{2}}|2,2\rangle
\end{equation}
The photon quadruplet detection probability is given as 
\begin{equation}\label{eq:four_detection}
P_{1234} = 4 a_{2} \eta_{1}\eta_{2}\eta_{3}\eta_{4}
\end{equation}
The photon triplet detection probabilities are  
\begin{equation}\label{eq:triplet_detection}
\begin{aligned}
P_{123} &= 2 a_{2} \eta_{1}\eta_{2}(2-\eta_{3})\eta_{3} \\
P_{124} &= 2 a_{2} \eta_{1}\eta_{2}(2-\eta_{4})\eta_{4} \\
P_{134} &= 2 a_{2} \eta_{3}\eta_{4}(2-\eta_{1})\eta_{1} \\
P_{234} &= 2 a_{2} \eta_{3}\eta_{4}(2-\eta_{2})\eta_{2} \\
\end{aligned}
\end{equation}
The correlated photon probabilities are 
\begin{equation}\label{eq:pair_detection}
\begin{aligned}
P_{13} &= \eta_{1}\eta_{3}(a_{1} + (2+\eta_{1})(2-\eta_{3})a_{2}) \\
P_{14} &=  \eta_{1}\eta_{4}(a_{1} + (2+\eta_{1})(2-\eta_{4})a_{2})\\
P_{23} &=  \eta_{2}\eta_{3}(a_{1} + (2+\eta_{2})(2-\eta_{3})a_{2})\\
P_{24} &=  \eta_{2}\eta_{4}(a_{1} + (2+\eta_{2})(2-\eta_{4})a_{2}) \\
\end{aligned}
\end{equation}

Taking into account the losses in the atomic ensemble, let the transmission efficiency of the Stokes mode be denoted as $\eta_{s}$ and the transmission in the anti-Stokes be denoted as $\eta_{as}$, Thus, the effective efficiency in each channel can be expressed as 
\begin{equation}
\eta_{i}^{\text{eff}} = \eta_{s} \eta_{i}\ . 
\end{equation}
This effective efficiency can replace $\eta_{i}$ in the preceding detection probability equations. 


\textbf{VP: To everyone: can we say something intelligent about the four photon correlation time from this measurement? If yes, any ideas how? Also can we say something about a Cauchy Swarz violation for 4 fold coincidences?}










\bibliography{reference}




\end{document}