# Bidimensional measurements of photon statistics within a multimodal temporal framework

C. Hainaut, K. Ouahrouche, A. Rancon, G. Patera, C. Ouarkoub, M. Le Parquier, P. Suret, A. Amo

► To cite this version:

C. Hainaut, K. Ouahrouche, A. Rancon, G. Patera, C. Ouarkoub, et al., Bidimensional measurements of photon statistics within a multimodal temporal framework. 2025. [hal-05397878](https://hal.science/hal-05397878v1)

HAL Id: hal-05397878

https://hal.science/hal-05397878v1

Preprint submitted on 4 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Bidimensional measurements of photon statistics within a multimodal temporal framework

C. Hainaut,¹,* K. Ouahrouche,¹ A. Rançon,¹,² G. Patera,¹ C. Ouarkoub,¹ M. Le Parquier,¹ P. Suret,¹ and A. Amo¹

¹Univ. Lille, CNRS, UMR 8523–PhLAM–Physique des Lasers, Atomes et Molécules, Lille, France

²Institut Universitaire de France

(Dated: November 17, 2025)

Ultrafast imaging of photon statistics in two dimensions is a powerful tool for probing non-equilibrium and transient optical phenomena, yet it remains experimentally challenging due to the simultaneous need for high temporal resolution and statistical fidelity. In this work, we demonstrate spatially resolved single-shot measurements of photon number distributions using difference-frequency generation (DFG) in a nonlinear BBO crystal. We show that our platform can discriminate between coherent and thermal photon statistics across two spatial dimensions with picosecond resolution. At the same time, we find that the retrieved distributions deviate from the ideal ones, a consequence of vacuum contamination and the multimodal response of the amplifier. To explain this, we develop a temporal mode decomposition framework that captures the essential physics of signal amplification and fluorescence, and quantitatively reproduces the experimental findings. This establishes a robust approach for measuring two-dimensional photon statistics while clarifying the fundamental factors that limit the fidelity of such measurements.

## I. INTRODUCTION

For decades, researchers have sought the ability to faithfully image spatial bi-dimensional information generated at incredibly fast timescales. Ultrafast chaotic behaviors, complex nonlinear optical phenomena, and rapid biological processes like protein rearrangement are just a few examples of events that occur on timescales ranging from picoseconds to femtoseconds. Additionally, transient events, specific to two-dimensional systems, reflect many major fundamental mechanisms occurring in biology, physics, or chemistry [1–3].

Traditionally, pump-probe methods [4–12] have been used to capture these fast dynamics through repeated measurements. However, many ultrafast phenomena are either non-repetitive or challenging to reproduce [13–18]. Significant efforts have been made to create platforms capable of performing ultrafast single-shot optical measurements [19–21], intending to be free of the paradigmatic repeatable measurement schemes. The main objective of these platforms is the possibility of measuring bi-dimensional photonic fields, at an ultrafast scale (picosecond to femtosecond information with high repetition rate) without alteration of the information of interest. Well-established strategies consist of using cutting-edge pulsed laser technologies to probe a system with strong optical pulses to have enough signal-to-noise ratio in a single-shot experiment and then retrieving the information of interest with various elaborated methods [22–27].

An important strategy to realize ultrafast sampling and perform single-shot measurements is based on the use of highly nonlinear optical processes such as Sum-Frequency Generation (SFG) or Difference-Frequency Generation (DFG), driven by a strong ultrashort pulse

in bulk nonlinear materials. Both processes convert the information contained in a 2D optical field to another wavelength and sample it on an ultrafast timescale determined by the pulse duration. The key distinction between the two methods lies in their respective properties: DFG can amplify the signal but introduces a fluorescence background, while SFG is background-free without vacuum contamination when reaching 100% conversion efficiency [28]. As a result, SFG is typically preferred for relatively strong input signals, where the signal significantly exceeds the dark count and readout noise of conventional detectors. In contrast, DFG is better suited for weak input signals due to its gain capabilities. While DFG-based bi-dimensional imaging has been demonstrated [29, 30], no focus has been set on the influence of the nonlinear effects and vacuum contamination on the spatial statistical distributions of photon counts in single-shot images. A system capable of extracting this information with high fidelity would enable the evaluation of bi-dimensional intensity correlation functions, providing additional information on ultrafast 2D phenomena.

In this study, we explore these questions by using parametric down-conversion to perform time sampling within a picosecond window of a two-dimensional image. Through a quantitative analysis of the system response, we demonstrate that the concept of temporal modes, originating from quantum optics and describing the response of nonlinear pulsed systems, serves as an effective framework to understand the system’s behavior. We further analyze the behavior of amplified vacuum fluctuations (i.e., fluorescence) and examine how known input 2D signals with distinct photon statistical distributions are altered by the DFG process. Finally, we provide detailed experimental analysis, theory and numerical investigation to explain quantitatively why the measurements deviate from the statistics of the input field due to fundamental vacuum contamination as well as due to the multimodal temporal response of the system. We obtain

<footnote>* clement.hainaut@univ-lille.fr</footnote>
predictions that match the experimental data remark-
ably. This work provides a comprehensive analysis of
the 2D photon distribution imaging framework using the
DFG scheme.

The paper is organized as follows. In the first part, we present two-dimensional single-shot measurements of photon statistics using coherent and thermal input states, showing that the system is capable of faithfully imaging complex spatial patterns while preserving qualitative differences between the two types of photon states. We then highlight the quantitative deviations of the measured statistics from the ideal photon distributions and identify their origin in vacuum contamination and the multi-modal response of the nonlinear amplifier. In the second part, we introduce a theoretical and numerical framework based on temporal mode decomposition, which provides a quantitative understanding of the system's response. Starting from a simplified plane-wave configuration, we validate the multimodal description by comparing its predictions with experiments performed with vacuum, coherent, and thermal input states. Finally, we extend this approach to the two-dimensional configuration, demonstrating that the model accurately accounts for the measured distributions, which establishes a comprehensive framework for DFG-based ultrafast-photon statistics imaging.

## II. 2D PHOTON STATISTICS
MEASUREMENTS

The system under study consists of a bulk material
with second-order nonlinearity $\chi^{(2)}$, namely, a Beta Bar-
ium Borate (BBO) crystal of length $l_c = 2$ mm, interact-
ing with two optical fields. A strong picosecond pulsed
field (≈ 0.1 mJ per pulse), referred to as the pump, with
a wavelength $\lambda_p = 400$ nm ($\omega_p = 2\pi \times 749.2$ THz) and a
continuous-wave signal at $\lambda_s = 840.1$ nm ($\omega_s = 2\pi \times 356.7$
THz) are spatially overlapped within the crystal. By
appropriately selecting the phase-matching angle of the
crystal, part of the pump field at 400 nm converts and
amplifies the signal field at $\lambda_s = 840.1$ nm and an idler
field centered at $\lambda_i = 763.3$ nm ($\omega_i = 2\pi \times 392.5$ THz) is
generated via Difference-Frequency Generation (DFG),
see level scheme in Fig. 1. The idler field is an amplified
replica of the signal field. The crystal and strong pump
beam act as a phase-insensitive amplifier, meaning that
the gain is independent of the phase of the input signal
field. The ultrafast nature of this scheme stems from the
ability to utilize a picosecond pulsed pump to achieve the
conversion, allowing for the sampling of picosecond-scale
temporal information from the signal field and transfer-
ring it into the idler field through stimulated parametric
amplification. However, in addition to the signal and
idler fields, the pump spontaneously amplifies vacuum
fluctuations at the wavelength of both signal and idler.
As we will see, these spontaneous fluorescence photons
will have an important effect in the measurement of pho-

<figure><img src="image_2.png" /><figcaption>FIG. 1. Experimental Setup. A 2 ps pump laser is sent to a BBO crystal where it is combined simultaneously with a coherent laser source (Ti:Sapph with typical power of 3-4 mW) and a thermal source (SLED with typical power of 40 µW after spectral filtering). Down-converted photons are filtered and recorded using a CMOS camera. The repetition rate of the camera, 1KHz, matches that of the pulsed laser, allowing for single shot image recording.</figcaption></figure>

ton distributions of the photon field.

A. Experimental setup and imaging

To demonstrate the capabilities of our set-up to mea-
sure photon statistics in space, we consider the case de-
picted in Fig. 1 where the input state sent into the system
is constituted of two input fields with different statistical
properties and spatial profiles. The first is a coherent
state originated from a Ti:Sapph laser, which has been
spatially tailored transversely using a Spatial Light Mod-
ulator (SLM) and a Gerchberg-Saxton-based Mixed Re-
gion Amplitude Freedom (GS-MRAF) algorithm to cre-
ate a spatial image in the form of the letter ”A” (See
appendix for more details about the SLM). The second
input field is generated from a thermal source (Exalos,
EXS210153-01) with a broad 20 nm spectrum, which we
spectrally filter with a spectrometer to retain a spectrum
of 0.2 nm width centered at 840 nm. This bandwidth
ensures a coherence time of 2 ps, suitable for observing
the thermal source with picosecond resolution. After-
wards, we employ a beam-shaping method to refocus the
thermal light into a small Gaussian beam with a waist
of 15 microns. Both input fields are combined with the
pulsed pump laser via a 10/90 beam splitter and im-
aged onto the BBO crystal. After DFG amplification
(type 1 (o) + (o) → (e) with a crystal angle of 29.14)
at the BBO crystal, a band-pass (FF01-747/33-25) and
<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) Measured average number of photons in a two-dimensional image using 10000 single-shot. (b) Variance of the photon count over the image.</figcaption></figure>

two high-pass and low-pass interference Semrock filters (TLP01-790.25x36 and TSP01-790-25x36) are added to select only the optical fields centered around 763.3 nm with roughly 1 nm bandwidth. The light emitted from the BBO crystal at these frequencies is then imaged on a Hamamatsu ORCA-Fusion CMOS camera. In this configuration, we collect 10000 single-shot 2D images at 1 kHz, the repetition frequency of the pulsed laser.

Figure 2(a) and (b) display, respectively, the average photon number and the variance on each pixel over the collection of 10000 images. A selection of the individual single shot images is shown in the appendix for completeness. Firstly, we observe that the output spatial photon distribution clearly exhibits an A-shaped pattern, which corresponds to the input image generated by the SLM using the coherent continuous-wave laser source. This result is achieved because the phase-matching tolerance was carefully chosen, ensuring that the transverse k components forming the image are all amplified with the same magnification. This measurement shows that the setup is capable of faithfully transforming complex images to another wavelength. Secondly, we observe a Gaussian shape spot at the center of the A-shaped distribution. This spot originates from the thermal light source and displays an average photon number similar to that of the coherent-source A shape. The variance at the small spot inside the A is much larger than that at the A shape itself. This important distinction indicates that the statistical distribution of photons forming the A shape is very different from that of the central small spot.

## B. 2D photon distribution analysis

We perform a quantitative analysis of the image by directly evaluating the probability of photon count number at two different positions within the image. To do so we use the 10000 single-shot images to construct two normalized histograms that define two photon distributions at two different positions marked in red and green squares in Fig. 3. The two distributions are shown in semi-logarithmic plots. The photon distribution measured at the green square location (right panel), originating from the thermal input field, shows after 30 photon counts per pixel, an exponential decay in photon number, a signature usually associated with thermal statistics. Nevertheless, if we compare with a pure thermal distribution that possesses the same average photon number as the experimental measurement (dashed green line) we see a quantitative discrepancy for photon count between 0 and 50. This indicates that the measured distribution is not purely thermal.

The left panel of Fig. 3 displays the measured photon distribution from the red square, originating from the continuous wave laser source. We can compare the distribution with that of a coherent state with average photon number corresponding to the measured one (dashed red line). We see a strong disagreement between the two, indicating that the measured photon distribution is not purely the one associated with a coherent state. Despite these disagreements, we can still conclude that this two-dimensional detection setup allows us to qualitatively distinguish between the coherent and thermal photon statistics present in the input state.

To be more quantitative, we can evaluate the intensity autocorrelation function

$$g_N^{(2)} = \frac{\langle N^2 \rangle}{\langle N \rangle^2} = 1 + \frac{\mathrm{Var}(N)}{\langle N \rangle^2}, \quad (1)$$

in the different regions of the image. In Eq. (1), N corresponds to the photon counts at a chosen location in the image and $\langle N \rangle$ its average value. However, as introduced above, the DFG set-up introduces amplified vacuum fluorescence at the wavelength of the measured idler field. This fluorescence also has an average photon number value $\langle N_{vac} \rangle$ and variance $\mathrm{Var}(N_{vac})$ that we denote as “vac” for vacuum for clarity, and it can be directly measured in an experiment with no input signals. To realize a precise computation of $g_N^{(2)}$ of the idler field, we need to remove the photon count contribution originating in the fluorescence from the measured variance and average photon count. This can be directly done because photons emitted by amplified spontaneous processes constitute independent processes compared to photons emitted via stimulated parametric processes. In the upcoming sections we will treat this in detail. Defining the corrected quantities:

$$\begin{aligned} \langle N_{\mathrm{corr}} \rangle &= \langle N \rangle - \langle N_{\mathrm{vac}} \rangle, \\ \mathrm{Var}(N_{\mathrm{corr}}) &= \mathrm{Var}(N) - \mathrm{Var}(N_{\mathrm{vac}}), \end{aligned} \quad (2)$$

we deduce a corrected autocorrelation $g_{N_{corr}}^{(2)}$:

$$g_{N_{corr}}^{(2)} = 1 + \frac{\mathrm{Var}(N_{corr})}{\langle N_{corr} \rangle^2}. \quad (3)$$

Applying these formulas to the photon counts measured at the position of the green square in Fig. 3 we obtain $g_{N_{corr}}^{(2),thermal} = 2.15$, which is close to the theoretical
<figure><img src="image_4.png" /><figcaption>FIG. 3. Center: Single shot of a two-dimensional spatial distribution and associated density distributions of the photon counts relative to two different positions in the image obtained from 10000 single-shot measurements. Dots are the experimental results, and in red and green lines are the theoretical curves calculated with $g_{2D} = 5.963 \times 10^{-13}$ s.m$^{-1}$ for both a coherent and a thermal state. The red and green dashed lines are coherent and thermal distribution, respectively, possessing the same average photon number as measured in the experiment.</figcaption></figure>

value of 2 expected for an ideal thermal source. Similarly, for photon counts measured at the position of the red square in Fig. 3 we obtain $g_{N_{corr}}^{(2),coherent} = 1.31$, which is 30% off compared to the ideal value of 1 for a perfect coherent source.

The two measurements of $g^{(2)}$ values and the full photon distributions at the green and red square locations differ both from the ideal thermal and coherent case. In the rest of this paper we will show that the difference between the measured distributions and those of the input signal are of fundamental origin. They arise from the amplification of spontaneous vacuum fluctuations and from the multimodal nature of the system “nonlinear BBO crystal + strong pump” acting as an amplifier. In the next sections, we propose to investigate quantitatively these aspects and show how we can describe comprehensively the measured photon distributions using a multimodal approach. The multimodal approach allows us to have predictions that are plotted as solid lines in Fig. 3. They capture remarkably well the measured distributions.

### III. TEMPORAL MULTIMODAL DESCRIPTION

In this section we present analytical, numerical, and experimental results that introduce and explain the temporal multimodal description of the DFG amplifier. In the following parts, for pedagogical purposes and to increase the complexity of the description progressively, we restrict the analysis to a “0D” description where we assume that all input fields are plane waves. We expect this approximation to hold as long as the transverse wavevector $k$ components of the image, which are converted through the DFG process, are all amplified with

the same gain owing to the phase-matching tolerance chosen for the BBO crystal. The experiments described in this section used a different BBO, compared to the case shown previously, with a size of $l_c = 3$ mm, optimized to work with $\lambda_s = 853$ nm and $\lambda_i = 753$ nm.

#### A. Singular Value Decomposition, Bogoliubov transformation and input-output relation

This part introduces the mathematical tools and conventions to perform the modal treatment of the DFG process and follows references [31–33]. We decompose the positive part of the electric field in terms of the Fourier components

$$ \hat{E}^{(+)}(t, z) = \sum_{n=s,i} \mathcal{E}_n \int d\Omega e^{-i(\omega_n + \Omega)t} e^{ik_z(\omega_n + \Omega)z} \hat{a}_n(\Omega, z), \quad (4) $$

where $\Omega$ are sidebands with respect to the carrier $\omega_n$ that can take the value $n = i$ for $\omega_i$ and $n = s$ for $\omega_s$ and with normalization factor

$$ \mathcal{E}_n = \sqrt{\frac{\hbar\omega_n^2}{2\epsilon_0 c^2 k_n}}, \quad (5) $$

where $\epsilon_0$ is the permittivity of vacuum, $c$ is the speed of light in vacuum, $\hbar$ is the reduced Planck constant, and $k_n$ the wave-vector of the optical field at frequency $\omega_n$. $\mathcal{E}_n$ has the meaning of a single photon field amplitude. The operators $\hat{a}_n^\dagger(\Omega, z)$ and $\hat{a}_n(\Omega, z)$ respectively create and annihilate photons at frequency $\Omega$ at position $z$ and obey the commutation relation

$$ [\hat{a}_n(\Omega, z), \hat{a}_{n'}^\dagger(\Omega', z')] = \delta_{n,n'} \delta(\Omega - \Omega') \delta(z - z'). \quad (6) $$

$\delta(\Omega - \Omega')$ and $\delta(z - z')$ are Dirac delta functions and $\delta_{n,n'}$ is the Kronecker delta. The index $n = \{s, i\}$ accounts for
the signal and idler fields collinearly propagating along the z-axis through the nonlinear medium. The interactions between these fields are mediated by a strong classical pump pulse of spectral amplitude $A_p$, and is such that we have perfect phase matching at carrier frequencies: $\omega_p = \omega_s + \omega_i$ and $k_p = k_s + k_i$. The wavevectors $k_p, k_s$, and $k_i$ are associated to the pump, signal, and idler fields, respectively. In the undepleted pump approximation (the pump amplitude does not vary over the propagation through the crystal i.e., it does not depend on $z$), the evolution of signal and idler is given by [33]

$$ \frac{\partial}{\partial z} \hat{a}_s(\Omega, z) = g \int d\Omega' K(\Omega, \Omega', z) \hat{a}_i^\dagger(\Omega', z), \quad (7) $$

$$ \frac{\partial}{\partial z} \hat{a}_i(\Omega, z) = g \int d\Omega' K(\Omega', \Omega, z) \hat{a}_s^\dagger(\Omega', z), \quad (8) $$

where $g$ is a parameter that depends on the crystal and field properties (see Appendix for more detail), $K(\Omega, \Omega', z) = A_p(\Omega+\Omega')e^{i\Delta(\Omega,\Omega')}z$ is the kernel of the two integro-differential equations, $\Delta(\Omega, \Omega') = k_{pz}(\Omega+\Omega') - k_{sz}(\Omega) - k_{iz}(\Omega')$ is the phase-mismatch angle that takes into account for the dispersion and the birefringent character of the nonlinear crystal, and $k_{nz}(\Omega) = k_z(\omega_n + \Omega)$ (for $n \in \{s,i\}$) is the $z$ component of the signal and idler wave-vectors. $A_p(\Omega+\Omega')$ is the unitless spectrum amplitude of the pump pulse (see Appendix for more details).

Since Eqs. (7) and (8) are linear in the field operators, a formal solution can be obtained in terms of a linear integral transformation (see for example [31–33]), generalizing the usual Bogoliubov transformation. Since $K$ is $z$-dependent, the exact solution can only be obtained by numerical integration. On the other hand, it has been shown [33] that the Magnus expansion truncated at first order provides an accurate analytical description up to a degree of amplification of about 12dB in fields, which corresponds to an amplification in photon number of a factor 16. In other words, this approximation is valid as long as the system gain is smaller than 16. The advantage of using truncated Magnus expansion, as we will show below, is that the solution can be cast in explicit form in terms of the normal modes of the evolution, offering a visualization and understanding of the modes in which photons are emitted. These modes are also known as “supermodes” or “temporal” modes; we will use the denomination “temporal modes” in this paper.

We start by rewriting Eqs. (7) and (8) in compact form, by defining:

$$ \hat{\xi}(\Omega, z) = \begin{pmatrix} \hat{a}_s(\Omega, z) \\ \hat{a}_i(\Omega, z) \\ \hat{a}_s^\dagger(\Omega, z) \\ \hat{a}_i^\dagger(\Omega, z) \end{pmatrix}, \quad (9) $$

so that the field evolution equation becomes

$$ \frac{\partial \hat{\xi}(\Omega, z)}{\partial z} = g \int d\Omega' \mathbb{K}(\Omega, \Omega', z) \hat{\xi}(\Omega', z), \quad (10) $$

with

$$ \mathbb{K}(\Omega, \Omega', z) = \begin{pmatrix} 0 & 0 & 0 & K(\Omega, \Omega', z) \\ 0 & 0 & K(\Omega', \Omega, z) & 0 \\ 0 & K^*(\Omega, \Omega', z) & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}, \quad (11) $$

As mentioned above, the solution of the $z$-dependent problem, Eq. (10), can be obtained in terms of the Magnus expansion truncated at the first order. This is equivalent to the infinite Dyson series when $z$-ordering is neglected. Then, associating an “out” label to the position after the crystal size at $z = l_c/2$ ($\hat{\xi}_{\text{out}}(\Omega) := \hat{\xi}(\Omega, l_c/2)$) and a “in” label for the position just before the crystal at $z = -l_c/2$ ($\hat{\xi}_{\text{in}}(\Omega) := \hat{\xi}(\Omega, -l_c/2)$), we do the Magnus expansion and integrate over the crystal size by considering the following linear and symplectic transformation:

$$ \hat{\xi}_{\text{out}}(\Omega) = \int d\Omega' \text{Exp}[gl_c M_1(\Omega, \Omega')] \hat{\xi}_{\text{in}}(\Omega'), \quad (12) $$

with

$$ M_1(\Omega, \Omega') = \begin{pmatrix} 0 & 0 & 0 & \tilde{M}_1(\Omega, \Omega') \\ 0 & 0 & \tilde{M}_1(\Omega', \Omega) & 0 \\ 0 & \tilde{M}_1^*(\Omega, \Omega') & 0 & 0 \\ \tilde{M}_1^*(\Omega', \Omega) & 0 & 0 & 0 \end{pmatrix}, \quad (13) $$

and

$$ \begin{aligned} \tilde{M}_1(\Omega, \Omega') &= \int_{-l_c/2}^{l_c/2} K^*(\Omega, \Omega', z) dz \\ &= A_p^*(\Omega + \Omega') \text{Sinc}\left(\frac{\Delta(\Omega, \Omega') l_c}{2}\right). \end{aligned} \quad (14) $$

$\tilde{M}_1(\Omega, \Omega')$ represents the generic Joint Spectral Amplitude (JSA) that describes the coupling between sidebands at frequency $\Omega$ and sidebands at $\Omega'$, and $\text{Sinc}(x) = \sin(x)/x$ is the sinc function.

From the structure of Eq. (14), it is clear that the basis of Fourier modes originally employed for describing the electric fields is not suitable for understanding the results of the dynamics of the pulsed DFG: i.e., plane-waves are not the eigenmodes of the problem. Instead, it is desirable to find the normal modes that decouple the solution of Eq. (12). This can be obtained by performing the Singular Value Decomposition (SVD) of the JSA (Eq. (14)) which leads to:

$$ \tilde{M}_1(\Omega, \Omega') = \sum_m \lambda_m \psi_m(\Omega) \phi_m^*(\Omega'), \quad (15) $$

where $\psi_m(\Omega)$ and $\phi_m(\Omega')$ are the temporal eigenmodes of the signal and the idler, respectively, and $\lambda_m \ge 0$ is the corresponding singular value. It is interesting to note that this decomposition also induces a SVD of Eq. (12) that preserves the symplectic structure of the problem, and it is known as Bloch-Messiah decomposition [31–36]. The latter can be written as

$$ \text{Exp}[gl_c M_1(\Omega, \Omega')] = \sum_m U_m(\Omega) D_m U_m^\dagger(\Omega'), \quad (16) $$
where

$$U_m(\Omega) = \begin{pmatrix} \psi_m(\Omega) & 0 & 0 & 0 \\ 0 & \phi_m^*(\Omega) & 0 & 0 \\ 0 & 0 & \psi_m^*(\Omega) & 0 \\ 0 & 0 & 0 & \phi_m(\Omega) \end{pmatrix}, \quad (17)$$

and

$$D_m = \begin{pmatrix} \cosh(g\lambda_m l_c) & 0 & 0 & \sinh(g\lambda_m l_c) \\ 0 & \cosh(g\lambda_m l_c) & \sinh(g\lambda_m l_c) & 0 \\ 0 & \sinh(g\lambda_m l_c) & \cosh(g\lambda_m l_c) & 0 \\ \sinh(g\lambda_m l_c) & 0 & 0 & \cosh(g\lambda_m l_c) \end{pmatrix}. \quad (18)$$

The two families of orthonormal functions {$\psi_m(\Omega)$}$_{m \in \mathbb{N}}$ and {$\phi_m(\Omega)$}$_{m \in \mathbb{N}}$ represent the components of the normal modes of the DFG in the Fourier basis. Thus, the normal modes in the temporal domain are given by

$$u_{s,m}(z,t) = \int d\Omega \, \psi_m(\Omega) e^{-i(\omega_s+\Omega)t} e^{ik_z(\omega_s+\Omega)z}, \quad (19)$$

$$v_{i,m}(z,t) = \int d\Omega \, \phi_m(\Omega) e^{-i(\omega_i+\Omega)t} e^{ik_z(\omega_i+\Omega)z}, \quad (20)$$

with the wavevector, $k_z$, having an explicit frequency dependence depending on the properties of the nonlinear medium. Most importantly, we can thus write explicitly the field (input/output) amplitude operators of these modes as

$$\hat{c}_s^m = \int d\Omega \, \psi_m(\Omega) \hat{a}_{s,in}(\Omega), \quad (21)$$

$$\hat{c}_i^m = \int d\Omega \, \phi_m(\Omega) \hat{a}_{i,in}(\Omega), \quad (22)$$

$$\hat{b}_s^m = \int d\Omega \, \psi_m(\Omega) \hat{a}_{s,out}(\Omega), \quad (23)$$

$$\hat{b}_i^m = \int d\Omega \, \phi_m(\Omega) \hat{a}_{i,out}(\Omega), \quad (24)$$

where we use the letter c for the input and the b for the output. In this basis, the problem is written as a set of decoupled two-mode Bogoliubov transformations

$$\begin{pmatrix} \hat{b}_s^m \\ \hat{b}_i^m \\ \hat{b}_s^{m,\dagger} \\ \hat{b}_i^{m,\dagger} \end{pmatrix} = D_m \begin{pmatrix} \hat{c}_s^m \\ \hat{c}_i^m \\ \hat{c}_s^{m\dagger} \\ \hat{c}_i^{m\dagger} \end{pmatrix}. \quad (25)$$

We can now focus on the possibilities of creating idler photons at the output of the crystal. From Eq. (25), we can identify two different independent processes. Firstly, we see that the lower left term of $D_m$ (Eq. (18)) allows for the conversion of an input signal photon in a given mode m into an output idler photon in the same mode. Secondly, the lower right term is responsible for, in a first place, stimulated emission of an idler field in the mode m at the output due to an existing idler field in the mode m at the input, and in the second place, spontaneous amplification of an idler vacuum field in a given mode m at

the output. Thus, we can summarize these different processes with the following input-output relation for each crystal eigenmode:

$$\hat{b}_i^{m\dagger} = \sinh(g\lambda_m l_c) \hat{c}_s^m + \cosh(g\lambda_m l_c) \hat{c}_i^{m\dagger}. \quad (26)$$

Defining the mode gain as

$$G_m := \cosh^2(g\lambda_m l_c), \quad (27)$$

we can rewrite Eq. (26) as

$$\hat{b}_i^{m\dagger} = \sqrt{G_m - 1} \hat{c}_s^m + \sqrt{G_m} \hat{c}_i^{m\dagger}. \quad (28)$$

For completeness, we also give the relation of the produced signal photons that can be obtained in the same manner:

$$\hat{b}_s^{m\dagger} = \sqrt{G_m - 1} \hat{c}_i^m + \sqrt{G_m} \hat{c}_s^{m\dagger}. \quad (29)$$

We have derived the fundamental equations that link the photon generated at the output of the crystal with respect to photons present at the input. As a result of the SVD, each mode m is independent and follows the same equation with a gain value $G_m$ that depends on the singular value $\lambda_m$ of the decomposition. A schematic representation of this process is given in Fig. 4. An input signal $S_{in}$ and an idler $I_{in}$ field are decomposed onto a set of orthonormal input modes {$\psi_m, \phi_m$}. Each mode pair m is independently amplified with a gain $G_m$. The total output fields ($S_{out}$ and $I_{out}$) are reconstructed from these amplified components. The dashed line in the input Idler field in Fig. 4 represents the process in which a vacuum idler field gets amplified to produce a nonzero average photon number at the output. These photons are usually referred to as fluorescence photons, which were already mentioned at the beginning of this paper.

## B. Analytical and numerical structure of the temporal modes

In the following, we aim to obtain the spectral form of the set of modes $\psi_m$ and $\phi_m$ and the singular values $\lambda_m$. Experimentally, we only measure the photons emitted at around $\lambda_i = 750$ nm, and signal photons emitted at around $\lambda_s = 840$ nm are filtered out via optical filters. We thus restrict the analysis on the emitted idler photons. Treating the emitted signal photons presents no difficulties. To make the connection between the experimental measurements and the theoretical description above explicit, we identify specific idler and signal side-bands in the Joint Spectral Amplitude (JSA) and define

$$M_1(\Omega_i, \Omega_s) := \tilde{M}_1^*(\Omega' = \Omega_i, \Omega = \Omega_s) \quad (30) \\ = A_p(\Omega_i + \Omega_s) \text{Sinc} \left( \Delta(\Omega_i, \Omega_s) \frac{l_c}{2} \right),$$

where $\Omega_s = \omega - \omega_s$ and $\Omega_i = \omega - \omega_i$. The phase-mismatch, responsible for the frequency-dependent conversion efficiency in the BBO crystal, is generically expressed as
<figure><img src="image_5.png" /><figcaption>FIG. 4. Schematic representation of the DFG process in the temporal eigenmode basis.</figcaption></figure>

$\Delta(\Omega_i, \Omega_s) = k_{pz}(\Omega_i + \Omega_s) - k_{sz}(\Omega_i) - k_{iz}(\Omega_s)$. In the plane-wave approximation, where the three fields only possess a k-component along $z$ (also called co-linear phase-matching), $\Delta(\Omega_i, \Omega_s)$ can be approximated in terms of a Taylor expansion truncated to the second order,

$$ \Delta(\Omega_i, \Omega_s) \approx (k_p - k_s - k_i) + [k'_p(\Omega_i + \Omega_s)k'_s\Omega_s - k'_i\Omega_i] + \frac{1}{2}[k''_p(\Omega_i + \Omega_s)^2 - k''_s\Omega_s^2 - k''_i\Omega_i^2]. \quad (31) $$

In this expression, $k_p$, $k_s$ and $k_i$ represent the wave-vectors at the carrier of pump, signal, and idler respectively. Similarly, $k'_p$, $k'_s$ and $k'_i$ represent the inverse of the group velocity $dk/d\omega$ at the pump, signal, and idler carrier, and $k''_p$, $k''_s$ and $k''_i$ represent the second order dispersion $d^2k/d\omega^2$ at the pump, signal, and idler carrier. We consider the experimental situation where the three waves are phase-matched ($k_p - k_s - k_i = 0$) according to the Type I configuration in which the pump at $\lambda_p = 400$ nm with polarization along the extraordinary axis produces a signal at $\lambda_s = 853$ nm and an idler at $\lambda_i = 753$ nm, both polarized along the ordinary axis. The frequency-dependent refraction index for a BBO crystal is given by the Sellmeier equations,

$$ n(\lambda) = \sqrt{A + \frac{B}{\lambda^2 - C} - D\lambda^2}. \quad (32) $$

A BBO crystal is an uniaxial nonlinear medium and, in its reference frame, the x and y axis are characterized by an ordinary refraction index $n_o(\lambda)$ with Sellmeier coefficients $A_o = 2.7359$, $B_o = 0.018782 \, \mu\text{m}^2$, $C_o = 0.01822 \, \mu\text{m}^2$ and $D_o = 0.01354 \, \mu\text{m}^{-2}$, while the z-axis is characterized by an extraordinary refraction index $n_e(\lambda)$ with Sellmeier coefficients $A_e = 2.3753$, $B_e = 0.01224 \, \mu\text{m}^2$, $C_e = 0.01667 \, \mu\text{m}^2$ and $D_e = 0.01516 \, \mu\text{m}^{-2}$. For the three selected wavelengths, the perfect phase-matching

<figure><img src="image_6.png" /><figcaption>FIG. 5. (a) Phase-matching and (b) the Joint Spectra Amplitude (JSA) profile in the $(\Omega_i, \Omega_s)$ space.</figcaption></figure>

condition ($k_p - k_s - k_i = 0$) is found when the corresponding waves propagate at an angle $\theta = 29.01$ deg with respect to the optical axis of the crystal. At this angle, the extraordinary refraction index seen by the pump is given by

$$ n_p = \left( \frac{\sin^2(\theta)}{n_e(\lambda_p)^2} + \frac{\cos^2(\theta)}{n_o(\lambda_p)^2} \right)^{-1/2} . \quad (33) $$

On the other hand, the signal and idler waves see the ordinary refraction index $n_o$ evaluated at their wavelengths. Knowing the dispersion parameters, see Table I, one can evaluate Eq. (31) and plot the phase-mismatch profile in Fig. 5(a) for a crystal length of $l_c = 3$ mm.

<table><caption>TABLE I. Dispersion parameters for critical phase-matching of BBO at room temperature.</caption><thead><tr><th></th><th>λ<sub>p</sub> = 400 nm, λ<sub>s</sub> = 853 nm, λ<sub>i</sub> = 753 nm</th></tr></thead><tbody><tr><td>k'<sub>p</sub> (×10<sup>-9</sup> s m<sup>-1</sup>)</td><td>5.81385</td></tr><tr><td>k'<sub>s</sub> (×10<sup>-9</sup> s m<sup>-1</sup>)</td><td>5.60459</td></tr><tr><td>k'<sub>i</sub> (×10<sup>-9</sup> s m<sup>-1</sup>)</td><td>5.62645</td></tr><tr><td>k''<sub>p</sub> (×10<sup>-25</sup> s<sup>2</sup> m<sup>-1</sup>)</td><td>1.96116</td></tr><tr><td>k''<sub>s</sub> (×10<sup>-26</sup> s<sup>2</sup> m<sup>-1</sup>)</td><td>6.71907</td></tr><tr><td>k''<sub>i</sub> (×10<sup>-26</sup> s<sup>2</sup> m<sup>-1</sup>)</td><td>8.20217</td></tr></tbody></table>
Multiplying the pump spectral function and the phase-matching functions gives the joint spectral amplitude $M_1(\Omega_i, \Omega_s)$ introduced earlier:

$$M_1(\Omega_i, \Omega_s) = A_p(\Omega_i + \Omega_s) \text{Sinc} \left( \Delta(\Omega_i, \Omega_s) \frac{l_c}{2} \right), \quad (34)$$

which is plotted in Fig. 5(b). It shows the strength of the frequency coupling in the nonlinear crystal, accounting for the spectral limitation of the pump. The pump's characterization is provided in the appendix. As explained in the previous part, temporal modes can generally be obtained through the numerical SVD of the JSA. However, apart from configurations where the pump and phase-matching are specifically engineered to tailor quantum states [32, 37], the JSA in most experimental situations can be well approximated by a bi-variate that we define as $M_1^{\text{ana}}(\Omega_i, \Omega_s)$:

$$M_1(\Omega_i, \Omega_s) \approx M_1^{\text{ana}}(\Omega_i, \Omega_s) \\ := e^{-(\Omega_i+\Omega_s)^2/2\Delta_1^2} e^{-(\Omega_i-\Omega_s)^2/2\Delta_2^2}, \quad (35)$$

where $\Delta_1$ accounts for the JSA bandwidth along the direction $\Omega_i - \Omega_s = 0$ and $\Delta_2$ accounts for the JSA bandwidth along the direction $\Omega_i + \Omega_s = 0$. This Gaussian approximation is valid as long as the frequency window of the emitted photons is small compared to the bending of the phase-matching due to second-order dispersion terms.

By using the approximation $\text{Sinc}(x) \approx \exp(-x^2/5)$, aiming at matching the FWHM of the Sinc function by an exponential, we identify Eq. (35) to Eq. (34) to obtain $\Delta_1^{-2} = 8 \ln(2)\Delta_\omega^{-2} + \Delta_{PM}^{-2}$ with $\Delta_{PM}^{-2} \approx (2k'_p - k'_s - k'_i)^2 l_c^2/40$. Since $\Delta_\omega \approx 1.38 \times 10^{12} \text{sec}^{-1}$ and $\Delta_{PM} \approx 5.31 \times 10^{12} \text{sec}^{-1}$, then $\Delta_1 \approx 0.584 \times 10^{12} \text{sec}^{-1}$; $\Delta_2^{-2} = (k'_i - k'_s)^2 l_c^2/40$, thus $\Delta_2 \approx 96.44 \times 10^{12} \text{sec}^{-1}$. These numbers correspond to the case of the experimental DFG with parameters shown in Table I and with the pump spectrum characterization presented in the appendix. Thus, replacing $M_1(\Omega_i, \Omega_s)$ by $M_1^{\text{ana}}(\Omega_i, \Omega_s)$ in Eq. (15) allows us to solve analytically the problem and obtain explicit solutions for the temporal modes in the form of Hermite-Gauss functions and an expression for the singular value in terms of a geometric progression. The obtained analytical results are

$$\lambda_m^{\text{ana}} = \frac{\sqrt{2\pi}\Delta_1\Delta_2}{\Delta_1 + \Delta_2} \left| \frac{\Delta_1 - \Delta_2}{\Delta_1 + \Delta_2} \right|^m, \quad (36)$$

$$\psi_m^{\text{ana}}(\Omega_i) = \frac{2^{1/4}\tau^{1/2}}{\sqrt{m!2^m\pi^{1/2}}} e^{-\tau^2\Omega_i^2} H_m(\sqrt{2}\tau\Omega_i), \quad (37)$$

$$\phi_m^{\text{ana}}(\Omega_s) = \frac{2^{1/4}\tau^{1/2}}{\sqrt{m!2^m\pi^{1/2}}} e^{-\tau^2\Omega_s^2} H_m(\sqrt{2}\tau\Omega_s), \quad (38)$$

where $\tau = 1/\sqrt{\Delta_1\Delta_2}$ and $H_n(x)$ are Hermite polynomials.

<figure><img src="image_7.png" /><figcaption>FIG. 6. (a) Signal temporal eigenmodes $\phi_m(\Omega_s)$ for mode $m=4$ and $m=1$ obtained analytically (in red and blue) using the formula of Eq. (38) and obtained by performing the numerical SVD of the JSA (black and green). (b) Analytical and numerical singular value $\lambda_m$ associated with each mode $m$.</figcaption></figure>

To illustrate the validity of this description, we compare the analytical modes calculated using the Gaussian approximation and the numerical ones obtained through numerical SVD of the JSA. Figure 6(a) shows representative temporal modes in the spectral domain. The analytical ones are Hermite polynomials and thus even or odd. The numerical ones present a slight asymmetry, as the Kernel is not symmetric due to the second-order dispersion in the phase-matching term. However, the analytical formula describes well the temporal modes and allows us to understand their oscillatory behavior originating from the Hermite polynomial functions. It is important to note that the representation in the time domain of such modes is also close to a Hermite polynomial, and thus strongly oscillating. The latter are represented on the right side of Fig. 4. Finally, Fig. 6(b) compares the numerical and the analytical singular values associated with each of the temporal modes. Again, there is an excellent agreement between the two, and the slight deviation observed is also due to the presence of second-order dispersion in the numerical calculation.

Although the analytical and numerical calculations agree well, we use the numerical results throughout this paper. The numerical approach incorporates second-order dispersion, which is present in the experiment, and is therefore used to derive key quantities such as photon statistics. The purpose of introducing the analytical modes here is to provide insights into the properties of the temporal modes and to understand their oscillatory nature in time and in the spectrum.

We have seen previously that the singular value decomposition of the joint spectral amplitude $M_1(\Omega_i, \Omega_s)$ can be written in the following form:

$$M_1(\Omega_i, \Omega_s) = \sum_m \lambda_m \psi_m(\Omega_i) \phi_m^*(\Omega_s), \quad (39)$$

where $\psi_m(\Omega_i)$ and $\phi_m(\Omega_s)$ are the respective crystal temporal eigenmodes of the idler and the signal associated with a common singular value $\lambda_m$. We recall that the
$\psi_m(\Omega_i)$ are orthogonal to each other; the same applies for $\phi_m(\Omega_s)$. These modes are of fundamental importance because they constitute the basis where two photons, simultaneously emitted in the same set of modes $\phi_m$ and $\psi_m$ through spontaneous or stimulated parametric emission, are fully correlated (entangled) and are statistically independent from photons emitted in modes with different $m$. This allows us to understand why this temporal mode basis is naturally suited to describe problems associated to photon statistics.

### C. Input-Output formalism and vacuum statistics

The evaluation of the photon number produced through the DFG process can be done considering the independent sets of eigenmodes obtained after singular value decomposition. From each set of eigenmode $m$ (for the idler and the signal) we have seen in section III A, Eqs.(28) and (29), that we can write a set of coupled equations in an input-output formalism,

$$ 
\begin{aligned}
b_i^m &= u_m c_i^m + v_m c_s^{m\dagger}, \\
b_s^m &= u_m c_s^m + v_m c_i^{m\dagger},
\end{aligned} 
\quad (40) 
$$

where $u_m = \sqrt{G_m}$ and $v_m = \sqrt{G_m - 1}$ are related to the mode gain of Eq. (27).

From Eq. (40), it is possible to compute the probability of having $n$ photons in the idler mode $m$ after the crystal (experimental situation of interest), starting from the vacuum. For a given $m$, we can call $|0\rangle_{\text{in}}$ the vacuum of the input modes (idler and signal), such that $c_s^m|0\rangle_{\text{in}} = c_i^m|0\rangle_{\text{in}} = 0$. It is convenient to introduce the state $|0\rangle_{\text{out}}$, the vacuum of the $b$ operators, $b_i^m|0\rangle_{\text{out}} = b_s^m|0\rangle_{\text{out}} = 0$. Following [38], one easily finds that

$$ |0\rangle_{\text{in}} = \frac{1}{u_m} e^{\frac{v_m}{u_m} b_i^{m\dagger} b_s^{m\dagger}} |0\rangle_{\text{out}} = \frac{1}{u_m} \sum_{n=0}^{\infty} \left( \frac{v_m}{u_m} \right)^n |n, n\rangle_{\text{out}}, \quad (41) $$

where $|n, n'\rangle_{\text{out}}$ corresponds to the state with $n$ photons of $b_i^m$ and $n'$ of $b_s^m$.

Therefore, the probability of finding exactly $n$ photons of $b_i^m$ starting from the vacuum of $c_i^m$ and $c_s^m$ is

$$ p_m(n) = \frac{1}{u_m^2} \left( \frac{v_m}{u_m} \right)^{2n} = p_m(1 - p_m)^n, \quad (42) $$

with $p_m = 1 - \left(\frac{v_m}{u_m}\right)^2 = \frac{1}{u_m^2} = 1/G_m$. From this expression, we can compute the average output photon number in mode $m$ of the idler,

$$ \langle n \rangle_m = G_m - 1, \quad (43) $$

and the corresponding variance,

$$ \langle \Delta n^2 \rangle_m = G_m(G_m - 1). \quad (44) $$

The probability distribution Eq. (42) can also be rewritten as

$$ p_m(n) = \frac{1}{1 + \langle n \rangle_m} \left( \frac{\langle n \rangle_m}{1 + \langle n \rangle_m} \right)^n, \quad (45) $$

where the only parameter is the average photon number in a mode $m$: $\langle n \rangle_m$. We see that the probability of emitting a photon into a mode $m$ knowing that we have the vacuum as an input state follows a thermal distribution with $\langle n \rangle_m = G_m - 1$ being the parameter of the law. Importantly, as we expect photons to be emitted in several modes, the resulting photon distribution is expected to be the convolution of several thermal distributions with different $\langle n \rangle_m$.

### D. Measurement of the multimodal vacuum photon distribution

To compare the previous results with experiments, we need, in a first place, to characterize the band-pass filter which is present on the experimental setup and which aims at selecting only a small frequency window around the central idler wavelength. To do so, we model the transmission function of the filter $T(\omega)$ with a high-pass and a low-pass sigmoid functions that we compare with experimental transmission measurements realized with a Ti:Sapp laser whose frequency is varied within a chosen scanning range; see Fig.7a). Using this characterization and the explicit form of the mode that we derived numerically, we can deduce the mode transmission by calculating

$$ T_m = \int_{-\infty}^{\infty} \phi_m^*(\omega) \phi_m(\omega) T(\omega) d\omega, \quad (46) $$

which is shown in Fig. 7b). From this transmission coefficient, using the calculation of the previous section, we can deduce the average photon number per mode that would go through the filter and thus will participate in the detection (hits the detector):

$$ \langle n \rangle_m^T = T_m(G_m - 1), \quad (47) $$

More importantly, we can use the theoretical calculation of the previous section, in association with the experimental characterization of the transmission, to verify if the predictions match experimental measurements of vacuum photon statistics distribution. The measurements are performed by recording 2000 single-shot measurements under pump only illumination of the BBO crystal (no input signal field). We do the statistical analysis on the pixels located at the center of the measured emission profile to remain in the zone of validity of a plane-wave, 0D approximation (see Fig. 8(a)). The histogram of the measured fluorescence after filtering with the use of the spectral filter described above is plotted in Fig. 8(b).
<figure><img src="image_8.png" /><figcaption>FIG. 7. (a) Pass band frequency filter of 13 nm centered at 749 nm. The curve is a fit of the experimental data fitted as a double sigmoid, see text. (b) Calculated transmission of the crystal eigenmode through the filters.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8. (a) Average spatial distribution of the measured fluorescence in a pump-only configuration. (b) Experimental photon distribution from pixels in the central region of the image in (a), and theoretical distribution computed with a fitted gain parameter of $g = g_0 = 8.4218 \times 10^{-13} \text{ s.m}^{-1}$ (red line).</figcaption></figure>

Using $g$ as a fitting parameter, we can calculate the joint probability of getting a photon on a detector knowing that different modes emit photons according to their specific law (see Eq. (45)). Looking at Fig. 9, we notice that nearly no photons are emitted in mode $m = 40$. In other words, modes with $m$ larger than 40 are not contributing to the photon statistics and can be disregarded. Thus, to obtain the red curve in the right panel of Fig. 8(b) we perform the numerical convolution of the first 40th thermal laws that follow Eq. (45) with $\langle n \rangle_m^T$ as a parameter. We obtain a remarkable agreement in the description of the photon statistics between the experiment and the theory, having only a unique fitting parameter $g = g_0 = 8.4218 \times 10^{-13} \text{ s.m}^{-1}$. This agreement validates the temporal multimodal description carried out in this work and shows that it provides a good description to address quantitatively the statistics of photon count emitted by amplified spontaneous emission processes.

<figure><img src="image_10.png" /><figcaption>FIG. 9. Calculated average photon number per mode emitted by the BBO crystal and passing through the filter.</figcaption></figure>

## E. Calculation of 0D Photon statistics with input coherent state

After having explained the behavior of the photon statistics emitted from the vacuum, we aim at describing a simple case of a pure input state which, specifically, in our experimental case, will be a coherent state originating from a Ti:Sapph laser. The first part of this section generalizes the calculation done in Sec. (III C) to the case of an arbitrary pure input state in the signal (and zero input photon in the idler).

As the temporal modes are independent, the calculation is performed for a single set of modes $m$, thus, the index $m$ is suppressed for clarity. We first compute the distribution of the number of output photons in the idler when starting from $n_0$ input signal photons. Similar results have been obtained in the context of Hawking radiation [39, 40]. We start by computing

$$ |0, n_0\rangle_{\text{in}} = \frac{(c_s^\dagger)^{n_0}}{\sqrt{n_0!}} |0\rangle_{\text{in}} \quad (48) $$

in terms of the out-states, that is (with $x := v/u$, $u = \sqrt{G}$ and $v = \sqrt{G-1}$)

$$ |0, n_0\rangle_{\text{in}} = \frac{(ub_s^\dagger - vb_i)^{n_0}}{u\sqrt{n_0!}} e^{xb_i^\dagger b_s^\dagger} |0\rangle_{\text{out}}. \quad (49) $$

Using

$$ \begin{aligned} (ub_s^\dagger - vb_i)e^{xb_i^\dagger b_s^\dagger}|0\rangle_{\text{out}} &= (ub_s^\dagger e^{xb_i^\dagger b_s^\dagger} - v[b_i, e^{xb_i^\dagger b_s^\dagger}])|0\rangle_{\text{out}}, \\ &= \frac{b_s^\dagger}{u}e^{xb_i^\dagger b_s^\dagger}|0\rangle_{\text{out}}, \end{aligned} \quad (50) $$

we have

$$ |0, n_0\rangle_{\text{in}} = \frac{1}{u^{n_0+1}\sqrt{n_0!}} (b_s^\dagger)^{n_0} e^{xb_i^\dagger b_s^\dagger} |0\rangle_{\text{out}}. \quad (51) $$

This implies that the probability of measuring the state $|n, n'\rangle_{\text{out}}$ is

$$ |\text{out}\langle n, n'|0, n_0\rangle_{\text{in}}|^2 = \delta_{n', n+n_0} \frac{1}{u^{2(n_0+1)}} \frac{(n+n_0)!}{n_0!n!} x^{2n}. \quad (52) $$
Considering an initial pure arbitrary state in the signal mode, i.e.

$$|\text{ini}\rangle = \sum_{n_0} c_{n_0} |0, n_0\rangle_{\text{in}}, \qquad (53)$$

where the $c_{n_0}$ are the superposition coefficients, we can use the results above, and the fact that the states with different $n'$ do not interfere with each other thanks to the Kronecker delta, to directly obtain

$$P(n) = \frac{1}{u^2 n!} x^{2n} \sum_{n_0} \frac{(n+n_0)!}{n_0! u^{2n_0}} |c_{n_0}|^2. \qquad (54)$$

If the initial state is a coherent state $|0, \alpha\rangle_{\text{in}}$, i.e. $c_{n_0} = e^{-|\alpha|^2/2} \frac{\alpha^{n_0}}{\sqrt{n_0!}}$, then the sum over $n_0$ reads

$$\sum_{n_0} \frac{(n+n_0)!}{n_0! u^{2n_0}} |c_{n_0}|^2 = e^{-|\alpha|^2} \sum_{n_0} \frac{(n+n_0)!}{n_0!^2} \left(\frac{|\alpha|^2}{u^2}\right)^{n_0}, \\ = e^{-(1-1/u^2)|\alpha|^2} n! L_n\left(-\frac{|\alpha|^2}{u^2}\right), \qquad (55)$$

where $L_n(z)$ is the n-th Laguerre polynomial, and we used

$$\sum_{n_0} \frac{(n+n_0)!}{n_0!^2} y^{n_0} = \partial_y^n \sum_{n_0} \frac{y^{n+n_0}}{n_0!}, \\ = \partial_y^n (y^n e^y), \\ = e^y e^{-y} \partial_y^n (y^n e^y), \qquad (56)$$

as well as Rodrigues' formula for the Laguerre polynomials $L_n(y) = \frac{e^y}{n!} \partial_y^n (y^n e^{-y})$.

Therefore, the number distribution for an initial coherent state reads

$$P^{(\alpha)}(n) = (1-x^2)e^{-x^2|\alpha|^2}x^{2n}L_n\left(-\frac{|\alpha|^2}{u^2}\right). \qquad (57)$$

Using the generating function of the Laguerre polynomials, $\sum_n t^n L_n(y) = \frac{e^{-\frac{ty}{1-t}}}{1-t}$, we obtain the corresponding generating function $Q^{(\alpha)}(z) = \sum_n z^n P^{(\alpha)}(n)$,

$$Q^{(\alpha)}(z) = \frac{1-x^2}{1-zx^2} \exp\left(-x^2|\alpha|^2 \frac{1-z}{1-zx^2}\right), \qquad (58)$$

from which we compute the photon number average and variance (e.g. $\langle n \rangle = \partial_z Q^{(\alpha)}|_{z=1}$)

$$\langle n \rangle = (|\alpha|^2 + 1) \frac{x^2}{1-x^2} = (|\alpha|^2 + 1)v^2, \\ \langle \Delta n^2 \rangle = (1 + |\alpha|^2(1+x^2)) \frac{x^2}{(1-x^2)^2}. \qquad (59)$$

For a given eigenmode $m$, this can be rewritten in terms of the mode gain $G_m$ and the decomposition of the initial state onto the eigenmode $m$, $\alpha_m$,

$$\langle n \rangle_m = (|\alpha_m|^2 + 1)(G_m - 1), \\ \langle \Delta n^2 \rangle_m = G_m(G_m - 1)[1 + (2 - 1/G_m)|\alpha_m|^2]. \qquad (60)$$

<figure><img src="image_11.png" /><figcaption>FIG. 10. (a) Square of the overlap integral between the input coherent state and the temporal eigenmodes (b) Normalized histogram of the output photon count measured with 2000 single-shot images. Inset: average image over 2000 single-shots of the emitted photons spatial profile. The histogram is obtained by selecting the value of the photon on the pixels located at the position of the blue disk.</figcaption></figure>

At this stage it is important to point that in an ideal single-mode and fluorescence-free response (removing the contribution when $|\alpha|^2 = 0$), the first two moments of the output photon distribution for large Gain ($G \gg 1$) are given by

$$\langle n \rangle_{\text{ideal}} = GN_0, \\ \langle \Delta n^2 \rangle_{\text{ideal}} = 2G^2 N_0, \qquad (61)$$

with $N_0 = \langle n \rangle_{\text{input}}$ the average photon number in the input coherent state. An input coherent state possesses the generic property that $\langle n \rangle_{\text{input}} / \sqrt{\langle \Delta n^2 \rangle_{\text{input}}} = \sqrt{N_0}$ that is, its signal-to-fluctuation ratio scales with the square root of the average photon number. In the ideal single-mode response discussed in Eq. (61), this ratio becomes $\langle n \rangle_{\text{ideal}} / \sqrt{\langle \Delta n^2 \rangle_{\text{ideal}}} = \sqrt{N_0}/2$. We thus see that the signal-to-fluctuation ratio scales less favorably with the average photon number. This indicates that even in an ideal single-mode scenario, there exists a fundamental limitation arising from vacuum contamination, which increases the relative fluctuations of the output photon distribution by a factor of $\sqrt{2}$. The amplifier does not preserve the coherent photon statistics through amplification processes.

In order to continue the description using a multimodal approach, it is important to find the $\alpha_m$ coefficients of Eq. (60). To do so, we need to compute the overlap integral between the initial signal state and the crystal temporal eigenmode $\phi_m$. For the sake of simplicity, we will assume that the initial state possesses the same spectral limitation as the one imposed by the pump pulse beam: in the experiment, the signal beam is a continuous wave source, and the response of the DFG processes is solely set by the pulsed pump beam. Thus, we will assume an input pulse having a spectrum

$$A_s(\Omega_s) = \exp\left(-4 \ln 2 \frac{\Omega_s^2}{\Delta_{\omega}^2}\right). \qquad (62)$$
<figure><img src="image_12.png" /><figcaption>FIG. 11. (a) Average output photon counts as a function of the average input photon number (b) Standard deviation ($\sqrt{\Delta N_{output}^2}$) of the output photon distribution as a function of the average input photon. The error on the experimental point is a systematic one, estimated to be ±10% due to optical waist and power measurements. The black line is the calculation done with $g_0 = 8.4218 \times 10^{-13} \text{ s.m}^{-1}$. The shadow zone around the theory curve shows that prediction if we allow for a variation of ±5% of the fitting parameter $g_0$.</figcaption></figure>

with $\Delta_{\omega} = 2.4$ THz as for the pump pulse. Thus, the overlap is calculated using

$$O_m = \int_{-\infty}^{\infty} A_s(\Omega_s) \phi_m(\Omega_s) d\Omega_s, \quad (63)$$

and shown in Fig. 10(a).

We see that the mode $m = 0$ contributes to nearly 30% of the decomposition. We notice that after $m = 20$ the different contributions lie in the percentage level. Finally, the $|\alpha_m|^2$ are just obtained from

$$|\alpha_m|^2 = N_{input} |O_m|^2, \quad (64)$$

where $N_{input}$ is the average number of photons in the input state.

### F. Measurements of 0D Photon statistics with input coherent state

To compare the 0D description carried-out previously with the experiments, we need to evaluate the experimental input photon number in the coherent state. We do so by evaluating the number of photons existing over an area corresponding to 1 pixel of the camera ($S = 66 \mu\text{m}^2$) over a time that include the duration of the pulse pump beam $\Delta_t = 2.96$ ps, knowing that the pixel of interest is chosen to be at the central position of the Gaussian intensity beam profile with a waist of $\omega_r = 140 \mu\text{m}$ and power $P_s$. Thus, the average input photon number at the central pixels of the measured spot is calculated as

$$N_{input} = \frac{2P_s \Delta_t S}{\hbar \omega_s \omega_r^2 \pi}. \quad (65)$$

We recorded 2000 single-shot images for each $N_{input}$ that we vary from 0 to 74. An example in the case of

$N_{input} = 74$ is given in Fig. 10(b) where the Region Of Interest (ROI) for the analysis is indicated as a blue disk in the inset. In Fig. 11 we show the measurement of the first two moments of the photon distributions realized with different power $P_s$ ranging from 0 mW to 2 mW corresponding to a variation of the average input photon number from 0 to 74. The point at $N_{input} = 0$ corresponds to the vacuum case presented in the previous section, where we extracted a fitted value of the $g$ parameter of $g_0 = 8.4218 \times 10^{-13} \text{ s.m}^{-1}$. Importantly, we keep this gain value estimated from vacuum statistics to compute the prediction with a non-zero signal, using Eq. (60) to compute the expected photon number emitted as well as the variance contribution of each mode $m$. As the modes are independent, a summation is done over the first 40 modes to get the total number of photons as well as the variance. Similarly to what has been done for the vacuum statistics, the transmission through the filter is also taken into account. Finally, we obtain the following formula for the prediction of the average and variance of the output photon number:

$$\begin{aligned} N_{output} &= \sum_m T_m (|\alpha_m|^2 + 1)(G_m - 1), \\ \Delta N_{output}^2 &= \sum_m T_m^2 G_m (G_m - 1) [1 + (2 - 1/G_m) |\alpha_m|^2]. \end{aligned} \quad (66)$$

$|\alpha_m|^2$ is given by Eq. (64), $G_m$ by Eq. (27), and $T_m$ by Eq. (47). The computed first two moments are plotted with solid lines in Fig. 11. Overall, we find a remarkable agreement between the prediction and experiment, confirming that the system in this configuration behaves effectively as a 0D amplifier whose photon distributions can be fully captured by the multimodal framework introduced above. An essential point is that the gain parameter extracted from vacuum statistics alone is sufficient to reproduce the results obtained with an input coherent state, without any additional fitting. This demonstrates both the predictive power and the internal consistency of the temporal mode approach. Moreover, the analysis highlights the crucial role of the large number of contributing modes –around forty in our case– in shaping the output photon distribution, underlying how multimodal effects intrinsically limit the preservation of input photon statistics.

### G. Calculation of 0D Photon statistics with input thermal state

For completeness, we provide here a generalized calculation that describes the case when the input state is a statistical mixture described by a thermal density matrix (keeping the input idler state as vacuum). Again, as the modes are independent, we perform the calculation on one set of modes and suppress the label $m$ for clarity.
The initial state is

$$ \rho_{\text{in}} = \sum_{n'} \frac{e^{-\beta\hbar\omega n'}}{Z} |0, n'\rangle_{\text{in}} \langle 0, n'|, \quad (67) $$

with $Z = \sum_{n'} e^{-\beta\hbar\omega n'}$, and $\beta$ the inverse temperature linked to the average photon number. The average number of photons in the input source is $\bar{n} = (e^{\beta\hbar\omega} - 1)^{-1}$, and the probability to measure $n$ photons in the input signal is the usual thermal distribution $\frac{1}{1+\bar{n}}(\frac{\bar{n}}{1+\bar{n}})^n$.

The probability of finding $n$ photons in the idler out-state is (using Eq. (52))

$$ \begin{aligned} P^{\text{th}}(n) &= \sum_{k} \text{out} \langle n, k | \rho_{\text{in}} | n, k \rangle_{\text{out}}, \\ &= \sum_{k, n'} \frac{e^{-\beta \hbar \omega n'}}{Z} |\text{out} \langle n, k | 0, n' \rangle_{\text{in}}|^2, \\ &= \sum_{k, n'} \frac{e^{-\beta \hbar \omega n'}}{Z} \delta_{k, n+n'} \frac{1}{u^{2(n'+1)}} \frac{(n+n')!}{n'!n!} x^{2n}, \\ &= \sum_{n'} \frac{e^{-\beta \hbar \omega n'}}{Z} \frac{1-x^2}{n'!} \frac{1}{u^{2n'}} \frac{(n+n')!}{n!} x^{2n}, \\ &= \frac{1-x^2}{Z} \frac{x^{2n}}{n!} \sum_{n'} \frac{(n+n')!}{n'!} \left( \frac{e^{-\beta \hbar \omega}}{u^2} \right)^{n'}. \end{aligned} \quad (68) $$

Using $\sum_{n'} \frac{(n+n')!}{n'!} y^{n'} = \frac{\partial^n}{\partial y^n} \sum_{n'} y^{n+n'} = \frac{\partial^n}{\partial y^n} \frac{y^n}{1-y} = \frac{n!}{(1-y)^{n+1}}$, we obtain

$$ \begin{aligned} P^{\text{th}}(n) &= \frac{1-x^2}{Z} \frac{x^{2n}}{(1-e^{-\beta\hbar\omega}/u^2)^{n+1}}, \\ &= \frac{1-e^{-\beta\hbar\omega}}{u^2-e^{-\beta\hbar\omega}} \left( \frac{v^2}{u^2-e^{-\beta\hbar\omega}} \right)^n, \quad (69) \\ &= p_{\text{th}}(1-p_{\text{th}})^n, \end{aligned} $$

defining $p_{\text{th}}$. We remark right away that this distribution is also thermal, with mean photon number $\langle n \rangle = \frac{1-p_{\text{th}}}{p_{\text{th}}} = \frac{v^2}{1-e^{-\beta\hbar\omega}} = (G-1)(\bar{n}+1)$ and

$$ P^{\text{th}}(n) = \frac{1}{1 + \langle n \rangle} \left( \frac{\langle n \rangle}{1 + \langle n \rangle} \right)^n. \quad (70) $$

Putting back the label $m$, we find that the output mean and variance of the idler photon number in mode $m$ reads

$$ \langle n \rangle_m = (G_m - 1)(\bar{n}_m + 1), \quad (71) $$

$$ \begin{aligned} \langle \Delta n^2 \rangle_m &= (G_m - 1)(\bar{n}_m + 1)[(G_m - 1)(\bar{n}_m + 1) + 1] \\ &= \langle n \rangle_m (\langle n \rangle_m + 1). \end{aligned} \quad (72) $$

The latter result is that of a thermal distribution. Importantly, we see that if the system is set in single mode configuration, a thermal input state will stay thermal (with a higher temperature) after the amplification. This is in stark contrast to the coherent case when, in a single-mode restriction, we have seen that the coherent photon statistic is not preserved after amplification.

<figure><img src="image_13.png" /><figcaption>FIG. 12. (a) Filter transmission in the configuration optimized for $\lambda_s = 840$nm. (b) Computed mode transmission $T_m$.</figcaption></figure>

## IV. EXPERIMENTAL SETTINGS FOR 2D MEASUREMENTS

At the beginning of the manuscript, we have presented two-dimensional measurements allowing us to distinguish qualitatively between two different photon number statistics. We have seen that the measurements do not show quantitatively the ideal behavior associated with thermal and coherent photon statistics. In this section, we will use the theoretical results displayed in Sec. III to derive predictions that capture the experimental behavior and, thus, explain the observed alteration of the photon statistics through the nonlinear amplifier.

We start by explicitly providing all the information about the parameters used to perform the measurements on the photon distributions of a coherent and thermal source, which were different compared to the 0D characterization developed in Sec. III. Indeed, to improve the imaging capability of the experimental setup, we changed the BBO crystal with a $l_c = 2$ mm length, allowing for a larger phase-matching bandwidth for the transverse $k$ component of the image of interest, and we doubled the spatial size of the pulse pump beam to increase the constant gain region in the transverse plane. In addition, as we want to obtain statistical information on the photons, we narrow down the frequency filter window to maximize the contribution originating mainly from one mode. The filter is configured as presented in Fig. 12(a), which leads to the mode transmission presented in Fig. 12(b). The phase-matching is optimized for a signal of 840 nm and an idler of 763.6 nm (this corresponds to an angle of 29.14 degrees).

We do not present here the full details of the calculation of the singular values and the temporal modes, as they are conceptually the same as presented in Sec III, with a crystal length of $l_c = 2$ mm and dispersion parameters presented in table II. We keep the same FWHM for the pump spectrum as in Sec III: $\Delta\omega = 2.4$ THz. From the new set of eigenmodes, we calculate the mode overlap with the input state accounting for the pump pulse duration, see Fig. 13(a).

To analyze the statistical distribution of the output
<table><caption>TABLE II. Dispersion parameters for critical phase-matching of BBO at room temperature optimized with an angle of 29.14 degree.</caption><thead><tr><th></th><th>λ<sub>p</sub> = 400 nm, λ<sub>s</sub> = 840 nm, λ<sub>i</sub> = 763.6 nm</th></tr></thead><tbody><tr><td>k'<sub>p</sub> (×10<sup>-9</sup> s m<sup>-1</sup>)</td><td>5.81385</td></tr><tr><td>k'<sub>s</sub> (×10<sup>-9</sup> s m<sup>-1</sup>)</td><td>5.61079</td></tr><tr><td>k''<sub>i</sub> (×10<sup>-9</sup> s m<sup>-1</sup>)</td><td>5.62757</td></tr><tr><td>k''<sub>p</sub> (×10<sup>-25</sup> s<sup>2</sup> m<sup>-1</sup>)</td><td>1.96116</td></tr><tr><td>k''<sub>s</sub> (×10<sup>-26</sup> s<sup>2</sup> m<sup>-1</sup>)</td><td>6.8984</td></tr><tr><td>k''<sub>i</sub> (×10<sup>-26</sup> s<sup>2</sup> m<sup>-1</sup>)</td><td>8.0428</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>FIG. 13. (a) Calculated mode overlap between the initial state and the temporal mode of the system when optimized for a signal wavelength of 840 nm and a filter in the configuration of Fig. 12. (b) Computed relative contribution of each mode to the statistical information at the output, Eq. 74.</figcaption></figure>

photons, it is useful to define a quantity that isolates the contribution of each temporal mode to the amplified signal of interest. We therefore introduce the average photon number per mode, corrected for fluorescence background, which reads

$$N_{\text{int}}^{m} = T_{m} N_{\text{input}} |O_{m}|^{2} (G_{m} - 1). \quad (73)$$

From this definition, one can directly evaluate the relative weight of each mode in carrying the statistical information of the input signal. This is quantified by

$$C_m = \frac{N_{\text{int}}^m}{\sum_m N_{\text{int}}^m}, \quad (74)$$

which informs us how much each temporal mode contributes to the transmission of the input photon statistics through the amplifier. A large value of $C_m$ indicates that most of the statistical information is carried by a single dominant mode, whereas smaller values distributed over many modes indicate that the information is spread across a multimodal response. Figure 13(b) displays the computed values of $C_m$. We see that the first mode contributes to 80% and the third one to 10%. The remaining contribution is split into the other modes. We designed the experiment to produce this response, highly peaked on the first mode, by selecting a small filtering window. The reason is that if we aim at preserving in the output field the statistical properties of the input field as much as possible, we would need to work only with 1 temporal

mode. Otherwise, if we have a large collection of modes that contribute significantly to the output, the measured statistics will be a collection of the emitted photons in all the modes, which will be described by the joint probability of all the independent mode distributions. The latter, in the limit of a large mode number, will finally result in a Gaussian distribution due to the central limit theorem and, thus, information on the input field statistics will be washed out in the amplification process.

We can now compute the expected photon distributions at the output of the BBO crystal using as input states the coherent or thermal sources used in the experiments and displayed in Figs. 2 and 3. It will allow us to compare the multimodal approach used in this article with the experimental observations. To do so, we compute the first 10 modal contributions to the idler output states using the parameters presented in this section and let the value of $g$ as a free parameter. We use the explicit formula for the photon distribution given by Eqs. (70) and (57) for the thermal and coherent statistics, respectively. In the case of the thermal distribution, the average number of photons per mode that serve as the parameter of the law is given by $\langle n \rangle_m = (G_m - 1)N_{\text{input}}^{\text{th}}|O_m|^2$ with $N_{\text{input}}^{\text{th}} = 35$. In the case of the coherent distribution, each mode has $|\alpha_m|^2 = N_{\text{input}}^{\text{ch}}|O_m|^2$ and $x_m^2 = (1 - 1/G_m)$ with $N_{\text{input}}^{\text{ch}} = 51.5$. The gains, $G_m$, calculated from Eq. (27) are given by the singular values of the SVD, the crystal length $l_c$ of 2mm, and a fitted value of $g = g_{2D} = 5.963 \times 10^{-13} \text{ s.m}^{-1}$. The smaller fitted value of $g$ compared to the experiments discussed in Sec. III is consistent with the fact of having an experimentally larger pump beam size in the case of the measurements of Figs. 2 and 3, while keeping the power constant. For each mode $m$, knowing $T_m$, we calculate the expected photon distribution after the filter. To be complete, we also consider the false count distribution coming from our detector. In our case, the ORCA fusion camera has readout noise distributed in a Gaussian manner with a standard deviation of $\sigma = 3.05$ photon count.

The calculated probability distributions for thermal and coherent input signals are plotted in Fig. 3. The predictions obtained from the multimodal framework reproduce the experimental photon distributions with excellent accuracy. This agreement confirms that the deviations observed from ideal coherent and thermal statistics originate naturally from the combined effects of vacuum contamination (fluorescence) and the multimode response of the amplifier, which offers a clear explanation of the experimental system's response.

## V. CONCLUSION AND OUTLOOK

In this work, we have demonstrated a method to measure spatially resolved photon statistics using difference-frequency generation (DFG) in a BBO crystal. By introducing a temporal multimode decomposition of the sys-
tem response, we have shown that the amplification pro-
cess can be described as a set of independent parametric
channels, each corresponding to a temporal eigenmode.
This framework provides an accurate and physically intu-
itive modeling of both fluorescence generation and signal
amplification, and explains how multimode effects shape
the measured photon statistics.

The central result of this study is twofold. First, we es-
tablish the capability of resolving photon statistics across
two spatial dimensions and demonstrate that our setup
can qualitatively discriminate between coherent and ther-
mal input states. Second, we show that the observed de-
viations from ideal distributions can be fully explained
by vacuum contamination and the multimodal response
of the amplifier. This dual outcome both validates the
experimental platform and highlights the importance of
temporal modes in determining the fidelity of photon
statistics measurements.

Beyond this proof of principle, our approach estab-
lishes a robust optical platform for bidimensional pho-
ton statistics measurements. Such a capability is par-
ticularly relevant for hybrid light-matter systems, where
photon correlations encode key signatures of many-body
physics. In exciton-polariton microcavities, for instance,
single-shot access to photon statistics in two dimensions
could provide crucial insights into the onset of superfluid-
ity [41–47], the role of fluctuations, and the emergence of
turbulent behavior [48, 49] driven by vortex interactions.

Looking ahead, several promising directions emerge
from this work. First, the integration of engineered
nonlinear materials, such as periodically poled crystals,
could offer increased flexibility in phase-matching and al-
low tailoring of the temporal mode structure. Consider-
ing Type-II phase-matching could offer twofold benefits:
the possibility of separating signal and idler based on
polarization, and the opportunity to engineer the pro-
cess to operate in a temporally single-mode regime. Fi-
nally, other nonlinear processes could also be considered
for one-shot two-dimensional photon statistics measure-
ments, such as sum-frequency generation or four-wave
mixing. The use of SFG could allow direct access to
the input statistics, without the need to infer them as
we have done here, provided that the 100% efficiency
limit is reached. While in single-pass configurations,
such as those in [28, 33, 50], the efficiency cannot ex-
ceed 80%, temporal-mode interferometry schemes, such
as those proposed by [51], could allow approaching near-
unity efficiency.

ACKNOWLEDGMENTS

This work was supported by the European Research
Council grant EmergenTopo (865151), the QuantERA
project MOLAR funded by the Agence National de
la Recherche (ANR-24-QUA2-0006), the projects ANR-
24-CE30-6695 FUSIoN, ANR-24-CE47-4949 UniQ-RingS
and ANR-23-PETQ-0001 Dyn1D at the title of France

2030, and by the European Union, the French govern-
ment through the Programme Investissement d'Avenir
(I-SITE ULNE /ANR16-IDEX-0004 ULNE) managed by
the Agence Nationale de la Recherche, the Labex CEMPI
(ANR-11LABX-0007) and the region Hauts-de-France
(CPER WAVETECH). This project has received fund-
ing from the European Union's Horizon 2020 research
and innovation programme under the Marie Sklodowska-
Curie grant agreement No 101108433.

DATA AVAILABILITY

The data that support the findings of this article are
openly available [52].
# APPENDIX: SINGLE-SHOT MEASUREMENTS

For completeness, in Fig.14) we present a few single shots from the pool of 10000 recorded images.

<figure><img src="image_15.png" /><figcaption>FIG. 14. Single-shots of a 2D spatial distributions. The A shape is made from a SLM and a Ti:Sapph laser, while the central spot originated from a filtered Exalos light source. The axis are in microns.</figcaption></figure>

# APPENDIX: SPECTRAL AND TEMPORAL PROPERTIES OF THE PUMP FIELD

The intense 400 nm pulse used for DFG is generated by frequency doubling an 800 nm pulse emitted from a Legend Elite Coherent laser system. The temporal width (FWHM) of the 800 nm pulse was measured with an autocorrelator (APE PulseCheck) and found to be $\Delta_t^{autoco} = 4.2$ ps. Assuming a Gaussian intensity profile, we apply a deconvolution factor of 0.707, yielding an intensity FWHM of $\Delta_t^{I,800nm} = 0.707 \times \Delta_t^{autoco} = 2.96$ ps. In the experiment, this pulse is frequency-doubled to produce the 400 nm beam. Ideally, the intensity FWHM of the resulting 400 nm pulse is expected to be reduced by a factor $\sqrt{2}$, giving $\Delta_t^{I,400nm} = \Delta_t^{I,800nm}/\sqrt{2} = 2.02$ ps. Correspondingly, the FWHM of the 400 nm field amplitude envelope is expected to be $\Delta_t = \sqrt{2} \times \Delta_t^{I,400nm} = 2.96$ ps.

Nevertheless, we measured the spectrum of the 400 nm pulse with an Optical Spectrum Analyzer (OSA) and found a FWHM of $\Delta_\omega = 2.4$ THz. For information, the Fourier limit due to the pulse duration for a Gaussian pulse is $\Delta_\omega^F = 0.44 \times 2\pi/\Delta_t = 0.93$ THz. The measured spectral width is thus approximately three times larger than this limit. This discrepancy likely arises from the frequency-doubling process used to generate the 400 nm pulse, which can introduce spectral phase distortions. As we do not have precise information about the resulting spectral phase, and since its impact is limited to minor modifications of the temporal mode shape and fitted gain values, we choose to neglect it. For simplicity, we will

therefore assume throughout the paper that the pump pulse spectrum is given by:

$$A_p^u(\Omega_p = \omega - \omega_p) = A_p^0 A_p(\Omega_p), \quad (75)$$

with $\Omega_p$ the sideband of the pump. By energy conservation ($\Omega_p = \Omega_s + \Omega_i$), we can rewrite the pump spectrum like: $A_p(\Omega_i + \Omega_s) = \exp[-4 \ln 2 \frac{(\Omega_i+\Omega_s)^2}{\Delta_\omega^2}]$ with $\Omega_i = \omega_1 - \omega_i$, $\Omega_s = \omega_2 - \omega_s$, $A_p^0$ the unitless peak spectral amplitude. $\omega_1$ and $\omega_2$ are two independent frequencies necessary to understand the nonlinear processes

# APPENDIX: SLM SETUP

This part provides a details about the creation of a complex input image using an SLM. The A input single image used in the second part of the manuscript to demonstrate faithful imaging was created by tailoring the coherent input beam from the Ti:Sapph laser with a phase-only Spatial Light Modulator (HOLOEYE Pluto 2.1 LCOS). Doing so we created a complex shape possessing a rich transverse $k_x$ and $k_y$ distribution. To realize this input state, we developed a homemade MRAF algorithm[53] allowing us to obtain a desired image with 20% overall error. The explicit Python code developed for creating the phase mask can be found in the corresponding GitHub repository [54]. The target intensity, the phase mask applied to the SLM, as well as the recorded image measured with a regular camera at the Fourier plane of the SLM, are represented in Fig. 15.

<figure><img src="image_16.png" /><figcaption>FIG. 15. (a) Target intensity pattern generated by the MRAF algorithm with mixing parameter m=0.4 after n= 100 iterations. (b) Corresponding phase mask generated and addressed on the as shown in Fig. (1). (c) Experimentally measured intensity distribution at the Fourier plane of an f = 150 mm lens.</figcaption></figure>

# APPENDIX: CRYSTAL GAIN AND RELATION WITH EXPERIMENTAL PARAMETERS

We have seen above that the mode gain is defined by

$$G_m = \cosh^2(g\lambda_m l_c). \quad (76)$$
with $l_c$ the crystal length over the propagation axis and $g$ is linked to the medium and fields through [55]

$$g = |A_p^0| \frac{\chi^{(2)}}{2} \sqrt{\frac{\omega_s \omega_i}{c^2 n(\omega_s) n(\omega_i)}}, \quad (77)$$

where $n(\omega_s)$ and $n(\omega_i)$ are the crystal refractive index evaluated at signal and idler frequencies respectively, and $\chi^{(2)} = 3.91 \times 10^{-12} \text{ m.V}^{-1}$ the nonlinear coefficient associated to the BBO crystal.

$A_p^0$ is the spectral amplitude of the pump. In the experiment, $A_p^0$ depends on the waist and the optical power of the pump beam, which can vary shot-to-shot. That is why in the main plot of the paper, we allowed for a variation of ±5% of $g$ to obtain the shadow area that illustrates the sensitivity to typical experimental variations of $g$.

[1] A. J. Bard and L. R. Faulkner, *Electrochemical Methods: Fundamentals and Applications*, 2nd ed. (Wiley, 2001).

[2] T. Gorkhover, S. Schorb, R. Coffee, M. Adolph, L. Foucar, D. Rupp, A. Aquila, J. D. Bozek, S. W. Epp, B. Erk, L. Gumprecht, L. Holmegaard, A. Hartmann, R. Hartmann, G. Hauser, P. Holl, A. Hömke, P. Johnson, N. Kimmel, K.-U. Kühnel, M. Messerschmidt, C. Reich, A. Rouzée, B. Rudek, C. Schmidt, J. Schulz, H. Soltau, S. Stern, G. Weidenspointner, B. White, J. Küpper, L. Strüder, I. Schlichting, J. Ullrich, D. Rolles, A. Rudenko, T. Möller, and C. Bostedt, Femtosecond and nanometre visualization of structural dynamics in superheated nanoparticles, Nature Photonics 10, 93–97 (2016).

[3] M. Imada, A. Fujimori, and Y. Tokura, Metal-insulator transitions, Reviews of Modern Physics 70, 1039–1263 (1998).

[4] P. Hockett, C. Z. Bisgaard, O. J. Clarkin, and A. Stolow, Time-resolved imaging of purely valence-electron dynamics during a chemical reaction, Nature Physics 7, 612–615 (2011).

[5] C. Y. Wong, R. M. Alvey, D. B. Turner, K. E. Wilk, D. A. Bryant, P. M. G. Curmi, R. J. Silbey, and G. D. Scholes, Electronic coherence lineshapes reveal hidden excitonic correlations in photosynthetic light harvesting, Nature Chemistry 4, 396–404 (2012).

[6] Y. Acremann, C. H. Back, M. Buess, O. Portmann, A. Vaterlaus, D. Pescia, and H. Melchior, Imaging precessional motion of the magnetization vector, Science 290, 492–495 (2000).

[7] T. Feurer, J. C. Vaughan, and K. A. Nelson, Spatiotemporal coherent control of lattice vibrational waves, Science 299, 374–377 (2003).

[8] R. Kodama, P. A. Norreys, K. Mima, A. E. Dangor, R. G. Evans, H. Fujita, Y. Kitagawa, K. Krushelnick, T. Miyakoshi, N. Miyanaga, T. Norimatsu, S. J. Rose, T. Shozaki, K. Shigemori, A. Sunahara, M. Tampo, K. A. Tanaka, Y. Toyama, T. Yamanaka, and M. Zepf, Fast heating of ultrahigh-density plasma as a step towards laser fusion ignition, Nature 412, 798–802 (2001).

[9] A. Velten, T. Willwacher, O. Gupta, A. Veeraraghavan, M. G. Bawendi, and R. Raskar, Recovering three-dimensional shape around a corner using ultrafast time-of-flight imaging, Nature Communications 3, 10.1038/ncomms1747 (2012).

[10] A. Barty, S. Boutet, M. J. Bogan, S. Hau-Riege, S. Marchesini, K. Sokolowski-Tinten, N. Stojanovic, R. Tobey, H. Ehrke, A. Cavalleri, S. Düsterer, M. Frank, S. Bajt, B. W. Woods, M. M. Seibert, J. Hajdu, R. Treusch, and H. N. Chapman, Ultrafast single-shot diffraction imaging of nanoscale dynamics, Nature Photonics 2, 415–419 (2008).

[11] J. Hajdu, C. M. Wilmot, R. Neutze, T. Sjögren, K. Edman, A. Szöke, and R. C. Wilmouth, Nature Structural Biology 7, 1006–1012 (2000).

[12] A. H. Zewail, Laser femtochemistry, Science 242, 1645–1653 (1988).

[13] B. Jalali, D. Solli, K. Goda, K. Tsia, and C. Ropers, Real-time measurements, rare events and photon economics, The European Physical Journal Special Topics 185, 145–157 (2010).

[14] D. R. Solli, C. Ropers, P. Koonath, and B. Jalali, Optical rogue waves, Nature 450, 1054–1057 (2007).

[15] P. R. Poulin and K. A. Nelson, Irreversible organic crystalline chemistry monitored in real time, Science 313, 1756–1760 (2006).

[16] V. Tuchin, *Tissue Optics*, (SPIE, 2007) pp. 303–304.

[17] N. Šiaulys, L. Gallais, and A. Melninkaitis, Direct holographic imaging of ultrafast laser damage process in thin films, Optics Letters 39, 2164 (2014).

[18] M. Sciamanna and K. A. Shore, Physics and applications of laser diode chaos, Nature Photonics 9, 151–162 (2015).

[19] J. Liang and L. V. Wang, Single-shot ultrafast optical imaging, Optica 5, 1113 (2018).

[20] L. Gao, J. Liang, and C. e. a. Li, Single-shot compressed ultrafast photography at one hundred billion frames per second (2014).

[21] C. Jin, Y. Xu, D. Qi, Y. Yao, Y. Shen, L. Deng, R. Han, Z. Pan, J. Yao, Y. He, Z. Huang, X. Pan, H. Tao, M. Sun, C. Liu, J. Shi, J. Liang, Z. Wang, J. Zhu, Z. Sun, and S. Zhang, Single-shot intensity- and phase-sensitive compressive sensing-based coherent modulation ultrafast imaging, Physical Review Letters 132, 10.1103/physrevlett.132.173801 (2024).

[22] X. Wang, L. Yan, J. Si, S. Matsuo, H. Xu, and X. Hou, High-frame-rate observation of single femtosecond laser pulse propagation in fused silica using an echelon and optical polarography technique, Appl. Opt. 53, 8395 (2014).

[23] N. Abramson, Light-in-flight recording by holography, Opt. Lett. 3, 121 (1978).

[24] Z. Li, R. Zgadzaj, X. Wang, Y.-Y. Chang, and M. C. Downer, Single-shot tomographic movies of evolving light-velocity objects, Nature Communications 5, 10.1038/ncomms4085 (2014).

[25] K. Nakagawa, A. Iwasaki, Y. Oishi, R. Horisaki, A. Tsukamoto, A. Nakamura, K. Hirosawa, H. Liao, T. Ushida, K. Goda, F. Kannari, and I. Sakuma, Sequentially timed all-optical mapping photography (stamp), Nature Photonics 8, 695–700 (2014).

[26] Q.-Y. Yue, Z.-J. Cheng, L. Han, Y. Yang, and C.-S. Guo, One-shot time-resolved holographic polarization
microscopy for imaging laser-induced ultrafast phenomena, Optics Express **25**, 14182 (2017).

[27] A. Ehn, J. Bood, Z. Li, E. Berrocal, M. Aldén, and E. Kristensson, Frame: femtosecond videography for atomic and molecular dynamics, Light Sci. Appl. **6**, e17045–e17045 (2017).

[28] A. Eckstein, B. Brecht, and C. Silberhorn, A quantum pulse gate based on spectrally engineered sum frequency generation, Opt. Express **19**, 13770 (2011).

[29] F. Devaux and E. Lantz, Ultrahigh-speed imaging by parametric image amplification, Optics Communications **118**, 25–27 (1995).

[30] P. M. Vaughan and R. Trebino, Optical-parametric-amplification imaging of complex objects, Opt. Express **19**, 8920 (2011).

[31] W. Wasilewski, A. I. Lvovsky, K. Banaszek, and C. Radzewicz, Pulsed squeezed light: Simultaneous squeezing of multiple modes, Physical Review A **73**, 10.1103/physreva.73.063819 (2006).

[32] G. Patera, C. Navarrete-Benloch, G. J. de Valcárcel, and C. Fabre, Quantum coherent control of highly multipartite continuous-variable entangled states by tailoring parametric interactions, The European Physical Journal D **66**, 241 (2012).

[33] A. Christ, B. Brecht, W. Mauerer, and C. Silberhorn, Theory of quantum frequency conversion and type-ii parametric down-conversion in the high-gain regime, New Journal of Physics **15**, 053038 (2013).

[34] R. S. Bennink and R. W. Boyd, Improved measurement of multimode squeezed light via an eigen-mode approach, Physical Review A **66**, 10.1103/physreva.66.053815 (2002).

[35] C. K. Law, I. A. Walmsley, and J. H. Eberly, Continuous frequency entanglement: Effective finite hilbert space and entropy control, Physical Review Letters **84**, 5304–5307 (2000).

[36] G. J. de Valcárcel, G. Patera, N. Treps, and C. Fabre, Multimode squeezing of frequency combs, Physical Review A **74**, 10.1103/physreva.74.061801 (2006).

[37] L. Serino, W. Ridder, A. Bhattacharjee, J. Gil-Lopez, B. Brecht, and C. Silberhorn, Orchestrating time and color: a programmable source of high-dimensional entanglement, Optica Quantum **2**, 339 (2024).

[38] S. Takagi, Vacuum Noise and Stress Induced by Uniform Acceleration: Hawking-Unruh Effect in Rindler Manifold of Arbitrary Dimension, Progress of Theoretical Physics Supplement **88**, 1 (1986).

[39] L. Parker, Probability distribution of particles created by a black hole, Phys. Rev. D **12**, 1519 (1975).

[40] M. Gasperini, Bogoliubov Transformations and Negative Binomial Distributions for Particle Creation by Gravitational Fields, Progress of Theoretical Physics **84**, 899 (1990).

[41] I. Carusotto and C. Ciuti, Quantum fluids of light, Reviews of Modern Physics **85**, 299–366 (2013).

[42] A. Amo, J. Lefrère, S. Pigeon, C. Adrados, C. Ciuti, I. Carusotto, R. Houdré, E. Giacobino, and A. Bramati, Superfluidity of polaritons in semiconductor microcavities, Nature Physics **5**, 805–810 (2009).

[43] A. Amo, D. Sanvitto, F. P. Laussy, D. Ballarini, E. d. Valle, M. D. Martin, A. Lemaître, J. Bloch, D. N. Krizhanovskii, M. S. Skolnick, C. Tejedor, and L. Viña, Collective fluid dynamics of a polariton condensate in a semiconductor microcavity, Nature **457**, 291–295 (2009).

[44] A. Amo, S. Pigeon, D. Sanvitto, V. G. Sala, R. Hivet, I. Carusotto, F. Pisanello, G. Leménager, R. Houdré, E. Giacobino, C. Ciuti, and A. Bramati, Polariton superfluids reveal quantum hydrodynamic solitons, Science **332**, 1167–1170 (2011).

[45] D. Sanvitto, F. M. Marchetti, M. H. Szymańska, G. Tosi, M. Baudisch, F. P. Laussy, D. N. Krizhanovskii, M. S. Skolnick, L. Marrucci, A. Lemaître, J. Bloch, C. Tejedor, and L. Viña, Persistent currents and quantized vortices in a polariton superfluid, Nature Physics **6**, 527–533 (2010).

[46] Q. Fontaine, D. Squizzato, F. Baboux, I. Amelio, A. Lemaître, M. Morassi, I. Sagnes, L. Le Gratiet, A. Harouri, M. Wouters, I. Carusotto, A. Amo, M. Richard, A. Minguzzi, L. Canet, S. Ravets, and J. Bloch, Kardar–parisi–zhang universality in a one-dimensional polariton condensate, Nature **608**, 687–691 (2022).

[47] D. Caputo, D. Ballarini, G. Dagvadorj, C. Sánchez Muñoz, M. De Giorgi, L. Dominici, K. West, L. N. Pfeiffer, G. Gigli, F. P. Laussy, M. H. Szymańska, and D. Sanvitto, Topological order and thermal equilibrium in polariton condensates, Nature Materials **17**, 145–151 (2017).

[48] F. Claude, S. V. Koniakhin, A. Maître, S. Pigeon, G. Lerario, D. D. Stupin, Q. Glorieux, E. Giacobino, D. Solnyshkov, G. Malpuech, and A. Bramati, Taming the snake instabilities in a polariton superfluid, Optica **7**, 1660 (2020).

[49] I. Amelio and I. Carusotto, Perspectives in superfluidity in resonantly driven polariton fluids, Physical Review B **101**, 10.1103/physrevb.101.064505 (2020).

[50] D. V. Reddy, M. G. Raymer, C. J. McKinstrie, L. Mejling, and K. Rottwitt, Temporal mode selectivity by frequency conversion in second-order nonlinear optical waveguides, Opt. Express **21**, 13840 (2013).

[51] D. V. Reddy, M. G. Raymer, and C. J. McKinstrie, Efficient sorting of quantum-optical wave packets by temporal-mode interferometry, Opt. Lett. **39**, 2924 (2014).

[52] The datasets corresponding to the figures of this article are available at DOI 10.57745/7MIDRA.

[53] M. Pasienski and B. DeMarco, A high-accuracy algorithm for designing arbitrary holographic atom traps, Opt. Express **16**, 2176 (2008).

[54] K. Ouahrouche, SLMUtils v1.0.0, https://doi.org/10.5281/zenodo.13867920 (2024).

[55] R. Loudon, *The Quantum Theory of Light*, 3rd ed. (Oxford University Press, 2000).