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

\begin{document}
\title{Multiphoton state from four-wave mixing}
\author{Yifan~Li}
\author{Chang Hoong Chow}
\author{Boon Long Ng}
\author{Vindhiya Prakash}
\author{Justin Yu Xiang Peh}
\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 investigate the multiple photon state in a bright narrowband correlated photon source.
This source is based on a double-$\Lambda$ four-wave mixing process in a cold atomic ensemble. Our study indicates a super-bunching phenomenon for both Stokes and anti-Stokes photons in the heralded second-order auto-correlation functions. Within a coincidence window of $ 20 \unit{\ ns}$,  we have observed a photon quadruplet rate of $ (1.0\pm 0.2) \times 10^{3} \unit{\ h^{-1}}$, implying an initial rate of approximately $1.0 \times 10^{6} \unit{\ s^{-1}}$, alongside an initial pair generation rate of around $3.7 \times 10^{6} \unit{\ s^{-1}}$. These findings signify an enhancement in the probability of multiple photon states compared to the thermal properties in one of the twin modes, underscoring the potential utility of preparing multiple-photon states from the atomic ensemble. 
\end{abstract}

%the photon pair rate is measured at $8.2 \pm 0.3 \times 10^{3}\unit{\ s^{-1}}$,
\maketitle
\section{Introduction}
Nonclassical photon states represent promising quantum resources, encompassing both continuous-variable states, and discrete-variable states.  The continuous-variable states are exemplified by single-mode squeezed states \cite{teichSqueezedStateLight1989} and two-mode squeezed states \cite{teichSqueezedStateLight1989, dengGeneralizedTwomodeSqueezed1993}, while discrete-variable states include  Fock states, correlated photon pairs \cite{PhysRevLett.75.4337}, and entangled multiple-photon state like GHZ states \cite{greenbergerGoingBellTheorem1989, pan2000experimental}, cluster states \cite{RevModPhys.79.135}. Among these, correlated photon pairs have become foundational in the advancement of quantum technologies including quantum entanglement \cite{PhysRevLett.47.460,PhysRevLett.111.130406}, quantum communications \cite{bennettQuantumCryptographyPublic2014, xuSecureQuantumKey2020}, and quantum metrology \cite{RevModPhys.90.035005}.  Conventionally, the generation of correlated photon pairs relies on spontaneous parametric down-conversion (SPDC) processes in nonlinear crystals \cite{PhysRevA.64.023802, doi:10.1126/science.1127676} or spontaneous four-wave mixing (FWM) in atomic media \cite{lukinQuantumNoiseCorrelations1999,shuSubnaturallinewidthBiphotonsDopplerbroadened2016,parkSpectralTemporalBiphoton2020,kuzmichGenerationNonclassicalPhoton2003, kolchinGenerationNarrowBandwidthPaired2006, srivathsanNarrowBandSource2013, cereCharacterizationPhotonpairSource2018}. The non-degenerate photon state in entangled twin modes can be effectively modeled as two-mode squeezing \cite{sangouardQuantumRepeatersBased2011, gutierrezExperimentalFockStateSuperradiance2018, spasibkoSpectralStatisticalProperties2020}. Particularly, in low squeezing or low gain conditions where higher-order terms can be neglected, such two-mode photon states can be approximated as correlated photon pairs, with non-detectable vacuum states being disregarded. While the generation of correlated photon pairs is feasible,  achieving high generation rates for entangled multiphoton states remains challenging. One existing approach involves interference multiplexing of correlated photon pairs associated with post-selection \cite{zhangExperimentalGenerationHighfidelity2016,dongExperimentalRealizationNarrowband2017,parkGenerationBrightFourphoton2022}. Alternative solutions include the cascaded SPDC process \cite{hubelDirectGenerationPhoton2010} or directly applying the multiple-photon nonlinear processes \cite{PhysRevA.57.2076}. However, these methods result in low generation rates. By implementing the PDC process in the high-gain regime \cite{spasibkoSpectralStatisticalProperties2020, PhysRevA.93.043849}, it is possible to achieve specific multiple-pair photon states. Nonetheless, attaining reasonable generation probabilities necessitates a substantial increase in instantaneous pump power, often by orders of magnitude, compared to the conventional continuous-wave pump used for generating correlated photon pairs. Consequently, pulsed pumping becomes crucial to generate multiple photon pair states within the pulse duration from the higher-order terms in the PDC process. 

%These correlated multiple photon-pair states from higher-order terms in the PDC process are valuable as they can serve as the initial photon basis for preparing deterministic entangled multiple photon states or improve precision in quantum metrology.

This paper investigates multiple photon states generated from a bright correlated photon source based on a double-$\Lambda$ energy level scheme in a cold $^{87}$Rb atomic ensemble. The characterization correlated photon pairs from this bright source have been reported in Ref.\cite{}. To analyze the photon statistics, we performed the Hanbury-Brown and Twiss (HBT) measurement on the twin modes, denoted by Stokes and anti-Stokes, respectively. Due to the limitations of single-photon detectors, which lack photon number-resolving capability, we cannot directly obtain the full photon probability distribution \cite{PhysRevA.97.053835}. Instead, we analyze the photon number fluctuations via heralded auto-correlation functions and demonstrate the enhanced generation probability of photon quadruplet states. These observations indicate that this four-wave mixing process operates in a regime similar to the high-gain PDC with a relatively large mean photon number $\langle n \rangle$. In this regime, photons generated at the beginning of the atomic ensemble seed the generation at the end \cite{fox2006quantum,spasibkoSpectralStatisticalProperties2020}. As a counterpart to the high-gain case, in the low-gain limit, the mean photon number is much smaller than the effective brightness of zero-point vacuum fluctuations, and the generated photons do not enhance the generation rate \cite{spasibkoSpectralStatisticalProperties2020}. Achieving such a high-gain regime for near-resonant four-wave mixing in a cold atomic ensemble is feasible with modest pump power under a continuous-wave configuration, typically at the order of several milliwatts. Furthermore, the observations suggest a signature of superradiant emission \cite{dickeCoherenceSpontaneousRadiation1954, jenSuperradiantCascadeEmissions2015}. This phenomenon can occur in nonlinear process \cite{kuzmichGenerationNonclassicalPhoton2003,yoshimura2012dynamics,jenSuperradiantCascadeEmissions2015} and is characterized by super-bunching features in the second-order correlation function\cite{jahnke2016giant,temnovPhotonStatisticsCooperative2009,bhattiSuperbunchingNonclassicalityNew2015, cygorekSignaturesCooperativeEmission2023}.

In this experiment, the FWM process involves double-$\Lambda$ atomic energy levels driven by two pump fields denoted as $\omega_{p}$ and $\omega_{c}$, as shown in Fig.\ref{fig:setup} (a). In the presence of the pump fields, an effective interaction emerges between twined Stokes ($\omega_{s}$) and anti-Stokes modes ($\omega_{as}$), which is described by an effective Hamiltonian characterizing the nonlinear parametric process as 
\begin{equation}\label{eq:effective_interaction_hamiltonian}
\hat{H}_{s, as} = \hbar \kappa (\hat{a}^{\dagger}_{as}\hat{a}^{\dagger}_{s} + \hat{a}_{as}\hat{a}_{s}) \ .
\end{equation}
Here, $\kappa$ denotes the nonlinear parametric coupling coefficient, given by $\kappa = -i\sqrt{\omega_{as}\omega_{s}}/2 \chi^{(3)}E_{p} E_{c}$ where $\chi^{(3)}$ is the third-order nonlinear susceptibility of this parametric process. Analogous to the non-degenerate SPDC, the evolution of this Hamiltonian results in a two-mode squeezing on the initial vacuum state in twined modes. In the low squeezing limit, the photon state \cite{sangouardQuantumRepeatersBased2011} 
 is expressed as
\begin{equation}\label{eq:two_mode_squeezing_state}
(1-\frac{p}{2})|0, 0\rangle - i\sqrt{p}|1, 1\rangle - p |2, 2\rangle +\mathcal{O}(p^{\frac{3}{2}}) \ ,
\end{equation}
where the state $|n,n\rangle$ represents a photon state in the Stokes mode with $n$ photons and simultaneously $n$ photons in the corresponding anti-Stokes mode, and $\sqrt{p}$ represents the photon pair probability amplitude. $|0,0\rangle$ signifies the vacuum state. $\mathcal{O}(p^{3/2})$ indicates the higher-order excitations, whose probabilities are equal or smaller than $p^{3}$. This photon state suggests that the probability of double excitations is quadratically lower than that of single excitations in the low squeezing approximation. Our findings demonstrate a deviation of multiphoton probability in photon statistics.

\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{setup.pdf}
    \caption{\label{fig:setup}(a) The double-$\Lambda$ atomic energy levels involved in the four-wave mixing process. In this experiment, the pumping field detuning is set at $\Delta_{p} = -50 \text{\ MHz}$ while the coupling field is resonant. (b) A schematic representation of the experimental setup. The pumping beam and coupling beam counter-propagate along the longest axis of a cigar-shaped atomic ensemble, with a collection path angled by approximately $1$ degree. The collected Stokes and anti-Stokes photons are split by a 50:50 fiber-based beam splitter and detected by two single-photon detectors, respectively.}
\end{figure}
 
\section{Auto-correlation}\label{section:auto_correlation}
The second-order auto-correlation function $g^{(2)}(\tau)$ is a measure of temporal correlations in single mode\cite{PhysRev.131.2766}. After tracing over one of the modes in the two-mode squeezing, the remaining mode is left in a thermal state \cite{mollowQuantumTheoryParametric1967b} with photon bunching behavior $g^{(2)}(0)=2$. In Fig.\ref{fig:auto_correlation} (a) and (b), the unheralded auto-correlation functions on twin modes give the zero-delay value $g^{(2)}_{\text{s-s}}(0) = 1.95 \pm 0.12$ and $g^{(2)}_{\text{as-as}}(0) = 1.99 \pm 0.11$, which verifies the thermal photon statistics and Schmidt mode purity \cite{zielnickiJointSpectralCharacterization2018} in each mode of the correlated photon pairs. The heralding efficiency of correlated photon pairs is approximately $1.9\%$ for Stokes mode and $2.9\%$ for anti-Stokes mode. This indicates that the unheralded auto-correlation derived from coincidence counting is predominantly attributed to unpaired photon events. If the initial state is thermal light, the attenuation does not change the value of $g^{(2)}(0)$ which remains equal to 2. If the initial photon state is nonclassical, optical losses can significantly impact $g^{(2)}(0)$, degrading its nonclassical photon statistics \cite{TEICH19881}. To address this issue, we measure the second-order correlation function conditioned on the detection of one photon in the twinned mode. The definition of heralded auto-correlation functions was proposed with the post-selected averaging on every term in the auto-correlation functions \cite{razaviCharacterizingHeraldedSinglephoton2009, bocquillonCoherenceMeasuresHeralded2009}, yielding a normalization value dominated by the cross-correlations between twin modes. However, some criticisms \cite{bettelliCommentCoherenceMeasures2010, bashkanskySignificanceHeraldingSpontaneous2014} have pointed out its limitations on revealing correlation properties in a single mode. Referring to the analysis detailed in Appendix.\ref{app:former_heralded}, this quantity stands for purity of heralded single-photon source \cite{srivathsanHeraldedSinglePhotons2015a,guoTIMEFREQUENCYENTANGLEMENTNARROWBAND2018}, representing the conditional probability of detecting another photon within the coincidence window once a correlated photon pair is confirmed. 

% For example, the second-order correlation function for $as$ mode conditioned on $s$ mode is defined as 
% \begin{equation}
% g^{(2)}_{\text{h-as}}(t_{1}, t_{2}|t_{s}) = \frac{\langle \hat{a}^{\dagger}(t_{1})\hat{a}^{\dagger}(t_{2})\hat{a}_{as}(t_{2})\hat{a}(t_{2}) \rangle_{\text{ps}}}{\langle \hat{a}^{\dagger}(t_{1}) \hat{a}(t_{1})\rangle_{\text{ps}}\langle \hat{a}^{\dagger}(t_{2}) \hat{a}(t_{2})\rangle_{\text{ps}}} \ . 
% \end{equation}
% Here, $\langle : \rangle_{\text{ps}}$ denotes an average over events post-selected on the detection of a Stokes photon. 
Conversely, we obtain the heralded normalized auto-correlation function by using paired photon events in twin modes as the trigger, instead of all photon events in one of the split channels for unheralded auto-correlation. The methodology for this post-selection process is detailed in Appendix. \ref{app:heralded_coincidence}. As a result, this heralded auto-correlation function explores the local correlation relationship between paired photons and photons in the other split channel. Theoretically, this coincidence counting measurement for anti-Stokes mode corresponds to 
\begin{equation}\label{eq:redefined_heralded_auto}
 g^{(2)}_{\text{h-as}}(\tau = t_{2}-t_{1}) =  \frac{G^{(2)}(t_{1}, t_{2}| t_{s})}{G^{(1)}(t_{1})G^{(1)}(t_{2})} \ , 
\end{equation}
in which the timing notations $t_{s}$ and $t_{1}$ represent a measured photon pair event in Stokes mode and split channel 1, $t_{2}$ denotes the photon event in the other split channel within the neighboring time of this photon pair. A similar definition can also be applied to Stokes photons to obtain the correlation function $g^{(2)}_{\text{h-s}}(\tau^{\prime})$. Since Stokes photons precede anti-Stokes photons in the correlated pair, post-selection involves Stokes photons in reverse time after detecting an anti-Stokes photon. 
% \begin{equation}
% \begin{aligned}
% & g^{(3)}(t_{1}, t_{2}, t_{s}) = \frac{G^{(3)}(t_{1}, t_{2}, t_{s})}{G^{(1)}(t_{1})G^{(1)}(t_{2})G^{(1)}(t_{s})} \\
% & = \frac{\langle \hat{a}^{\dagger}(t_{s})\hat{a}^{\dagger}(t_{1})\hat{a}^{\dagger}(t_{2})\hat{a}(t_{2})\hat{a}(t_{1})\hat{a}(t_{s})\rangle}{\langle\hat{a}^{\dagger}(t_{s})\hat{a}(t_{s})\rangle\langle\hat{a}^{\dagger}(t_{1})\hat{a}(t_{1})\rangle\langle\hat{a}^{\dagger}(t_{2})\hat{a}(t_{2})\rangle} \ , 
% \end{aligned}
% \end{equation}
In these heralded auto-correlation functions, we observe the super-bunching features. As depicted in Fig.\ref{fig:heralded_auto} (a) and (b), Using the paired photon events within $20 \mathrm{\ ns}$ coherence time, the heralded normalized auto-correlation function for Stokes has a maximum $g^{(2)}_{\text{h-s}} = 3.41 \pm 0.32$ while the heralded auto-correlation function for anti-Stokes shows a maximum  $g^{(2)}_{\text{h-s}} = 3.38 \pm 0.33$. Fig.\ref{fig:heralded_auto} (c) and (d) show the heralded auto-correlations based on $20\text{-}40 \mathrm{\ ns}$ pair window. This indicates that paired photons within the second oscillation interval have a larger correlation with photons located in the first oscillation interval. According to the relationship between auto-correlation and photon fluctuation, expressed as $g^{(2)}(0) = 1 + (\text{Var}(n)-\langle n\rangle)/\langle n\rangle^2$, these zero-delay values larger than those of thermal light indicate higher photon number fluctuations within the coherence time for correlated photon pairs. 

Superradiance, which has been reported in the parametric process \cite{chaneliereQuantumTelecommunicationBased2006, jenPositivePhasespacemethodSimulation2012, srivathsanNarrowBandSource2013, gulatiPolarizationEntanglementQuantum2015, jenSuperradiantCascadeEmissions2015, PhysRevLett.122.143601}, predict the super-bunching feature with $g^{(2)}(0)>2$ \cite{temnovPhotonStatisticsCooperative2009,auffevesFewEmittersCavity2011,hoiGenerationNonclassicalMicrowave2012,grujicRepulsivelyInducedPhoton2013,bhattiSuperbunchingNonclassicalityNew2015,pagelNonclassicalLightFew2015,PhysRevApplied.4.044018}. From the perspective of single-mode squeezing,  this also indicates a nonclassical photon distribution with large photon number fluctuations, distinct from the Bose-Einstein distribution of thermal light \cite{kimPropertiesSqueezedNumber1989,2000LoudonQuantum}. However, unlike anti-bunching ($g^{(2)}(\tau)<0.5$), which is a clear criterion for nonclassical nature, super-bunching behavior can potentially be generated from classical mechanisms, like the gain competition in bimodal lasers \cite{PhysRevA.87.053819} and the constructive interference of multiple indistinguishable two-photon paths with thermal light \cite{jenPositivePhasespacemethodSimulation2012}. Therefore, photon coincidences at zero time delay alone are not sufficient to demonstrate the presence of superradiance\cite{cygorekSignaturesCooperativeEmission2023}. It is also necessary to examine the enhancement of the spontaneous radiative decay rate in cross-correlation \cite{}. 

% Unlike anti-bunching ($g^{(2)}(\tau)<0.5$), which is a criterion for nonclassical nature, super-bunching that photons tend to bunch together can potentially be explained by several mechanisms. Superradiant coupling in an ensemble of independent emitters results in super-bunching with $g^{(2)}(0)>2$ \cite{temnovPhotonStatisticsCooperative2009,auffevesFewEmittersCavity2011,hoiGenerationNonclassicalMicrowave2012,grujicRepulsivelyInducedPhoton2013,bhattiSuperbunchingNonclassicalityNew2015,pagelNonclassicalLightFew2015,PhysRevApplied.4.044018}. Single-mode squeezing also provides a nonclassical photon distribution distinct from the Bose-Einstein distribution of thermal light \cite{kimPropertiesSqueezedNumber1989,2000LoudonQuantum}. Additionally, the gain competition in bimodal lasers can lead to super-bunching light \cite{PhysRevA.87.053819}. Classical mechanisms, such as constructive interference of multiple indistinguishable two-photon paths with thermal light, can also generate super-bunching \cite{jenPositivePhasespacemethodSimulation2012}.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{auto_correlation_wo_fit.pdf}
    \caption{\label{fig:auto_correlation}Figures (a) and (b) illustrate the auto-correlation functions associated with the unheralded Stokes photon and anti-Stokes photon, respectively. This data originates from the identical dataset presented in Fig.\ref{fig:heralded_auto}. The mechanism responsible for the oscillation feature in (a) remains unidentified.}
\end{figure}
%In (a), the dashed grey line represents the Sinc function fit for the Stokes photons, yielding a coherence time of $8.41 \pm 0.54$ ns, whereas in (b), the dashed grey lines represent an exponential decay fit for the anti-Stokes photons with a coherence time of $8.98 \pm 0.36\mathrm{\ ns}$.
\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{heraldauto_correlation_wo_fit.pdf}
    \caption{\label{fig:heralded_auto}The heralded auto-correlation functions defined as Eq. \ref{eq:redefined_heralded_auto} for Stokes photons and Anti-Stokes photons are evaluated. In panels (a) and (b), we observe the heralded auto-correlation for Stokes and anti-Stokes photons, respectively, with selected pairs falling within a $20\mathrm{\ ns}$ coincidence window, corresponding to the first peak in the cross-correlation. In panels (c) and (d), the auto-correlation functions for Stokes and anti-Stokes photons are depicted for pairs occurring within a $20\text{-}40 \mathrm{\ ns}$ coincidence window, corresponding to the second peak in the cross-correlation. The delayed peak observed in Fig.(c) suggests a higher probability of detecting another Stokes photon after a certain delay relative to the Stokes photon paired in a $20\text{-}40 \mathrm{\ ns}$ coincidence window. Conversely, the early peak observed in Fig.(d) indicates an increased probability of detecting another anti-Stokes photon preceding the anti-Stokes photons paired within a $20\text{-}40 \mathrm{\ ns}$ coincidence window.}
\end{figure}

\section{Photon statistics}\label{section:photon_statistics}
\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{probability_w_fit.pdf}
    \caption{\label{fig:probability_vs_count_rate} (a) illustrates the ratios of photon pair rate, photon triplet rate, and photon quadruplet rate relative to the Stokes photon count rate. In contrast, (b) depicts the same parameters relative to the anti-Stokes photon count rate. Both axes are plotted on a logarithmic scale.}
\end{figure}
Referring to the photon-resolving analysis method of superradiant photon statistics as described in Ref.\cite{gutierrezExperimentalFockStateSuperradiance2018}, we investigate the rate of photon pair ($R_{2}$), photon triplet ($R_{3}$), and photon quadruplet ($R_{4}$) with respect to the single-photon rate. According to the model detailed in Appendix.\ref{app:detection_model}, these photon detection rates can be expressed as a linear combination of products of photon state probabilities and corresponding effective efficiency functions which is determined by the transmission, detection efficiency, and integration time. In Fig.\ref{fig:probability_vs_count_rate} (a) and (b), the ratios between multiple-photon detection rates and single-photon detection rates are plotted on a logarithmic scale relative to the detection rates of Stokes photons $R_{s}$ and anti-Stokes photons $R_{as}$, respectively. We observe that as the single photon detection rate increases, $R_{2}/R_{s}$ and $R_{2}/R_{as}$ exhibit two distinct plateaus. On a logarithmic scale, $R_{3}/R_{s}$ demonstrates an approximately linear relationship with $R_{s}$, characterized by a fitted slope parameter of $0.98 \pm 0.01$. Similarly, $R_{4}/R_{s}$ is proportional to $R_{s}$ with a fitted slope parameter of $1.06 \pm 0.06$. Compared to the anti-Stokes photon rate $R_{as}$, these fitted slopes are given as $0.81\pm 0.02$ for $R_{3}/R_{as}$ and $0.89 \pm 0.06$ for $R_{3}/R_{as}$. These observations confirm that the probability $p_{2}$ of the initial photon state $|2,2\rangle$ scales roughly quadratically with the probability $p_{1}$ of the initial photon states $|1,1\rangle$ as the single-photon probability. However, this method cannot resolve the gain factor $\beta$ in the relation $p_{2}=\beta p_{1}^{2}$. Identifying this gain factor necessitates the knowledge of effective functions. In the subsequent analysis, we incorporated the transmission ratio and detection efficiency to deduce the initial photon state generation rate and estimate this gain factor. 

\begin{figure}
\includegraphics[width=\linewidth]{triplet_coind40.pdf}
\caption{\label{fig:triplet_photon_distribution} Triplet photon coincidence within a coincidence window $0\text{-}80 \mathrm{\ ns}$. (a) The histogram shows the detected triplet events consisting of one Stokes photon and two anti-Stokes photons in Ch3 and Ch4 respectively. The x-axis denotes the time delay of photon events between Ch3 and Stokes mode, while the y-axis represents the time delay between CH4 and Stokes mode. The theoretical accidental count per bin is $R_{s}R_{\text{ch3}}R_{\text{ch4}}t_{\text{coi}}^{2}t_{\text{meas}}/ n^{2}\simeq418$ (b) The histogram exhibits the triplet events consisting of two Stokes photons and one anti-Stokes photon in Ch3 and Ch4 respectively. The theoretical accidental count per bin is $R_{as}R_{\text{ch1}}R_{\text{ch2}}t_{\text{coi}}^{2}t_{\text{meas}}/n^{2}\simeq 623$.}
\end{figure}

For an effective measurement duration of approximately 7 hours, with experimental settings of a pumping laser power of about $ 800 \mathrm{\ \mu W}$ and coupling laser power of roughly $10 \mathrm{\ mW}$, we acquired an average Stokes count rate of $(3.32 \pm 0.07)\times 10^{5} \mathrm{\ s^{-1}}$ and an average anti-Stokes photon count rate of $ (2.21 \pm 0.08) \times 10^{5} \mathrm{\ s^{-1}}$. Within a $20\mathrm{\ ns}$ coincidence window,  the photon pair counting rate was measured at $ (8138 \pm 290) \mathrm{\ s^{-1}}$ with heralding efficiencies $2\%$ for Stokes photons and  $2.9\%$ for anti-Stokes photons.  Fig.\ref{fig:triplet_photon_distribution} illustrates the triplet temporal distribution relative to their relative delays within a large $80\mathrm{\ ns}$ coincidence window. The majority of photon triplets are concentrated within a $20\mathrm{\ ns}$ coincidence window. With this coincidence window, the rate of photon triplets consisting of one Stokes photon and two anti-Stokes photons was recorded as $(28\pm 2)\mathrm{\ s^{-1}}$, while the rate for triplets consisting of two Stokes photons and one anti-Stokes photon was $ (41\pm 3)\mathrm{\ s^{-1}}$. Furthermore, a four-photon event was recorded at a rate of $(1044 \pm 216)\mathrm{\ h^{-1}}$. Fig.\ref{fig:fourphoton_dis} shows the temporal distribution of detected four-photon events for a $80\mathrm{\ ns}$ coincidence window, revealing an oscillation pattern inherited from second-order cross-correlation. The total transmission ratio for Stokes channels is approximately $0.13$ with quantum efficiencies $70\%$ for two Stokes single-photon detectors.  For the anti-Stokes channels, the transmission ratio is $0.12$, with quantum efficiencies of $50 \%$ for anti-Stokes single-photon detectors. Resorting to this knowledge of transmission ratios in collection paths and quantum efficiencies, the analysis detailed in Appendix.\ref{app:initial_state} reconstructed the effective probability functions for photon event detection. This enables us to infer the initial photon state from the atomic ensemble within the collection solid angle in the experimental setup. The analysis indicates that, in our collection setup, the photon pair generation rate out of the atomic cloud reaches $1 \times 10^{6}\mathrm{\ s^{-1}}$ with expected heralding efficiencies around $27\%$ for both Stokes channel and anti-Stokes channels. Simultaneously, the photon quadruplet generation rate was determined to be about $4 \times 10^{4} \mathrm{\ s^{-1}}$, accounting for roughly $1\%$ of all Stokes photon events and anti-Stokes photon events. The residual three-photon event generation rates are calculated as approximately $9 \times10^{4}\mathrm{\ s^{-1}}$ for two Stokes photon and one anti-Stokes and $1 \times10^{5}\mathrm{\ s^{-1}}$ for one Stokes photon and two anti-Stokes photons. We further assume that all photon triplet states originate from the photon quadruplet state and are caused by photon loss in the atomic ensemble and the collection path. Under these assumptions, we estimate the initial photon quadruplet generation rate to be up to around $9.5 \times 10^{5} \mathrm{\ s^{-1}}$, alongside a pure photon pair generation rate of $3.7 \times 10^{6} \mathrm{\ s^{-1}}$, with estimated loss coefficients of approximately $0.56$ for Stokes photons and $0.54$ for anti-Stokes photons. Thus, the initial photon quadruplet generation rate is nearly of the same order of magnitude as the photon pair generation rate. 

\begin{figure}
\includegraphics[width=\linewidth]{four_photon_coind40.pdf}
\caption{\label{fig:fourphoton_dis}(a) The green histogram illustrates the four-fold photon event counts within an $80\mathrm{\ ns}$ coincidence window over an observation period of $7$ hours. The dashed line represents a theoretical estimation of approximately $110$ accidental four-fold photon events per $2\mathrm{\ ns}$ bin during the same measurement duration. The accidental four-fold photon detection ratio ($R_{\mathrm{acc4}}$) is determined by $R_{1}R_{2}R_{3}R_{4} t_{\mathrm{coi}}^{3}$, where $R_{i}$ denotes the detection rate in each single-photon detector, and $t_{\mathrm{coi}}$ is the coincidence window. The grey histogram presents an example of accidental four-fold photon event counts from an uncorrelated coincidence window. (b) This 2D histogram illustrates the spatial distribution of four-photon event counts against the relative delay between two Stokes photons and the delay between two anti-Stokes photons in a single four-fold photon event. (c) This displays the four-fold event counts versus the delay time of the anti-Stokes photon after the Stokes photon. Ch1 and Ch2 are two detector channels for Stokes photons, while Ch3 and Ch4 are two channels for anti-Stokes photons.}
\end{figure}

Given the assumption that the generation of correlated photons in the atomic medium follows a Poisson process, the probabilities associated with specific photon states, such as $|1,1\rangle$ and $|2,2\rangle$, can be formulated with respect to a given time interval $\Delta t$. The estimated probability $\bar{p}_{1}$ of observing the state $|1,1\rangle$ within $\Delta t$ is expressed as $\bar{p}_{1} = R_{2}\Delta te^{-R_{2}\Delta t}$, while the  probability $\bar{p}_{2}$ of detecting the state $|2,2\rangle$ in the same interval is given by $\bar{p}_{2} = R_{4}\Delta te^{-R_{4}\Delta t}$. These probabilities are influenced by the selection of the time window. Within a defined $20 \mathrm{\ ns}$ coincidence window, the probability model yields a photon pair rate of $\bar{p}_{1}\simeq 0.07$ and photon quadruplet rate of $\bar{p}_{2}\simeq 0.019 $. The photon quadruplet probability, $\bar{p}_{2}$, is enhanced by a factor of approximately $3.85$ compared to the squared one-photon probability, $\bar{p}_{1}^{2}$. Furthermore, it exceeds the predicted value of $P(2)\simeq 0.0044$ for a photon state with a mean photon number $\langle n\rangle \simeq 0.074$ based on photon number distribution for thermal light. According to the Bose-Einstein distribution, the probability of observing n photons is given by $P(n)=1/(\langle n \rangle + 1)(\langle n \rangle /(\langle n \rangle + 1))^{n}$.  Although the complete photon number distribution cannot be determined due to the lack of photon-number-resolving detectors, the observed ratio between the photon pair and photon quadruplet probabilities indicates a higher likelihood of multiple photon events than predicted by the Bose-Einstein distribution. This deviation suggests an enhancement in the occurrence of multi-photon events.

\section{Conclusion}\label{section:conclusion}
In conclusion, besides generating correlated photon pairs, the FWM process in a cold atomic ensemble can also produce correlated four-photon states with a decent generation rate. Our investigation focused on these four-photon states. The correlated photons were collected from the phase-matching mode in the paraxial regime of the counterpropagating pump fields. Using the HBT setup, we examined the auto-correlation, heralded auto-correlation from the photon coincidence counting measurement. The enhanced probability of two same-mode photons within a zero-delay time bin can be verified from heralded auto-correlation. For the photon quadruplet state $|2,2\rangle$, in addition to the quadratic relation with the single-pair probability, a linear gain is also observed. This photon statistics deduction reveals that the initial photon states exhibit a high generation rate of four-photon states, making this approach competitive with previous high-photon-number state generation protocols. Given the potential for optimizing the collection and filter system, as well as enhancing the optical depth of the atomic cloud, a higher detection rate of multiple-photon states can be anticipated. Such our-photon states hold significant potential for a variety of applications, These include the generation of heralding higher-order Fock states\cite{PhysRevA.100.04180}, implementation of multiparty teleportation \cite{doi:10.1080/09500340008244039}, improvement of the phase sensitivity in interferometry experiments \cite{HashemiRafsanjani:17}, detecting sub-wavelength interference \cite{zhai2014direct} and enhanced visibility in ghost imaging \cite{alleviSuperthermalLightImaging2017}. 

\bibliography{reference}

\appendix
\section{Coincidence histogram}
The unheralded auto-correlation function is obtained from a coincidence counting histogram in the Hanbury-Brown and Twiss (HBT) measurement, The collected photons pass through a 50:50 beam splitter, which directs them into two separate arms, each equipped with a single-photon detector. The detection rates in these arms are denoted as $R_{1}$ for arm 1 and $R_{2}$ for arm 2. Coincidence counting uses arm 1 as the trigger, recording photon events in arm 2 before and after the trigger within a relatively long time window. This implies that the coincidence histogram is a superposition of $N_{1} = R_{1}T_{m}$ copies of individual histograms conditioned on a photon event in arm 1. 

In the far delay regime $\tau = t_{2} - t_{1} \gg \tau_{c}$ (where $\tau$ is the relative delay between two-photon events in the two arms and $\tau_{c}$ is the coherence time), the probability of another photon event existing within a time bin $\Delta t$ is independent of $R_{1}$, originating from accidental background counts in arm 2. Therefore, assuming that the photon emission follows a Poisson process, the probability per bin $p_{b}$ is given by:
\begin{equation}
p_{b} = R_{2}\Delta t e^{- R_{2}\Delta t} \approx R_{2}\Delta t
\end{equation}
This approximation holds as $\Delta t = 2$ ns for the minimal time bin in our timestamp and $R_{2}\Delta t \ll 1$. Hence, the average background count per bin in the far delay regime is:
\begin{equation}
\bar{N}_{b} = p_{b}N_{1} \approx R_{1}R_{2}\Delta t T_{m}
\end{equation}
Experimentally, the background count serves as the normalization factor to obtain the auto-correlation function $g^{(2)}(\tau)$, corresponding to the theoretical normalization value $\langle a^{\dagger}(t+\tau)a(t+\tau)\rangle \langle a^{\dagger}(t)a(t)\rangle$ with a constant time value $\Delta t T_{m}$. This method allows us to infer the auto-correlation function by analyzing the coincidence histogram in the Hanbury Brown and Twiss (HBT) setup.

The zero-delay count $N_{\tau=0}$ within the same measurement time $T_{m}$ and time bin $\Delta t$ can be expressed by the conditional coincidence probability $p_{\tau=0}$ as:
\begin{equation}
N_{\tau=0} = p_{\tau=0} N_{1}
\end{equation}
The thermal properties yield $g^{(2)}(\tau =0) = 2$, hence the zero-delay count $N_{\tau=0}$ within the zero-delay bin is written in terms of the average background count per bin in the far delay regime as ${N_{\tau=0}}/{\bar{N}_{b}} = 2$. Therefore, the zero-delay conditional coincidence probability is given by:
\begin{equation}
p_{\tau =0} = 2 R_{2}\Delta t
\end{equation}
\section{Heralded coincidence histogram}\label{app:heralded_coincidence}
Considering the heralding case involving the photon rate in the Stokes channel denoted as $R_{s}$, we can record the pair rate $C_{s1}$ between this channel and arm 1 for the anti-Stokes channel. Post-selection procedures account only for photon events in arm 1 paired with the Stokes mode. Once the pair is confirmed, the photon events in arm $2$ before and after this photon event in arm $1$ will accumulate to the whole histogram. The coincidence histogram is the sum of $N_{1}^{h} = C_{s1}T_{m}$ individual counting histograms. Even when arm 1 photon events are post-selected, the probability per bin in the far delay regime in arm 2 is independent of $R_{s}$, $R_{1}$, or $C_{si}$, with $p_{b} \approx R_{2}\Delta t$ due to the absence of correlation between them. Moreover, it is not necessary to post-select photon events in arm 2, as the photon counting histogram inherently records only those photon events in arm 2 that are close to the paired photon events in arm 1. Consequently, the average background count per bin in the far delay regime is
\begin{equation}
\bar{N}_{b}^{h} = p_{b}N_{1}^{h} = p_{b}C_{s1}T_{m}
\end{equation}
At the zero-delay bin, if the photon event number $N_{\tau=0}^{h} = p_{\tau=0}^{h}N_{1}^{h}$ is measured in relation to the far-delay average background count as 
\begin{equation}
N_{\tau=0}^{h} / \bar{N}_{b}^{h} = k 
\end{equation}
Then, we obtain 
\begin{equation}
p_{\tau=0}^{h} = k R_{2}\Delta t =\frac{k}{2} p_{\tau=0}
\end{equation}
Here, $p_{\tau=0}$ refers to the probability of two photons within zero-delay time bin $\Delta t$ for all photons in the same mode, while $p_{\tau=0}^{h}$ refers to the probability of two photons existing within zero-delay time bin $\Delta t$ for paired photons. This analysis indicates that super-bunching with $g^{(2)}(0)$ represents the enhancement of conditional coincidence probability. 

\section{Former definition of heralded auto-correlation function}\label{app:former_heralded}
Previously, the second-order correlation function for one arm conditioned on the observation of one photon in the other arm is proposed with an averaging procedure over post-selection, as described in previous works \cite{razaviCharacterizingHeraldedSinglephoton2009, bocquillonCoherenceMeasuresHeralded2009,hockelDirectMeasurementHeralded2011}.  For example, the second-order correlation function for ``$as$'' mode conditioned on ``$s$'' mode is defined as 
\begin{equation}
g^{(2)}_{\text{ps}}(t_{1}, t_{2}|t_{s}) = \frac{\langle \hat{a}_{s}^{\dagger}(t_{1})\hat{a}_{as}^{\dagger}(t_{2})\hat{a}_{as}(t_{2})\hat{a}_{s}(t_{2}) \rangle_{\text{ps}}}{\langle \hat{a}_{s}^{\dagger}(t_{1}) \hat{a}_{s}(t_{1})\rangle_{\text{ps}}\langle \hat{a}_{as}^{\dagger}(t_{2}) \hat{a}_{as}(t_{2})\rangle_{\text{ps}}} \ . 
\end{equation}
Here, $\langle : \rangle_{\text{ps}}$ denotes an average over events post-selected on the detection of a Stokes photon. Analog to the unheralded auto-correlation, $g^{(2)}(0)=2p(2_{as})/p^{2}(1_{as})$, this quantity can be expressed by the conditional probability \cite{Panyukov:22} as 
\begin{equation}
g^{(2)}_{\text{ps}}(0)= \frac{2 p(2_{as}|1_{s})}{p^{2}(1_{as}|1_{s})}
\end{equation}
where $p(n_{as}|1_{s})$ denotes the probability of detecting one Stokes photon and at the subsequent coincidence window the Stokes mode containing $n_{as}$ photons. The conditional probabilities are related to the photon state probabilities using Bayes's theorems as 
\begin{equation}
p(n_{as}|1_{s}) = \frac{p(n_{as}, 1_{s})}{p_{s}} = \frac{p(n_{as}, n_{s})}{\sum_{n} p(n_{s}, n_{as})} \approx \frac{p(n_{as}, n_{s})}{ p(1_{as}, 1_{s})} 
\end{equation}
in which $p(n_{as}, n_{s})$ is the joint probabilities of $n_{as}$ photons in anti-Stokes mode and $n_{s}$ photon in the Stokes mode, also it represents the probabilities of photon state $|n_{s}, n_{as}\rangle$. The last approximation refers to the case where the single-pair probability is dominant compared to multiple-pair states. 

If we consider a two-mode squeezing state $\sqrt{p}|1,1\rangle + p|2,2\rangle$, neglecting the vacuum and the $n>2$ higher terms, this post-selected normalized $g_{\text{ps}}^{(2)}$ is given as 
\begin{equation}
\begin{aligned}
p(2_{as}|1_{s}) &\approx p^{2}/(p+2p^{2}) \\
p(1_{as}|1_{s}) &\approx p/(p+2p^{2}) \\
g^{(2)}_{\text{ps}}(0) &\approx 2p + 4p^{2}    
\end{aligned}
\end{equation}
It is easy to observe that this post-selected quantity is determined by the single-pair probability $p$. In the low squeezing limit where $p\ll1$, it is clear that $g^{(2)}_{\text{ps}}$ consistently exhibits $g^{(2)}_{\text{ps}}\ll 1$. This example underscores that $g^{(2)}_{\text{ps}}$ does not function as a correlation function, instead, it signifies the purity of the heralded single photon. Specifically, if a pair of photons in the twin modes is detected, $g^{(2)}_{\text{ps}}$ denotes the conditional probability of another photon's existence within the same coincidence window in the anti-Stokes mode. Therefore, we can predict that the property $g^{(2)}_{\text{ps}}(0)<1$ applies to the photon pairs generated from SPDC and FWM, given that the condition $p\ll 1$ is satisfied for the minimal time bin. This implies that with a lower photon generation rate, a deeper dip in $g^{(2)}_{\text{ps}}(\tau)$ can be achieved. This dip indicates strong nonclassical cross-correlation between two modes rather than sequential photons within a single mode \cite{bettelliCommentCoherenceMeasures2010, hockelDirectMeasurementHeralded2011,bashkanskySignificanceHeraldingSpontaneous2014}. Consequently, any correlated photon pairs exhibiting strong cross-correlation and low generation probability can be considered as a heralded single-photon source. This post-selected quantity $g^{(2)}_{\text{ps}}$ can be used to evaluate the quality of the heralded single-photon. 

\section{Multiple photon detection model}\label{app:detection_model}
Assuming that all detected single photons, correlated photon pairs, photon triplets, and photon quadruplets originate from a quantum photon state represented by $\sqrt{p_{1}}|1,1\rangle + \sqrt{p_{2}}|2,2\rangle$, resembling the form of a two-mode squeezed state neglecting vacuum and higher-order terms. Here, $p_{1}$ denotes the probability of the state $|1,1\rangle$, and $p_{2}$ represents the probability of $|2,2\rangle$. In the approximation of a two-mode squeezed state, we can anticipate the relationship $p_{2} = p_{1}^{2}$.

Furthermore, we assume that all photon triplet events originate solely from the state $|2,2\rangle$, rather than from a distinct process generating three photons, and the detectors do not have photon-number resolving ability. The rates of single-photon generation $R_{s}$ and $R_{as}$, photon pair generation $R_{2}$, photon triplet generation $R_{3}$, and photon quadruple generation $R_{4}$ are proportional to the corresponding state probabilities as follows:
\begin{equation}\label{eq:event_rate}
\begin{aligned}
& R_{s} \propto \eta_{s} (p_{1} + p_{2})  \\
& R_{as} \propto \eta_{as} (p_{1} + p_{2})  \\
& R_{2} \propto \eta_{2} (\eta_{s}, \eta_{as}) (p_{1} + p_{2})\\
& R_{3} \propto \eta_{3} (\eta_{s}, \eta_{as}) p_{2}   \\
& R_{4} \propto \eta_{4} (\eta_{s}, \eta_{as}) p_{2}  \\
\end{aligned}
\end{equation}
Under the approximation of $p_{2} \ll p_{1}$, we have 
\begin{equation}\label{eq:approx_event_rate}
\begin{aligned}
& R_{s} \propto  \eta_{s} p_{1} \\
& R_{as} \propto \eta_{s} p_{1} \\
& R_{2} \propto \eta_{2} (\eta_{s}, \eta_{as}) p_{1} \\
\end{aligned}
\end{equation}
where $\eta_{i}$ (for $i = s, as$) represents an effective function determined by transmission, coupling efficiency, integration time, and detection efficiency in the corresponding channels. Similarly, $\eta_{i}$ (for $i = 2, 3, 4$) is a combination of $\eta_{s}$ and $\eta_{as}$, dependent on the detection configuration. These effective probability functions depend solely on the experimental setup and detection configuration, independent of the initial photon state probabilities $p_{1}$ and $p_{2}$.

Therefore, the ratio between the correlated photon pair rate and the Stokes single-photon rate $R_{2}/R_{s}$ is a constant relative to the $R_{s}$ as 
\begin{equation}\label{eq:theo_pair_rate_vs_single_rate}
\frac{R_{2}}{R_{s}} = \frac{\eta_{2}(\eta_{s}, \eta_{as})}{\eta_{s}} = \text{const}
\end{equation}
The ratio of photon triplet rate and Stokes single-photon rate is given as 
\begin{equation}
\frac{R_{3}}{R_{s}} = \frac{\eta_{3} p_{2}}{\eta_{s}(p_{1} + p_{2})}
\end{equation}
If $p_{2}$ has a quadratic relation with $p_{1}$ as $p_{2} = \beta p_{1}^{2}$ with a gain factor $\beta$ and we obtain 
\begin{equation}
\begin{aligned}
\frac{R_{3}}{R_{s}} & = \frac{\eta_{3} p_{2}}{\eta_{s}(p_{1} + p_{2})} \\
&\simeq \frac{\eta_{3}\beta p_{1}}{\eta_{s}} \
\end{aligned}
\end{equation}
Therefore $R_{3}/R_{s}$ is proportional to $R_{s}$, and in logarithmic scale, it will show a linear relation with slope $1$.  

The photon quadruple rate, defined by Eq. \ref{eq:event_rate}, shares the same structure as the photon triplet rate but with a different efficiency function $\eta_{4}$. Consequently, the conclusions drawn from the photon triplet analysis can also apply to the photon quadruple rate. Furthermore, extending this analysis to ratios relative to the anti-Stokes single-photon rate is feasible by substituting $R_{s}$ with $R_{as}$ and replacing $\eta_{s}$ with the corresponding detection efficiency parameter $\eta_{as}$.

\section{Collection transmission and quantum efficiencies}
The total collection efficiencies (total transmission) are determined by losses in the collection channels, primarily due to imperfect coupling between the collected Gaussian spatial mode and the single-mode fiber, transmission loss in the Etalon systems, and attenuations in the single-mode fibers and fiber adaptors. These losses are assessed by measuring the attenuation of a laser beam in these systems. Therefore, the total transmissions used in the initial state deduction represent upper bounds for the real transmission condition. Due to the small collection solid angle, a portion of correlated photon pairs with opposite wave vectors may not overlap well with the collected spatial mode, resulting in a significant chance that only one of the correlated photon pairs is collected. These unpaired single photons in the collection spatial mode contribute substantially to the noise single-photon count rates $S_{s}^{(i)}$ and $S_{as}^{(i)}$ for each mode. This also indicates that the model in the following appendix provides an estimate for the worst case. The actual initial photon quadruplet rate under wider phase-matching conditions is expected to be higher than the evaluation.

The two single-photon detectors in Stokes channels are fiber-coupled Excelitas detectors with a nominal quantum efficiency of approximately $70 \%$, while the two single-photon detectors in anti-Stokes channels are free-space detectors with homemade coupling systems, resulting in estimated quantum efficiencies for anti-Stokes photons of $50 \%$. 

\section{Initial state deduction}\label{app:initial_state}
Firstly, we model the loss and absorption of photons in the atomic ensemble.  Here, we designate the loss probability for Stokes photons as $A_{s}$ and the loss probability for anti-Stokes photons as $A_{a}$.
The residual four-photon rate is given as
\begin{equation}
R_{4}^{(i)} = R_{4}(1-A_{s})^{2}(1-A_{a})^{2}
\end{equation}
Owing to the loss of single photons, a fraction of the four-photon states undergo the loss of a single photon into a three-photon state with the following rates:
\begin{equation}
\begin{split}
&R_{3a}^{(i)} = 2R_{4}(1-A_{s})A_{s}(1-A_{a})^2 \\
&R_{3s}^{(i)} = 2R_{4}(1-A_{s})^{2}A_{a}(1-A_{a})
\end{split}
\end{equation}
Here, the subscript $3a$ signifies one Stokes and two anti-Stokes photons, while $3s$ represents two Stokes and one anti-Stokes photon. The Stokes anti-Stokes pair originates from two sources: one from the original two-photon pairs, and the other contributed by the loss of photon in photon quadruples.
\begin{equation}
R_{2}^{(i)} = R_{2}(1-A_{s})(1-A_{a}) + 4R_{4}(1-A_{s})A_{s}(1-A_{a})A_{a}
\end{equation}
% The two photons remain in the same wavelength generation rates are given as 
% \begin{equation}
%   \begin{split}
%       R_{2s}^{(i)} = & R_{4}(1-A_{s})^{2} A_{a}^{2} \\ 
%       R_{2a}^{(i)} = & R_{4}A_{s}^{2}(1-A_{a})^{2}
%   \end{split}
% \end{equation}
Incoherent scattering rates for Stokes and anti-Stokes are expressed as follows:
\begin{equation}
\begin{split}
S_{s}^{(i)} &= S_{s}(1-A_{s}) + R_{2}(1-A_{s})A_{a} + 2R_{4}(1-A_{s})A_{s}A_{a}^{2} \\
S_{a}^{(i)} &= S_{a}(1-A_{a}) + R_{2}A_{s}(1-A_{a}) + 2R_{4}A_{s}^{2}A_{a}(1-A_{a})
\end{split}
\end{equation}
The total loss in the collection channels encompasses spatial collection loss, total fiber coupling loss, Etalon filter loss, and fiber-based beam splitter loss. Correspondingly, the transmission rate in the Stokes channel is denoted as $T_{s}$, while the transmission in the anti-Stokes channel is denoted as $T_{a}$. We utilize the superscript $(t)$ to indicate the transmitted event rate. The rates for different photon states are expressed as follows:
\begin{enumerate}
    \item The photon quadruplets
\begin{equation}
R_{4}^{(t)} = R_{4}^{(i)}T_{s}^{2}T_{a}^{2}
\end{equation}
\item The photon triplets($3a$ for two anti-Stokes photons and $3s$ for two Stokes photons)
\begin{equation}
    \begin{split}
        R_{3s}^{(t)} &= R_{3s}^{(i)}T_{s}^{2}T_{a} + 2R_{4}^{(i)}T_{s}^{2}T_{a}(1-T_{a}) \\
        R_{3a}^{(t)} &= R_{3a}^{(i)}T_{s}T_{a}^{2} + 2R_{4}^{(i)}T_{s}(1-T_{s})T_{a}^{2}
    \end{split}
\end{equation}
\item The photon pairs 
\begin{equation}
\begin{split}
    R_{2}^{(t)} =  & R_{2}^{(i)}T_{s}T_{a}  
     + 2R_{3a}^{(i)}T_{s}T_{a}(1-T_{a}) \\
    & + 2R_{3s}^{(i)}T_{s}(1-T_{s})T_{a} \\
    & + 4R_{4}^{(i)}T_{s}(1-T_{s})T_{a}(1-T_{a})
\end{split}
\end{equation}
\item The two-photon both in Stokes or anti-Stokes 
\begin{equation}
\begin{split}
    R_{2s}^{(t)} = & R_{3s}^{(i)}T_{s}^{2}(1-T_{a})  + R_{4}^{(i)}T_{s}^{2}(1-T_{a})^{2} \\ 
    R_{2a}^{(t)} = & R_{3a}^{(i)}(1-T_{s})T_{a}^{2}  + R_{4}^{(i)}(1-T_{s})^{2}T_{a}^{2} \\ 
\end{split}
\end{equation}
\item Single photon scattering 
\begin{equation}
\begin{split}
S_{s}^{(t)} = & R_{2}^{(i)}T_{s}(1-T_{a}) \\ 
                  & + R_{3a}^{(i)}T_{s}(1-T_{a})^{2} + 2R_{3s}^{(i)}T_{s}(1-T_{s})(1-T_{a}) \\ 
                  & + 2R_{4}^{(i)}T_{s}(1-T_{s})(1-T_{a})^{2} \\ 
                  & + S_{s}^{(i)}T_{s} \\
S_{a}^{(t)} = & R_{2}^{(i)}(1-T_{s})T_{a} \\ 
                  & + 2R_{3a}^{(i)}(1-T_{s})T_{a}(1-T_{a}) + R_{3s}^{(i)}(1-T_{s})^{2}T_{a} \\ 
                  & + 2R_{4}^{(i)}(1-T_{s})^{2}T_{a}(1-T_{a})\\ 
                  & + S_{a}^{(i)}T_{a}
\end{split}
\end{equation}
\end{enumerate}
Having accounted for the entire loss in both $L_{s}$ and $L_{a}$, we describe the likelihood of a photon traversing a specific arm of the beam splitter and reaching the corresponding detection channels as $P_{1}, P_{2}, P_{3}, P_{4}$, respectively. These probabilities adhere to the conditions $P_{1} + P_{2} = 1$ and $P_{3} + P_{4} = 1$, where channel 1 and channel 2 represent the two arms for Stokes photons, while channel 3 and channel 4 denote the two outputs for anti-Stokes photons. The probability distribution for a photon cluster state across four channels is expressed as follows:
\begin{equation}
\begin{aligned}
P(n, a-n, m, b-m) =  \\ C(a, n)C(b, m)P_{1}^{n}P_{2}^{a-n}P_{3}^{m}P_{4}^{a-m} 
\end{aligned}
\end{equation}
Here, the variable ``a" denotes the number of Stokes photons within a given photon state, with possible values of 0, 1, or 2. Similarly, "b" represents the count of anti-Stokes photons in the same photon state, with potential values of 0, 1, or 2. 

In the HBT setup, photons in four channels are detected using single photon detectors with detection efficiencies  $Q_{s}$ for Stokes photon and $Q_{a}$ for anti-Stokes photon. Therefore, we can express the photon counting measurement rates in terms of the initial photon state rates. 
\begin{widetext}
\begin{equation}  R^{(m)}_{\mathrm{4\text{-}photon}} = R_{4}^{(t)}P_{1111}Q_{s}^{2}Q_{a}^{2}
\end{equation}
\begin{equation} 
\begin{aligned}
  R^{(m)}_{\mathrm{1S,2aS}} = 
  R_{4}^{(t)}[ 2P_{1111}Q_{s}(1-Q_{s})Q_{a}^{2} +  (P_{2011} + P_{0211})Q_{s}Q_{a}^{2}] 
  + R_{3a}^{(t)}(P_{1011} + P_{0111})Q_{s}Q_{a}^{2}  
\end{aligned}
\end{equation}

\begin{equation}
\begin{aligned}
 R^{(m)}_{\mathrm{2S,1aS}} = 
 R_{4}^{(t)} [2P_{1111}Q_{s}^{2}(1-Q_{a})Q_{a} + (P_{1120} + P_{1102})Q_{s}^{2}Q_{a}] 
+ R_{3s}^{(t)}(P_{1110} + P_{1101})Q_{s}^{2}Q_{a} 
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
R^{(m)}_{\mathrm{2}}=  & R_{2}^{(t)}Q_{s}Q_{a} \\
& + R_{3a}^{(t)}[2(P_{1011}+P_{0111})Q_{s}Q_{a}(1-Q_{a}) + (P_{1020}+P_{1002}+P_{0120}+P_{0102})Q_{s}Q_{a}]\\ 
& + R_{3s}^{(t)}[2(P_{1110}+P_{1101})Q_{s}(1-Q_{s})Q_{a} + (P_{2010}+P_{0210}+P_{2001}+P_{0201})Q_{s}Q_{a}]\\
& + R_{4}^{(t)}[4P_{1111}Q_{s}Q_{a}(1-Q_{s})(1-Q_{a}) + 2(P_{2011} + P_{0211})Q_{s}Q_{a}(1-Q_{a})  \\ 
& + 2(P_{1120} + P_{1102})Q_{s}(1-Q_{s})Q_{a} + (P_{2020}+P_{2002}+P_{0220}+P_{0202})Q_{s}Q_{a}]\\
\end{aligned}
\end{equation}
\end{widetext}
\end{document}
