[1]Federico Presutti [1,2]Logan G. Wright Department of Applied Physics, Yale
University, CT, USA [1,4]Peter L. McMahon 1]School of Applied and Engineering
Physics, Cornell University, NY, USA 2]NTT Physics and Informatics
Laboratories, NTT Research, Inc., CA, USA 3]Intelligence Community
Postdoctoral Research Fellowship Program, Cornell University, NY, USA 4]Kavli
Institute at Cornell for Nanoscale Science, Cornell University, NY, USA

#  Highly multimode visible squeezed light with programmable spectral
correlations through broadband up-conversion

<EMAIL_ADDRESS><EMAIL_ADDRESS>Shi-Yuan Ma  Tianyu Wang  Benjamin K.
Malia  Tatsuhiro Onodera<EMAIL_ADDRESS>[  [  [  [

###### Abstract

Multimode squeezed states of light have been proposed as a resource for
achieving quantum advantage in computing and sensing. Recent experiments that
demonstrate multimode Gaussian states to this end have most commonly opted for
spatial or temporal modes, whereas a complete system based on frequency modes
has yet to be realized. Instead, we show how to use the frequency modes
simultaneously squeezed in a conventional, single-spatial-mode, optical
parametric amplifier when pumped by ultrashort pulses. Specifically, we show
how adiabatic frequency conversion can be used not only to convert the quantum
state from infrared to visible wavelengths, but to concurrently manipulate the
joint spectrum. This near unity-efficiency quantum frequency conversion, over
a bandwidth >45 THz and, to our knowledge, the broadest to date, allows us to
measure the state with an electron-multiplying CCD (EMCCD) camera-based
spectrometer, at non-cryogenic temperatures. We demonstrate the squeezing of
>400 frequency modes, with a mean of approximately 700 visible photons per
shot. Our work shows how many-mode quantum states of light can be generated,
manipulated, and measured with efficient use of hardware resources – in our
case, using one pulsed laser, two nonlinear crystals, and one camera. This
ability to produce, with modest hardware resources, large multimode squeezed
states with partial programmability motivates the use of frequency encoding
for photonics-based quantum information processing.

Figure 1:  Frequency domain, multimode, visible squeezed state preparation and
detection. a. Highly multimode squeezed vacuum is generated in the degenerate
optical parametric amplifier (DOPA), at near infrared wavelengths. The
squeezed modes occupy orthogonal frequency spectra. This may be represented by
a set of squeezing operators $\hat{S}_{i}$ acting the vacuum state $\lvert
0\rangle$. Each operator squeezes a distinct frequency mode with some
squeezing parameter $r_{i}$. The squeezing occurs around a central frequency
$\omega_{0}$ equal to half the pump central frequency $2\omega_{0}$.
Concretely, these correspond to the wavelengths 1550 nm and 775 nm,
respectively. b. Adiabatic frequency conversion (AFC) efficiently converts the
squeezed light to visible wavelengths. The pump, a broadband pulse centered at
1033 nm ($3\omega_{0}/2$), combines with the broadband 1550 nm signal to yield
620 nm ($5\omega_{0}/2$) light. This transformation between the two sets of
frequencies can be changed by varying the spectro-temporal profile of the
pump. This operation may be represented as a linear unitary transformation $U$
acting on the infrared and visible fields, represented by the operators
$\hat{a}_{\omega_{\text{ir}}}$ and $\hat{a}_{\omega_{\text{vis}}}$. c. The
final state, incident on a diffraction grating (DG), is split into frequency
modes, and frequency-resolved photon counting is performed with an EMCCD
camera. Each measurement yields some photon-number sequence, or vector,
$\accentset{\rightharpoonup}{n}$, whose probability distribution depends on
the state (determined by $\hat{S}_{i}$ and $U$). The camera measures the
spectrum in a discrete manner, as the pixels capture the photons within some
frequency “bins,” denoted $\Delta\omega_{i}$. Hence, the state’s overlap with
these bin-basis Fock states, $\lvert n_{\Delta\omega_{i}}\rangle$, determines
the probability distribution.

The generation, control, and measurement of entangled multimode Gaussian
states of light are crucial elements of continuous-variable (CV) quantum
information processing [43, 62, 67, 6]. Most quantum technologies based on
multimode quantum optics benefit from being able to use as many modes as
possible. As an example, a Gaussian boson sampler (GBS) [52, 21] is a special-
purpose quantum computer that can – at least in the ideal case – perform
certain calculations that are believed to be intractable on a classical
computer when the number of modes and number of photons in the GBS are
sufficiently large [21]. The recent demonstration of Gaussian boson sampling
in the regime of quantum computational supremacy, with tens to hundreds of
squeezed modes and detected photons per shot [73, 72, 13, 42], was a milestone
in the development of CV-based quantum systems that was achieved because of
the success in pushing to large numbers of modes and photons. A GBS executes a
sequence of three steps, which are common to many CV quantum-information-
processing protocols: (1) generate squeezed states, (2) apply a unitary
transformation to entangle them, and (3) measure the final state (by photon
counting).

Optics gives us the choice of several possible degrees of freedom within which
to encode information – most importantly: space, time and frequency (or any
combination thereof). While large-scale GBS experiments have been realized
using space [73, 72, 13] and time [42] encodings, the frequency domain remains
to be explored. Frequency encoding offers important potential advantages over
space or time encoding for both the generation and the manipulation (unitary
control) of multimode squeezed light: reduced hardware resources and
complexity, and reduced loss. The extremely broad bandwidth of light enables
frequency-encoded systems to operate on many frequency modes in a compact
system [39, 27, 29, 22, 53, 7, 2]. Many demonstrations of large-scale
multimode squeezing, for example, use the frequency domain, e.g., in broadband
frequency combs [53, 7, 35, 49, 69, 55, 8, 54, 70, 25, 51, 10, 68, 9, 16].
Linear unitary operations in the frequency domain, i.e., acting on the
frequency modes, can be implemented in a hardware-efficient way, operating on
all frequency modes in parallel [37, 38]. One approach is to use one or more
electro-optic modulators [26, 71, 37, 38] (although there are limitations on
the unitary from the driving microwave bandwidth being $\lesssim$100 GHz);
another is to use nonlinear wave mixing to convert photons in each frequency
mode to photons in a combination of other frequency modes [36]. Unitaries
based on nonlinear wave mixing mediated by a classical field, such as four-
wave mixing [28], provide a route to realizing programmable unitaries that can
operate over wide bandwidths, in compact hardware with low loss. The
programming of the unitary can be done by shaping the classical field(s) used
to control the wave mixing of the modes containing quantum light, and the wave
mixing can be implemented compactly in a single-spatial-mode waveguide.
However, a programmable, frequency-domain unitary working over optical ($>$1
THz) bandwidths for quantum light has yet to be realized at even moderate
scale (more than 2 modes [28]).

The final step, state measurement by photon counting, is a challenge for many
multimode architectures. Since the preferred nonlinear optical materials for
generating squeezed vacuum work best at longer wavelengths, squeezing at
wavelengths centered near 1550 nm is typical [73, 72, 13, 42, 22, 53, 7, 2,
35, 49, 69, 55, 8, 54, 70, 25, 16], also in part due to the convenience of
being able to use optical components from telecommunications technologies.
However, good (high quantum-efficiency, low dark-count) single-photon
detectors at 1550 nm, namely superconducting nanowire detectors [46, 17], are
very expensive and require cryogenic cooling. Silicon-based camera sensor
technologies – both charge-coupled device (CCD) and complementary
metal–oxide–semiconductor (CMOS) detectors – are well-established,
comparatively inexpensive and compact, and each camera comprises
$10^{5}$–$10^{6}$ individual pixels, i.e., detectors. Cameras capable of
detecting single photons with low noise have recently become available, and
there is a growing literature of quantum-optics and sensing experiments that
were performed with visible wavelengths and these cameras [11, 4, 12, 30, 34,
61]. Electron-multiplying CCD (EMCCD) cameras are arguably the current state-
of-the-art, and low-noise CMOS photon-number-resolving cameras are also a
promising tool within this domain [14, 41].

Here, we demonstrate how to use frequency conversion [31] to enable the use of
these visible-light cameras in combination with techniques for strong
squeezing possible at longer wavelengths. In addition, our demonstration is an
improvement over existing methods of quantum frequency conversion: previous
demonstrations are limited to either modest bandwidths or efficiencies [65,
64, 3, 57, 5, 1, 59]. However, adiabatic frequency conversion (AFC)
essentially eliminates this trade-off [60, 45]. We show how this method allows
us to obtain robust, efficient and broadband conversion over >45 THz
(1390–1750 to 590–650 nm) and near-unity efficiency. Furthermore, it allows
unitary control of the multimode entanglement (with, in principle, no
additional loss) by manipulation of the complex profile of the broadband pump
used to drive the conversion. This architecture provides the best of both
worlds: squeezing at telecommunications wavelengths, and photon detection at
visible wavelengths. We will show that we are able to generate strong
squeezing in over 400 frequency supermodes, resulting in states having a
measured mean photon number of nearly 700. By using AFC to efficiently convert
the squeezed light to visible wavelengths, and using the highly parallel
photon counting made possible by a modern EMCCD, we can directly measure these
states. We will also show that we can control the entanglement between
different modes by using different spectrally shaped classical fields as the
pump of the AFC process, resulting in different measured correlations between
photon detections across the frequency modes.

## Results

The experimental setup is illustrated in Fig. 1; an overview is as follows. We
use a waveguided degenerate optical parametric amplifier (DOPA) pumped with a
pulsed laser: this provides squeezing in a single spatial mode and over many
frequency modes (known as the “supermodes”). An adiabatic frequency conversion
(AFC) crystal subsequently converts this near-infrared squeezed light to the
visible. The temporal profile of the AFC pump pulse is shaped, which controls
the conversion process – the linear transformation between infrared and
visible light frequencies. Finally we detect the visible squeezed light with
an electron-multiplying CCD (EMCCD) camera, serving as an array of high-
quality photodetectors.

In the following section, we present characterizations of each of these
stages. We first measure the properties of the squeezed light in the infrared,
directly out of the DOPA. Our characterizations enable us to measure the
bandwidth of the squeezed light, estimate the number of modes, and verify
squeezing. We then demonstrate efficient conversion using AFC. In addition,
with the camera, we are able to observe spectral photon-number correlations
throughout the whole bandwidth at high resolution. With this ability, we
generate qualitatively different joint spectra as a proof of concept of
frequency-domain unitary transformations by pulse shaping. Finally, we discuss
detection and the photon numbers generated by this architecture. Refer to the
Methods section for details on the DOPA, AFC, pulse shaper, the photon-
counting spectrometer, as well as specific experiments.

### Highly multimode squeezing in the frequency domain

Figure 2:  Frequency-multimode squeezing in a lithium niobate waveguide
optical parametric amplifier. a. Simulated photon-number covariance matrix of
the state produced in the degenerate optical parametric amplifier (DOPA). This
is made up of classical “thermal” correlations (due to photon-number variance)
and non-classical entanglement correlations. b and c show the decomposition of
the simulated state into single-mode squeezed states, referred to as the
supermodes. The plots share the vertical axis, which indexes these modes. b.
Simulated average photon number ($\langle n\rangle$) and squeezing values of
the first 1000 supermodes. c. Simulated supermode basis of the DOPA. Each row
represents a mode and its spectrum is represented by the colormap. The high
gain supermodes are approximately Hermite–Gauss functions while the weaker
supermodes resemble two-mode squeezed vacuum (TMSV) states with two distinct
frequency peaks. For a-c, the pump power is approximately equivalent to the
highest used in the experiment. d. Measurement of the DOPA output spectrum and
comparison to simulation. As the detector quantum efficiency (QE) vanishes at
longer wavelengths, the simulation spectrum is also plotted multiplied by the
nominal relative QE for comparison. e. Parametric gain of the DOPA, fit to the
number of modes, and comparison to simulation. The photon number $\langle
n\rangle$ is obtained by dividing the average power measurement by the
repetition rate and the photon energy at the central wavelength. In the fit:
$M$ represents the number of modes; $\eta$ the detection efficiency; $P$ the
pump power; $P_{0}$ the fitted scaling factor. f. Coincidence detection,
comparing squeezed light, attenuated laser light (coherent) and broadband
incandescent lamp (thermal). The prediction model is derived in Appendix C.
Sketches of the experimental setups corresponding to c–e are shown below their
respective plots. The DOPA output is measured with either a single-photon
avalanche diode (SPAD) or a power meter (PM), after isolating the signal with
a long-pass filter (LPF).

We begin by studying the degenerate optical parametric amplification (DOPA)
process, and looking for behavior consistent with highly multimode squeezing.
Three experiments were performed to verify numerical predictions: a spectrum
measurement, a parametric gain measurement and a coincidence detection
measurement. The purpose of the first is to verify that the bandwidth is
indeed as broad as predicted. The second is to infer the number of modes,
based on how the output energy scales with pump energy. The third is to verify
that this light is made up of squeezed states, through photon statistics.

The DOPA was first simulated, and the photon-number covariance is plotted in
Fig. 2a, which has a notable “X” shape. The anti-diagonal is due to the non-
classical frequency correlations of entangled photons (the joint spectral
intensity, or “entanglement correlations”) and energy conservation dictates
its shape: down-converted photon pair energies must add up to the pump photon
energy. The diagonal represents the photon-number variance, and the
correlations are classical (“thermal correlations,” as these statistics have
the same properties as those of thermal states [40]): this has a finite width,
given by the phase-matching function of the OPA (typically a sinc function for
periodic poling). The Bloch–Messiah decomposition helps us interpret how the
DOPA generates squeezing [66, 23], by reducing the output to a set of
supermodes and their corresponding squeezing or anti-squeezing values. These
are shown in Fig. 2b and c. For a 4 cm LN waveguide, and 775 nm 200 fs pulse
pump, we predict $\sim$600 supermodes with substantial average photon numbers,
spanning over 50 THz, or over 400 nm. The weakly squeezed supermodes (indices
above 200) resemble two-mode squeezed vacuum (TSMV) states, occupying two sinc
peaks on opposite sides of the central wavelength. The highly squeezed
supermodes (indices below 200) appear to approximately be the typical
Hermite–Gauss functions [66]. The supermodes are not constant with respect to
the pump power, and tend to transition from TMSV to Hermite–Gauss as the
parametric gain increases.

The spectrum measurement, whose outcome is shown in Fig. 2d, confirms that the
spectrum is broad, spanning hundreds of nanometers. The spectrum could not be
measured past about 1700 nm due to the detector response cut-off. The spectrum
may be thought of as the incoherent sum of the supermodes, scaled by average
photon number.

The parametric gain experiment (Fig. 2e) measures the output power of the DOPA
at varying pump powers. The number of down-converted photons scales linearly
with the number of supermodes $M$ and the detection efficiency $\eta$, while
each mode responds nonlinearly, containing on average $\langle
n\rangle=\sinh^{2}\sqrt{P/P_{0}}$ photons, for pump power $P$ and scaling
factor $P_{0}$. We may estimate the number of modes by fitting to $\eta
M\sinh^{2}\sqrt{P/P_{0}}$, which assumes equally squeezed modes (not truly the
case, as shown in Fig. 2b, hence this is a lower bound). From this measurement
we infer the presence of at least $M\geq\eta M=431$ squeezed modes – on the
order of 650 squeezed modes accounting for average detector responsivity,
$\eta\approx 71\%$. In order to achieve a sufficently high – measurable –
average output power, this experiment requires different operating conditions
(a higher repetition rate pump laser), and we estimate that we would measure a
slightly lower number of modes, 386 (465 corrected), under the usual
conditions (Appendix I).

The coincidence experiment (inspired by Ref. [63]) is shown in Fig. 2f. The
DOPA output beam is evenly split and the photon-number covariance between the
two beams is measured. The coincidences depend on the photon statistics of the
state, and at low photon numbers, only the coincidence rate of highly-
multimode squeezed light is expected to scale linearly with respect to the
average photon number. We performed coincidence detection with the DOPA
output, attenuated 1550 nm laser light, and broadband incandescent lamp light;
indeed only the former was nonzero.

These three experiments verify our predictions that the DOPA is a source of
highly multimode squeezed light, spanning hundreds of frequency modes. A more
thorough discussion of these results, including explanations of the
experiment-theory discrepancies, is presented in the Methods.

### Adiabatic frequency conversion (AFC) of squeezed light

Figure 3:  Broadband, adiabatic frequency conversion (AFC) of squeezed light.
Experiments corresponding to plots a and b use a diffraction grating (DG) to
use the EMCCD camera as a spectrometer. Experiments corresponding to plots c–e
use a Wollaston prism configured as a 50:50 beamsplitter (BS). a. Converted
spectrum measurement. This is compared to simulation under similar conditions
and with 92% conversion within the conversion bandwidth. The shading in the
simulation curve indicates the unconverted energy, to show that the conversion
profile should be approximately flat within the target bandwidth. b. Measured,
wavelength-resolved photon-number covariance matrix. The diagonal and anti-
diagonal features can be recognized as the thermal (classical) and
entanglement (non-classical) correlations. c. Coincidence detection after
conversion. This experiment measures the transmission of the state (the
effective fidelity or purity) as 40%. d. Parametric gain measured after
conversion. The fit to this curve estimates the number of squeezed modes
present in the converted state as 433, accounting for loss. e. Conversion as a
function of pump power. A fit to the saturation curve estimates the conversion
efficiency as 92.5%.

The next part of the experiment involves frequency conversion of the squeezed
light. The infrared squeezed light, centered around 1550 nm, is up-converted
with a pulsed 1033 nm laser, to a central wavelength of 620 nm, by sum
frequency generation. This ultra-broadband nonlinear process is possible
through AFC [60, 45], which acts as a unitary transformation in the frequency
domain, and most importantly allows us to measure highly squeezed light with a
sensitive, silicon-based camera. (See Appendix G for the crystal quasi-phase-
matching design used in this experiment.)

In this section we demonstrate frequency conversion of squeezed light, show
that it remains squeezed and measure the overall system fidelity, and provide
estimates of the conversion efficiency and the number of converted modes. To
this end, two experiments are presented in this section: one based on our
visible light, single photon spectrometer (Fig. 3a-b); and one based on two-
mode coincidence detection (Fig. 3c-d). The spectrometer uses a diffraction
grating to separate the squeezed light into its wavelength components, which
are imaged onto the EMCCD camera. The coincidence experiment also uses the
camera, but a with prism, functioning as a 50:50 beamsplitter, replacing the
grating. The two beams are imaged onto the camera, which emulates two
detectors via two small regions of interest on the CCD.

The spectrum measured with the EMCCD camera visible light spectrometer is
shown in Fig. 3a. The converted spectrum matches the designed conversion band.
The shoulders in the spectrum (roughly between 590–600 and 640–650 nm), which
differ from the simulation, are consistent with Fig. 2c (around 1400 nm):
these features are from the original spectrum, not the conversion. The
intensity implies that, on average, the most brightly-illuminated pixels
receive almost 2 photons per pulse.

The photon-number covariance matrix is shown in Fig. 3b. Shot-by-shot
intensities are obtained by triggering the camera to capture individual
pulses. The camera is used in “analog” mode (not thresholded) due to the high
photon numbers per pixel (see Appendix B). For this experiment, the AFC is
pumped in such a way as to maximize the conversion and to perform a
transformation as close as possible to the identity, although the pump
bandwidth causes some broadening (see Appendix G; an identity would require a
monochromatic pump). The photon covariance matrix has the “X” shape we expect
from Fig. 2a, with both classical and non-classical correlations.

The results of the coincidence experiment, similar to the one presented in the
previous section, are shown in Fig. 3c. We can now use this experiment to
measure the effective fidelity purity or of the state: where the photon
statistics lie between those of pure squeezed vacuum and thermal noise, as
parametrized by an effective loss or transmission. For few photons in each
supermode, $\langle n\rangle\ll 1$, we expect linear scaling; at intermediate
photon numbers, where $\langle n\rangle=\sinh^{2}\sqrt{P/P_{0}}\approx
P/P_{0}$, we expect the covariance to scale quadratically; and at greater
photon number we cannot predict the scaling without precise knowledge of the
supermodes (see Appendix C). In the first two cases, the coefficient of the
linear component gives us an estimate of the overall transmission (a weighted
average over all modes). Hence, we perform a quadratic fit to the lowest
photon number portion of the data. The fit implies a 40% overall effective
transmission for all the modes.

The same data is used for Fig. 3d, to quantify the parametric gain. By
assuming the 40% overall (effective) transmission, the scaling implies the
presence of over 430 squeezed modes. Some parasitic second harmonic generation
is observed in the DOPA at higher pump energies, which could reduce the pump
energy available for squeezing and possibly cause this value to be slightly
underestimated.

The conversion as a function of AFC pump power is used as an estimate of the
conversion efficiency, which is shown in Fig. 3e. Conversion efficiency
(saturation) is described by the Landau–Zener formula: unconverted signal
follows an exponential decay with respect to the pump power. The fit implies
92.5% conversion; this is a reasonable upper-bound estimate on the conversion
efficiency. We have been careful to make the pump beam much larger than the
signal, and have verified that it has minimal spatial chirp, as a fit to the
Landau–Zener formula may overestimate the efficiency if there is substantial
mismatch in the spatio-temporal overlap in the pump and signal. Nonetheless,
we can obtain a lower bound estimate from the 40% total fidelity or effective
transmission. By measuring the transmission of each optical element in the
squeezed light beam path and taking into account the camera quantum efficiency
(QE), we estimate the end-to-end passive optical transmission to be at least
60% (excluding AFC). This bounds the conversion efficiency from below to at
least $\sim$66%. This is a fairly loose bound, as the fidelity figure includes
non-loss decoherence (e.g. detector excess noise; weakly converted photons at
wavelengths outside the conversion band) and the combined linear loss in the
beam path probably exceeds 40%. Finally, we point out that the conversion
efficiency is pump-power limited, meaning that achieving close to 100% over
these bandwidths is realistic without requiring phase-locking.

### Upconversion as a frequency-domain unitary transformation

Figure 4:  Preparing the joint spectrum via frequency conversion. By pulse
shaping the AFC pump $A(\omega_{\text{p}})$, we can modify the linear unitary
transformation $U$ performed by the AFC process, hence the state and the
measured spectral photon-number correlations. The pump spectral intensity and
phase both play a role in the transformation. a. Pulse shaping involves two
operations: intensity and phase modulation. The former
($\mu(\omega_{\text{p}})$) changes the spectral intensity of the pump, which
in turn mainly affects the magnitude of elements of $U$, and the latter
($\phi(\omega_{\text{p}})$) changes the spectral phase of the pump, which
mainly affects the complex phase of the elements of $U$. b. The phase-averaged
linear transformation performed by the AFC. The rows of the unitary strongly
resemble the pump spectrum, as shown in the subplots: the first pair shows how
it compares to the original spectrum, and the second is an example of how it
changes given an intensity-modulated pump. The frequency-difference axis is
shared among the pump and unitary-row subplots. c. These plots show how a
phase-modulated pump affects the conversion process and produces more
complicated correlation structure. Example spectral intensity correlations are
plotted on the right, and the corresponding (inferred) pump intensity profiles
on the left. Specifically, we plot the entanglement contribution to photon
correlation matrices (correlation matrices with thermal-like components
subtracted). The first shows the same data in Fig. 3b, for reference, with no
pulse shaping. The next three are the result of random pulse shaping. The
corresponding pump intensities, plotted on the right, are not measured
directly, but estimated based on the phase modulation applied and the pulse
shaper calibration.

Most quantum sensing or computing protocols require unitary operations to be
performed on the overall state, because each application requires some
specific entanglement structure. Here we demonstrate how the covariance matrix
can be transformed by implementing frequency unitaries through the broadband
sum frequency generation process. The unitary is applied during the
conversion; shaping the AFC pump pulse affects the nonlinear dynamics,
changing which wavelengths convert to which. This means that, within certain
constraints, we can program the pulse to achieve a desired transformation.

Intuitively, the relative phase of a pump frequency affects the phase of a
signal it converts, thereby constructive and destructive interference may
promote or suppress the conversion of a signal from one frequency to another,
effectively forming a frequency-domain interferometer network. This is
illustrated in Fig. 4a-b. We discuss this further in Appendix D.

Some examples of this process and how it influences the covariance matrix are
shown in Fig. 4c. Phase modulations are applied with the pulse shaper, on top
of a fixed quadratic phase (chirp) to guarantee a certain pump pulse duration.
The predicted pulse shapes (inferred from the applied phase modulation) are
shown for each example. For clarity, we show the correlation rather than
covariance matrices, and we subtract a fit to the thermal (classical) part in
order to show only the entanglement contribution (the joint spectrum). For
reference, Fig. 4a is derived from the same data as Fig. 3b. As we do not
measure the relative phases in the covariance, we focus on transformations in
the joint intensity.

In our experiment the bandwidth of the pump (<5 THz) is much smaller than the
bandwidth of the squeezed light ($\sim$47 THz), which prevents all-to-all
coupling (the pump frequencies mediate the change in signal frequency).
Additionally, the AFC pump peak intensity is reduced by pulse-shaping,
therefore reducing the efficiency. Despite these limitations, it is possible
to achieve qualitatively different joint spectra by simply changing the pulse
shape.

### Multimode quantum state sampling at visible wavelengths

Figure 5:  Parallel single-photon detection for multimode-quantum-state
measurement: sampling the many-mode, many-photon distribution. a. Receiver
operating characteristic (ROC) curve for the EMCCD camera’s detectors. This
quantifies the trade-off between false click rate and photon-detection
efficiency (PDE). b. Histogram of the total photon number per event: analog
(noisy) and thresholded. We use the term photoelectrons to refer to the
amplified charge on the CCD divided by the gain. Some example, high-photon-
number events, are labeled on the histogram with callouts i–iv; these are
referenced in the subsequent subfigures. c. Example unprocessed camera image
from one sampled event. d. Samples integrated vertically (left) and
thresholded then summed vertically (right). The experimental schematic
(bottom) refers to the experimental configuration for collecting the data
shown in b–d. The data in a consists of dark frames collected with a closed
shutter (see Appendix B and H.)

The final ingredient in any quantum-optical computing or sensing protocol is
measurement, and in this work we focus on photon-counting measurements. Any
quantum advantage for computing or sensing typically scales with the number of
photons detected so high average photon numbers and high detection efficiency
are crucial. Here we demonstrate sampling of the photon-count distribution of
our generated states using an EMCCD camera.

We used 5 rows and all 512 columns of a 512$\times$512 pixel EMCCD camera,
whose CCD sensor has a QE of $\sim$95% at 620 nm [48]. The receiver operating
characteristic (ROC) curve for thresholding for this camera is shown in Fig.
5a: this characterizes the false click rate against the photon-detection
efficiency (PDE), parametrized by the threshold value (Appendix B). The green
marker indicates the threshold used to generate the subsequent plots,
resulting in a competitive photon-detection efficiency of $\sim$80%.

Fig. 5b shows a histogram of the number of photons detected per shot. Blue
(left) is prior to thresholding, green (right) is with thresholding. The gain
noise in EMCCDs is large, which makes it practically impossible to distinguish
between photon numbers $\geq$1 incident on one pixel. However, statistical
averages tend to be accurate, as this is zero mean noise [32] (see Appendix
B). Henceforth, we use the term photoelectrons to refer to the amplified
charge on the CCD divided by the gain: since the gain is a noisy process, this
is not quantized. The first histogram implies an average of almost 700 photons
per shot; with thresholding the average number of clicks is just under 500 per
shot. The discrepancy is due to both the effective QE and the high rate of
multiple photons incident on one pixel. In Appendix G we analyze the
contribution of dark counts and photons from parasitic processes: these are
not significant in number compared to the squeezed light photons.

For these frequency-resolving experiments we use more than one pixel-rows of
the camera. This is done in order to capture all the photons, due to the
spectrometer’s point spread function occupying a space larger than one pixel
(Gaussian width of $\sigma_{\text{PSF}}\approx 0.6$ pixels). An example of a
raw sample is shown in Fig. 5c, in units of photoelectrons. With this
configuration we could sample at a rate of just over 800 Hz. Using a single
row would allow an average sampling rate of 15.7 kHz, which is limited by data
readout times. Fig. 5d shows samples integrated vertically, for each frequency
bin mode. The point spread function in the vertical direction can allow a
thresholding camera to act as a pseudo-photon-number resolving detector, by
illuminating more pixels (via astigmatic focus) (an established method [18,
50, 58] used e.g. in Refs [13]). The point spread function in the horizontal
dimension, however, acts as an effective decoherence and must be carefully
engineered to fit within the dimensions of a pixel (see Appendix E).

## Discussion

The creation, manipulation, and detection of highly squeezed, highly multimode
entangled states are important ingredients of many continuous-variable (CV)
schemes for quantum computing, sensing, and communication using photonics.
When compared to space or time encoding, frequency encoding has significant
advantages for integration and scalability, but it is not always practical to
perform unitary transformations on frequency modes without substantial loss.
We have experimentally demonstrated efficient, broadband, quantum frequency
conversion of highly multimode squeezed light generated in the near infrared
into the visible using adiabatic frequency conversion (AFC). This
simultaneously allows efficient, parallel photon counting across over 400
squeezed and 500 detection modes using CCD-based photon-detector arrays, and
the programming of frequency correlations in the multimode squeezed light. Our
approach requires no active phase locking, uses a single optical beam path for
the quantum light, and the number of modes and shot rates should be scalable
well beyond what we have demonstrated.

To the best of our knowledge, the quantum-optical states we produced are the
largest partially programmable, photon-counted, multimode squeezed states
produced by a moderate factor (about 2–4 times larger than previous results,
albeit with different limitations and caveats), and by an order of magnitude
the largest in the frequency domain (see Appendix J). Our work provides a path
to constructing large-scale Gaussian boson samplers using frequency encoding.
111The experiments we report contain the key parts of a Gaussian boson sampler
(GBS): squeezing of multiple modes, a unitary, and photon detection. Why then
do we not declare our experiments a realization of the world’s largest GBS
already? The issue is we have imperfect calibration of the squeezed modes and
of the unitary transformation realized by AFC, making it unreasonable for us
to compare the experimental sampling results we obtained with samples one
could obtain from a simulator of our experiment running on a classical
computer. Since the simulator would not have a sufficiently accurate
description of the modes and the unitary, one wouldn’t expect the samples – or
the statistics of the samples – to match those from the experiments, which
would prevent us from verifying the operation of the GBS in the way that has
been done by prior studies [73, 72, 13, 42]. In essence, we have demonstrated
a large-scale but poorly calibrated GBS.

There are currently limitations to our source (the combination of our
generator of multimode squeezed light, followed by the AFC, which performs a
unitary), particularly with respect to control of the programmable unitary.
The class of unitaries we could experimentally realize were limited to fairly
local connectivity (approximately to 10% of the nearest frequency modes):
hence we could programmably entangle each mode to at most $\sim$20% of the
rest (Fig. 4b-c). Ideally we would like to be able to realize unitaries with
all-to-all connectivity. This may be possible with broader AFC pump bandwidths
and higher pump intensities (as shown in simulations presented in Appendix K).
In addition, although we have shown control of the joint spectral probability
distributions of the state, we currently lack a prescription for what pattern
we should program on the pulse shaper for the AFC pump pulse to realize a
specific unitary. In principle, this is achievable if our experimental setup
were to have a sufficiently fast means to characterize the realized unitary,
allowing us to optimize the pulse-shape patterns to realize a given unitary.
This can be resolved in the future by building a different characterization
setup, for example involving two in-phase, fast-wavelength-sweeping
continuous-wave lasers: this would allows to perform fast, phase-sensitive
measurements of our source. There are also imperfections with the detection
setup, as discussed in Appendix H. Most importantly, the frequency binning is
suboptimal, which may be a cause of decoherence. A disadvantage of using an
EMCCD camera instead of, for example, an array of superconducting nanowire
single-photon detectors for measurement is speed: the camera’s frame rate of
800 Hz is slow. However, the EMCCD camera we used has a linerate of 2 MHz (up
to 3.33 MHz with higher noise), which is in principle the relevant speed when
only using a single row (or few rows) of pixels, as is the case for our
measurements. However, there is currently a practical bottleneck to reading
out the data, which is limited to line readouts at up to 15.7 kHz on our
camera model. Several modern EMCCD cameras, including ours, allow reading out
bursts of frames (or regions of whole frames, including lines) at up to MHz
rates for short time windows; this could allow fast detection for a limited
number of shots in the future. Finally, the photon-number sampling we
demonstrated used threshold detection or pseudo-photon-number-resolving
measurements; ideally, future implementations could be capable of true photon-
number-resolving measurements up to high photon numbers per mode (e.g., using
ultra-low-noise CMOS cameras [14, 41]).

Despite these current limitations, the experiment we reported already has some
distinct advantages over other platforms. First of all, to the best of our
knowledge, our experiment marks the first instance of simultaneous sampling
with 500 detectors with a similar number of squeezed modes. In comparison to
spatial-domain architectures that reach a count of 100s of modes with
superconducting detectors (Appendix J) our approach is far less complicated
and less expensive to realize for the same number of modes. Second, a crucial
consideration in building large quantum-optical systems is loss: while, as far
as we are aware, this has not been rigorously studied, we speculate that in
practice frequency unitaries may be far less prone to loss (as do the authors
of Refs. [39, 27, 29, 22, 53, 7]). Performing arbitrary frequency unitaries
requires hardware resources that scale as $M^{2}$ (in the number of modes
$M$), just as is the case for spatial unitaries. However, frequency unitaries
realized using our AFC method place the burden of scaling on the classical
driving pump, in energy and spectro-temporal complexity requirements (as well
as on the interaction length in the AFC crystal, which may also need to
increase as a function of $M$). Much of the $M^{2}$ resources can implemented
using pixels of a spatial light modulator in the pulse shaper of the classical
pump – which has the advantage both that spatial light modulators with
millions of pixels are commercially available, and that any losses incurred in
shaping the pump pulse are far less deleterious than losses in a typical
spatial-unitary implementation, where the quantum light passes through $M^{2}$
controllable elements, and the loss of each will directly degrade the quantum
state. It is important to note, however, that arbitrary frequency
transformations may require full control of the AFC crystal dispersion
(Appendix D). It appears possible to substantially increase the number of
modes and detected photons from what we have already demonstrated, as well as
the complexity of the programmed unitaries; crucially, for all these
enhancements, the total system transmission is not expected to change
significantly, which would be in marked contrast to most space- and time-
domain implementations.

We will conclude with some thoughts on future research that could build on
what we have already shown. First, engineering the bandwidth of the squeezed
light to be compatible with the bandwidth of the pump will enable all-to-all
correlations. Second, there are many opportunities for improvements using
integrated, on-chip, devices [47, 33]. For example, one could implement the
combined DOPA and AFC source on a single chip, with single-spatial-mode
waveguiding; this could eliminate detrimental spatio-temporal effects in AFC.
It would allow for longer interaction lengths than in free space (limited by
Rayleigh length) and allow for high intensities with less power, both
resulting in more powerful unitary control. Combining the DOPA and AFC
processes into one device would also eliminate most interface losses of the
quantum light prior to detection. Using integrated photonics would also allow
the engineering of dispersion to optimize the squeezing bandwidth, the number
of modes in the DOPA, and the interaction complexity imparted by AFC pump.
Frequency-encoded squeezed light is not restricted to using supermodes of a
continuous basis: quantum states can be engineered to inhabit a (more
conventional) basis of (ultrafast) frequency combs lines or discrete frequency
bins [49, 24, 15, 44, 19], through dispersion and pump engineering. Similarly,
shaping the DOPA pump (and possibly the domain poling) to engineer the DOPA
joint spectrum [49, 24, 15, 44, 56] provides additional avenues for
programmability, including the squeezing level of each individual supermode.
Finally, understanding how to program the AFC unitary, e.g. through the
development of more quantitative models and better calibration procedures, is
a valuable avenue of research; it is also important to understand what parts
of the state (Hilbert) space can be reached through DOPA and AFC pulse
shaping, as well as potentially incorporating additional frequency-domain
operations (Appendix D discusses the theory of AFC transformations). Overall
we hope our work will enable more widespread study of many-mode, many-photon
entangled quantum states, and provide a useful building block for large-scale
frequency-encoded CV quantum technologies.

Acknowledgments

We acknowledge helpful discussions with Noah Flemens, Alexey Gorshkov and
Mandar Sohoni. We thank Valeria Cimini, Alen Senanian, James Williams and
Ryotatsu Yanagimoto for helpful comments regarding the manuscript.

Funding We thank NTT Research for their financial and technical support.
Portions of this work were supported by the National Science Foundation (award
CCF-1918549) and a David and Lucile Packard Foundation Fellowship. P.L.M.
acknowledges membership of the CIFAR Quantum Information Science Program as an
Azrieli Global Scholar. This research was supported in part by an appointment
to the Intelligence Community Postdoctoral Research Fellowship Program at
Cornell University, administered by Oak Ridge Institute for Science and
Education through an interagency agreement between the U.S. Department of
Energy and the Office of the Director of National Intelligence.

Data availability Experimental data and scripts to replicate the figures in
this paper based on this data, are available at
http://doi.org/10.5281/zenodo.10492956.

Code availability All simulations were performed using our software library
for nonlinear-optical design, simulation, and Gaussian quantum optics, which
is found at github.com/mcmahon-lab/MultimodeNonlinearOptics.

Author contributions F.P. and L.G.W. designed the experimental setup. F.P.
performed the simulations; designed the nonlinear optics; and built and ran
the experiment with assistance from L.G.W. B.K.M. and S.-Y.M. S.-Y.M.
evaluated and calibrated photon-counting cameras for multimode squeezed light
detection. T.W. designed and tested the EMCCD-camera spectrometer. F.P.,
L.G.W., T.O. and S.-Y.M. worked on the quantum-optical modeling of the
experiment. P.L.M. and L.G.W. conceived and supervised the project. All
authors contributed to preparing the manuscript.

## Methods

### Highly multimode squeezing in the frequency domain: details

To study the squeezing process, the DOPA is first simulated with a split step-
based partial differential equation solver, based on material properties of
LN, and the pump laser (such as dispersion and nonlinearity; pulse duration
and energy). The Green’s function, which describes how the DOPA acts on the
vacuum, is computed with this solver. The Bloch–Messiah decomposition is then
applied: this helps us interpret how the DOPA generates squeezing [66, 23], by
reducing the Green’s function to a set of supermodes and corresponding
squeezing or anti-squeezing values. These are shown in Fig. 2a and b. For a 4
cm LN waveguide, and 775 nm 200 fs pulse pump, we predict $\sim$600 supermodes
with substantial average photon numbers, spanning over 50 THz, or over 400 nm.
Refer to Appendix A for the DOPA equations of motion and the relationship
between the Green’s function and the Bloch–Messiah supermodes.

Three experiments were subsequently performed to verify the numerical
predictions: a spectrum measurement, a parametric gain measurement and a
coincidence detection measurement. The purpose of the first is to verify that
the bandwidth is indeed as broad as predicted. The second is to infer the
number of modes, based on how the output energy scales. The third is to verify
that this light is made up of squeezed states, through photon statistics. The
spectrum may be thought of as the incoherent sum of the supermodes, scaled by
average photon number, while for the parametric gain, we measure the total
photon number as a function of pump power.

The spectrum is measured by monochromation of the DOPA output beam to 1 nm
resolution and detecting with a single-photon avalanche diode (SPAD). The
spectrum could not be measured past about 1700 nm, due to the SPAD’s InGaAs
detector response cutoff. As predicted by simulation, the output is broad,
spanning hundreds of nanometers (see Fig. 2c). The spectrum is in fact broader
than simulation, with a shoulder appearing and extending beyond the prediction
around 1400 nm: we believe this might be due to pump higher order spatial
modes phase-matching the OPA process at the outer wavelengths. The 775 nm
pump, unless otherwise specified, is generated from a 200 fs, 1033 nm laser by
nonlinear frequency conversion.

The parametric gain experiment involves measuring the output power of the DOPA
with a power meter at varying pump powers (see Fig. 2d). A high repetition
rate laser is necessary to obtain a sufficiently high average output power,
hence an 80 MHz Ti:S laser is used to pump the DOPA. The data is fit to the
function $\eta M\sinh^{2}\sqrt{P/P_{0}}$, where $\eta$ represents overall
transmission, $M$ number of modes, $P$ pump power, and $P_{0}$ power scaling
factor (a single squeezed state has $\langle n\rangle=\sinh^{2}r$, with
squeezing parameter $r\propto\sqrt{P}$). This fit assumes $M$ equally squeezed
modes each contributing the same number of photons (thus, a lower bound).
Remarkably, the simulation, which includes the responsivity of the detector as
a function of wavelength, scales ostensibly the same (albeit with a different
nonlinear efficiency $P_{0}$). From this measurement we can infer the presence
of at least 430 squeezed modes – on the order of 650 modes if accounting for
detector responsivity. The Ti:S pulses, 100 fs, have more bandwidth than the
200 fs upconverted 775 nm. This is accounted for in simulation, which predicts
that the squeezing bandwidth will be moderately broader, hence a marginally
larger $M$ than in the 200 fs case, but the central squeezed modes are not
significantly affected. If a 200 fs experiment were possible at a lower
repetition rate, we predict a slightly lower measurement of $\eta M$ (see
Appendix I for comparisons).

The coincidence experiment (inspired by Ref. [63]) involves splitting the DOPA
output beam with a broadband 50:50 beamsplitter. Each beam is detected by a
gated SPAD, with the goal of measuring the photon-number covariance between
the two beams. The coincidences depend on the photon statistics of the state
(namely the bunching, which may also be thought of in terms of g(2) or Mandel
Q). Fig. 2e shows the covariance between SPAD clicks against the average click
rate (as these are threshold detectors). At low photon numbers, only the
coincidence rate of highly-multimode squeezed light is expected to scale
linearly with respect to the average photon number (see Appendix C for
details). To verify this, we performed coincidence detection with the DOPA
squeezed light, attenuated 1550 nm laser light, and broadband incandescent
lamp light. Indeed, only the squeezed light has nonzero, linear covariance
where $\langle n\rangle\ll 1$. As the thresholding detectors yield Bernoulli
random variables, the covariance between the two is restricted to a parabola,
hence the covariance curves back down around $\langle n\rangle=1/2$. The slope
of the function depends on the transmission (loss). In the next section we use
this method as one way to estimate the purity of the detected visible light,
but due to the significantly inhomogeneous responsivity of the detectors, this
is not possible in the infrared. We can, however, compare to theory
(“Predicted” in Fig. 2e): a simple model (see Appendix C) can approximately
predict the output, based on the spectrum in Fig. 2c and the detector QE,
assuming biphotons centered around 1550 nm. It predicts a similar scaling to
the one measured. As this coincidence measurement is quite sensitive to
wavelength-dependent loss, the prediction estimates a decoherence amounting to
98% effective loss while we measured 96%, hence the factor of two between
measurement and prediction. This deviation is due to inaccuracies in the QE
and spectrum measurement. We emphasize that this decoherence is due to the
limited spectral range and QE of the detector, and not due to the purity of
the squeezed states.

These three experiments verify our predictions that the DOPA is a source of
highly multimode squeezed light, spanning hundreds of frequency modes.

### DOPA Waveguide

The DOPA waveguide consists of 4 cm MgO-doped LN ridge waveguide (Covesion
WG-1550-40 WGCK40), poled for second harmonic generation of 1550 to 775 nm and
designed to be single-mode at near-IR wavelengths. The LN device is mounted in
an oven (Covesion PV40). The oven is mounted with a V-shaped mount (Thorlabs
VC3C) and a 3D-printed sleeve to ensure horizontal orientation. The oven and
objective lenses (Newport 5726-B-H, 5726-C-H) are mounted on linear axes
stages (Newport 562-XYZ, 561D-XYZ-LH, 561D-XYZ respectively). Prior to the
waveguide input, a dichroic mirror is used to merge the 775/1550 nm pump and
alignment beams (Layertec 109068). A custom mirror (Layertec) is used
following the output to separate the pump and signal. The output pump beam is
imaged onto a camera (Basler acA800-510um) to view and optimize the pump’s
spatial profile in the waveguide.

### AFC

The crystal used for AFC is a 3 cm poled KTP crystal (Raicol Crystals). The
poling profile is designed such that the spatial frequency (inverse of the
poling period) varies linearly over the length of the crystal, except at the
facets of the crystal. The poling spatial frequency is designed to phase match
the conversion of all the target wavelengths. The front and back 0.25 mm of
the crystal are poled with a rapidly varying tanh profile. This improves the
rapid passage of the wavelengths that convert near the ends of the crystal.
See Appendix G for more details. The crystal is designed for use at 48°C and
held in an oven (Eksma Optics HP30), mounted on a 7-axis mount (Thorlabs K6XS,
SM1P1).

To optimize the conversion efficiency, bright CW laser light is sent through
the DOPA: it thus has same spatial mode as the squeezed light. This alignment
beam is used to optimize the spatial overlap of the pump with the signal. The
conversion efficiency of the squeezed light is then optimized by maximizing
the temporal overlap. This consists of a sweep of the time delay, followed by
a sweep of the pump pulse chirp (i.e. duration). As AFC is phase-insensitive,
there is no need for phase stabilization. The delay line, a motorized stage
with a retro-reflecting right-angle mirror, adjusts the timing of 775 nm pump
prior to the DOPA waveguide. The chirp is programmed by adjusting the spectral
phase profile of the pulse shaper. See the spectrometer section below for more
details.

To mitigate spatio-temporal effects in the AFC, we designed the pump to have a
beam waist of 300 µm, and the signal beam to have a waist of approximately 70
µm when focused in the crystal.

To measure the phase-averaged response of the AFC process, monochromated
super-continuum (NKT Origami, Thorlabs HN1550, and the infrared monochromator
described below) is sent through the AFC, one wavelength at a time, and the
resulting visible spectrum measured on the EMCCD camera spectrometer
(described below).

### 775 nm generation from 1033 nm

Pulsed 775 nm light is generated by cascading a second harmonic generation
(SHG) and OPA process. Both processes use lithium triborate (LBO; Newlight
Photonics). Both crystals are temperature tuned and non-critically phase
matched (the optical polarization are aligned to the crystal axes, to avoid
spatial walk-off). The lengths are designed to match the temporal walk-off of
the 200 fs pulses. The first crystal, 6 mm long, held at 182.6°C, generates
517 nm ($y$-axis polarized) from 1033 nm ($z$-axis). The second, 10 mm long,
held at 145.4°C, carries out the 517 nm ($y$) $-$ 1548.5 nm ($z$)
$\rightarrow$ 775 nm ($z$) process. The 1548.5 nm seed is generated by a diode
laser (ILX Lightwave 79800C), amplitude-modulated by a stabilized (Oz Optics
MBC-SUPER-PD-3A) electro-optic modulator (EOM) (Eospace AZ-DS5-10-PFA-PFA-LV-
LR), to generate 10 ns square pulses. These seed pulses are optically
amplified (Pritel PMFA-20), following an isolator. The EOM pulsing is
triggered by the Amplitude Satsuma 1033 nm laser, with an appropriate time
delay (IDQuantique ID900). The EOM is driven by a pulse/function generator
(HP-8116A).

### Parametric gain measurements

For this measurement, the DOPA is pumped with a 775 nm, 80 MHz Ti:S laser
(Spectra-Physics Tsunami). (All other experiments use the 775 nm pump
generated from a 1033 nm laser.) The high repetition rate ensures a
sufficiently high average power of squeezed light. The 1550 nm output beam
power is measured with a pW-sensitive power meter (Thorlabs S150C), after
coupling into multimode fiber (Thorlabs FG200LEA), which is about 95%
efficient. The DOPA pump power is removed with a dichroic mirror (Layertec)
and measured concurrently (Thorlabs S130VC/S130C). The squeezed light is
filtered further with long-pass filters (Thorlabs FELH1150).

At visible wavelengths, the 620 nm power is measured on the camera by
integrating the signal on all illuminated pixels, after appropriate filtering
of the pump (Thorlabs DMLP900; Semrock FF01-632/148).

### Single photon spectrometers

The infrared wavelength spectrometer is based on a diffraction grating on
motorized rotation stage (Thorlabs GR25-0616; K10CR1). The first order
reflection couples into SMF-28 single mode fiber (Thorlabs F260APC-1550),
which monochromates the input. The fiber coupling efficiency is approximately
35%. For single photon detection, this coupled into a InGaAs SPAD (IDQuantique
ID Qube NIR Gated), set to the nominal 15% QE. The wavelength-angle
correspondence is calibrated using tunable lasers (JDSU mTLG-C1C1L1) between
1527 nm and 1609 nm and fitting to the grating equation,
$m\lambda=d(\sin(\theta_{i0}+\Delta)-\sin(\theta_{ir}-\Delta)=2d\cos[(\theta_{i0}+\theta_{r0})/2]\sin[(\theta_{i0}-\theta_{r0})/2+\Delta]$,
where $d$ is the grating constant, $\theta_{i0},\theta_{r0}$ are some
reference incidence and reflection angles, and $\Delta$ is the rotation angle
of the grating. The extrapolation to uncalibrated wavelengths is deemed
correct as the spectrum stops sharply at 1150 nm, which matches the long pass
filter cutoff (Thorlabs FELH1150). However, the coupling efficiency as a
function of wavelength uncalibrated.

The visible light spectrometer is based on a diffraction grating (Ibsen
PCG-1908/675-972) imaged by an objective lens (Olympus UPLFLN4x) onto the NüVü
HNü 512 IS EMCCD camera. See Appendix H for more details. The wavelength-pixel
correspondence is calibrated by monochromating supercontinuum (NKT Origami,
Thorlabs HN1550) with the infrared monochromator, and converting this light
through AFC. The AFC pump is amplitude-modulated to a narrow bandwidth, thus
also effectively monochromated. The monochromated supercontinuum wavelength-
to-angle is calibrated with an optical spectrum analyzer (Ando AQ6317B). An EM
gain of 3000 is used on the camera.

The SPAD and camera are triggered by the Amplitude Satsuma laser, with an
appropriate time delay (IDQuantique ID900).

### Coincidence detection

In the infrared, the beam out of the DOPA is split with a broadband
beamsplitter and achromatic waveplate (Thorlabs UFBS50502, AHWP10M-1600). The
two beam paths are coupled into two SPADs (IDQuantique ID Qube NIR Gated)
using fixed collimators (Thorlabs F260APC-1550). The SPADs are set to the
nominal 15% QE. The fiber coupling efficiency is approximately 85%. Data
acquisition and triggering is configured with an IDQuantique ID900 time
controller. This includes feedback for enforcing a global dead time after
either of the SPADs fired, and reducing the triggering rate to 2 MHz (maximum
SPAD response rate) when using the 80 MHz Spectra-Physics Tsunami Ti:S.

The coincidence detection of visible light is performed with the NüVü EMCCD
camera. The diffraction grating is replaced with a Wollaston prism (Thorlabs
WP10-A) rotated to achieve 50:50 in the two output beams, at around 45°. The
beams are focused with a 30 mm lens (Thorlabs AC254-030-A-ML). The camera is
configured to detect in two separate regions of interest (opposite corners of
the CCD). The signal from these two 16$\times$16 pixel regions are integrated
to obtain the overall power. An EM gain of 500 is used on the camera, which is
used in analog mode.

The SPADs and camera are triggered by the Amplitude Satsuma 1033 nm laser,
with an appropriate time delay (IDQuantique ID900).

Refer to Appendix C for a theoretical motivation for these experiment.

### Temporal pulse shaper

Light diffracted by a transmission grating (Ibsen Photonics PCG-1765-808-981)
is focused by a 150 mm cylindrical lens (Thorlabs LJ1629L2-B) onto an spatial
light modulator (SLM; Meadowlark P1920-0600-1300-PCIe). A vertically-oriented
blazed grating is written on the SLM (orthogonal to the direction that
frequency components are dispersed). The SLM is oriented such that the light
reflected in the 1st order travels back through the cylindrical lens and
grating, following the same path except for a slight downward angle, such that
backward-traveling beam is separated by a pick-off mirror. The spatial phases
(vertical translations) of the blazed grating along the wavelength axis impart
spectral phases to the pulses. Amplitude modulation is also possible by
reducing the grating amplitude at a given wavelength, thereby reducing the
diffraction efficiency.

The pixel-column-to-wavelength calibration is performed by writing narrow
Gaussian-shaped amplitude modulation patterns along the wavelength axis of the
SLM. The beam is then measured with an optical spectrum analyzer (Ando
AQ6317B) to determine the transmitted wavelength.

The distance between the three optics is optimized by rounding the ellipticity
in the beam shape, as viewed on a camera, and minimizing the pulse duration,
measured with an autocorrelator (APE Pulsecheck). The orientations of the SLM,
grating, and input mirrors are tuned to eliminate the spatial chirp after the
initial alignment. The spatial chirp as a function of the transverse
directions of the beam is measured by translating the end of a multimode fiber
connected to an optical spectrum analyzer (Ando AQ6317B).

The pulse shaper is about 60% efficient without any phase modulations. Phase
modulations can reduce the efficiency if the features are smaller than the
imaging resolution of one wavelength onto the SLM.

This pulse shaper is based on the design in Ref. [20].

### EMCCD camera operation

The NüVü HNü EMCCD camera is water-cooled by a thermo-electric chiller (Solid
State TCube Edge), and the camera operates with a CCD temperature of -60°C,
cooled by the built-in TEC. The camera is mounted on one rotation and two
linear stages (Newport RS65; Thorlabs PT1A). One linear stage is used to
center the camera with respect to the beam. The other two are actuated by
motorized micrometers (Thorlabs Z812B; Newport TRA25CC), in order to align the
camera to the objective’s focal plane, such that the entire spectrum is in
focus when used as a spectrometer. These two degrees of freedom are optimized
by minimizing the point spread function of two different wavelengths,
alternating until convergence.

The laser repetition rate of the Satsuma is chosen to match the maximum frame-
rate of the camera (under a given experimental configuration), and the camera
is triggered by the laser (IDQuantique ID900 is used to adjust the delay). The
blanking and exposure times are generally set to 0.1 ms, as these values do
not adversely affect the frame-rates within the experimental configurations
reported.

The camera digital readout (“pixel value”, $p$) is converted to photo-
electrons, $\langle n_{e}\rangle$, by subtracting the bias $b$ and dividing by
the total gain. The latter is comprised of the electron multiplication (EM)
gain and the analog-to-digital conversion factor $k$ (here 21.43
photoelectrons per pixel unit). Hence $\langle n_{e}\rangle=kp/g-b$. Bias
subtraction is calibrated per pixel where possible, as these exhibited small
variations. See Appendix B for an explanation of the measurements performed
with an EMCCD.

See Appendix F for photographs.

Figure M1:  Simplified experimental diagram. EDFA: erbium doped fiber
amplifier. WG: waveguide. NLF: nonlinear fiber. CW: continuous wave.

## References

  * [1] M. Allgaier, V. Ansari, L. Sansoni, C. Eigner, V. Quiring, R. Ricken, G. Harder, B. Brecht, and C. Silberhorn (2017-01) Highly efficient frequency conversion with bandwidth compression of quantum light.  Nature Communications 8 (1).  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [2] J. R. Basani, M. Heuck, D. R. Englund, and S. Krastanov (2022) All-photonic artificial neural network processor via non-linear optics.  External Links: 2205.08608 Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [3] C. Baune, J. Gniesmer, A. Schönbeck, C. E. Vollmer, J. Fiurášek, and R. Schnabel (2015-06) Strongly squeezed states at 532 nm based on frequency up-conversion.  Opt. Express 23 (12),  pp. 16035–16041.  External Links: Link, Document Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [4] E. Bolduc, D. Faccio, and J. Leach (2017-04) Acquisition of multiple photon pairs with an EMCCD camera.  Journal of Optics 19 (5),  pp. 054006.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [5] K. A. G. Bonsma-Fisher, P. J. Bustard, C. Parry, T. A. Wright, D. G. England, B. J. Sussman, and P. J. Mosley (2022-11) Ultratunable quantum frequency conversion in photonic crystal fiber.  Phys. Rev. Lett. 129,  pp. 203603.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [6] J. E. Bourassa, R. N. Alexander, M. Vasmer, A. Patil, I. Tzitrin, T. Matsuura, D. Su, B. Q. Baragiola, S. Guha, G. Dauphinais, K. K. Sabapathy, N. C. Menicucci, and I. Dhand (2021-02) Blueprint for a scalable photonic fault-tolerant quantum computer.  Quantum 5,  pp. 392.  External Links: ISSN 2521-327X, Document Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [7] M. Cabrejo-Ponce, A. L. M. Muniz, M. Huber, and F. Steinlechner (2023-07) High-dimensional entanglement for quantum communication in the frequency domain.  Laser & Photonics Reviews,  pp. 2201010.  External Links: Document, Link Cited by: Discussion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [8] Y. Cai, J. Roslund, G. Ferrini, F. Arzani, X. Xu, C. Fabre, and N. Treps (2017-06) Multimode entanglement in reconfigurable graph states using optical frequency combs.  Nature Communications 8 (1),  pp. 15645.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [9] K. Chang, X. Cheng, M. C. Sarihan, A. K. Vinod, Y. S. Lee, T. Zhong, Y. Gong, Z. Xie, J. H. Shapiro, F. N. C. Wong, and C. W. Wong (2021-03) 648 Hilbert-space dimensionality in a biphoton frequency comb: entanglement of formation and Schmidt mode decomposition.  npj Quantum Information 7 (1),  pp. 48.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [10] M. Chen, N. C. Menicucci, and O. Pfister (2014-03) Experimental realization of multipartite entanglement of 60 modes of a quantum optical frequency comb.  Phys. Rev. Lett. 112,  pp. 120505.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [11] H. Defienne, B. Ndagano, A. Lyons, and D. Faccio (2021-02) Polarization entanglement-enabled quantum holography.  Nature Physics 17 (5),  pp. 591–597.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [12] H. Defienne, J. Zhao, E. Charbon, and D. Faccio (2021-04) Full-field quantum imaging with a single-photon avalanche diode camera.  Phys. Rev. A 103,  pp. 042608.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [13] Y. Deng, Y. Gu, H. Liu, S. Gong, H. Su, Z. Zhang, H. Tang, M. Jia, J. Xu, M. Chen, J. Qin, L. Peng, J. Yan, Y. Hu, J. Huang, H. Li, Y. Li, Y. Chen, X. Jiang, L. Gan, G. Yang, L. You, L. Li, H. Zhong, H. Wang, N. Liu, J. J. Renema, C. Lu, and J. Pan (2023-10) Gaussian boson sampling with pseudo-photon-number-resolving detectors and quantum computational advantage.  Phys. Rev. Lett. 131,  pp. 150601.  External Links: Document, Link Cited by: Multimode quantum state sampling at visible wavelengths, footnote 1, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [14] K. Dhimitri, S. M. Fullerton, B. Coyle, K. E. Bennett, T. Miura, T. Higuchi, and T. Maruno (2022) Scientific CMOS (sCMOS) camera capabilities with a focus on quantum applications.  In Photonics for Quantum 2022, D. F. Figer (Ed.),  PC122430L, Vol. PC12243, Rochester, NY, United States.  External Links: Document, Link Cited by: Discussion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [15] C. Drago and A. M. Brańczyk (2022-10) Tunable frequency-bin multimode squeezed vacuum states of light.  Phys. Rev. A 106,  pp. 043714.  External Links: Document, Link Cited by: Discussion.
  * [16] A. Eldan, O. Gilon, A. Lagimi, E. Forman, and A. Pe’er (2023) Multiplexed processing of quantum information across an ultra-wide optical bandwidth.  External Links: 2310.17819 Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [17] I. Esmaeil Zadeh, J. Chang, J. W. N. Los, S. Gyger, A. W. Elshaari, S. Steinhauer, S. N. Dorenbos, and V. Zwiller (2021-05) Superconducting nanowire single-photon detectors: a perspective on evolution, state-of-the-art, future developments, and applications.  Applied Physics Letters 118 (19),  pp. 190502.  External Links: ISSN 0003-6951, Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [18] M. J. Fitch, B. C. Jacobs, T. B. Pittman, and J. D. Franson (2003-10) Photon-number resolution using time-multiplexed single-photon detectors.  Phys. Rev. A 68,  pp. 043814.  External Links: Document, Link Cited by: Multimode quantum state sampling at visible wavelengths.
  * [19] P. Folge, M. Santandrea, M. Stefszky, B. Brecht, and C. Silberhorn (2022) Multimode squeezed states in frequency conversion based time-frequency networks.  In Frontiers in Optics + Laser Science 2022 (FIO, LS),  FM5B.5, Rochester, NY, United States.  External Links: Link, Document Cited by: Discussion.
  * [20] E. Frumker and Y. Silberberg (2007-12) Phase and amplitude pulse shaping with two-dimensional phase-only spatial light modulators.  J. Opt. Soc. Am. B 24 (12),  pp. 2940–2947.  External Links: Document Cited by: Temporal pulse shaper.
  * [21] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex (2017-10) Gaussian boson sampling.  Phys. Rev. Lett. 119,  pp. 170501.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [22] T. Hiemstra, T.F. Parker, P. Humphreys, J. Tiedau, M. Beck, M. Karpiński, B.J. Smith, A. Eckstein, W.S. Kolthammer, and I.A. Walmsley (2020-07) Pure single photons from scalable frequency multiplexing.  Phys. Rev. Appl. 14,  pp. 014052.  External Links: Document, Link Cited by: Discussion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [23] A. Hosaka, T. Kawamori, and F. Kannari (2016-11) Multimode quantum theory of nonlinear propagation in optical fibers.  Phys. Rev. A 94,  pp. 053833.  External Links: Document, Link Cited by: Highly multimode squeezing in the frequency domain, Highly multimode squeezing in the frequency domain: details.
  * [24] I. Hurvitz, A. Karnieli, and A. Arie (2023-06) Frequency-domain engineering of bright squeezed vacuum for continuous-variable quantum information.  Optics Express 31 (12),  pp. 20387.  External Links: Document, Link Cited by: Discussion.
  * [25] M. Jahanbozorgi, Z. Yang, S. Sun, H. Chen, R. Liu, B. Wang, and X. Yi (2023-08) Generation of squeezed quantum microcombs with silicon nitride integrated photonic circuits.  Optica 10 (8),  pp. 1100–1101.  External Links: Link, Document Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [26] U. A. Javid, R. Lopez-Rios, J. Ling, A. Graf, J. Staffa, and Q. Lin (2023-06) Chip-scale simulations in a quantum-correlated synthetic space.  Nature Photonics 17 (10),  pp. 883–890.  External Links: Document, ISSN 1749-4893 Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [27] C. Joshi, A. Farsi, S. Clemmen, S. Ramelow, and A. L. Gaeta (2018-02) Frequency multiplexing for quasi-deterministic heralded single-photon sources.  Nature Communications 9 (1).  External Links: Document, Link Cited by: Discussion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [28] C. Joshi, A. Farsi, A. Dutt, B. Y. Kim, X. Ji, Y. Zhao, A. M. Bishop, M. Lipson, and A. L. Gaeta (2020-04) Frequency-domain quantum interference with correlated photons from an integrated microresonator.  Phys. Rev. Lett. 124,  pp. 143601.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [29] C. Joshi, A. Farsi, and A. Gaeta (2017) Frequency-domain boson sampling.  In Conference on Lasers and Electro-Optics,  FTu1F.1, San Jose, CA, United States.  External Links: Link, Document Cited by: Discussion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [30] A. Kumar and A. M. Marino (2019-12) Spatial squeezing in bright twin beams generated with four-wave mixing: constraints on characterization with an electron-multiplying charge-coupled-device camera.  Phys. Rev. A 100,  pp. 063828.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [31] P. Kumar (1990-12) Quantum frequency conversion.  Opt. Lett. 15 (24),  pp. 1476–1478.  External Links: Link, Document Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [32] E. Lantz, J. Blanchet, L. Furfaro, and F. Devaux (2008-04) Multi-imaging and Bayesian estimation for photon counting with EMCCDs.  Monthly Notices of the Royal Astronomical Society 386 (4),  pp. 2262–2270.  External Links: ISSN 0035-8711, Document Cited by: Multimode quantum state sampling at visible wavelengths.
  * [33] L. Ledezma, R. Sekine, Q. Guo, R. Nehra, S. Jahani, and A. Marandi (2022-03) Intense optical parametric amplification in dispersion-engineered nanophotonic lithium niobate waveguides.  Optica 9 (3),  pp. 303–308.  External Links: Link, Document Cited by: Discussion.
  * [34] F. Li, T. Li, M. O. Scully, and G. S. Agarwal (2021-04) Quantum advantage with seeded squeezed light for absorption measurement.  Phys. Rev. Appl. 15,  pp. 044030.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [35] N. Liu, Y. Liu, J. Li, L. Yang, and X. Li (2016-01) Generation of multi-mode squeezed vacuum using pulse pumped fiber optical parametric amplifiers.  Optics Express 24 (3),  pp. 2125–2133.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [36] H. Lu, M. Liscidini, A. L. Gaeta, A. M. Weiner, and J. M. Lukens (2023-12) Frequency-bin photonic quantum information.  Optica 10 (12),  pp. 1655–1671.  External Links: Document Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [37] H. Lu, J. M. Lukens, N. A. Peters, O. D. Odele, D. E. Leaird, A. M. Weiner, and P. Lougovski (2018-01) Electro-optic frequency beam splitters and tritters for high-fidelity photonic quantum information processing.  Phys. Rev. Lett. 120,  pp. 030502.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [38] H. Lu, E. M. Simmerman, P. Lougovski, A. M. Weiner, and J. M. Lukens (2020-09) Fully arbitrary control of frequency-bin qubits.  Phys. Rev. Lett. 125,  pp. 120503.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [39] J. M. Lukens and P. Lougovski (2017-01) Frequency-encoded photonic qubits for scalable quantum information processing.  Optica 4 (1),  pp. 8–16.  External Links: Link, Document Cited by: Discussion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [40] A. I. Lvovsky (2015-01) Squeezed light.  Photonics: Scientific Foundations, Technology and Applications 1,  pp. 121–163.  External Links: ISBN 9781119009719, Document Cited by: Highly multimode squeezing in the frequency domain.
  * [41] J. Ma, S. Masoodian, D. A. Starkey, and E. R. Fossum (2017-12) Photon-number-resolving megapixel image sensor at room temperature without avalanche gain.  Optica 4 (12),  pp. 1474–1481.  External Links: Link, Document Cited by: Discussion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [42] L. S. Madsen, F. Laudenbach, M. Falamarzi. Askarani, F. Rortais, T. Vincent, J. F. F. Bulmer, F. M. Miatto, L. Neuhaus, L. G. Helt, M. J. Collins, A. E. Lita, T. Gerrits, S. W. Nam, V. D. Vaidya, M. Menotti, I. Dhand, Z. Vernon, N. Quesada, and J. Lavoie (2022-06) Quantum computational advantage with a programmable photonic processor.  Nature 606 (7912),  pp. 75–81.  External Links: Document, Link Cited by: footnote 1, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [43] N. C. Menicucci, P. van Loock, M. Gu, C. Weedbrook, T. C. Ralph, and M. A. Nielsen (2006-09) Universal quantum computation with continuous-variable cluster states.  Phys. Rev. Lett. 97,  pp. 110501.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [44] C. L. Morrison, F. Graffitti, P. Barrow, A. Pickston, J. Ho, and A. Fedrizzi (2022-06) Frequency-bin entanglement from domain-engineered down-conversion.  APL Photonics 7 (6),  pp. 066102.  External Links: ISSN 2378-0967, Document, Link Cited by: Discussion.
  * [45] J. Moses, H. Suchowski, and F. X. Kärtner (2012-05) Fully efficient adiabatic frequency conversion of broadband Ti:sapphire oscillator pulses.  Opt. Lett. 37 (9),  pp. 1589–1591.  External Links: Link, Document Cited by: Adiabatic frequency conversion (AFC) of squeezed light, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [46] C. M. Natarajan, M. G. Tanner, and R. H. Hadfield (2012-04) Superconducting nanowire single-photon detectors: physics and applications.  Superconductor Science and Technology 25 (6),  pp. 063001.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [47] R. Nehra, R. Sekine, L. Ledezma, Q. Guo, R. M. Gray, A. Roy, and A. Marandi (2022) Few-cycle vacuum squeezing in nanophotonics.  Science 377 (6612),  pp. 1333–1337.  External Links: Document, Link Cited by: Discussion.
  * [48] Nuvu Cameras HNü 512 – 512 x 512 EMCCD.  Note: https://www.nuvucameras.com/products/hnu-512/Accessed: 2024-01-07. External Links: Link Cited by: Multimode quantum state sampling at visible wavelengths.
  * [49] A. Omi, A. Hosaka, M. Tomita, Y. Yamagishi, K. Wakui, S. Niimura, K. Takahashi, M. Takeoka, and F. Kannari (2021-07) Independently programmable frequency-multiplexed phase-sensitive optical parametric amplification in the optical telecommunication band.  Opt. Express 29 (14),  pp. 21683–21697.  External Links: Link, Document Cited by: Discussion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [50] H. Paul, P. Törmä, T. Kiss, and I. Jex (1996-04) Photon chopping: new way to measure the quantum state of light.  Phys. Rev. Lett. 76,  pp. 2464–2467.  External Links: Document, Link Cited by: Multimode quantum state sampling at visible wavelengths.
  * [51] M. Pysher, Y. Miwa, R. Shahrokhshahi, R. Bloomer, and O. Pfister (2011-07) Parallel generation of quadripartite cluster entanglement in the optical frequency comb.  Phys. Rev. Lett. 107,  pp. 030505.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [52] S. Rahimi-Keshari, A. P. Lund, and T. C. Ralph (2015-02) What can quantum optics say about computational complexity theory?.  Phys. Rev. Lett. 114,  pp. 060501.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [53] C. Reimer, M. Kues, P. Roztocki, B. Wetzel, F. Grazioso, B. E. Little, S. T. Chu, T. Johnston, Y. Bromberg, L. Caspani, D. J. Moss, and R. Morandotti (2016) Generation of multiphoton entangled quantum states by means of integrated frequency combs.  Science 351 (6278),  pp. 1176–1180.  External Links: Document, Link Cited by: Discussion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [54] V. Roman-Rodriguez, D. Fainsin, G. L. Zanin, N. Treps, E. Diamanti, and V. Parigi (2023) Spectrally multimode squeezed states generation at telecom wavelengths.  External Links: 2306.07267 Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [55] J. Roslund, R. M. de Araújo, S. Jiang, C. Fabre, and N. Treps (2014-12) Wavelength-multiplexed quantum networks with ultrafast frequency combs.  Nature Photonics 8 (2),  pp. 109–112.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [56] E. Rozenberg, A. Karnieli, O. Yesharim, J. Foley-Comer, S. Trajtenberg-Mills, D. Freedman, A. M. Bronstein, and A. Arie (2022-06) Inverse design of spontaneous parametric downconversion for generation of high-dimensional qudits.  Optica 9 (6),  pp. 602–615.  External Links: Link, Document Cited by: Discussion.
  * [57] A. Samblowski, C. E. Vollmer, C. Baune, J. Fiurášek, and R. Schnabel (2014-05) Weak-signal conversion from 1550 to 532 nm with 84% efficiency.  Opt. Lett. 39 (10),  pp. 2979–2981.  External Links: Link, Document Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [58] J. Sperling, W. Vogel, and G. S. Agarwal (2012-02) True photocounting statistics of multiple on-off detectors.  Phys. Rev. A 85,  pp. 023820.  External Links: Document, Link Cited by: Multimode quantum state sampling at visible wavelengths.
  * [59] H. Suchowski, B. D. Bruner, Y. Israel, A. Ganany-Padowicz, A. Arie, and Y. Silberberg (2016) Broadband photon pair generation at 3$\omega$/2.  Applied Physics B 122,  pp. 1–5.  External Links: Document Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [60] H. Suchowski, V. Prabhudesai, D. Oron, A. Arie, and Y. Silberberg (2009-07) Robust adiabatic sum frequency conversion.  Opt. Express 17 (15),  pp. 12731–12740.  External Links: Link, Document Cited by: Adiabatic frequency conversion (AFC) of squeezed light, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [61] P. Svihra, Y. Zhang, P. Hockett, S. Ferrante, B. Sussman, D. England, and A. Nomerotski (2020-07) Multivariate discrimination in quantum target detection.  Applied Physics Letters 117 (4),  pp. 044001.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [62] S. Tan, B. I. Erkmen, V. Giovannetti, S. Guha, S. Lloyd, L. Maccone, S. Pirandola, and J. H. Shapiro (2008-12) Quantum illumination with Gaussian states.  Phys. Rev. Lett. 101,  pp. 253601.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [63] V. D. Vaidya, B. Morrison, L. G. Helt, R. Shahrokshahi, D. H. Mahler, M. J. Collins, K. Tan, J. Lavoie, A. Repingon, M. Menotti, N. Quesada, R. C. Pooser, A. E. Lita, T. Gerrits, S. W. Nam, and Z. Vernon (2020) Broadband quadrature-squeezed vacuum and nonclassical photon number correlations from a nanophotonic device.  Science Advances 6 (39),  pp. eaba9186.  External Links: Document, Link Cited by: Highly multimode squeezing in the frequency domain, Highly multimode squeezing in the frequency domain: details.
  * [64] C. E. Vollmer, C. Baune, A. Samblowski, T. Eberle, V. Händchen, J. Fiurášek, and R. Schnabel (2014-02) Quantum up-conversion of squeezed vacuum states from 1550 to 532 nm.  Phys. Rev. Lett. 112,  pp. 073602.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [65] X. Wang, X. Jiao, B. Wang, Y. Liu, X. Xie, M. Zheng, Q. Zhang, and J. Pan (2023-04) Quantum frequency conversion and single-photon detection with lithium niobate nanophotonic chips.  npj Quantum Information 9 (1).  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [66] W. Wasilewski, A. I. Lvovsky, K. Banaszek, and C. Radzewicz (2006-06) Pulsed squeezed light: simultaneous squeezing of multiple modes.  Phys. Rev. A 73,  pp. 063819.  External Links: Document, Link Cited by: Highly multimode squeezing in the frequency domain, Highly multimode squeezing in the frequency domain: details.
  * [67] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd (2012-05) Gaussian quantum information.  Rev. Mod. Phys. 84,  pp. 621–669.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [68] Z. Xie, T. Zhong, S. Shrestha, X. Xu, J. Liang, Y. Gong, J. C. Bienfang, A. Restelli, J. H. Shapiro, F. N. C. Wong, and C. W. Wong (2015-06) Harnessing high-dimensional hyperentanglement through a biphoton frequency comb.  Nature Photonics 9 (8),  pp. 536–542.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [69] Y. Yamagishi, A. Hosaka, K. Tanji, S. Kurimura, and F. Kannari (2021) Arbitrary mixing of frequency-range multimode quantum states using nonlinear waveguide crystals based on dispersion engineering.  In OSA Nonlinear Optics 2021,  NM2B.4, Washington, DC, United States.  External Links: Link, Document Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [70] Z. Yang, M. Jahanbozorgi, D. Jeong, S. Sun, O. Pfister, H. Lee, and X. Yi (2021-08) A squeezed quantum microcomb on a chip.  Nature Communications 12 (1),  pp. 4781.  External Links: Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [71] L. Yuan, A. Dutt, and S. Fan (2021-07) Synthetic frequency dimensions in dynamically modulated ring resonators.  APL Photonics 6 (7),  pp. 071102.  External Links: ISSN 2378-0967, Document, Link Cited by: Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [72] H. Zhong, Y. Deng, J. Qin, H. Wang, M. Chen, L. Peng, Y. Luo, D. Wu, S. Gong, H. Su, Y. Hu, P. Hu, X. Yang, W. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N. Liu, J. J. Renema, C. Lu, and J. Pan (2021-10) Phase-programmable Gaussian boson sampling using stimulated squeezed light.  Phys. Rev. Lett. 127,  pp. 180502.  External Links: Document, Link Cited by: footnote 1, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.
  * [73] H. Zhong, H. Wang, Y. Deng, M. Chen, L. Peng, Y. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, P. Hu, X. Yang, W. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N. Liu, C. Lu, and J. Pan (2020) Quantum computational advantage using photons.  Science 370 (6523),  pp. 1460–1463.  External Links: Document, Link Cited by: footnote 1, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion, Highly multimode visible squeezed light with programmable spectral correlations through broadband up-conversion.

[sections]

## Appendix Contents

[sections]0

##  Appendix A Photon statistics in multimode Gaussian quantum optics

The aim of this section is to provide a working knowledge of the physics of
multimode Gaussian quantum optics, especially as it relates to our
experiments. We provide an explanation of the underlying theory in our
preferred formalism, show how to apply the nonlinear optics equations, and how
to predict the supermodes and squeezing values. Finally, we describe the
photon number properties of Gaussian states.

We refer to as “Gaussian” the class of states that are fully described by the
first two moments of its field operators, or, equivalently, whose phase space
distribution is a Gaussian. These states are the most typical in experimental
quantum optics (this is a consequence of weak optical nonlinearities and the
large – classical – driving fields that are required to achieve them).
Furthermore, while multimode or multivariable quantum-mechanical states have
generally exponentially large representations, Gaussian states have efficient
phase-space representations. This makes working with large multimode Gaussian
states tractable.

### The bosonic covariance matrix

To fully characterize a zero-mean Gaussian state, we must know the squeezing,
the thermal noise, the loss and the correlations between all modes; all this
information is encoded in the covariance matrix. In typical convention
[Olivares, FerraroOlivaresParis]:

$\sigma=\frac{1}{2}\left\langle\left\\{\hat{\xi},\hat{\xi}^{\dagger}\right\\}\right\rangle-\left\langle\hat{\xi}\right\rangle\left\langle\hat{\xi}^{\dagger}\right\rangle,~{}\hat{\xi}^{\intercal}=\left[\hat{a}_{1},\ldots,\hat{a}_{M},\hat{a}_{1}^{\dagger},\ldots,\hat{a}_{M}^{\dagger}\right],$

where $\hat{a}_{i}^{\dagger}$ and $\hat{a}_{i}$ are the bosonic creation and
annihilation operators for some mode $i$, defined in the usual manner.
$\sigma$ is Hermitian and, for later convenience, we can write it in terms of
the submatrices:

$\sigma=\left(\begin{matrix}V+I_{M}/2&U\\\
U^{*}&V^{\intercal}+I_{M}/2\end{matrix}\right)$

such that $V$ is Hermitian and $U$ is symmetric; $I_{M}$ is the $M\times M$
identity. As we shall see,
$V\sim\langle\hat{a}_{i}^{\dagger}\hat{a}_{j}+h.c.\rangle$ contains the
information pertaining to the state’s classical properties and thermal
correlations, and
$U\sim\langle\hat{a}_{i}^{\dagger}\hat{a}_{j}^{\dagger}+h.c.\rangle$ encodes
the entanglement and higher-order-correlation physics.

Note that it is more common in the quantum optics and information literature
to see the quadrature covariance matrix, which use the real-valued canonical
variables $\hat{x}\propto\hat{a}+\hat{a}^{\dagger}$ and
$\hat{p}\propto\hat{a}-\hat{a}^{\dagger}$. The quadrature basis has the
advantage of being a real-valued symplectic space. However, the
$\hat{a}$-basis is a more natural convention for photon-number properties of
the field.

Overall, the Gaussian states defined by $\sigma$ have $2M^{2}$ free parameters
(real numbers), plus an additional $M$ local phase degree of freedom that may
be ignored, as they have no effect on the photon number statistics.

### Constructing and decomposing the covariance matrix

Here we summarize some of the body of work pertaining to the matrix
representations of multimode Gaussian systems. For more depth, see Refs. [28,
OpatrnyLeuchs, 15, Olivares, AdessoLee, 29, FerraroOlivaresParis].

In this formalism, the space of Gaussian states is closed under Gaussian
operations, some of which are represented as linear transformations. We call
the transformation a Green’s function if it represents the outcome of a
process described by a linear differential equation. In a discrete basis of
$\hat{a}$-operators, any lossless Gaussian operation can be represented as:

$\displaystyle\hat{a}_{\text{out}}$
$\displaystyle=C\hat{a}_{\text{in}}+S\hat{a}_{\text{in}}^{\dagger}$
$\displaystyle\hat{\xi}_{\text{out}}$
$\displaystyle=\left(\begin{matrix}C&S\\\
S^{*}&C^{*}\end{matrix}\right)\hat{\xi}_{\text{in}}=G\hat{\xi}_{\text{in}}$

or in a continuous basis parametrized by $\omega$:

$\displaystyle\hat{a}_{\text{out}}=\int
C(\omega,\omega^{\prime})\hat{a}_{\text{in}}(\omega^{\prime})+S(\omega,\omega^{\prime})\hat{a}_{\text{in}}(\omega^{\prime})^{\dagger}d\omega^{\prime}.$

$S$ may only be nonzero only if there is squeezing (if the Hamiltonian
contains an $\hat{a}^{\dagger}\hat{a}^{\dagger}+c.c.$ term).

In our experiment, specifically, we are interested in the Green’s function of
the OPA process, a combination of parametric amplification and dispersion
described by the differential equation:

$\frac{d\hat{a}}{dz}(\Delta\omega)=iD(\Delta\omega)\hat{a}+iA(z,\Delta\omega)*\hat{a}^{\dagger}$

where $D$ represents the dispersion (phase matching, group velocity
difference, and higher order dispersion), and the coupling term $A$ represents
the classical pump field, convolved with the quantum field. This equation is
the same as you would find in classical nonlinear optics, but with the
operator $\hat{a}$ replacing the classical field term. Classically, we could
interpret this as the expectation value of the field operator and its
evolution; quantum mechanically, we should think of this as the evolution of
$C,S$ that act on the operator. In our experiment, we expect $S$ to be
predominantly anti-diagonal due to energy conservation, and $C$ to depend on
the phase-matching function.

Such a Green’s function can be decomposed by the Bloch–Messiah decomposition,
which informs us of how much squeezing there is and over what modes. To take
advantage of standard computational methods, this is performed in the
quadrature basis of $\hat{x},\hat{p}$ (denote this Green’s function by
$G^{\prime}$). This matrix decomposition returns canonically conjugate
squeezed and anti-squeezed “supermodes.” Concretely: 222Note: this assumes the
convention $\hat{x}=\hat{a}^{\dagger}+\hat{a}$ and
$\hat{p}=i(\hat{a}^{\dagger}-\hat{a})$.

$\displaystyle G^{\prime}$ $\displaystyle=O_{\text{out}}\Sigma
O_{\text{in}}^{\intercal}$ $\displaystyle\Sigma$
$\displaystyle=\tfrac{1}{2}\text{diag}(\\{s_{i}\\}_{i=1}^{M},\\{s_{i}^{-1}\\}_{i=1}^{M})=\tfrac{1}{2}\text{diag}(\\{\mathrm{e}^{r_{i}}\\}_{i=1}^{M},\\{\mathrm{e}^{-r_{i}}\\}_{i=1}^{M})$

where $s_{i}$ are the symplectic eigenvalues, related to the squeezing
parameters $r_{i}$, and orthogonal matrices $O$ contain the “input” and
“output” symplectic eigenvectors or supermodes, which come in pairs (squeezed
and anti-squeezed). The input supermodes are not important if the initial
state is vacuum, as the squeezed vacuum only depends on the output supermodes.
Intuitively, Bloch–Messiah may be thought of as a singular value decomposition
that preserves commutation relations. These bases diagonalize the Green’s
function into single mode squeezing operations. Fig. A1 illustrates the
process. Each $s_{i}$ represents an independent squeezed mode, and a source of
photons with $\langle n\rangle=\sinh^{2}r_{i}$ distributed over the mode
$O_{\text{out},i}$, or from a quantum noise perspective, a $20\log_{10}s_{i}$
dB noise reduction in the mode $O_{\text{out},i}$. A change of basis from
$\hat{x},\hat{p}$ back to $\hat{a}$ transforms the orthogonal matrices into
unitary, and the diagonal matrix into a matrix where the quadrants are
diagonal:

$\sigma=\left(\begin{matrix}\text{diag}(\\{\cosh
r_{i}\\}_{i=1}^{M})&\text{diag}(\\{\sinh r_{i}\\}_{i=1}^{M})\\\
\text{diag}(\\{\sinh r_{i}\\}_{i=1}^{M})&\text{diag}(\\{\cosh
r_{i}\\}_{i=1}^{M})\end{matrix}\right).$

Figure A1:  Parametric amplification and the Bloch–Messiah decomposition. The
Bloch–Messiah decomposition reduces the parametric amplification to a phase-
dependent amplification of an orthogonal set of modes. Each input mode is
transformed to an output mode during this process: these modes are generally
not the equal due to other effects in the OPA that may alter the signal (e.g.
dispersion). a. In the case of vacuum input, squeezed vacuum is generated in
each output mode independently. b. In the case of bright input, a signal may
be decomposed into input modes, mapped to output modes, and each component is
(de)amplified depending on the phase.

With vacuum input, a covariance matrix can be generated as follows:

$\sigma_{ij}=\tfrac{1}{2}\\{G_{ik}\hat{\xi}_{k},G_{jk}\hat{\xi}_{k}\\}\Rightarrow\sigma=\tfrac{1}{2}GG^{\dagger}.$

Additional operations can be similarly applied to the covariance matrix:

$\sigma\rightarrow G\sigma G^{\dagger}.$

Note that vacuum $\sigma=I_{2M}/2$.

As previously mentioned, (going back to the $\hat{x},\hat{p}$ basis) the input
supermodes cancel out in the vacuum case:

$\tfrac{1}{2}G^{\prime}G^{\prime\intercal}=\tfrac{1}{2}O_{\text{out}}\Sigma
O_{\text{in}}^{\intercal}O_{\text{in}}\Sigma
O_{\text{out}}^{\intercal}=\tfrac{1}{2}O_{\text{out}}\Sigma^{2}O_{\text{out}}^{\intercal}$

and we see that the Bloch–Messiah decomposition is also valid on the the
covariance matrix, although one must account for squared diagonal matrix
(however, we note that the Williamson decomposition must be used when the
state is not pure: in the case of loss, thermal noise, or if the state has
been partially traced out).

We can observe that:

$\sigma=\tfrac{1}{2}GG^{\dagger}=\frac{1}{2}\left(\begin{matrix}SS^{\dagger}+I&CS+SC\\\
C^{\dagger}S^{\dagger}+S^{\dagger}C^{\dagger}&S^{\dagger}S+I\end{matrix}\right)=\left(\begin{matrix}V+I/2&U\\\
U^{*}&V^{\intercal}+I/2\end{matrix}\right).$

Note that for the $U$ quadrants to be nonzero, a squeezing operation must be
involved through a nonzero $S$. No unitary $C$ applied to any thermal state
can populate $U$. On the other hand, the squeezing contributes to “thermal”
components $V$, as it adds photons to the field, thus increasing the photon-
number variance.

The formalism introduced so far – where all the operations are unitary or
symplectic and preserve commutation relations – cannot account for losses or
inefficiencies that must be considered during the frequency conversion and
detection steps of our experiment. In both cases, the state occupies
unobserved modes, that are traced out. Tracing out modes in a covariance
matrix is simply equivalent to removing the corresponding rows and columns, in
other words, taking a principal submatrix that omits the traced out modes. For
example, this is used to derive the action of a lump loss or noise on a
covariance matrix (modeled as passing through a fictitious beamsplitter and
tracing out the second port):

$\sigma\rightarrow\sqrt{\eta\eta^{\intercal}}\circ\sigma+(1-\eta)\circ\nu\circ
I$

where $\eta$ is the transmission, $\nu={\bar{n}}+1/2$ represents any thermal
noise added ($1/2$ for vacuum). Both are vectors in the general case. The
$\circ$ operator is the element-wise or Hadamard product.

We now use these results to describe the physics relevant to our experiment:
we convert the infrared squeezed light to visible, but do not detect the
remaining infrared light. For notational simplicity let
$\hat{\xi}=[\hat{\accentset{\rightharpoonup}{a}}_{\text{vis}},\hat{\accentset{\rightharpoonup}{a}}_{\text{vis}}^{\dagger},\hat{\accentset{\rightharpoonup}{a}}_{\text{ir}},\hat{\accentset{\rightharpoonup}{a}}_{\text{ir}}^{\dagger}]$.
The initial state is:

$\sigma_{\text{tot}}(z=0)=\left(\begin{matrix}\sigma_{\text{vis}}&\sigma_{\text{vis,ir}}\\\
\sigma_{\text{vis,ir}}^{\dagger}&\sigma_{\text{ir}}\end{matrix}\right)=\left(\begin{matrix}I_{2M}/2&0\\\
0&\sigma_{\text{ir}}(0)\end{matrix}\right).$

The AFC (over crystal length $L$) acts as a unitary on the entire state:

$\displaystyle\sigma_{\text{tot}}(z=L)$
$\displaystyle=\left(\begin{matrix}G_{\text{vis,vis}}&G_{\text{vis,ir}}\\\
G_{\text{ir,vis}}&G_{\text{ir,ir}}\end{matrix}\right)\left(\begin{matrix}I_{2M}/2&0\\\
0&\sigma_{\text{ir}}\end{matrix}\right)\left(\begin{matrix}G_{\text{vis,vis}}&G_{\text{vis,ir}}\\\
G_{\text{ir,vis}}&G_{\text{ir,ir}}\end{matrix}\right)^{\dagger}.$

Expanding and tracing out the infrared modes we obtain the covariance matrix
of the observed visible modes:

$\sigma_{\text{vis}}(L)=G_{\text{vis,ir}}\sigma_{\text{ir}}(0)G_{\text{vis,ir}}^{\dagger}+G_{\text{vis,vis}}G_{\text{vis,vis}}^{\dagger}/2.$

As the conversion tends to unity, $G_{\text{vis,ir}}$ becomes unitary and
$G_{\text{vis,vis}}$ vanishes.

The sum frequency generation equations that yield the Green’s function in this
case are:

$\displaystyle\frac{d\hat{a}_{\text{vis}}}{dz}(\Delta\omega)$
$\displaystyle=iD_{\text{vis}}(z,\Delta\omega)\hat{a}_{\text{vis}}+iA(z,\Delta\omega)*\hat{a}_{\text{ir}}$
$\displaystyle\frac{d\hat{a}_{\text{ir}}}{dz}(\Delta\omega)$
$\displaystyle=iD_{\text{ir}}(z,\Delta\omega)\hat{a}_{\text{ir}}+iA^{*}(z,\Delta\omega)*\hat{a}_{\text{vis}}$

where, again, $A$ is the pump field and $D$ is the dispersion, which is
notably a function of $z$, due to aperiodic poling.

### Photon number properties of zero-mean Gaussian states

For convenience, in the context of photon number statistics, we can make the
covariance matrix symmetric and complex-valued, by defining the transformation
$X$, as first introduced in [13]:

$\begin{gathered}\sigma^{\prime}=\sigma-I_{2M}/2,\quad
X\sigma^{\prime}=\left(\begin{matrix}U&V^{\intercal}\\\
V&U^{*}\end{matrix}\right),\\\ X=\left(\begin{matrix}0&I_{M}\\\
I_{M}&0\end{matrix}\right).\end{gathered}$

$\sigma^{\prime}$ defined in this manner is useful, because we can use it to
find any of the following photon number expectation values:

$\left\langle\prod_{i}\hat{n}_{i}^{m_{i}}\right\rangle=\left\langle\prod_{i}\left(\hat{a}_{i}^{\dagger}\hat{a}_{i}\right)^{m_{i}}\right\rangle=\text{Haf}(X\sigma_{\accentset{\rightharpoonup}{m}}^{\prime}),~{}\accentset{\rightharpoonup}{m}\in\\{{0,1\\}}^{M}.$

using the conventional definition of the Hafnian (Haf).
$\sigma_{\accentset{\rightharpoonup}{m}}^{\prime}$ indicates the principal
submatrix with indices given by $\accentset{\rightharpoonup}{m}$. Suffice to
say, the Hafnian function is central to the (Gaussian or perturbative)
multimode physics of bosons, since it can be thought of as an implementation
of Wick’s theorem for Gaussian integrals; for zero-mean Gaussian variables
$x_{i_{j}}$:

$\langle x_{i_{1}}\ldots x_{i_{2m}}\rangle=\sum_{\mathcal{P}}\langle
x_{k_{1}}x_{k_{2}}\rangle\ldots\langle x_{k_{2m-1}}x_{k_{2m}}\rangle$

where the sum is over all pairings $\mathcal{P}$ – all possible ways to group
the $i$-indices into $m$ pairs of $k$-indices – hence the Hafnian, a function
defined for this very purpose [Caianiello]. And of course the operators
$\hat{a}_{i}$ are jointly Gaussian (in their quasi-probability distributions).

The general case with $\accentset{\rightharpoonup}{m}\in\mathbb{N}^{M}$ is
more complicated, as it requires more careful consideration of operator
ordering. Indeed, we were being hasty: $\hat{a}$ is not a random variable, it
is an operator. However, Wick’s theorem tells us that this was allowed as long
as the state and operator have compatible representations. For example, one
way to “treat our operators as though random variables” is by converting our
expression to Weyl ordering (denoted by $:\mathrel{\ldots}:_{W}$). This
requires a prescribed order of $\hat{a}$’s and $\hat{a}^{\dagger}$’s, and once
we have this expression, we may move the operators within without incurring a
commutation relation, e.g. $:\mathrel{\hat{a}^{\dagger}\hat{a}}:_{W}\ =\
:\mathrel{\hat{a}\hat{a}^{\dagger}}:_{W}$. The correct procedure is to expand
out the operator expression into Weyl-ordered expressions, and replace each
one with the corresponding Hafnian, i.e.:

$\displaystyle\langle\hat{n}_{1}^{m_{1}}\hat{n}_{2}^{m_{2}}\ldots\rangle$
$\displaystyle=\langle(\hat{a}_{1}^{\dagger}\hat{a}_{1})^{m_{1}}(\hat{a}_{2}^{\dagger}\hat{a}_{2})^{m_{2}}\ldots\rangle$
$\displaystyle=\sum_{\text{orderings}}c_{\text{order}}\langle:\mathrel{{\hat{a}_{1}^{\dagger}}^{m_{\text{order}1}}\hat{a}_{1}^{m_{\text{order}2}}\ldots}:_{W}\rangle$
$\displaystyle=\sum_{\text{orderings}}c_{\text{order}}\left[\sum_{\mathcal{P}(\text{order})}\prod_{(k_{1},k_{2})}\langle:\mathrel{\hat{a}_{k_{1}}\hat{a}_{k_{2}}}:_{W}\rangle\right]$
$\displaystyle=\sum_{\text{orderings}}c_{\text{order}}\text{Haf}(X\sigma_{\bar{m}_{\text{order}}}).$

It is well known how to relate normal (operator expressions where all the
$\hat{a}^{\dagger}$ precede the $\hat{a}$, denoted by $:\mathrel{\ldots}:$)
and Weyl forms [CahillGlauber, FerraroOlivaresParis], in our case:

$:\mathrel{{\hat{a}^{\dagger}}^{m}\hat{a}^{m}}:\ =\sum_{k=0}^{m}k!{m\choose
k}^{2}\left(-\frac{1}{2}\right)^{k}:\mathrel{{\hat{a}^{\dagger}}^{m-k}\hat{a}^{m-k}}:_{W}$

and all modes are treated separately. Additionally, the normal form of
$(\hat{a}^{\dagger}\hat{a})^{m}$ can be obtained by repeatedly applying
commutation relations. This can be generalized to a sum represented by
Stirling numbers $S$ [Katriel]:

$(\hat{a}^{\dagger}\hat{a})^{m}=\sum_{i=1}^{m}S(m,i){\hat{a}^{\dagger}}^{i}\hat{a}^{i}.$

In the case of $\accentset{\rightharpoonup}{m}\in\\{{0,1\\}}^{M}$, as above,
$\hat{a}^{\dagger}\hat{a}=\ :\mathrel{\hat{a}^{\dagger}\hat{a}}:_{W}-1/2$, and
we can subtract the constant directly from the covariance matrix as we did for
$\sigma^{\prime}$ ($\sigma$ is conventionally defined as above, with the
anticommutator, such that it is conveniently in Weyl form).333Note, however,
that because in this case $\hat{n}^{1}=\hat{a}^{\dagger}\hat{a}$ is already
normal ordered, we end up implicitly converting $\sigma$ to normal ordering by
subtracting $I/2$ (i.e.
$\sigma^{\prime}_{ii}=\langle:\mathrel{\hat{a}_{i}^{\dagger}\hat{a}_{i}}:\rangle=\langle\hat{a}_{i}^{\dagger}\hat{a}_{i}\rangle$
from
$\sigma_{ii}=\langle:\mathrel{\hat{a}_{i}^{\dagger}\hat{a}_{i}}:_{W}\rangle=\langle\\{\hat{a}_{i}^{\dagger},\hat{a}_{i}\\}\rangle$).
Similarly, for the second order case,

$\langle\hat{n}^{2}\rangle=\langle(\hat{a}^{\dagger}\hat{a})^{2}\rangle=\langle:\mathrel{{\hat{a}^{\dagger}}^{2}\hat{a}^{2}+\hat{a}^{\dagger}\hat{a}}:\rangle=\langle:\mathrel{{\hat{a}^{\dagger}}^{2}\hat{a}^{2}-\hat{a}^{\dagger}\hat{a}}:_{W}\rangle\Rightarrow\text{Haf}(X\sigma^{(2)})-\text{Haf}(X\sigma).$

In this example we use $\sigma^{(2)}$ to denote a larger covariance matrix
where we repeat the rows and columns of $\sigma$ for this mode.

These equations involving Hafnians will look familiar to the reader familiar
with the recent Gaussian Boson Sampling literature. Indeed, for multimode
Gaussian states, calculating the probabilities is largely the dual of
calculating the expectation values. The photon number distribution is
essentially a generalized Bose–Einstein distribution that incorporates modes,
interference and entanglement. To quote the result proven in [22, 13], the
probability of a given measurement,
$\accentset{\rightharpoonup}{n}\in\mathbb{N}^{M}$, has the form:

$\displaystyle P(\accentset{\rightharpoonup}{n})$
$\displaystyle=\frac{1}{\accentset{\rightharpoonup}{n}!\lvert\sigma+I_{2M}/2\rvert^{1/2}}\text{Haf}\left(XA_{\accentset{\rightharpoonup}{n}}\right)$
$\displaystyle A$
$\displaystyle=I_{2M}-(\sigma+I_{2M}/2)^{-1},\quad\accentset{\rightharpoonup}{n}!={\textstyle\prod_{i}}n_{i}!$

where the Hafnian’s argument refers to the principal submatrix with indices
given by $\accentset{\rightharpoonup}{n}$. To gain some intuition, one may
notice the Bose–Einstein resemblance via
$\sigma_{ii}=\langle\hat{a}_{i}^{\dagger}\hat{a}_{i}+h.c.\rangle/2=\langle
n_{i}\rangle+1/2$, so $(\sigma+I/2)_{ii}=\langle n_{i}\rangle+1$, thus in some
limiting cases:

$\displaystyle\lvert\sigma+I/2\rvert^{-1/2}\rightarrow\frac{1}{\langle
n\rangle+1},\quad I-(\sigma+I/2)^{-1}\rightarrow\frac{\langle
n\rangle}{\langle n\rangle+1}$

More precisely, the Bose–Einstein distribution, which is simply geometric, can
be thought of as stemming from the recurrence relation:

$P_{\textrm{BE}}(n)=P_{\textrm{BE}}(n-1)a=P_{\textrm{BE}}(0)a^{n}=\frac{1}{\langle
n\rangle+1}\left(\frac{\langle n\rangle}{\langle n\rangle+1}\right)^{n}$

while this general version follows the similar, yet more complex, Hafnian
(Wick’s theorem) recursive property:

$\displaystyle\forall i\text{ s.t.\ }n_{i}>0,P(n_{1},\ldots,n_{M})$
$\displaystyle=\frac{P(0,\ldots,0)}{\accentset{\rightharpoonup}{n}!}\sum_{{j=1}\atop{j\neq
i}}^{2N}a_{ij}\text{Haf}(XA_{\accentset{\rightharpoonup}{n}-\\{i,j\\}}),\quad
N=\sum_{i}{\accentset{\rightharpoonup}{n}}_{i}$

where the notation $A_{\accentset{\rightharpoonup}{n}-\\{i,j\\}}$ denotes the
subtraction of rows and columns $i,j$ from the matrix
$A_{\accentset{\rightharpoonup}{n}}$ (note e.g. with $j=i+N$,
$\text{Haf}(A_{\accentset{\rightharpoonup}{n}-\\{ij\\}})\propto
P(n_{1},\ldots,n_{i}-1,\ldots,n_{M})$, $a_{ij}=P(0,\ldots,n_{i}=1,\ldots,0)$).

From the quantum computing and complexity theory perspective, the Hafnian is
of interest as it belongs to #P, a class of classically intractable functions
[ArkhipovAaronson, 13].

A brief sketch of a derivation of the probabilities is as follows; it is
perhaps a bit more cumbersome than using phase space formalism (which
abstracts away the use of operators and orders) but it is hopefully more
transparent. We use the normal form, as this allows us to use a convenient
representation of the photon-number Fock state projection operator: the vacuum
projection is known to be

$\lvert 0\rangle\langle 0\rvert=\ :\mathrel{\exp(-\hat{a}^{\dagger}\hat{a})}:$

(see, for example, [Fan] Eqns. 17–19 for a short proof), and it follows that
[FerraroOlivaresParis]:

$\displaystyle\lvert\accentset{\rightharpoonup}{n}\rangle\langle\accentset{\rightharpoonup}{n}\rvert$
$\displaystyle=\lvert n_{1},\ldots,n_{m}\rangle\langle
n_{1},\ldots,n_{m}\rvert=\prod_{i}(\hat{a}_{i}^{\dagger})^{n_{i}}\lvert
0\rangle\langle 0\rvert\prod_{i}\hat{a}_{i}^{n_{i}}/n_{i}!=\
:\mathrel{\prod_{i=1}^{M}\exp(-\hat{a}_{i}^{\dagger}\hat{a}_{i})\frac{\left(\hat{a}_{i}^{\dagger}\hat{a}_{i}\right)^{n_{i}}}{n_{i}!}}:$
$\displaystyle P(\accentset{\rightharpoonup}{n})$
$\displaystyle=\left\langle\lvert\accentset{\rightharpoonup}{n}\rangle\langle\accentset{\rightharpoonup}{n}\rvert\right\rangle=\frac{1}{\accentset{\rightharpoonup}{n}!}\left\langle:\mathrel{\prod_{i=1}^{M}\exp(-\hat{a}_{i}^{\dagger}\hat{a}_{i})\left(\hat{a}^{\dagger}\hat{a}\right)^{n_{i}}}:\right\rangle.$

Since the operator is normal-ordered, we must also use the normal-ordered
covariance matrix, which is $\sigma^{\prime}$ as defined above (since, e.g.
for the first matrix quadrant,
$\sigma^{\prime}_{ij}=\langle:\mathrel{\hat{a}_{i}^{\dagger}\hat{a}_{j}}:\rangle=\langle\hat{a}_{i}^{\dagger}\hat{a}_{j}\rangle=\langle\\{\hat{a}_{i}^{\dagger},\hat{a}_{j}\\}\rangle-\delta_{ij}=\langle:\mathrel{\hat{a}_{i}^{\dagger}\hat{a}_{j}}:_{W}\rangle-\delta_{ij}=\sigma_{ij}-\delta_{ij}$).
Lastly, notice that we can set

$\lvert 0\rangle\langle 0\rvert=\
:\mathrel{\exp\left(-\tfrac{1}{2}{\hat{\xi}}^{\dagger}\hat{\xi}\right)}:$

so that the notation is compatible with that of the state’s covariance matrix.

Hence, the problem is reduced once more to calculus; multivariate Gaussian
integrals. We first find $P(0)$, and the rest of the probabilities follow from
Wick’s theorem, with a modified covariance matrix due to the above $\lvert
0\rangle\langle 0\rvert$ Gaussian. The former can be solved with the usual
tricks. Let $Z\sim\mathcal{N}(\accentset{\rightharpoonup}{0},I)$ be a vector
of standard normal random variables. By affine transformation, we can
transform it to a vector of arbitrary jointly Gaussian variables:
$\sigma^{\prime 1/2}Z+\mu$. Since we are only interested in the case of zero
displacement, let the means $\mu=0$. We thus simplify our calculation by
computing the expectation over the multivariate standard normal distribution:

$\left\langle\exp\left(-\tfrac{1}{2}(\sigma^{\prime
1/2}Z)^{\dagger}\sigma^{\prime
1/2}Z\right)\right\rangle=(2\pi)^{-M}\int_{\mathbb{R}^{2M}}dZ\exp\left(-\tfrac{1}{2}Z^{\dagger}(\sigma^{\prime}+I)Z\right)=|\sigma^{\prime}+I|^{-1/2}.$

This is the zero-photon probability. We note that
$\sigma^{\prime}+I=\sigma+I/2$ corresponds to the anti-normal covariance
matrix (all the $\hat{a}$ precede the $\hat{a}^{\dagger}$).

Finally, we can include monomial terms in the integral to calculate other
probabilities, using Wick’s theorem once more, after solving for the effective
covariance matrix, $A$ – since the distribution and operator Gaussians
combine:

$\exp\left(-\tfrac{1}{2}Z^{\dagger}Z-\tfrac{1}{2}Z^{\dagger}\sigma^{\prime-1}Z\right)=\exp\left(-\tfrac{1}{2}Z^{\dagger}A^{-1}Z\right).$

$A$ can be found by invoking the matrix inversion lemma:

$\displaystyle A^{-1}$ $\displaystyle=\sigma^{\prime-1}+I$ $\displaystyle A$
$\displaystyle=(\sigma^{\prime-1}+I)^{-1}$
$\displaystyle=I-\sigma^{\prime-1}(\sigma^{-1}+I)^{-1}$
$\displaystyle=I-(\sigma^{\prime}+I)^{-1},$

as per convention. We have not been completely rigorous, however, since
$\sigma^{\prime}$ may not be invertible (if, in a diagonal basis, sigma has
some $\langle:\mathrel{\hat{a}_{i}^{\dagger}\hat{a}_{i}}:\rangle=0$, 444The
field can never be a pure vacuum at finite temperature, although for all
intents and purposes it is at optical frequencies, since the photon energy
$\hbar\omega\gg k_{B}T$ (the Boltzmann constant temperature product) at room
temperature. in which case the Gaussian tends towards a delta function);
hence, besides for the last line, the above equalities may not be strictly
valid. Nonetheless, the moments exist, which is what we need for the moment
expansion in Wick’s theorem.

### Photon number expectation values

We can solve for some lower order photon number moments in terms of $\sigma$’s
elements. To eliminate local phase degrees of freedom, we set all diagonal
$U_{ii}$ terms to be real and positive. Here, wherever relevant, we have
$i<j<k$ (switching order requires taking complex conjugates of some terms and
can therefore introduce inconsistencies).

$\displaystyle\langle n_{i}\rangle$
$\displaystyle=\text{Haf}(X\sigma_{i}^{\prime})=V_{ii}$ $\displaystyle\langle
n_{i}n_{j}\rangle$
$\displaystyle=\text{Haf}(X\sigma_{ij}^{\prime})=|U_{ij}|^{2}+|V_{ij}|^{2}+\langle
n_{i}\rangle\langle n_{j}\rangle$ $\displaystyle\langle n_{i}^{2}\rangle$
$\displaystyle=\text{Haf}(X\sigma_{ii})-\text{Haf}(X\sigma_{i})=|U_{ii}|^{2}+2\langle
n_{i}\rangle^{2}+\langle n_{i}\rangle$ $\displaystyle\langle
n_{i}n_{j}n_{k}\rangle$ $\displaystyle=\text{Haf}(X\sigma_{ijk}^{\prime})$
$\displaystyle=\langle n_{i}\rangle\langle n_{j}n_{k}\rangle+\langle
n_{j}\rangle\langle n_{i}n_{k}\rangle+\langle n_{k}\rangle\langle
n_{i}n_{j}\rangle-2\langle n_{i}\rangle\langle n_{j}\rangle\langle
n_{k}\rangle$
$\displaystyle\quad+2\mathcal{R}\left(U_{ij}^{*}(V_{ik}U_{jk}+V_{jk}U_{ik})+V_{ij}(U_{ik}^{*}U_{jk}+V_{ik}^{*}V_{jk})\right)$
$\displaystyle\langle n_{i}^{2}n_{j}\rangle$
$\displaystyle=\text{Haf}(X(\sigma_{iij}-\tfrac{1}{2}\delta_{jj})-\text{Haf}(X(\sigma_{ij}-\tfrac{1}{2}\delta_{jj}))$
$\displaystyle=\langle n_{i}n_{j}\rangle(4\langle n_{i}\rangle+1)+\langle
n_{j}\rangle(|U_{ii}|^{2}-2\langle
n_{i}\rangle^{2})+2U_{ii}(U_{ij}^{*}V_{ij}+U_{ij}V_{ij}^{*})$
$\displaystyle\langle n_{i}n_{j}^{2}\rangle$
$\displaystyle=\text{Haf}(X(\sigma_{ijj}-\tfrac{1}{2}\delta_{ii})-\text{Haf}(X(\sigma_{ij}-\tfrac{1}{2}\delta_{ii}))$
$\displaystyle=\langle n_{i}n_{j}\rangle(4\langle n_{j}\rangle+1)+\langle
n_{i}\rangle(|U_{jj}|^{2}-2\langle
n_{j}\rangle^{2})+2U_{jj}(U_{ij}V_{ij}+U_{ij}^{*}V_{ij}^{*})$
$\displaystyle\langle n_{i}^{2}n_{j}^{2}\rangle$
$\displaystyle=\text{Haf}(X\sigma_{iijj})+\text{Haf}(X\sigma_{ij})-\text{Haf}(X\sigma_{iij})-\text{Haf}(X\sigma_{ijj})$
$\displaystyle=\langle n_{i}^{2}n_{j}\rangle(4\langle n_{j}\rangle+1)+\langle
n_{i}n_{j}^{2}\rangle(4\langle n_{i}\rangle+1)-\langle
n_{i}n_{j}\rangle(4\langle n_{i}\rangle+1)(4\langle n_{j}\rangle+1)$
$\displaystyle\quad+4(\langle n_{i}n_{j}\rangle-\langle n_{i}\rangle\langle
n_{j}\rangle)^{2}+(2\langle n_{i}\rangle^{2}-|U_{ii}|^{2})(2\langle
n_{j}\rangle^{2}-|U_{jj}|^{2})$
$\displaystyle\quad+2U_{ii}U_{jj}(V_{ij}^{2}+U_{ij}^{2}+{V_{ij}^{*}}^{2}+{U_{ij}^{*}}^{2})+8|U_{ij}|^{2}|V_{ij}|^{2}.$

A few important observations. The equations for the statistics can be
separated into a trivial component (composed of lower order statistics) and an
interference term, due to complex-valued $U$ and $V$ terms. If the $U$ terms
are zero, there are no nontrivial higher order correlations: everything can be
described by second- or first-order statistics. Hence why thermal states do
not have interesting higher order correlations.

Relatedly, where $U=0$, there are a vast number of states that have the same
photon number distribution, as there is no information beyond photon number
mean and covariance. Measurements carry no information about the underlying
supermodes.

The value of $U_{ii}\in[0,\langle n_{i}\rangle(\langle n_{i}\rangle+1)]$, and
the two extremes correspond to two-mode and single-mode squeezing
respectively. States that have intermediate values can be thought of as a
multimode generalization of these two concepts, where the value of $U_{ii}$
represents the locality of the entanglement. When $U_{ii}=0$, certain higher
order correlations also lose interference terms. Therefore, measurements of
two-mode squeezed vacuum also lack certain information about the underlying
state.

Finally, for a state with sufficient information encoded in $U$, the inverse
problem of retrieving the state and modes from observations is possible,
within some limitations. This is within the realm of possibility and future
experiments should seek to achieve this.

Returning to the experiment, the photon covariance matrix in Fig. 3b (and the
simulation in Fig. 2a) consists of contributions from both the $V$ (diagonal)
and $U$ (anti-diagonal) components.

##  Appendix B EMCCD camera statistics

The electron-multiplying (EM) gain process used in EMCCD cameras is what makes
these instruments highly sensitive: by magnifying the charge of a few
photoelectrons to macroscopic levels, the camera is capable of detecting
single photons. However, the nature of the gain process is stochastic, which
introduces noise. Here we introduce the model for the gain process, and
discuss how it influences the measurements we are interested in.

### Photon statistics after EM gain

The EM gain process adds a high amount of noise, which typically precludes the
possibility of resolving the exact photon numbers in EMCCD cameras. When $n$
photons are captured as photoelectrons on the CCD sensor and amplified with EM
gain $g$, the conditional probability of yielding $x$ amplified electrons
follows a gamma (Erlang) distribution [16]:

$\displaystyle P(x;n)$ $\displaystyle=\frac{x^{n-1}\exp(-x/g)}{g^{n}(n-1)!}.$

To estimate the statistics of the EM gain output $x$, let $P_{n}$ be the
photon-number distribution incident on a pixel. We are first able to calculate
the conditional moments of $x$ with respect to $n$, use these to calculate
(unconditional) moments of $x$, and finally calculate the photon number
moments:

$\displaystyle\langle x^{k}\rangle_{n}$
$\displaystyle=\int_{0}^{\infty}x^{k}P(x;n)dx=g^{k}\frac{(n+k-1)!}{(n-1)!},$
$\displaystyle\langle x^{k}\rangle$ $\displaystyle=\sum_{n\geq 0}P_{n}\langle
x^{k}\rangle_{n},$ $\displaystyle g^{-k}\langle x^{k}\rangle$
$\displaystyle=\sum_{n\geq
0}P_{n}\prod_{i=0}^{k-1}n+i=\left\langle\prod_{i=0}^{k-1}n+i\right\rangle.$

This can be extended to correlations between multiple variables, e.g.:

$\displaystyle\langle x_{i}^{k}x_{j}^{l}\rangle$
$\displaystyle=\sum_{n_{i},n_{j}\geq 0}P_{n_{i}n_{j}}\langle
x_{i}^{k}\rangle_{n_{i}}\langle
x_{j}^{l}\rangle_{n_{j}}=g^{k+l}\left\langle\left(\prod_{m=0}^{k-1}n_{i}+m\right)\left(\prod_{m=0}^{l-1}n_{j}+m\right)\right\rangle.$

We may then solve for the photon number moments by inverting the above
equations. Hence we can write down the formulae for our photon statistics of
interest (omitting the gain term $g$):

$\displaystyle\langle n\rangle$ $\displaystyle=\langle x\rangle$
$\displaystyle\langle n_{i}n_{j}\rangle$ $\displaystyle=\langle
x_{i}x_{j}\rangle$ $\displaystyle\langle n_{i}n_{j}n_{k}\rangle$
$\displaystyle=\langle x_{i}x_{j}x_{k}\rangle$ $\displaystyle\langle
n^{2}\rangle$ $\displaystyle=\langle x^{2}\rangle-\langle x\rangle$
$\displaystyle\langle n_{i}^{2}n_{j}\rangle$ $\displaystyle=\langle
x_{i}^{2}x_{j}\rangle-\langle x_{i}x_{j}\rangle$ $\displaystyle\langle
n_{i}^{2}n_{j}^{2}\rangle$ $\displaystyle=\langle
x_{i}^{2}x_{j}^{2}\rangle-\langle x_{i}^{2}x_{j}\rangle-\langle
x_{i}x_{j}^{2}\rangle+\langle x_{i}x_{j}\rangle.$

Finally, the electronic stages between the EM gain process and a digitized
pixel value will introduce additional noise, such as readout noise. However
this noise is independent of the signal $x$, hence it does not affect the
above our ability to estimate the photon number.

Figure A2:  Electron multiplying gain model and thresholding. a. The
probability distribution of an amplified signal for a given number of
photoelectrons, prior to additional noise. We model this with the Erlang
distribution. The signal $x$ is measured in electrons ($e^{-1}$) divided by
the gain $g$ (the expected value of this is the original number of
photelectrons). b. The probability distribution of the signal with Gaussian
readout noise. Although the two probability distributions begin to overlap, in
this example, thresholding can be used to distinguish between 0 and 1 (or
more) photons with some degree confidence, incurring a tradeoff between the
photon-detection efficiency and false count rate. The readout noise has
standard deviation $\sigma$. The inset shows the analytical (ideal) ROC curve
for this example.

### Thresholded operation

As discussed above, the electronic signal $x$, ensuing an EM process with a
gain of $g$, follows a random distribution parametrized by the number of
photons captured $n$ – this is shown in Fig. A2. Due to the stochastic nature
of the EM gain process, it is impossible to resolve the original number of
photoelectrons $n$ on the CCD sensor from the number of amplified electrons
$x$. However, if there are no photons absorbed by the CCD sensor during the
detection window, there is no signal (no electrons) to amplify in the gain
process. Hence, there is no excess noise from the EM process, and the variance
of the output signal depends solely on the readout noise $\sigma$. With a high
EM gain, the effective readout noise, $\sigma\ll 1\,e^{-}$, we can set a
threshold on the output signal to discriminate between the absence or presence
of photoelectrons with a high degree of confidence.

We evaluate the EMCCD camera in the single-photon detection mode using a
threshold to detect the absence or presence of a photon. The threshold value
will determine the photon-detection efficiency (PDE) and the false click rate.
In general, a lower threshold increases the effective PDE but increases the
false click rate, and vice versa. See Fig. A2b. By varying the threshold
value, we obtain a receiver operating characteristic (ROC) curve to
characterize the performance of the photon counter (Fig. A2b inset). The ROC
curve quantifies the trade-off in the effective quantum efficiency and dark
count rate. In practice, we obtain the false click rate as a function of the
threshold by obtaining dark frames – images from the sensor with no
illumination – and count how many times a given threshold is exceeded. We
obtain the PDE by determining the readout noise from this data, and computing
the PDE through the model by comparing $g$ and $e$, multiplied by the quantum
efficiency (QE) of the CCD sensor.

##  Appendix C Coincidence detection

In this section we explain how coincidence detection may be used as a tool to
distinguish between different pulsed photonic states, as well as measure the
purity of the state.

### Derivation

Figure A3:  Toy model for the coincidence detection experiment.

Consider a beam of squeezed light incident on a balanced beamsplitter, which
is followed by two detectors on each side of the output (Fig. A3). The first
arm of the beamsplitter, with operator $\hat{a}_{1}$, has some
$\langle\hat{n}\rangle$ and
$\operatorname{Var}(\hat{n})=2\langle\hat{n}\rangle(\langle\hat{n}\rangle+1)$,
in the lossless case ($\eta=1$). The second input port, with operator
$\hat{a}_{2}$ has vacuum input.

The third and fourth ports, defined by operators $\hat{a}_{3},\hat{a}_{4}$,
have