%\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}.
\textcolor{red}{There are too many citations using this scheme, and these two papers are not done by our group, it is confusing to mention "similar to"}
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. \textcolor{red}{When talking about the third-order nonlinear interaction, this term involves the generated field, it is odd to express  ``third-order nonlinear interaction... generates...''} This process produces a two-mode squeezed state of the correlated Stokes and anti-Stokes fields described by \cite{Leonhardt, Lvovsky, sangouardQuantumRepeatersBased2011}
 \textcolor{red}{CK has asked us not to use the term ``two-mode squeezed state''.}
\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, \textcolor{red}{``the integral'' definition is ambiguous. the effective interaction Hamiltonian is not given.} $n$ is the number of photons and $|n,n\rangle$ indicates Fock states with $n$ photons \textcolor{blue}{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} which analyse the nonlinear interaction in great detail. \textcolor{red}{The other representations of Hamiltonian are irrelevant to our paper}. 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 \textcolor{blue}{Stokes-anti-Stokes} pair ($p_2$), within a single SFWM event (i.e $p_4=(p_2)^2 = \alpha^2 $). \textcolor{red}{According to the equations, it seems that $\alpha$ is defined as the probability amplitude. However, the probability amplitude should be expressed in its full form, including the normalization factor, because we can always factor out a constant from the bracket.}. From Eq. \ref{eq:two_mode_squeezing_state2} In this case the 4 photons in the double pair are correlated in time and frequency. 

  \textcolor{red}{The term ``FWM events'' seems odd. Typically, the two-mode photon state is used as the effective macroscopic representation. If the photon state is expressed as a two-mode squeezed state, the following scenario will not occur. In the reference paper, SPDC is done within a pulse. there is no correlation between pulses, in this case the definition of SPDC events is reasonable. }
  
 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 ...

 

  \textcolor{red}{British spelling or American spelling?}
\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 at the start of an experimental sequence   \textcolor{red}{ a parametric process is a process that remains the population of atoms unchanged. "at the start of" could be misleading.} (b) Schematic of the 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. \textcolor{red}{these two sentences are not clear about whether the MOT is continuously running or if the atoms are released from the MOT within the measurement window.} 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) \textcolor{red}{300?} 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}
\textcolor{red}{Since the topic of this paper is about correlated photon pairs and photon quadruplets, I recommend using creation (annihilation) operators to define the correlation functions, rather than field operators.}

\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 \textbf{input value}. 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)$). The accidental coincidence rate is \textbf{input value} }
    \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$, detected at time $t_i$, and $\hat{E}_j$ detected at time $t_j = t_i+\tau$ is \cite{glauberQuantumTheoryOptical1963}  
\textcolor{red}{'Stationary' is not necessary. Reasons:The definition of Glauber correlation functions does not require this assumption.
In the case of a pulse-pump configuration, we can still obtain the normalized autocorrelation with the relative delay time as the parameter (pulsed coincidence), if the correlation function is defined using photon creation and annihilation operators, and $\tau$ is the relative delay between single-photon events.
For coherent or thermal light, we can define the correlation function with field operators and 'stationary' fields to perform the parameter transformation for $\tau$}
\begin{equation}\label{eq:normg2}
    g^{(2)}_{i,j}= \frac{\langle\hat{E}^{\dagger}_{i}(t_{i})\hat{E}^{\dagger}_{j}(t_{i}+\tau)\hat{E}_{j}(t_{i}+\tau)\hat{E}_{i}(t_{i})\rangle}{\langle \hat{E}^{\dagger}_{j}(t_{i}+\tau)\hat{E}_{j}(t_{i}+\tau)\rangle \langle \hat{E}^{\dagger}_{i}(t_{i})\hat{E}_{i}(t_{i})\rangle }\ . 
\end{equation}

In our case, $ i,j \in\{s, as\}$ for the Sokes($s$) and anti-Stokes ($as$) mode. Since the joint \textcolor{blue}{Stokes-anti-Stokes} output state of SFWM is a zero-mean Gaussian state, the Gaussian moment factoring theorem (Isserli's theorem) can be used to reduce Eq.\ref{eq:normg2} in terms of the phase sensitive first order cross-correlation $  C(\tau) = \langle \hat{E}^\dagger_s(t+\tau) \hat{E}_{as}(t)\rangle,$ and the autocorrelation $R(\tau) = \langle \hat{E}_i(t+\tau) \hat{E}_i(t)\rangle$ \cite{shapiroSemiclassicalQuantumBehavior1994}

(Note $R(0)$ is the pair generation rate).  The intensity cross-correlation is
\begin{equation}\label{eq:normcross}
    g^{(2)}_{s,as}(\tau)=1+\frac{|C(\tau)|^2}{|R(0)|^2},
\end{equation}
and the intensity auto-correlation is
\begin{equation}\label{eq:normauto}
    g^{(2)}_{i,i}(\tau)=1+\frac{|R(\tau)|^2}{|R(0)|^2}.
\end{equation}
\textcolor{red}{CK does not like ``zero-mean Gaussian state''}

\textcolor{red}{The theoretical expression of second-order correlation is not our topic}
\begin{comment}

\begin{equation}\label{eq:normg2}
    g^{(2)}_{ij}= \frac{\langle\hat{E}^{\dagger}_{i}(t_{i})\hat{E}^{\dagger}_{j}(t_{i}+\tau)\hat{E}_{j}(t_{i}+\tau)\hat{E}_{i}(t_{i})\rangle}{\langle \hat{E}^{\dagger}_{j}(t_{i}+\tau)\hat{E}_{j}(t_{i}+\tau)\rangle \langle \hat{E}^{\dagger}_{i}(t_{i})\hat{E}_{i}(t_{i})\rangle }\ . 
\end{equation}
and is proportional to the rate of coincidences detected between modes $i$ and $j$. Here $ i,j \in\{s, as\}$ for the Sokes($s$) and anti-Stokes ($as$) modes. 


\textbf{VP: Do we need this expression below? Also, here the second order correlation results are normalised whereas later the third and fourth order correlation results are not normalised. Should everything be consistent? We can decide on this later.}
\begin{equation}
g^{(2)}_{{i\text{-}j}}(\tau) = \frac{\langle \hat{a}^{\dagger}_{i}(t_{i})\hat{a}^{\dagger}_{j}(t_{i}+\tau)\hat{a}_{j}(t_{i}+\tau)\hat{a}_{i}(t_{i})\rangle
}{\langle \hat{a}^{\dagger}_{j}(t_{i}+\tau)\hat{a}_{j}(t_{i}+\tau)\rangle \langle \hat{a}^{\dagger}_{i}(t_{i})\hat{a}_{i}(t_{i})\rangle }\ . 
\end{equation}

Here $ i,j \in\{s, as\}$, $\hat{a}_{s(as)}$ is the annihilation operator in the Stokes (anti-Stokes) mode and $\langle \dots \rangle$ is the expectation value over the generated two-mode output state. 
\end{comment}

The second-order correlations $g^{(2)}_{{s,s}}(\tau)$, $g^{(2)}_{{as,as}}(\tau)$ and $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. 
\textcolor{red}{The waveform is not exactly the exponential decay. the definition of correlation time is not clear}
The oscillations in the coincidences are due to the strong coupling beam \cite{kolchinElectromagneticallyinducedtransparencybasedPairedPhoton2007} \textcolor{red}{``the strong coupling beam'' this reason is too trivial. at least we need to use some descriptions in the references}
which has an effective Rabi frequency of $2\pi \times 55\, \text{MHz}$. The pair generation rates inferred from these measurements are given in \textbf{include section reference}
\textcolor{red}{It is abrupt to mention the pair rate directly here, as all discussions regarding event rates are in the next section.}
In accordance with Eq.\ref{eq:normauto} 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)}_{\text{s-s}}(0) = 2.07 \pm 0.02$ for the Stokes mode and $g^{(2)}_{\text{as-as}}(0) = 2.02 \pm 0.07$ for the anti-Stokes mode. \textcolor{red}{The importance is that autocorrelation aligns with thermal photon statistics, not by equations}



\subsection{ Triple coincidences }
\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{triplet_count.pdf}
    \caption{\textbf{VP: I have not edited this yet}(a)~Photon triplet coincidence counts $N(\tau_{3\text{-}s}, \tau_{4\text{-}s})$, with paired events in Ch3 and Ch4, heralded by a photon event in Stoke mode. (b) Solid line: normalized cross-correlation function $g^{(2)}_{\text{as-s}}(\tau_{3\text{-}s})$ between Stokes and Ch3. Red dots: mean trace averaged over $\tau_{4\text{-}s}$ from 20 ns to 60 ns, normalized by the accidental triplet counts $R_{s}R_{3}R_{4}\delta t^{2}T_{m}$, where $\delta t =2 \ \text{ns}$ is the time bin and $T_{m}$ is the total measurement time. The offset above the cross-correlation function is due to diagonal contributions, indicating the auto-correlation of the anti-Stokes mode. (c) Vertical trace along vertical ridge $\tau_{3\text{-}s}$, normalized by the far-delay background. The maximal triplet coincidence count is approximately four times higher than the vertical ridge value.}
    \label{fig:triplet_w_corr}
\end{figure}

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. \textcolor{red}{the rate is discussed in next section. the coincidence window is not defined.}

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 \textcolor{red}{t3 t4 have exact meaning in our case: they refer to the photon detection timing in ch3 and ch4. ``anti-Stokes photons at times $t_3$ and $t_4$'' is too general to point out any photon in anti-Stokes mode}
\begin{equation}
    P(t_3,t_4,t_s) \equiv \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}

\textcolor{red}{The $P$ notation is meaningless. I suggest using Glauber correlation notations $G^{(3)}$ to express the joint third-order correlation function through triplet coincidence.}
This expression can be further reduced with the Gaussian moment factoring theorem and equations \ref{eq:normauto} and \ref{eq:normcross} to \cite{razaviCharacterizingHeraldedSinglephoton2009,bettelliCommentCoherenceMeasures2010}
 \textcolor{red}{``reduce'' Gaussian moment factoring theorem does not reduce the terms. it just allows to calculate higher-order moments using lower-order moments}

\begin{equation}
    \begin{aligned}
        P(\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^3(0)g^{(2)}_{as, as}(\tau_{34})+ R^3(0)g^{(2)}_{as, s}(\tau_{3s})\\ 
        &+R^3(0)g^{(2)}_{as, s}(\tau_{4s}) - 2 R^3(0) \\
        &  +2 R^3(0)\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}
 \textcolor{red}{This equation is too long}
\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}
Here we have replaced absolute detection times with relative delays $\tau_{3s}\equiv t_3-t_s, \tau_{4s}\equiv t_4-t_s$ and $\tau_{34}=t_3-t_4$. 
\textcolor{red}{The relative delay times are simply the time difference between two photon events. There is no need to use the equivalence sign to emphasize their equivalence.}
A measurement of the triple coincidence rate for an anti-Stokes detection each in channels 3 and 4 and a Stokes detection in either of the channels 2 or 1 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}. \textbf{VP: Yifan please convert that plot into a rate by dividing the counts by the measurement duration}.
\textcolor{red}{This expression duplicates the definition of triple coincidence.}

\textcolor{red}{The triple coincidence is shown as the event count rather than the count rate. There are several reasons for this:
1. Since we do not perform normalization for the triple coincidence, I prefer to show the event count directly as the result of the measurement.
2. The triplet coincidence aligns with the subsequent quadruplet coincidence.
3. This figure is plotted using Python's histogram function, which counts the number of triplet events in each bin.
4. This section focuses solely on temporal correlation. The rate would complicate the following discussion on the triplet rate. }

\textcolor{red}{It is too detailed to describe how the triple coincidence appears. According to BL's suggestions, we should represent the relationship using autocorrelation and cross-correlation}

The features in Fig. \ref{fig:triplet_w_corr} are understood from analyzing Eq.\ref{eq:TripletReduced} over various delays. For a coherence time $\Delta t$ for the Stokes/anti-Stokes photon, when $\tau_{3s}, \tau_{4s}, \tau_{34} \gg \Delta t$, $P(\tau_{3s},\tau_{4s},\tau_{34})$ reduces to the background accidental rate $R^3(0)$. 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 $\tau_{3s}, \tau_{4s}, \gg \Delta t$ but $\tau_{34} \lesssim \Delta t$. For such delays Eq.\ref{eq:TripletReduced} $P(\tau_{3s},\tau_{4s},\tau_{34}) \rightarrow R^3(0) g^{(2)}_{as,as}(\tau_{34})$. Thus the diagonal in Fig. \ref{fig:triplet_w_corr} is proportional to the intensity auto-correlation in the anti-Stokes mode. For $\tau_{3s}, \tau_{34} \gg \Delta t$ but $\tau_{4s} \lesssim \Delta t$, $P(\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 Stokes-anti-Stokes intensity cross-correlation. Similarly the vertical ridge is $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$. 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)$. \textcolor{red}{not correct}
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)$. \textcolor{red}{no need for a low accidental triplet rate. it is not exactly 4 times} This is a signature that the SFWM output is dominated by correlated double pairs generated from the same SFWM event. \textcolor{red}{The measured triplet consists of both correlated and accidental triplets. The descriptions above and below are worth examining.} 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$ 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. \textcolor{red}{In this paper, within a single pulse, the SPDC process (event) is well-defined, with no correlation between different pulses. In our case, however, the situation is different.} 

Due to the long coherence time of the Stokes and anti-Stokes photons, these measurements are not limited by averaging effects due to detector resolution. Thus we are able to better resolve the rich structures in the triple coincidences unlike previously reported measurements performed with broadband SPDC outputs \cite{razaviCharacterizingHeraldedSinglephoton2009}. Furthermore, we are also able to 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 measure triplet rate as compares to the peak outside the 20 ns Stokes-anti-Stokes correlation window. 

 \textcolor{red}{No observation supports ``strong''.``even with cw-pumping'' implies the cw-pumping is a special scenario. ``from the four-fold increase in the peak measure triplet rate as compares to the peak outside the 20 ns Stokes-anti-Stokes correlation window.''  The description of a four-fold increase in the peak measured triplet rate, compared to the peak outside the 20 ns Stokes-anti-Stokes correlation window'' is inaccurate. The maximum value is only four times the ridge value, and this ratio results from the autocorrelation and cross-correlation} 
 

 \textcolor{red}{The naming of the subsections is inconsistent}
\subsection{Four-fold coincidences}
\begin{figure}
\includegraphics[width=\linewidth]{four_photon_coind_dis.pdf}
 \caption{\textbf{VP:Caption to be worked on}(a) Temporal distribution of photon quadruplet events $(t_{1},t_{2},t_{3},t_{4})$ within a $60 \ \text{ns}$ time window, projected to the parameters delay time $\tau_{\text{3-1}}$ (between Ch3 and Ch1) and the delay time $\tau_{\text{4-1}}$ (between Ch4 and Ch2), over a measurement period of $0.7$ hour. Red square represents the $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}
 \textcolor{red}{``two photons in the anti-Stokes mode'' implies an additional case in which two photons are detected in the same detectors. But this is not involved in this paper}
 
We look for four-fold coincidences between detections of two photons in the anti-Stokes mode and two photons in the Stokes mode to understand the temporal structure of the double pairs, obtain the double pair generation rate and as an additional check to verify the direct generation of correlated double-pairs from our SFWM source.  \textcolor{red}{additional check?}
  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}
\begin{aligned}
 & Q(t_1, t_2, t_3, t_4) = \\
 & \langle \hat{a}^{\dagger}_{as}(t_{\text{4}})
 \hat{a}^{\dagger}_{as}(t_{\text{3}}) \hat{a}^{\dagger}_{s}(t_{2})\hat{a}^{\dagger}_{s}(t_{1})\hat{a}_{s}(t_{1})\hat{a}_{s}(t_{2})\hat{a}_{as}(t_{\text{3}})\hat{a}_{as}(t_{\text{4}}) \rangle.
\end{aligned}
\end{equation}
 \textcolor{red}{Q?}
 \textcolor{red}{The time definition is ambiguous with the following}
By replacing absolute detection times with relative delays, $Q$ can be reduced to a function of three parameters. 

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 coincidence window of 100 ns.  \textcolor{red}{Coincidence window is ambiguous with others. Probably use the 'analysis window.'}
For easier visualisation, we represent the results as a projection onto the plane \textcolor{red}{two-dimensional plot?} 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). We see an increase in quadruplets by several-folds for relative delays within 20 ns.
\textcolor{red}{This sentence is confusing because there is no clear counterpart, and it is unclear whether 'this bunching' refers to the maximal counts or the average counts.}

\textcolor{red}{The one-dimensional plot does not provide any correlation information; the only useful insight is that the count numbers are high, which can easily be seen in Figures 4a and 4b}

\textbf{VP: Yifan please include the 2-D representation of the histogram so readers can better see the shape}. \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?}






\appendix
\section{Double, Triple and Quadruple coincidences}

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}^\dagger_i(t+\tau_{ij}) \hat{E}_{j}(t)\rangle
\end{equation}
\textcolor{red}{Cross correlation is not defined with $ \langle \hat{E}^\dagger_i(t+\tau_{ij}) \hat{E}_{j}(t)\rangle$}
and the first-order autocorrelation
\begin{equation}
    R(\tau_{ij}) = \langle \hat{E}_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}
    \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})\}. \\
    \end{aligned}
\end{equation}
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}
    \begin{aligned}
        P_Q=&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}

\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 channels 1-2, 1-3, 1-4, 2-3, 2-4 and 3-4. \textcolor{red}{pair coincidences are only performed for 1-3,1-4, 2-3,2-4} We note the respective photon arrival timestamp data. \textcolor{red}{this sentence is not clear} The triplet coincidences are then identified based on the pair detections that share a photon arrival timestamp and quadruplet events are identified from triplet detections that share a common timestamp with a pair detection.

For instance, to identify triplet events involving channels 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 Ch1 (at $t_1^*$) and Ch2 (at time $t_2$), between Ch2 (at $t_2$) and Ch3 (at time $t_3^*$) and between between Ch2 (at $t_2$) and Ch3 (at time $t_4^*$)  with the previously identified triplet event with times $t_{1}, t_{3}, t_{4}$. If $t_{1}=t_{1}*$ or $t_{3}=t_{3}*$ or $t_{4}=t_{4}*$, the events are combined to form a quadruplet detection at times $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.
















\bibliography{reference}



\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{triplet_count.pdf}
    \caption{(a)~Photon triplet coincidence counts $N(\tau_{\text{3-s}}, \tau_{\text{4-s}})$, with paired events in Ch3 and Ch4, heralded by a photon event in Stoke mode. (b) Solid line: normalized cross-correlation function $g^{(2)}_{\text{as-s}}(\tau_{\text{3-s}})$ between Stokes and Ch3. Red dots: mean trace averaged over $\tau_{\text{4-s}}$ from 20 ns to 60 ns, normalized by the accidental triplet counts $R_{s}R_{3}R_{4}\delta t^{2}T_{m}$, where $\delta t =2 \ \text{ns}$ is the time bin and $T_{m}$ is the total measurement time. The offset above the cross-correlation function is due to diagonal contributions, indicating the auto-correlation of the anti-Stokes mode. (c) Vertical trace along vertical ridge $\tau_{\text{3-s}}= 8\ \text{ns}$, normalized by the far-detuned background. The maximal triplet coincidence count is approximately four times higher than the vertical ridge value.}
    \label{fig:triplet_w_corr}
\end{figure}


\begin{figure}
\includegraphics[width=\linewidth]{four_photon_coind_dis.pdf}
 \caption{(a) Temporal distribution of photon quadruplet events $(t_{1},t_{2},t_{3},t_{4})$ within a $60$ ns time window, projected to the delay times $\tau_{\text{3-1}}$ (between Ch3 and Ch1) and the delay time $\tau_{\text{4-1}}$ (between Ch4 and Ch2), over a measurement period of $0.7$ hour. The red square represents the selected $20$ ns coincidence window. (b) Same data plotted against the Stokes delay $\tau_{\text{1-2}}$ (between Ch1 and Ch2) and anti-Stokes  $\tau_{\text{3-4}}$ (between Ch3 and Ch4).}
    \label{fig:fourphoton_dis}
\end{figure}

\end{document}