# Photon number distribution of squeezed light from a silicon nitride
microresonator measured without photon number resolving detectors

Emanuele Brusaschi  Massimo Borghi  corresponding author:
<EMAIL_ADDRESS>Marcello Bacchi  Dipartimento di Fisica, Università di
Pavia, Via Bassi 6, 27100 Pavia, Italy.  Marco Liscidini  Matteo Galli
Dipartimento di Fisica, Università di Pavia, Via Bassi 6, 27100 Pavia, Italy.
Daniele Bajoni  Dipartimento di Ingegneria Industriale e dell’Informazione,
Università di Pavia, Via Ferrata 5, 27100 Pavia, Italy.

###### Abstract

The measurement of the photon number distribution (PND) allows one to extract
metrics of non-classicality of fundamental and technological relevance, but in
principle it requires the use of detectors with photon number resolving (PNR)
capabilities.

In this work we reconstruct the PND of two-mode pulsed squeezed light
generated from a silicon nitride microresonator using threshold detectors and
variable optical attenuations. The PNDs are characterized up to $\sim 1.2$
photons/pulse, through which we extracted an on-chip squeezing level of
$6.2(2)$ dB and a noise reduction factor of $-3.8(2)$ dB. The PNDs are
successfully reconstructed up to an Hilbert space dimension of $6\times 6$.
The analysis performed on the photon-number basis allows us to characterize
the influence of a spurious thermal background field that spoils the photon
number correlations. We evaluate the impact of self and cross phase modulation
on the generation efficiency in case of a pulsed pump, and validate the
results through numerical simulations of the master equation of the system.

††preprint: APS/123-QED

##  I  Introduction

The characterization of quantum optical states in the photon number basis is
of great relevance in fundamental science, quantum communication, and quantum
information processing. For example, non-classical optical states are
characterized by a negative Mandel-Q parameter, measuring the departure of
their PND from that of a coherent state [1]. Similarly, from the PND we can
determine the sign of the Wigner function at the origin [2] and of the
eigenvalues of the matrix containing the high-order correlation moments [3].
In both cases, a negative value represents a criterion of non-classicality. In
decoy-state quantum key distribution (QKD), one can estimate the amount of
information an eavesdropper could get by sending pulses with different mean
photon numbers and of which the PND is known [4]. The latter requirement stems
from the fact that even if mean photon number is known, the photon number
distribution still affects the performance of the QKD system. Therefore, it is
important to characterize it to select the states which are more robust to
photon number splitting attacks [5].
The PND of bipartite optical states carries information on both inter-beam
photon number correlations and intra-beam photon statistics. Each of the
conjugate beams generated by spontaneous parametric down conversion (SPDC) or
spontaneous four wave mixing (SFWM), is usually described by a multi-mode
thermal state. The number of these modes is called the Schmidt number and
determines the amount of entanglement of the bipartite state [6]. The Schmidt
number is tailored on the requirements of the application. For example, in
quantum computing and simulation, single photons should be heralded in a well
defined spectro-temporal mode to achieve high visibility quantum interference
[7]. Inter-beam correlations in the photon number are instead quantified by a
reduced variance in the photon number difference between the two beams with
respect to that of a coherent state [8]. Finally, the shrinkage of the
variance is called noise reduction factor (NRF). A value below unity indicates
photon number difference squeezing, which is a signature of non-classicality.
A prominent example of a quantum optical state with strong photon number
correlation is squeezed light, which is a fundamental resource for quantum
information processing and simulation. In particular, squeezed light can be
efficiently generated from integrated microresonators. Broadband single [9]
and two-mode [8, 10] quadrature squeezing have been already reported for these
types of sources and characterized by coherent photodetection. While homodyne
tomography allows us to fully reconstruct the complex Wigner function of the
optical state, and thus also the coherences that may exist between different
Fock states, its implementation can be very challenging for two-mode squeezed
light because of the need of phase locking two local oscillators to the pump
beam with hundreds of GHz of separation. [8, 9]. Alternatively, the PND of
squeezed light can be measured by photon number resolving detectors (PNR). For
example, transition edge sensors [11] have been used in [8] to measure the NRF
of squeezed light generated by a silicon nitride resonator, and in [12, 13] to
detect patterns of bosonic coalescence at the output of sampling machines.
However, PNR detectors have long dead-times ($>$ tens of $\mu s$), they are
not commercial and they work at temperatures below 100 mK.
The large-scale deployment of integrated sources of squeezed light on
different technology platforms and in multidisciplinary research fields drive
the demand for characterization procedures which are faster, cheaper and
simpler than those which implement PNR detectors. An ideal solution is the use
of standard superconducting nanowire single photon detectors (SNSPD) or off-
the-shelf avalanche photodiodes. However, their PNR capability is so far
limited to a small number of photons ($<3-4$) and they heavily rely on
electronic post processing [14]. Nevertheless, they can provide PNR
capabilities when employed in temporal or spatial multiplexing schemes. These
strategies fall under the umbrella of pseudo-PNR, which include fiber loop
detectors [15], beamsplitter trees [16], cascaded delay-lines [2] .
A parallel strategy for retrieving the PND that is efficient in hardware
resources exploits the nonlinear photodection probability to different photon
number states [17] and in presence of losses [18]. The nonlinearity introduces
photon number sensitivity to the detected on-off statistics, from which the
photon number distribution can be inferred by statistical methods, such as the
Expectation-Maximization algorithm (EM) [19]. The PND of a broad class of
quantum optical states have been retrieved using this strategy, including
heralded photons from SPDC, superposition of single photon states, coherent
and thermal states [20].
In this work, we extend this approach to pulsed two-mode squeezed vacuum (TMS)
generated by SFWM in a silicon nitride microresonator. We show that from the
reconstructed PND we can detect the presence of parasitic thermal noise,
mainly originating from spontaneous Raman scattering in the resonator, which
reduces the photon number correlation between the squeezed modes. We recover
the average photon number per pulse, the squeezing parameter, and the NRF.
Finally, we compare our results with numerical simulations of the master
equation of the system, including the influence of self and cross-phase
modulation (SPM/XPM), and discuss the limits of application of the method.

##  II  Theory of photon number reconstruction using on-off detectors

Here we briefly review the method presented in [21] for the reconstruction of
the PND using on-off detectors and variable attenuation. The PND $\rho_{n}$ of
a single mode quantum optical state is given by
$\rho_{n}=\textrm{Tr}(\hat{\rho}\ket{n}\bra{n})=\left[\hat{\rho}\right]_{nn}$,
where $\left[\hat{\rho}\right]_{nn}$ are the diagonal elements of the density
matrix $\hat{\rho}$ on the photon number basis $\ket{n}$. Here we define
$\ket{n}=(n!)^{-\frac{1}{2}}\left(A^{\dagger}\right)^{n}\ket{\textrm{vac}}$,
where the operator $A^{\dagger}$ creates a single photon in a well defined
spectral-temporal mode. The operator $A^{\dagger}$ will be associated in
Section II.1 to the single Schmidt mode in which the signal and the idler
photon will be generated. The measurement process of an on/off detector with
detection efficiency $\eta$ is described by a two-value positive operator-
valued measure (POVM) $\Pi_{\textrm{on}/\textrm{off}}(\eta)$. The off-event in
which no photons are detected is represented by

$\Pi_{\textrm{off}}(\eta)=\sum_{n=0}^{\infty}(1-\eta)^{n}\ket{n}\bra{n},$ (1)

and the related probability outcome
$p_{0}(\eta)=\textrm{Tr}(\Pi_{\textrm{off}(\eta)}\hat{\rho})$ is

$p_{0}(\eta)=\sum_{n=0}^{\infty}(1-\eta)^{n}\rho_{n}.$ (2)

By performing $M$ measurements with different quantum efficiencies
$\eta_{\mu}$ ($\mu=1,...,M)$ we construct a vector $\bm{p_{0}}$ of off-events
probabilities

$\bm{p_{0}}=\big{[}p_{0}(\eta_{\mu}):p_{0}(\eta_{\mu})=\sum_{n=0}^{\infty}B_{\mu
n}\rho_{n}\quad\mu=1,2,..,M\big{]},$ (3)

where $B_{\mu n}=(1-\eta_{\mu})^{n}$. Equation (3) corresponds to a linear
system that can be written in compact form as $\bm{p_{0}}=B\bm{\rho}$. It can
be interpreted as a linear positive statistical model for the unknown
parameters $\rho_{n}$ that can be solved through the EM algorithm [21]. In
this iterative approach, the estimate of the PND element $\rho_{n}^{(i+1)}$ at
the step $i+1$ is calculated as

$\rho_{n}^{(i+1)}=\rho_{n}^{(i)}\sum_{\mu=1}^{M}\frac{B_{\mu
n}f_{0}(\eta_{\mu})}{(\sum_{\nu}B_{\nu
n})p_{0}(\eta_{\mu})\\{\bm{\rho^{(i)}}\\}},$ (4)

where $f_{0}(\eta_{\mu})$ is the experimental off-probability measured at the
detection efficiency $\eta_{\mu}$, while
$p_{0}(\eta_{\mu})\\{\bm{\rho}^{(i)}\\}$ is calculated by Eq.(3) through the
estimate of $\bm{\rho}$ at step $i$. The convergence of the algorithm is
checked by evaluating the error parameter $\epsilon$, which is defined as the
geometric distance between $f_{0}({\eta_{\mu}})$ and the reconstructed off-
probabilities at the $i^{\textrm{th}}$ step

$\epsilon^{(i)}=\frac{1}{M}\sum_{\mu=1}^{M}|f_{0}(\eta_{\mu})-p_{0}(\eta_{\mu})\\{\bm{\rho^{(i)}}\\}|.$
(5)

The similarity between $\bm{\rho}$ and a target theoretical distribution
$\bm{\rho}^{\textrm{theo}}$ is evaluated with the fidelity $\mathcal{F}$,
defined as

$\mathcal{F}^{(i)}=\sum_{n=0}^{N}\sqrt{\rho_{n}^{(\textrm{theo})}\rho_{n}^{(i)}}.$
(6)

The method can be generalized to the bipartite case (see Appendix C) to
reconstruct a two-dimensional PND $\rho_{nn^{\prime}}$ from the on-off click
probabilities and the coincidences between the two channels $x$ and $y$. In
the bipartite case there are four possible outcomes $p_{xy}(\eta)$ for each
measurement, where $x(y)=0$ labels a no-click event in the detector $x(y)$,
while $x(y)=1$ denotes a click.
In general, the efficacy of the method depends on the maximum efficiency
$\eta_{\textrm{max}}$ and on the level of noise in the data. First, if
$n\eta_{\textrm{max}}\ll 1$, then $(1-\eta_{\textrm{max}})^{n}\sim
1-\eta_{\textrm{max}}n+\mathcal{O}((n\eta_{\textrm{max}})^{2})$, and Eq.(3)
reduces to $p_{0}(\eta_{\mu})=\sum_{n}\eta n\rho_{n}=\eta\langle n\rangle$,
i.e, the method can only solve the average number of photons but not the
individual $\rho_{n}$. Hence, we require that
$\eta_{\textrm{max}}\sim\frac{1}{\langle n\rangle}$ for optimal performance,
which means that losses should be kept as low as possible in order to
reconstruct the PND of states with small $\langle n\rangle$. Second, if we
apply the binomial expansion in Eq.(3) we obtain
$p_{0}(\eta)=\sum_{k=0}^{\infty}\left(\sum_{n=0}^{k}\rho_{n}C_{nk}\right)\eta^{k}=\sum_{k=0}^{\infty}Q_{k}\eta^{k}$,
where $C_{nk}=(-1)^{k}{n\choose k}$ and $Q_{k}=\sum_{n=0}^{k}\rho_{n}C_{nk}$.
Then, $\frac{\partial^{m}p_{0}}{\partial\eta^{m}}=m!Q_{m}$, which implies that
only the elements of $\rho_{n}$ from $n=0$ to $n=m$ contribute to the
$m^{\textup{th}}$ derivative of $p_{0}$. Thus, the next leading order term
$\rho_{m+1}$ in the PND can be reliably estimated only if
$\frac{\partial^{m+1}p_{0}}{\partial\eta^{m+1}}$ is accurately determined from
the experiment, which imposes stringent demands on the precision and accuracy
with which $p_{0}(\eta)$is measured, especially for large $n$. Rigorous
Cramér-Rao bounds on the estimate of $\rho_{n}$ are derived in [21]. In
Appendix B we experimentally validate the method by reconstructing the PND of
known single-mode optical thermal and the coherent states. In both cases, we
obtain fidelities higher than $0.99$ with the theoretical target
distributions.

###  II.1 Photon number statistics of squeezed light from a lossy
microresonator

Figure 1: (a) Sketch of a lossy microresonator in which photons can decay in
the waveguide at a rate $\gamma_{e}$ and being lost with a rate $\gamma_{i}$.
(b) Equivalent model of the lossy resonator shown in panel (a), in which the
losses are accounted by coupling the resonator with the auxiliary Loss
waveguide. (c) Series of transformations that the TMS $\rho$ undergoes from
generation to detection. (b). The three panels above the beamsplitters with
transmittivity $\eta_{e}$ and $\eta$ sketch the joint probability distribution
of the signal and the idler photon numbers $n_{s}$ and $n_{i}$ at the points
(1), (2) and (3).

We aim to reconstruct the PND of two-mode squeezed light at the output of a
lossy silicon nitride microresonator [8]. The state is produced by pumping the
resonator with a strong pulsed pump with central frequency $\omega_{p}$, which
generates signal and idler conjugate beams centered at frequencies
$\omega_{s}$ and $\omega_{i}$, respectively, with
$2\omega_{p}=\omega_{s}+\omega_{i}$. In order to understand, interpret and
compare the experimental PND to that expected from theory, it is important to
develop a model of SFWM in the resonator which takes into account both the
pulsed nature of the pump, the presence of losses, and of the background
noise. First, the presence of losses and spurious noise alters the photon
number correlation between the squeezed modes. Second, we can relate Eq.(1) to
the Mandel operator describing a missed photodetection only in the case where
the signal and idler photons are confined to a single and well defined
spectral-temporal mode [22]. The model should then be used to identify these
modes, and validate to which extent SPM and XPM may alter their total number
and distribution [6].
We start to build our model from the geometry shown in Fig.1(a). The resonator
is coupled to a bus waveguide with a rate $\gamma_{e}$, and at the same time
leaks photons into the external environment with a rate $\gamma_{i}$. This
system is equivalent to a lossless resonator in which an additional waveguide,
called _phantom channel_ , is coupled to the cavity with a decay rate
$\gamma_{i}$ [23], which is schematically shown in Fig.1(b). One way to solve
the spontaneous problem is to write the hamiltonian as
$H=H_{\textrm{wg}}+H_{\textrm{ph}}+H_{\textrm{r}}+H_{\textrm{c}}+H_{\textrm{NL}}$,
where the first three terms describe the free-evolution of the
waveguide(phantom) and resonator modes, $H_{c}$ mediates the coupling between
the ring and the two channels and $H_{\textrm{NL}}$ is responsible for SFWM,
self and cross-phase modulation (SPM/XPM) [24]. Alternatively, the linear part
of $H$ can be diagonalized to find the eigenmodes of the coupled ring-
waveguide(phantom) system, and $H_{\textrm{NL}}$ can be written on this bases.
These modes are called asymptotic fields and they span a continuum of
frequencies [25]. On these basis, the hamiltonian $H_{\textup{NL}}$ in the
interaction picture can be written as [26]

$\begin{gathered}H_{\textup{NL}}(t)=\hbar\sum_{m,n}\left\\{\int\Delta_{s}^{(m,n)}(\omega_{s},\omega_{s}^{\prime},t)a_{m}^{\dagger}(\omega_{s})a_{n}(\omega_{s}^{\prime})d\omega_{s}d\omega_{s}^{\prime}+\right.\\\
+\int\Delta_{i}^{(m,n)}(\omega_{i},\omega_{i}^{\prime},t)a_{m}^{\dagger}(\omega_{i})a_{n}(\omega_{i}^{\prime})d\omega_{i}d\omega_{i}^{\prime}+\\\
+\left.\left(\int\xi^{(m,n)}(\omega_{s},\omega_{i},t)a_{m}^{\dagger}(\omega_{s})a_{n}^{\dagger}(\omega_{i})d\omega_{s}d\omega_{i}+\textrm{H.c.}\right)\right\\},\end{gathered}$
(7)

where the terms $\Delta^{(m,n)}_{s(i)}$ describe the frequency shifts induced
by XPM (diagonal terms) and possible mode couplings (off-diagonal terms),
while $\xi^{(m,n)}$ accounts for SFWM and depends on the strength of the
nonlinearity, the internal pump power, the field enhancements at the pump,
signal and idler wavelengths and on the mode volume of the resonator [26]. The
indices $(m,n)=\\{\textrm{wg},\textrm{ph}\\}$ indicate the channel in which
photons are scattered in the asymptotic limit of $t\rightarrow\infty$, when
all nonlinear interactions have already occurred. In writing Eq.(7), the pump
is treated as a classical field and pump-depletion is neglected. It is shown
in Appendix D that the following relations hold

$\frac{\xi^{\textrm{wg,ph}}}{\xi^{\textrm{wg,wg}}}=\frac{\xi^{\textrm{ph,wg}}}{\xi^{\textrm{wg,wg}}}=\sqrt{\frac{\gamma_{i}}{\gamma_{e}}},\quad\frac{\xi^{\textrm{ph,ph}}}{\xi^{\textrm{wg,wg}}}=\frac{\gamma_{i}}{\gamma_{e}},$
(8)

and identical properties relate the quantities $\Delta^{(m,n)}_{s(i)}$. In
Eq.(8) we assumed that the decay rates $\gamma_{e}$ and $\gamma_{i}$ are the
same at the pump, signal and idler frequencies. One can show (see Appendix D)
from Eq.(8) and Eq.(7) that $H$ describes a TMS impinging on a beamsplitter
with transmittivity $\frac{\gamma_{e}}{\gamma_{e}+\gamma_{i}}=\eta_{e}$. This
equivalent picture is shown in Fig.1(c). We can embed the loss $\eta_{e}$ on
the transmitted photons in the detection efficiency
$\eta^{\prime}=\eta_{e}\eta$, and focus on the PND of the state immediately at
the output of the lossless resonator (point (1) in Fig.1(c)). The quadratic
hamiltonian in Eq.(7) generates a multimode squeezed vacuum state over a set
of $\mathcal{K}$ orthogonal modes $A^{(s)}_{k}$ and $A^{(i)}_{k}$, called the
Schmidt modes [6, 26]. We can then write the state $\ket{\Psi}$ at the output
of the resonator (point (1) in Fig.1(c)) as
$\ket{\Psi}=\otimes_{k=1}^{\mathcal{K}}S_{k}(r_{k})\ket{\textrm{vac}}$, where
$S_{k}=\exp\left(r_{k}{A_{k}^{(s)}}^{\dagger}{A_{k}^{(i)}}^{\dagger}-\textrm{H.c.}\right)$
is the two-mode squeezing operator associated to the $k^{\textup{th}}$ Schmidt
mode and $r_{k}$ is the corresponding squeezing parameter. If the pump, the
signal and the idler quality factors
$Q=\frac{\omega}{(\gamma_{e}+\gamma_{i})}$ are equal, and if the pump pulse
duration is shorter than the dwelling time of photons within the cavity, the
Schmidt number $K$ is $K\sim 1.075$ [27]. These conditions are satisfied in
our experiment, so we assume to work in the approximation of one Schmidt mode.
Under this hypothesis, the PND is given by
$\rho^{\textrm{SFWM}}_{nn^{\prime}}=\left(\frac{\tanh^{n}(r)}{\cosh(r)}\right)$
[26]. It follows that the marginal PND of the signal and the idler beam is
thermal with a mean photon number of $\langle n_{s}\rangle=\langle
n_{i}\rangle=\langle n\rangle=\sinh^{2}(r)$. At low powers, the relation
between the squeezing parameter and the input power $P$ is linear, i.e.,
$r=aP$. Furthermore, when $aP\ll 1$, we have that $\langle n\rangle\sim
aP^{2}$, which implies that the rate of signal and idler photons generated by
SFWM quadratically grows with the pump power. This breaks at higher power,
where the spectral distribution of the signal and the idler photons modifies
due to SPM/XPM shifts and to time ordering effects in the state evolution
governed by the hamiltonian in Eq.(7) [28, 26, 29]. In this strongly nonlinear
regime, no simple relation holds between $r$ and $P$. The real PND should take
into account the contribution of spurious noise photons that can be generated
in the resonator due to spontaneous Raman scattering or fluorescence, whose
rate scales as $\sim b_{s(i)}P$, in which $b_{s(i)}$ describes the efficiency
of the parasitic process with the pump power. This contribution is
uncorrelated from that of SFWM, and modifies the density matrix
$\rho=\ket{\Psi}\bra{\Psi}$ as
$\rho\rightarrow\rho\otimes\rho^{\textrm{th,s}}\otimes\rho^{\textrm{th,i}}$,
where $\rho^{\textrm{th,s(i)}}$ is the density matrix of a thermal state with
mean photon number $\langle n_{\textrm{s(i),th}}\rangle=b_{s(i)}P$. The
modified PND $\rho_{nn^{\prime}}$ is the convolution of the three uncorrelated
PNDs

$\rho_{nn^{\prime}}=\sum_{k,q,p}\rho^{\textrm{SFWM}}_{kk}\rho^{\textrm{th,s}}_{p}\rho^{\textrm{th,i}}_{q}\delta_{n,k+p}\delta_{n^{\prime},k+q}.$
(9)

The PND at the output of the lossy resonator (point (2) in Fig.1(b,c))
$\rho_{nn^{\prime}}^{\textrm{out}}$ is given by
$\textrm{Tr}\big{[}U_{\textrm{BS}}(\rho\otimes\rho^{\textrm{th,s}}\otimes\rho^{\textrm{th,i}})U_{\textrm{BS}}^{\dagger}\ket{nn^{\prime}}\bra{nn^{\prime}}\big{]}$,
where $U_{\textrm{BS}}$ is the generator of the beamsplitter transformation
with transmittivity $\eta_{e}$. With reference to Fig.1(b,c) and by using the
method described in Section II, $\rho_{nn^{\prime}}^{\textrm{out}}$ can be
reconstructed by setting the losses to $\eta$. Indeed, these are the losses
experienced by the signal and the idler beam from point (2) in Fig.1(c,b) to
the detectors. On the other hand, $\rho_{nn^{\prime}}$ is retrieved by setting
the losses to $\eta\eta_{e}$ (point (1) in Fig.1(c)). The expression in Eq.(9)
strictly holds only for $\rho_{nn^{\prime}}$, because it implies a perfect
photon number correlation between the signal and the idler in the limit of
vanishing noise.

##  III  Device and experimental set-up

Figure 2: (a) Sketch of the experimental setup. The pump (green) pulse
duration is $\tau$ and the repetition rate $1/T$. The signal and idler beams
are shown with red and blue arrows. VOA: variable optical attenuator, TIA:
time-interval analyzer. (b) Transmission spectra of the microresonator. The
pump, signal and idler resonances involved in the SFWM process are highlighted
in green, red and blue respectively. The bottom right inset shows an optical
microscope image of the device.

Our device is a silicon nitride microresonator of cross section that is
coupled to a single bus waveguide. Details of sample fabrication and device
geometry are provided in Appendix A. We use a tunable laser to probe the
transmission spectrum of the resonator, which is reported in Fig.2(b). The
free spectral range (FSR) is about $200$ GHz. Two-mode squeezed light is
generated at the signal and idler resonance wavelengths $\lambda_{s}=1552.3$
nm and $\lambda_{i}=1536.7$ nm by resonantly pumping the device at a
wavelength of $\lambda_{p}=1544.5$ nm. The loaded Q of these resonances are
$Q_{s}\sim Q_{i}=2.6\times 10^{5}$ and $Q_{p}=3.3\times 10^{5}$. From the
extinction of the resonances, we estimated an intrinsic quality factor of
$Q_{\textrm{int}}\sim 3.5\times 10^{6}$, and an escape efficiency
$\eta_{e}=(Q_{\textrm{int}}-Q)/Q_{\textrm{int}}\sim 0.926$. The generation of
squeezed light from the microresonator is investigated using the setup
sketched in Fig.2(a). A continuous wave pump laser is passed through an
electro-optic amplitude modulator(EOM) which carves pulses of duration
$\tau=300$ ps at a repetition rate of $2.5$ MHz. After the EOM, the pump is
amplified by an Erbium-Doped Fiber Amplifier (EDFA). We stabilize the pump
power with high precision ($<0.01\%$ rms of relative power fluctuation) by
using a variable optical attenuator and a closed-loop optical feedback
controlled by a field-programmable gate array (FPGA). A fiber polarization
controller is used to set the polarization to TE, which maximizes the SFWM
efficiency. We suppress the laser background noise at the signal and idler
wavelengths with a passband filter. Light is coupled in and out from the chip
using two UHNA4 fibers and an index matching gel for improved coupling
efficiency and mechanical stability ($1.0(5)$ dB loss/facet). The output is
passed through a VOA which controls the overall loss $\eta$ of the system. The
VOA attenuation $\eta_{\textrm{VOA}}$ (expressed net of the insertion loss) is
calibrated with an accuracy below $2\%$ (including the polarization dependent
loss) in the spectral range $[1500-1600]$ nm. The signal and the idler modes
are separated by two Dense Wavelength Division Multiplexing (DWDM) modules
with a bandwidth of $100$ GHz. Each filter is followed by a tunable fiber
Bragg grating with a passband of $12.5$ GHz to reduce the in-band noise.
Photons are detected by two superconducting nanowire single-photon detectors
(SNSPD, from PhotonSpot). At $\eta_{\textrm{VOA}}=1$, the losses from the chip
to the detector are $3.5(5)$ dB for the signal and the idler channels,
including the fiber-to-chip coupling losses, the detector efficiency of $85\%$
and the escape efficiency $\eta_{e}$.

##  IV Model validation and limits of application

Figure 3: (a) Coincidence click probability (black dots) between the signal
and the idler detector as a function of the input pump input power. The red
dashed line is calculated from model in Eq.(9) and fits the data for $P<2.5$
mW. (b) Average number of photons per pulse from squeezing ($\langle
n\rangle$) and from the thermal background in the signal ($\langle
n_{\textrm{s,th}}\rangle$) and the idler mode ($\langle
n_{\textrm{i,th}}\rangle$). Confidence bounds covering one standard deviation
are shown with colored regions. These are calculated from the uncertainty
$\Delta\eta$ in the measurement of the losses, which is $\Delta\eta=0.5$ dB.

The model of the PND in Eq.(9) holds as long as the intra-cavity mean photon
number $\langle n_{s(i)}\rangle$ is much lower than the number of pump
photons, i.e., in the undepleted pump approximation. In our work, this
approximation is well justified since, as it will be shown later, $\langle
n_{s(i)}\rangle\sim 1$. However, the linear relation between the squeezing
parameter $r$ and the input pump power $P$ holds only when SPM and XPM effects
are negligible. This is validated by measuring the joint click probabilities
$p_{11}$ as a function of the average power coupled to the chip. For each
power, we maximize the count rate by adjusting the wavelength of the pump
laser to compensate the SPM and XPM induced shifts on the cold resonance
frequencies. In Fig.3(a) we show the measured coincidence probability $p_{11}$
as a function of the pump power $P$. We fit $p_{11}$ using Eq.(10) with the
PND given by Eq.(9), leaving $a$ and $b_{s(i)}$ as free parameters. We found
that the $R^{2}$ of the fit decreases when click-probabilities at input powers
higher than $2.5$ mW are included in the fit. We then limit the fit to the
range $1$ mW $\leq P\leq 2.5$ mW, and show the result in Fig.3(a). We see that
the model correctly predicts the experimental data for input powers $<2.5$ mW,
while at higher powers it overestimates the generation rate. It is clear that
for $P>3$ mW, the linear relationship between the squeezing parameter and the
input power no longer holds. The behavior is different from that of a CW pump,
in which the SPM and XPM induced resonance shifts can be fully compensated by
adjusting the pump frequency, and the generation rate scales with the square
of the input power [30]. Using the parameters $a$ and $b_{s(i)}$, we
calculated the mean number of photons per pulse from squeezing ($\langle
n\rangle=\sinh^{2}(aP)$) and noise ($\langle
n_{\textrm{s(i),th}}\rangle=b_{s(i)}P$), which are shown in Fig.3(b).

##  V Reconstruction of the photon number distribution

Figure 4: (a) Single ($p_{01}$, $p_{10}$) and joint ($p_{11}$) click
probabilities (dots) for different VOA attenuations $\eta_{\textrm{VOA}}$ at a
pump power of $2.2$ mW. The black line is a fit of the data using the model in
Eq.(9) and Eq.(10). (b) Reconstructed photon number distribution
$\rho_{\textrm{EM}}$ at the output of the resonator for $2.2$ mW of input pump
power, corresponding to an estimated on chip squeezing level of $6.2(2)$ dB.
(c,d) Same as in (b), but reconstructed at the output of the chip (c) and at
the detectors (squeezing level of $2.2(2)$ dB) (d). (e) Comparison between the
experimental PND shown in panel (b) and the simulation with highest fidelity.
(f) Comparison between the experimental marginal photon number distribution of
the signal and the simulated marginal distribution with highest fidelity.

The PND is reconstructed using the method presented in Sec.II. We implemented
$M=50$ measurement settings by varying the VOA transmission
$\eta_{\textrm{VOA}}$ from a minimum of $0.05$ to a maximum of $0.95$ in equal
steps. The total transmission is then
$\eta=\eta_{\textrm{exp}}\eta_{\textrm{VOA}}$, where $\eta_{\textrm{exp}}$ is
transmission $\eta$ measured with the VOA set at the minimum of attenuation.
For each setting $\mu$, we extracted the four click-probabilities
$p_{xy}(\eta_{\mu})$ by collecting $12.5$ millions of events. The
probabilities are input to the EM algorithm and a uniform PND is used as an
initial uninformative ansatz. The algorithm is stopped when the variation in
relative error $\frac{\Delta\epsilon}{\epsilon}$ in Eq.(5) is below $10^{-3}$.
The dimension of the Hilbert space is truncated to a dimension of $10$. In
Fig.4(a) we show the measured click probabilities $p_{xy}$ as a function of
the VOA transmission for an input power of $2.5$ mW. We also show the
predicted click probabilities $p_{xy}(\eta)\\{\rho_{\textrm{EM}}\\}$ given the
PND $\rho_{\textrm{EM}}$ at the last iteration step of the EM. These are in a
very good agreement with the experiment. The reconstructed PND is reported in
Fig.4(b) and shows strong photon number correlations between the signal and
the idler mode, as expected for a TMS. The off-diagonal elements in the PND
matrix $\rho_{\textrm{EM}}$ arise from two contributions. The first is the
presence of uncorrelated noise photons in the signal and idler arms, mainly
due to spontaneous Raman scattering. This represents the dominant
contribution. The second could be due to an imperfect estimation of the losses
$\eta_{\textrm{exp}}$. If the latter are higher than those used to retrieve
$\rho_{\textrm{EM}}$, i.e. $\eta_{\textrm{exp}}^{\prime}<\eta_{\textrm{exp}}$,
then the reconstructed PND will describe the statistics of the state generated
by the (lossless) resonator and subjected to an extra-attenuation of
$\frac{\eta_{\textrm{exp}}^{\prime}}{\eta_{\textrm{exp}}}$. This loss spoils
the perfect photon number correlation in absence of noise.
With our method, we can easily track the changes of the PND from the chip to
the detectors. In Fig.4(c,d) we report respectively the reconstructed PND
immediately out of the chip and at the detectors. The two PNDs are obtained by
substituting
$\eta_{\textrm{exp}}\rightarrow\frac{\eta_{\textrm{exp}}}{\eta_{\textrm{chip}}}$
and $\eta_{\textrm{exp}}\rightarrow 1$ in the EM algorithm, where
$\eta_{\textrm{chip}}=0.8$ is the chip-to-fiber coupling efficiency. Compared
to the on-chip PND, the photon number correlation is gradually lost by adding
losses, and tends to that of the product of two uncorrelated thermal states.
This is consistent with the fact that for a TMS the variance $V_{\Delta n}$ of
the photon number difference $\Delta_{n}=n_{s}-n_{i}$ decreases with the
losses as $V_{\Delta n}=\eta\langle n_{s}+n_{i}\rangle$ [8]. Since the EM
reconstructs $\rho_{\textrm{EM}}$ without relying on any model but solely from
the raw click probabilities, it is important to calculate the fidelity
$\mathcal{F}$ between $\rho_{\textrm{EM}}$ and the physics-driven model
distribution in Eq.(9). We obtained a maximum fidelity of $\mathcal{F}=0.98$
when $r=0.63$, $\langle n_{s}\rangle=0.11$ and $\langle n_{i}\rangle=0.10$.
The comparison between $\rho_{\textrm{EM}}$ and the theoretical distribution
is shown in Fig.4(e). The remarkable agreement between the two indicates that
$\rho_{\textrm{EM}}$ is a reliable reconstruction of the PND of the state, and
that we can use Eq.(9) to infer the squeezing parameter $r$. We also apply the
EM to reconstruct the marginal photon number distribution of the signal(idler)
arm, and evaluated the fidelity with the distribution obtained by
marginalizing Eq.(9) with respect to the idler(signal) photon number. The
comparison between the two PND is reported in Fig.4(f) and indicates a very
good agreement, with a fidelity of $\mathcal{F}=0.999$. As a side remark, we
noticed the emergence of artifacts in the reconstructed PND for input powers
higher than $3.25$ mW. We traced this fact back to an overfitting problem,
which becomes relevant at high powers due to the increased number of not-
vanishing elements in the PND and to the high sensitivity of those elements to
experimental noise.

###  V.1 Calculation of the mean photon number, squeezing parameter and NRF
from the PND

The measurement of the PND allows us to estimate important metrics. The first
one we calculate is the squeezing parameter $r$. We reconstructed the PND for
different input pump powers, and for each we find the values of $r$ and
$\langle n_{\textrm{s(i),th}}\rangle$ that maximize the fidelity with the
distribution in Eq.(9). These are shown in Fig.5(a) for $P<2.5$ mW, for which
Eq.(9) is a faithful model.

Figure 5: (a) Squeezing parameter $r$ extracted from the PND as a function of
the input pump power (black dots). The red line is a linear fit of the data.
(b) Simulated average number of photons per pulse $\langle n\rangle$ for
different input powers. The blue dots consider the effects of SPM and XPM,
while the red dots are obtained with SPM and XPM switched off. The red dashed
line is a fit of the data without SPM and XPM using the function
$\sinh^{2}(aP)$. The black dots are the values of $\langle n\rangle$ extracted
from the PND in the experiment.

The result shows that $r$ scales linearly with $P$, which agrees with the
theory and is consistent with the power scaling of the click probabilities
shown in Fig.3(a), with the important difference that in the former case no
prior relations between the two variable is assumed. We can exploit the PND to
calculate the mean photon number $\langle n\rangle$ per pulse. This quantity
had already been calculated with the model in Eq.(9) (see Fig.3(b)), but the
results were justified only in the low power regime. In Fig.5(b) we show
$\langle n\rangle$ as a function of the input power up to $P=3.25$ mW, i.e.,
the highest power for which no artifacts were observed in the reconstructed
PNDs. To validate the reliability of the reconstructed PND and of the value of
$\langle n\rangle$ in this extended power range, we numerically solved the
stochastic master equation for the signal and idler modes under the action of
the hamiltonian in Eq.(7) and subjected to continuous photodetection of both
the signal and idler modes [31]. The simulation details are reported in
Appendix E. From the simulated PND, we calculated the mean number of
signal(idler) photons scattered into the bus waveguide for different input
powers $P$, which we shown in Fig.5(b). The experimental values of $\langle
n\rangle$ are in very good agreement with the simulation. The full list of the
fidelities between the simulated and the experimental PND is reported in Table
2 of Appendix E. We also calculated $\langle n\rangle$ by switching off SPM
and XPM, and overlaied the results in Fig.5(b). The effect of SPM and XPM is
to decrease the average number of photons generated by SFWM, which explains
why the relation $r=aP$ overestimates the coincidence click probabilities
shown in Fig.3(a).

Figure 6: (a) Variance $V_{\Delta n}$ of the photon number difference between
the signal and the idler arm as a function of the average photon number
$\langle n_{\textrm{tot}}\rangle$. The blue dots refer to a bipartite coherent
state while the black dots to the squeezed light from the resonator. Solid
lines are linear fit of the data. Confidence bounds covering one standard
deviation are shown with colored regions. These are calculated from the
uncertainty $\Delta\eta$ in the measurement of the losses, which is
$\Delta\eta=0.5$ dB. (b) Reconstructed PND for a bipartite coherent state with
an average photon number of $\langle n\rangle=1$.

We noticed that, similarly to the CW case, there is an optimal and power
dependent pump frequency detuning from resonance for which $\langle n\rangle$
is maximized. Intuitively, the nonlinear resonance shifts can not be totally
compensated because the pump detuning is kept fixed, while the SPM and XPM-
induced resonance shifts follow the temporal variations of the pump intensity
in the resonator. As shown in Fig.5(b), the mean number of generated photons
per pulse scales as $\langle n\rangle=\sinh^{2}(aP)$ if SPM and XPM effects
are not cosidered. We used the numerical simulation to calculate the time-
integrated unheralded $\bar{g}_{2}$ (see Appendix E) of the signal/idler
modes, which is connected to the Schmidt number $K$ as $K=g_{2}-1$ [6]. We
found that up to $4$ mW of input power, its value changes by less than $1\%$
from that at low power ($\bar{g}_{2}=1.89$), implying that SPM and XPM impact
the squeezing strength more than the number of Schmidt modes.
As a last step, from the PND we calculated the variance $V_{\Delta n}$ of
difference between the signal and the idler photon numbers, which is given by
$V_{\Delta n}=V_{n_{s}}+V_{n_{i}}-2\big{(}\langle n_{s}n_{i}\rangle-\langle
n_{s}\rangle\langle n_{i}\rangle\big{)}$, where $V_{n_{s(i)}}$ is the variance
of the signal(idler) photon number. For an ideal TMS subjected to a loss
$\eta$, this quantity is equal to $V_{\Delta n}=\eta(\langle
n_{s}\rangle+\langle n_{i})\rangle=\eta\langle n_{\textrm{tot}}\rangle$, while
for a coherent and a thermal state is respectively equal to $V_{\Delta
n}=\langle n_{\textrm{tot}}\rangle$ and $V_{\Delta n}=\langle
n_{\textrm{tot}}\rangle+\langle n_{s}\rangle^{2}+\langle n_{i}\rangle^{2}$ [1,
8]. A reduction of the slope between $V_{\Delta n}$ and $\langle
n_{\textrm{tot}}\rangle$ from unity indicates correlations between the number
of signal and idler photons. The ratio $\frac{V_{\Delta n}}{\langle
n_{\textrm{tot}}\rangle}$ is the NRF, or photon number difference squeezing. A
value below unity is considered a metric of non-classicality of the state
[32]. The variance $V_{\Delta n}$ as a function of $\langle
n_{\textrm{tot}}\rangle$ is calculated from the PND reconstructed immediately
at the output of the lossy resonator, and is reported in Fig.6(a). The slope
of the linear fit is $0.42(2)$, corresponding to a photon number difference
squeezing of $-3.8(2)$ dB. The expected value of the NRF calculated from the
escape efficiency of the resonator is $10\log(1-\eta_{e})\sim-11$ dB. However,
the correlation between the signal/idler photon number difference is reduced
by the background noise. Indeed, by calculating $\frac{V_{\Delta n}}{\langle
n_{\textrm{tot}}\rangle}$ using the the PND in Eq.(9) we obtain
$\frac{V_{\Delta n}}{\langle n_{\textrm{tot}}\rangle}=0.4(1)$ ($-4.0(1)$ dB of
photon number difference squeezing), which is consistent with the NRF value in
Fig.6(a). We complement the analysis by comparing the NRF of the TMS from that
of a coherent state. The latter is obtained by passing the attenuated pulsed
pump laser into a 50/50 fiber beamsplitter and by reconstructing its PND at
different powers. An example of such reconstruction is shown in Fig. 6(b). The
variance $V_{\Delta n}$ as a function of $\langle n_{\textrm{tot}}\rangle$ is
shown in Fig.6(a). The calculated NRF is $0.96(2)$ and lies well above that of
the TMS.

##  VI Conclusions

In this work we reconstructed the PND of pulsed squeezed light generated from
a silicon nitride microresonator without the use of photon number resolving
detectors. The method exploits the nonlinear response of on-off photodetectors
to photon number states undergoing variable attenuation. When SPM and XPM
effects have negligible impact on the rate of photon pairs generated by SFWM,
the reconstructed photon statistics is well described by the convolution of a
TMS and two uncorrelated thermal fields. The latter account for spurious noise
photons generated by the resonator, and mainly arise from spontaneous Raman
scattering. At low power, we found that the squeezing parameter scales
linearly with the pump power. At higher input powers, the squeezing strength
is limited by SPM and XPM. In contrast to the CW case, it is not possible to
fully compensate the nonlinear resonance shifts by adjusting the pump
frequency. The PND are used to calculate the mean photon number per pulse, the
squeezing parameter and the noise reduction factor. The PNDs are characterized
up to $\sim 1.2$ photons/pulse, through which we extracted an on-chip
squeezing level of $6.2(2)$ dB and a noise reduction factor of $-3.8(2)$ dB.
This approach gives reliable results when the number of parameters in the PND
is not too large, which may depend on several factors. The most relevant are
the complexity of the PND, the accuracy and reproducibility of the VOA
settings, the stability of the measuring apparatus over time, the overall
losses and the finite sampling size. In our case, we found that reliable
results are obtained for a Hilbert space dimension $<6\times 6$. We expect
that larger dimensions can be reached by complementing this method with the
spatial and temporal multiplexing strategies used in pseudo-PNR.

### Acknowledgements

D.B acknowledges the support of Italian MUR and the European Union - Next
Generation EU through the PRIN project number F53D23000550006 - SIGNED. E.B.,
M.B., M.B, M.G and M.L. acknowledge the PNRR MUR project PE0000023-NQSTI. All
the authors acknowledge the support of Xanadu Quantum Technology for providing
the samples and the useful discussions.

## Appendix A: Device geometry and fabrication

The device is fabricated on a silicon nitride photonic chip in a standard
multi-project wafer run by Ligentech. All the waveguides have width of 1.75
$\mu$m and a thickness of 800 nm, and are covered by a silicon dioxide
cladding of thickness $4.1\,\mu$m. The gap between the resonator and the bus-
waveguide is 0.4 $\mu$m, and the corresponding coupling coefficient is
$\kappa^{2}=0.016$. The resonator is made by two U-shaped Euler bends of
length $360\,\mu$m, which are designed to minimize the bending loss, thus
maximizing the intrinsic quality factor of the devices. The linear propagation
losses of the Euler have been found comparable to those of the straigth
waveguide sections, which are $\sim 0.11$ dB/cm. A metallic heater of Titanium
Nitride of $0.4\,\mu$m thickness and $2\,\mu$m width is placed $1.7\,\mu$m
above the waveguides. It allows to change the resonance wavelength through the
thermo-optic effect.

Figure 7: Optical microscope image of the silicon nitride microresonator used
in this work. Indicated are the folded TiN heater used to locally heat the
waveguide underneath, the Aluminium routing layer from the TiN heaters to the
metal pads, the input and the output optical ports.

##  Appendix B: reconstruction of single mode coherent and thermal states

We first validated the method presented in Sec.2 by reconstructing some known
optical states, which we chose to be a single mode coherent state and a
thermal state. To highlight that the measurement is sensitive to the full PND
and not just to its first order moment, i.e. the average photon number
$\langle n\rangle$, we aimed to reconstruct states with the same $\langle
n\rangle\sim 1.2$. The coherent state is obtained by attenuating the pump
laser, while the second from the unheralded signal beam generated by SFWM in
the microresonator. Figure 8(a) reports the probability of missed detections
$p_{0}(\eta)$ of the two states as a function of the variable detection
efficiency $\eta$. We notice that up to $\eta\sim 0.15$, the click
probabilities are almost identical, which agrees with the fact that for
$\eta_{\textrm{max}}\ll 1$ the probability of an missed detection is
$p_{0}=\eta\langle n\rangle$. As $\eta$ increases, the photodetection
statistics of the two states become distinguishable. This can be appreciated
also from the reconstructed PND of the two states, which are shown in
Fig.8(b). The experimental PNDs are then compared to their theoretical
distributions, obtaining a fidelity of $\mathcal{F}=0.997$ for the thermal
state, and $\mathcal{F}=0.998$ for the coherent state.

##  Appendix C: The bipartite case

Figure 8: (a) Measured off-detection probabilities as a function of the system
loss $\eta$ for a thermal and a coherent state with the same average photon
number of $\langle n\rangle\sim 1.2$. (b) Reconstructed PND of the two states.

The method presented in Sec.2 generalizes to the bipartite case, in which a
single detector is used to sense each mode. The four possible photodetection
events $p_{xy}$ occur with probabilities

$\begin{split}&p_{00}(\eta)=\sum_{n,k}B_{\eta n}B_{\eta k}\rho_{nk},\\\
&p_{10}(\eta)=\sum_{n,k}(1-B_{\eta n})B_{\eta k}\rho_{nk},\\\
&p_{01}(\eta)=\sum_{n,k}B_{\eta n}(1-B_{\eta k})\rho_{nk},\\\
&p_{11}(\eta)=1-p_{00}(\eta)-p_{10}(\eta)-p_{01}(\eta),\end{split}$ (10)

in which the diagonal elements of the two-mode density matrix $\hat{\rho}$ are

$\rho_{nk}=\bra{nk}\hat{\rho}\ket{nk},\quad\ket{nk}=\ket{n}\otimes\ket{k}.$
(11)

The indices $(n,k)$ in Eq.(10) range from $0$ to $N$, where $N$ is the
dimension of the truncated Hilbert space. We now define the vectorized PND
$\bm{\rho}$ as

$\bm{\rho}=(\rho_{00},\rho_{01},\dots,\rho_{0N},\rho_{10},\rho_{11},\dots\rho_{N0},\rho_{N1},\dots\rho_{NN})^{T},$
(12)

and the vector $\bm{p}$ collecting the $3M$ click-probabilities from the $M$
measurement settings as

$\bm{P}=(p_{00}^{\eta_{1}},p_{00}^{\eta_{2}},\dots,p_{00}^{\eta_{M}},p_{01}^{\eta_{1}},p_{01}^{\eta_{2}},\dots,p_{01}^{\eta_{M}},p_{10}^{\eta_{1}},p_{10}^{\eta_{2}},\dots
p_{10}^{\eta_{M}})^{T}.$ (13)

The system of equations in Eq.(10) can be written in compact form as
$\bm{p}=B\bm{\rho}$, in which the $3M\times(N+1)^{2}$ matrix $B$ has elements

$B=\begin{cases}B_{\mu n}B_{\mu k}\qquad\quad\,\,\,\mu=1,\dots,M\\\ B_{\mu
n}(1-B_{\mu k})\quad\mu=M+1,\dots,2M\\\ (1-B_{\mu n})B_{\mu
k}\quad\mu=2M+1,\dots,3M\end{cases}$ (14)

The EM step at the iteration $i$ becomes

$\rho_{p}^{(i+1)}=\rho_{p}^{(i)}\sum_{\mu=1}^{3M}\frac{B_{\mu
p}f_{\mu}}{(\sum_{\nu}B_{\nu p})p_{\mu}\\{\bm{\rho^{(i)}}\\}},$ (15)

where we have defined the experimental frequencies as

$f_{\mu}=\begin{cases}f_{00}=n_{00\mu}/n_{\mu}\quad\mu=1,\dots,M\\\
f_{01}=n_{01\mu}/n_{\mu}\quad\mu=M+1,\dots,2M\\\
f_{10}=n_{10\mu}/n_{\mu}\quad\mu=2M+1,\dots,3M\\\ \end{cases}$ (16)

where $n_{\mu}=R_{\mu}T_{\mu}$ is the total number of events, which is given
by the repetition rate of the pump laser $R_{\mu}$ and the integration time
$T_{\mu}$ of the $\mu$ measurement setting, while $n_{xy\nu}$ is the number of
recorded events. In Eq.(15), we implemented the following mapping of the
indices

$\rho_{nk}\rightarrow,\rho_{p}\qquad p=1+k+n(1+N),$ (17)

##  Appendix D: Relation between the $\xi$ coefficients in the hamiltonian

The function $\xi^{(m,n)}(\omega_{s},\omega_{i},t)$ in Eq.(7) of the main
manuscript is equal to [33]

$\xi^{(m,n)}(\omega_{s},\omega_{i},t)=\mathcal{N}_{p}(\omega_{s},\omega_{i},t)\phi_{s,m}(\omega_{s})\phi_{i,n}(\omega_{i}),$
(18)

in which

$\begin{gathered}\mathcal{N}_{p}(\omega_{s},\omega_{i},t)=\mathcal{N}^{\prime}\int\phi_{p}(\omega_{p_{1}})\phi_{p}(\omega_{p_{2}})\\\
\beta(\omega_{p_{1}})\beta(\omega_{p_{2}})e^{i\Delta\omega
t}d\omega_{p_{1}}d\omega_{p_{2}},\end{gathered}$ (19)

and $\mathcal{N}^{\prime}$ is a constant. In Eqs.(18,19), $\phi_{p/s/i,m}$ are
the field enhancements near the pump, signal and idler frequencies when the
light is input to port $m=\\{\textrm{wg,ph}\\}$, $\beta$ is the spectral
profile of the pump and
$\Delta\omega=\omega_{s}+\omega_{i}-\omega_{p_{1}}-\omega_{p_{2}}$. Equation
(19) assumes that the relation between the incident pump field and the (slowly
varying) internal field $E_{p}$ (normalized such that
$[|E_{p}|^{2}]=\textrm{Watts}$) in the resonator is linear, i.e,
$E_{p}\propto\beta\phi_{p}$. When SPM effects are not negligible, we have to
replace $\beta\phi_{p}$ with the Fourier transform of the internal field
$E_{p}(t)$. The latter satisfies the differential equation

$\frac{dE_{p}}{dt}=[-\gamma_{\textrm{tot,p}}+i(\Lambda|E_{p}|^{2}+\Delta_{p})]E_{p}+i\sqrt{\frac{2\gamma_{ep}}{\tau_{\textrm{rt}}}}\beta(t),$
(20)

in which $\gamma_{\textrm{tot,p}}$ is the total decay rate of the field,
$\gamma_{ep}$ is the coupling rate between the resonator and the bus
waveguide, $\Delta_{p}=\omega_{c}-\omega_{p0}$ is the detuning of the central
frequency $\omega_{c}$ of the pump laser from the cold resonance frequency of
the pump $\omega_{p0}$ and $\tau_{\textrm{rt}}$ is the round-trip time of
light in the resonator. The nonlinear coefficient $\Lambda$ accounts for SPM
and is defined as

$\Lambda=\frac{\omega_{p}cn_{2}}{V_{\textrm{eff}}n_{0}^{2}}\tau_{\textrm{rt}},$
(21)

where $n_{2}$ is the nonlinear refractive index of silicon nitride, $n_{0}$ is
the refractive index of the core material and $V_{\textrm{eff}}$ is the
effective volume of the resonator. We can obtain $\phi_{p}$ in the low power
regime by eliminating the SPM term in Eq.(20). This gives

$\phi_{p}(\omega)=\frac{i\sqrt{\frac{2\gamma_{ep}}{\tau_{\textrm{rt}}}}}{\gamma_{\textrm{tot}}-i(\Delta_{p}+\omega)}.$
(22)

Similarly, $\phi_{s(i),m(n)}$ are given by

$\phi_{s(i),m(n)}(\omega)=\frac{i\sqrt{\frac{2\gamma_{s(i)}^{m(n)}}{\tau_{\textrm{rt}}}}}{\gamma_{\textrm{tot,s(i)}}-i(\omega_{s(i)}-\omega_{s(i)0})},$
(23)

where $\gamma_{s(i)}^{m}=\gamma_{es(i)}$ if $m=\textrm{wg}$ and
$\gamma_{s(i)}^{m}=\gamma_{i}$ if $m=\textrm{ph}$. From the definition in
Eq.(23) we have that
$\phi_{\textrm{s(i),wg}}=\sqrt{\frac{\gamma_{es(i)}}{\gamma_{i}}}\phi_{\textrm{s(i),ph}}$,
which combined with Eq.(18) gives the relations in Eq.(8) in the main text,
where $\gamma_{es}=\gamma_{ei}=\gamma_{e}$ has been used for simplicity.
The XPM induced by the pump to the signal and idler fields is described by the
term

$\Delta_{s(i)}^{(m,n)}=2\Lambda|E_{p}(t)|^{2}\phi_{s(i),m}(\omega_{s(i)})\phi_{s(i),n}^{*}(\omega_{s(i)}^{\prime})e^{i\delta_{s(i)}t},$
(24)

where $\delta_{s(i)}=\omega_{s(i)}-\omega_{s(i)}^{\prime}$. Armed with these
tools, we can show that the hamiltonian in Eq.(7) of the main text is
equivalent to that of a lossless resonator coupled to a single bus waveguide
with a rate $\gamma_{e}+\gamma_{i}$, in which the outgoing photons impinge on
a beamsplitter with transmittivity
$\eta_{e}=\frac{\gamma_{e}}{\gamma_{e}+\gamma_{i}}$. The nonlinear hamiltonian
of the lossless resonator, in the interaction picture, is given by

$\begin{gathered}H_{\textup{NL}}(t)=\hbar\left\\{\int\Delta_{s}^{\textup{(wg,wg)}}(\omega_{s},\omega_{s}^{\prime},t)a_{\textup{wg}}^{\dagger}(\omega_{s})a_{\textup{wg}}(\omega_{s}^{\prime})d\omega_{s}d\omega_{s}^{\prime}+\right.\\\
+\int\Delta_{i}^{\textup{(wg,wg)}}(\omega_{i},\omega_{i}^{\prime},t)a_{\textup{wg}}^{\dagger}(\omega_{i})a_{\textup{wg}}(\omega_{i}^{\prime})d\omega_{i}d\omega_{i}^{\prime}+\\\
+\left.\left(\int\xi^{\textup{(wg,wg)}}(\omega_{s},\omega_{i},t)a_{\textup{wg}}^{\dagger}(\omega_{s})a_{\textup{wg}}^{\dagger}(\omega_{i})d\omega_{s}d\omega_{i}+\textrm{H.c.}\right)\right\\},\end{gathered}$
(25)

because once they are generated by SFWM, photons can only leave resonator from
the bus waveguide. Those photons impinge on a beasmplitter with transmittivity
$\eta_{e}$, which realize the transformation
$a_{\textup{wg}}\rightarrow\eta_{e}a_{\textup{wg}}+\sqrt{1-\eta_{e}}a_{\textup{ph}}$,
where $a^{\dagger}_{\textup{ph}}$ creates a photon in the phantom loss channel
in Fig.1(b) of the main manuscript. By inserting this operator transformation
into Eq.(25) and by using Eq.(8) we obtain the hamiltonian in Eq.(7) of the
main manuscript, concluding our proof.

##  Appendix E: Stochastic master equation for the signal and idler fields

The use of asymptotic fields to represent the Hamiltonian in Eq.(7) of the
main text is useful when SPM and XPM effects are negligible. In this case, the
time dependence of the hamiltonian is solely contained in the frequency
mismatch term in Eq.(19), which is approximated by
$\delta(\omega_{s}+\omega_{i}-\omega_{p_{1}}-\omega_{p_{2}})$ upon time
integration [26]. As a consequence, the frequency domain description of the
generated quantum state is naturally provided by the continuum of asymptotic
states. On the contrary, in the high gain regime, it is more practical to
solve the system dynamics entirely in the time domain [26, 30, 24]. We use
this approach to write the master equation for the signal/idler density
matrix. Following [26], we write the full hamiltonian of the resonator coupled
to the bus waveguide and the _phantom_ loss channel as
$H=H_{\textrm{wg}}+H_{\textrm{ph}}+H_{\textrm{r}}+H_{\textrm{c}}+H_{\textrm{NL}}$,
where the first three terms describe the free-evolution of the waveguide-
phantom channels $a_{\textrm{wg-ph}}(k)$ (forming a continuum) and of the
discrete resonator modes $a_{s(i)}$. The cavity-channel coupling is given by

$H_{c}=\hbar\sqrt{v_{g}}\sum_{J=(i,s)}\left(\sqrt{2\gamma_{eJ}}\Psi_{\textrm{wg,J}}(0)+\sqrt{2\gamma_{i}}\Psi_{\textrm{ph,J}}(0)\right)a_{J}^{\dagger}+\textrm{H.c.},$
(26)

where $v_{g}$ is the group velocity and

$\Psi_{\textrm{wg,ph}}(z)=\int\frac{dk}{\sqrt{2\pi}}a_{\textrm{wg,ph}}(k)e^{i(k-k_{J})z}.$
(27)

The integration is performed in the neighborhood of the wavevector
$k_{J=(s,i)}$ of the signal(idler) resonance. The nonlinear hamiltonian
$H_{\textrm{NL}}$ describes the effects of XPM and SFWM, and is given by

$H_{\textrm{NL}}=-[\bar{\Lambda}(\langle
a_{p}(t)\rangle^{*})^{2}a_{s}^{\dagger}a_{i}^{\dagger}+\textrm{h.c.}]-2\bar{\Lambda}|\langle
a_{p}(t)\rangle|^{2}(a_{s}^{\dagger}a_{s}+a_{i}^{\dagger}a_{i})$ (28)

where

$\bar{\Lambda}=\frac{\hbar\omega_{p}^{2}cn_{2}}{n_{0}^{2}V_{\textrm{eff}}}$
(29)

and $\langle a_{p}(t)\rangle$ is the solution of the differential equation

Figure 9: Comparison between the experimental PND of the TMS at the output of
the resonator and the one predicted by the stochastic master equation Eq.(32).
Panels (a) and (b) refer respectively to an input pump power of $1$ mW and
$2.5$ mW.

$\begin{gathered}\left[\frac{d}{dt}+\gamma_{\textrm{tot}}-2i\bar{\Lambda}|\langle
a_{p}\rangle|^{2}\right]\langle
a_{p}\rangle=-i\sqrt{2\gamma_{ep}}\beta(t)e^{i\Delta_{p}t}\end{gathered}$ (30)

In writing Eq.(28) we replaced the pump operators $a_{p}$ and
$a_{p}^{\dagger}$ with the complex numbers $a_{p}\rightarrow\langle
a_{p}\rangle$ and $a_{p}^{\dagger}\rightarrow\langle a_{p}\rangle^{*}$, which
is the usual mean field approximation for a strong pump in a coherent state
[30]. We also assume that the mean number of pump photons in the cavity
$\int|\langle a_{p}(t)\rangle|^{2}dt\ll\langle
c_{s(i)}^{\dagger}c_{s(i)}\rangle$, which justifies the absence of pump
depletion terms in Eq.(30). Within these approximations, the pump equation is
self-contained. We choose the driving field $\beta(t)$ to be a top-hat
function of duration $T$

$\beta(t)=\sqrt{\frac{P}{R\hbar\omega_{p}}}\left(\Theta(t)-\Theta(t-T)\right),$
(31)

where $P$ is the average power and $\Theta(t)$ is a step function defined as
$\Theta(t)=1$ for $(t>=0)$ and $\Theta(t)=0$ for $(t<0)$.

Parameter | Value
---|---
$\gamma_{\textrm{tot,i}}$ | $2.25\times 10^{-3}\,\textrm{ps}^{-1}$
$\gamma_{\textrm{tot,p}}$ | $1.84\times 10^{-3}\,\textrm{ps}^{-1}$
$\gamma_{\textrm{tot,s}}$ | $2.25\times 10^{-3}\,\textrm{ps}^{-1}$
$\eta_{es}$ | $0.926$
$\eta_{ei}$ | $0.926$
$\tau_{\textrm{rt}}$ |  $5$ ps
$n_{g}$ | $2.09$
$\bar{\Lambda}$ | $1.72\,\textrm{s}^{-1}$
$\Lambda$ | $6.72\times 10^{7}\,\textrm{J}^{-1}$
Table 1:  List of the parameters used in the simulation of the master equation
of the signal and idler modes.

The master equation for the signal and idler density matrix is obtained by
tracing out the continuum of modes of the waveguide and the phantom channel.
Due to the colorless model of the minimal coupling between the ring and the
channel modes in Eq.(26), we can use the standard second-order Born-Markov
approximation to write the master equation for the signal and the idler modes
$a_{s(i)}$ in Lindblad form

$\dot{\rho}=-\frac{i}{\hbar}[H_{\textrm{NL}},\rho]+\left(\mathcal{D}[\sqrt{2\gamma_{\textrm{tot,s}}}a_{s}]+\mathcal{D}[\sqrt{2\gamma_{\textrm{tot,i}}}a_{i}]\right)\rho,$
(32)

where $\mathcal{D}[a]=a\rho a^{\dagger}-\frac{1}{2}(a^{\dagger}a\rho+\rho
a^{\dagger}a)$. We numerically solve Eq.(32) using the open access Python
library Qutip [34]. The Fock space dimension for the signal and the idler
photon is truncated to $12$ excitations, without any noticeable change as a
result of moving into a larger space. The list of simulation parameters is
shown in Table 1. The density matrix $\rho(t)$ is used to calculate the
average photon number in the cavity
$\textrm{Tr}(\rho(t)a_{s(i)}^{\dagger}a_{s(i)})$, which is related to the
number of scattered photons in the channel waveguide per pulse $\langle
n_{s(i)}\rangle$ as

$\langle
n_{s(i)}\rangle=2\gamma_{es(i)}\int\textrm{Tr}(\hat{\rho}(t)a_{s(i)}^{\dagger}a_{s(i)})dt$
(33)

The Schmidt number $K$ can be inferred from the time-integrated $\bar{g}_{2}$

$\bar{g}_{2}=\frac{\int\langle
a^{\dagger}(t)a^{\dagger}(t+\tau)a(t+\tau)a(t)\rangle
dtd\tau}{\left(\int\langle a^{\dagger}(t)a(t)\rangle dt\right)^{2}}$ (34)

as $K=\bar{g}_{2}-1$ [6]. We calculated the numerator using the Qutip method
`qutip.correlation_3op_2t`, which calculates the two-time correlation function
by applying the quantum regression theorem, while the denominator is
calcuilated using Eq.(33).
The PND is simulated using a stochastic master equation describing the signal
and the idler modes subjected to continuous photodetection [31]. Since only
those photons which are scattered into the waveguide can be detected, we first
split the dissipation term in Eq.(32) as

$\mathcal{D}\left[2\sqrt{\eta_{es(i)}\gamma_{\textrm{tot,s(i)}}}c_{s(i)}\right]+\mathcal{D}\left[\sqrt{2(1-\eta_{es(i)})\gamma_{\textrm{tot,s(i)}}}c_{s(i)}\right]$
(35)

and let the first and the second term to govern respectively the stochastic
and the deterministic part of the evolution of each quantum trajectory [31,
35]. The stochastic master equation is solved by using the method
`qutip.stochastic.photocurrent_mesolve` [34]. To calculate the PND, we
simulated $N_{\textrm{tot}}=3000$ trajectories and for each of them we
recorded the number of joint signal/idler detections $N_{n_{s},n_{i}}$. The
dataset is then binned on a two-dimensional grid representing the signal and
the idler photon numbers $n_{s}$ and $n_{i}$, and the probabilities are
approximated as $p(n_{s},n_{i})\sim N_{n_{s},n_{i}}/N_{\textrm{tot}}$. Figure
9 shows some examples of simulated PND at different pump powers compared to
experimental ones. The fidelities between the simulation and the experiment
are reported in Table 2.

Pump power (mW) | $\mathcal{F}$
---|---
$1.00$ | $0.985$
$1.26$ | $0.984$
$1.41$ | $0.986$
$1.58$ | $0.982$
$1.78$ | $0.981$
$2.00$ | $0.978$
$2.23$ | $0.976$
$2.51$ | $0.956$
$2.82$ | $0.941.$
Table 2:  Fidelity between the simulated PND and the experimental ones for
different pump powers.

## References

  * Loudon [2000] R. Loudon,  _The quantum theory of light_ (OUP Oxford, 2000).
  * Harder  _et al._ [2014] G. Harder, C. Silberhorn, J. Rehacek, Z. Hradil, L. Motka, B. Stoklasa, and L. L. Sánchez-Soto, Time-multiplexed measurements of nonclassical light at telecom wavelengths, Physical review A 90, 042105 (2014).
  * Laiho  _et al._ [2022] K. Laiho, T. Dirmeier, M. Schmidt, S. Reitzenstein, and C. Marquardt, Measuring higher-order photon correlations of faint quantum light: a short review, Physics Letters A 435, 128059 (2022).
  * Shan  _et al._ [2021] Y.-G. Shan, G.-J. Fan-Yuan, S. Wang, W. Chen, D.-Y. He, Z.-Q. Yin, G.-C. Guo, and Z.-F. Han, Characterizing photon-number distribution of pulsed sources using non–photon-number-resolving detectors for quantum key distribution, Physical Review A 104, 032406 (2021).
  * Dynes  _et al._ [2018] J. Dynes, M. Lucamarini, K. Patel, A. Sharpe, M. Ward, Z. Yuan, and A. Shields, Testing the photon-number statistics of a quantum key distribution light source, Optics Express 26, 22733 (2018).
  * Christ  _et al._ [2011] A. Christ, K. Laiho, A. Eckstein, K. N. Cassemiro, and C. Silberhorn, Probing multimode squeezing with correlation functions, New Journal of Physics 13, 033027 (2011).
  * Shchesnovich [2014] V. Shchesnovich, Sufficient condition for the mode mismatch of single photons for scalability of the boson-sampling computer, Physical Review A 89, 022333 (2014).
  * Vaidya  _et al._ [2020] V. D. Vaidya, B. Morrison, L. Helt, R. Shahrokshahi, D. Mahler, M. Collins, K. Tan, J. Lavoie, A. Repingon, M. Menotti, _et al._ , Broadband quadrature-squeezed vacuum and nonclassical photon number correlations from a nanophotonic device, Science advances 6, eaba9186 (2020).
  * Zhang  _et al._ [2021] Y. Zhang, M. Menotti, K. Tan, V. Vaidya, D. Mahler, L. Helt, L. Zatti, M. Liscidini, B. Morrison, and Z. Vernon, Squeezed light from a nanophotonic molecule, Nature communications 12, 2233 (2021).
  * Zhao  _et al._ [2020] Y. Zhao, Y. Okawachi, J. K. Jang, X. Ji, M. Lipson, and A. L. Gaeta, Near-degenerate quadrature-squeezed vacuum generation on a silicon-nitride chip, Physical Review Letters 124, 193601 (2020).
  * Schmidt  _et al._ [2018] M. Schmidt, M. Von Helversen, M. López, F. Gericke, E. Schlottmann, T. Heindel, S. Kück, S. Reitzenstein, and J. Beyer, Photon-number-resolving transition-edge sensors for the metrology of quantum light sources, Journal of Low Temperature Physics 193, 1243 (2018).
  * Arrazola  _et al._ [2021] J. M. Arrazola, V. Bergholm, K. Brádler, T. R. Bromley, M. J. Collins, I. Dhand, A. Fumagalli, T. Gerrits, A. Goussev, L. G. Helt, _et al._ , Quantum circuits with many photons on a programmable nanophotonic chip, Nature 591, 54 (2021).
  * Madsen  _et al._ [2022] L. S. Madsen, F. Laudenbach, M. F. Askarani, F. Rortais, T. Vincent, J. F. Bulmer, F. M. Miatto, L. Neuhaus, L. G. Helt, M. J. Collins, _et al._ , Quantum computational advantage with a programmable photonic processor, Nature 606, 75 (2022).
  * Cahall  _et al._ [2017] C. Cahall, K. L. Nicolich, N. T. Islam, G. P. Lafyatis, A. J. Miller, D. J. Gauthier, and J. Kim, Multi-photon detection using a conventional superconducting nanowire single-photon detector, Optica 4, 1534 (2017).
  * Banaszek and Walmsley [2003] K. Banaszek and I. A. Walmsley, Photon counting with a loop detector, Optics letters 28, 52 (2003).
  * Piacentini  _et al._ [2015] F. Piacentini, M. Levi, A. Avella, M. López, S. Kück, S. Polyakov, I. P. Degiovanni, G. Brida, and M. Genovese, Positive operator-valued measure reconstruction of a beam-splitter tree-based photon-number-resolving detector, Optics letters 40, 1548 (2015).
  * Banner  _et al._ [2024] P. R. Banner, D. Kurdak, Y. Li, A. Migdall, J. V. Porto, and S. L. Rolston, Number-state reconstruction with a single single-photon avalanche detector, Optica Quantum 2, 110 (2024).
  * Mogilevtsev [1999] D. Mogilevtsev, Reconstruction of quantum states with binary detectors, acta physica slovaca 49, 743 (1999).
  * Vardi and Lee [1993] Y. Vardi and D. Lee, From image deblurring to optimal investments: Maximum likelihood solutions for positive linear inverse problems, Journal of the Royal Statistical Society Series B: Statistical Methodology 55, 569 (1993).
  * Zambra  _et al._ [2005] G. Zambra, A. Andreoni, M. Bondani, M. Gramegna, M. Genovese, G. Brida, A. Rossi, and M. G. Paris, Experimental reconstruction of photon statistics without photon counting, Physical review letters 95, 063602 (2005).
  * Brida  _et al._ [2011] G. Brida, M. Genovese, M. Gramegna, A. Meda, F. Piacentini, P. Traina, E. Predazzi, S. Olivares, and M. G. Paris, Quantum state reconstruction using binary data from on/off photodetection, Advanced Science Letters 4, 1 (2011).
  * Walls and Harvey [2012] D. F. Walls and J. D. Harvey,  _Quantum Optics VI: Proceedings of the Sixth International Symposium on Quantum Optics, Rotorua, New Zealand, January 24–28, 1994_ , Vol. 77 (Springer Science & Business Media, 2012).
  * Vernon and Sipe [2015a] Z. Vernon and J. E. Sipe, Spontaneous four-wave mixing in lossy microring resonators, Physical Review A 91, 053802 (2015a).
  * Vernon  _et al._ [2019] Z. Vernon, N. Quesada, M. Liscidini, B. Morrison, M. Menotti, K. Tan, and J. Sipe, Scalable squeezed-light source for continuous-variable quantum sampling, Physical Review Applied 12, 064024 (2019).
  * Liscidini  _et al._ [2012] M. Liscidini, L. Helt, and J. Sipe, Asymptotic fields for a hamiltonian treatment of nonlinear electromagnetic phenomena, Physical Review A 85, 013833 (2012).
  * Quesada  _et al._ [2022] N. Quesada, L. Helt, M. Menotti, M. Liscidini, and J. Sipe, Beyond photon pairs—nonlinear quantum photonics in the high-gain regime: a tutorial, Advances in Optics and Photonics 14, 291 (2022).
  * Vernon  _et al._ [2017] Z. Vernon, M. Menotti, C. Tison, J. Steidle, M. Fanto, P. Thomas, S. Preble, A. Smith, P. Alsing, M. Liscidini, _et al._ , Truly unentangled photon pairs without spectral filtering, Optics letters 42, 3638 (2017).
  * Triginer  _et al._ [2020] G. Triginer, M. D. Vidrighin, N. Quesada, A. Eckstein, M. Moore, W. S. Kolthammer, J. Sipe, and I. A. Walmsley, Understanding high-gain twin-beam sources using cascaded stimulated emission, Physical Review X 10, 031063 (2020).
  * Vendromin  _et al._ [2024] C. Vendromin, Y. Liu, Z. Yang, and J. E. Sipe, Highly squeezed states in ring resonators: Beyond the undepleted pump approximation, arXiv preprint arXiv:2404.15563 (2024).
  * Vernon and Sipe [2015b] Z. Vernon and J. Sipe, Strongly driven nonlinear quantum optics in microring resonators, Physical Review A 92, 033840 (2015b).
  * Wiseman and Milburn [2009] H. M. Wiseman and G. J. Milburn,  _Quantum measurement and control_ (Cambridge university press, 2009).
  * Aytür and Kumar [1990] O. Aytür and P. Kumar, Pulsed twin beams of light, Physical review letters 65, 1551 (1990).
  * Borghi [2020] M. Borghi, Phase-resolved joint spectra tomography of a ring resonator photon pair source using a silicon photonic chip, Optics Express 28, 7442 (2020).
  * Johansson  _et al._ [2012] J. R. Johansson, P. D. Nation, and F. Nori, Qutip: An open-source python framework for the dynamics of open quantum systems, Computer Physics Communications 183, 1760 (2012).
  * Onodera  _et al._ [2022] T. Onodera, E. Ng, C. Gustin, N. Lörch, A. Yamamura, R. Hamerly, P. L. McMahon, A. Marandi, and H. Mabuchi, Nonlinear quantum behavior of ultrashort-pulse optical parametric oscillators, Physical Review A 105, 033508 (2022).