Also at ]Fraunhofer Institute for Applied Optics and Precision Engineering,
Albert-Einstein-Straße 7, 07745 Jena, Germany

#  Non-perturbative theory of spontaneous parametric down-conversion in open
and dispersive optical systems

Aleksa Krstić<EMAIL_ADDRESS>Institute of Applied Physics, Abbe
Center of Photonics, Friedrich-Schiller University Jena, Albert-Einstein-
Straße 15, 07745 Jena, Germany  Frank Setzpfandt  [  Institute of Applied
Physics, Abbe Center of Photonics, Friedrich-Schiller University Jena, Albert-
Einstein-Straße 15, 07745 Jena, Germany  Sina Saravi  Institute of Applied
Physics, Abbe Center of Photonics, Friedrich-Schiller University Jena, Albert-
Einstein-Straße 15, 07745 Jena, Germany

###### Abstract

We develop a non-perturbative formulation based on the Green-function
quantisation method, that can describe spontaneous parametric down-conversion
in the high-gain regime in nonlinear optical structures with arbitrary amount
of loss and dispersion. This formalism opens the way for description and
design of arbitrary complex and/or open nanostructured nonlinear optical
systems in quantum technology applications, such as squeezed-light generation,
nonlinearity-based quantum sensing, and hybrid quantum systems mediated by
nonlinear interactions. As an example case, we numerically investigate the
scenario of integrated quantum spectroscopy with undetected photons, in the
high-gain regime, and uncover novel gain-dependent effects in the performance
of the system.

††preprint: APS/123-QED

##  I Introduction

With the ever growing importance of optical quantum technologies, new ways to
leverage nonclassical properties of light are in high demand [1, 2]. Some of
the most prevalent methods for generating non-classical light are based on
nonlinear sources of photon pairs in which nonlinear optical processes, such
as spontaneous parametric down-conversion (SPDC) or spontaneous four-wave
mixing (SFWM), are used to convert input light into pairs of photons which
exhibit non-classical correlations [3, 4], highly relevant for many
applications in quantum technologies. Such sources also offer the advantage
that they can be implemented in integrated platforms [2, 4], as well as
utilised in hybrid quantum photonic systems, to overcome the limitations of
monolithic integrated systems [5].

Apart from generating entangled photon pairs, nonlinear photon-pair sources
are also commonly used as a means to obtain heralded single photons [6, 7, 8].
In both cases the sources are operated in the so-called low-gain regime, where
the dominant contribution to the output state (apart from the vacuum) is a
single photon-pair state [7]. This is achieved by keeping the input pump beam
of the nonlinear source at a power that is low enough to keep the multiple
photon-pair contributions in the output state negligible. Although such low-
gain probabilistic sources of photon pairs are of importance in many areas of
quantum technologies, their development in the high-gain regime is also of
interest [9]. In the high-gain regime, the SPDC and SFWM sources can be used
as sources of squeezed light, where squeezed light has wide applications in
continuous variable (CV) quantum computation [10, 11, 12, 13], quantum
communication [14], as well as quantum sensing [15, 16, 17]. Another
application of high-gain operation of such sources can be for generating
multi-photon Fock states [18, 19], which have applications in many branches of
quantum technology as well as metrology and fundamental tests of quantum
entanglement [20, 21, 22].

The rapidly growing interest in the utilisation of high-gain nonlinear sources
of quantum light in nanostructured and hybrid platforms [13, 23, 24, 25, 26,
27, 28] necessitates the development of new theoretical formalisms to describe
high-gain SPDC and SFWM, as such systems can be generally open/lossy with
highly complex spatial and spectral properties [29, 30]. Formalisms for
describing SPDC and SFWM in the high-gain regime in closed systems are well-
established, but either neglect loss altogether [31, 32, 33] or are only able
to account for weak losses that have negligible effect on the modal structure
of the system [34, 9], but are not capable of treating systems where loss is
an inherent part of the system’s optical properties, such as nanoresonators
[35], plasmonic systems [30], or structures with inherently leaky guided modes
[36, 37], or strongly lossy systems that can be encountered in quantum sensing
applications [38]. Finally, many established high-gain formalisms use a weak-
dispersion approximation for the involved modes or rely on approximate
expansions to accommodate higher-order dispersions [32, 9, 39], and do not
offer a straightforward path to inclusion of complex dispersion relations,
such as those that can appear in photonic crystals with internal loss or gain,
where the combination of loss/gain and evanescent modes create complex
dispersion diagrams [40, 41].

The Green’s function (GF) quantisation method [42, 43, 44] offers a way to
describe light quantisation in optical systems with arbitrary dispersion,
loss, and nanostructuring, as the classical GF of the system incorporated into
this quantisation approach naturally takes all these effects into account. In
contrast to closed systems, where the normal modes of the field are quantised
to obtain photon creation/annihilation operators, the GF method quantises the
local bosonic excitations of the combined field and medium system, which
allows it to naturally describe any form of loss as a coupling between
different field-matter modes. A detailed overview of the GF quantisation
procedure can be found in Ref. [45]. The GF formalism has already been used to
describe photon-pair generation in plasmonic structures [46] and dielectric
nanoresonators and nanoparticles [47, 48, 49], as well as in the presence of
quantum emitters and significant loss [50, 38, 51], but in all these cases, it
has only been used do describe pair generation in the low-gain regime.

In this work, we present a formalism, based on the GF quantisation method, to
describe high-gain SPDC in arbitrary open and dispersive optical systems,
within the undepleted pump approximation. Our developed non-perturbative
formalism allows the exact calculation of the field operators, while
intrinsically taking into account arbitrary dispersion and losses of
nanostructured systems through the classical electromagnetic GF of the optical
system. This method can be of special interest for the description of
engineered squeezed-light generation in nanostructured systems in CV quantum
computing applications [13] or systems where losses cannot be treated
perturbatively (e.g. nonlinear metasurfaces for quantum light generation [52,
53, 54]), as well as high-gain quantum sensing and imaging applications, where
loss is not necessarily a weak effect [38, 51, 15, 16]. Additionally, our
formalism may open the path for the high-gain description of hybrid nonlinear
systems, which involve direct interfacing of quantum emitters with nonlinear
systems [50].

The paper is organised as follows: In Sec. II we introduce the theoretical
framework for describing SPDC in arbitrarily dispersive and lossy systems
using fields quantised using the GF quantisation method. We derive coupled-
mode equations that allow the calculation of the full electric field operator
at arbitrary times and amounts of gain (within the undepleted pump
approximation). Then, we introduce a modified formalism for the calculation of
frequency-domain field components, which is extremely useful in evaluating the
spectral properties of the output quantum state. In Sec. III, as an example
application scenario, we apply the formalism to the problem of integrated
quantum spectroscopy with undetected photons (QSUP) and investigate the
effects of a spectrally localised loss on the spectrum of the output photons
of a waveguide SPDC source, operating in the high-gain regime. We show that
the performance of the sensing system is gain-dependent and its signatures
become more prominent as gain is increased, but also that they saturate at
even higher gain values. Finally, in Sec. IV we discuss our results and review
the advantages and use cases of our non-perturbative formalism.

##  II Non-perturbative description of SPDC in open and dispersive systems

We are considering a generally inhomogeneous, dispersive and lossy dielectric
system, characterised by the linear permittivity
$\varepsilon(\mathbf{r},\omega)=\varepsilon^{\prime}(\mathbf{r},\omega)+i\varepsilon^{\prime\prime}(\mathbf{r},\omega)$,
where $\varepsilon^{\prime}$ and $\varepsilon^{\prime\prime}$ are its real and
imaginary parts, respectively. To make our expressions less cumbersome, we
restrict our present analysis to an isotropic system and note that a
generalisation to anisotropic systems is possible by using an appropriately
generalised GF quantisation procedure, such as the one derived in Ref. [55].

We consider the system to have a second-order nonlinearity, characterised by
the second-order nonlinear susceptibility tensor
$\chi^{(2)}_{ijk}(\mathbf{r})$. In the Heisenberg picture, the SPDC
interaction Hamiltonian has the following form [56]:

$\displaystyle\hat{H}_{\mathrm{SPDC}}(t)$ $\displaystyle=$
$\displaystyle-\varepsilon_{0}\int\differential{\mathbf{r}}\chi^{(2)}_{ijk}(\mathbf{r})E^{(-)}_{P,i}(\mathbf{r},t)$
(1)
$\displaystyle\times\hat{E}^{(+)}_{j}(\mathbf{r},t)\hat{E}^{(+)}_{k}(\mathbf{r},t)$
$\displaystyle+H.c.,$

where $\hat{E}^{(+)}_{j,k}(\mathbf{r},t)$ are Cartesian components of the
positive-frequency part of the generated quantum fields,
$E^{(-)}_{P,i}(\mathbf{r},t)$ are the Cartesian components of the negative-
frequency part of the pump field and $\varepsilon_{0}$ is the dielectric
permittivity of vacuum. We emphasise that we are working in the undepleted
pump approximation and thus assume the pump field to be an undepleted
classical pulse, which results in the Hamiltonian being quadratic in the field
operators and is a crucial assumption for finding the input-output operator
relations later in this work. Treating the extremely high-gain scenarios, in
which the pump can also deplete substantially or ones where the pump is in a
few-photon state at the input, results in more complex dynamics which cannot
be described by the Hamiltonian of the quadratic form given in Eq. (1).
However, such scenarios can be used for generation of non-Gaussian quantum
states [57, 58] and advancing the GF formalism to go beyond the undepleted
pump approximation could be an interesting direction for future works.

In addition, we note that for our calculations, we assume a non-dispersive and
real-valued $\chi^{(2)}$ nonlinearity, which corresponds to operating the
nonlinear material in spectral regions that are far from its material
resonances. Such an approximation is widely used in treatment of high-gain
squeezed-light generation in most scenarios of practical interest [9, 31, 32,
33, 34, 56, 59]. In principle, our approach could be generalized to treat
dispersive nonlinearities, by starting the calculations from a more complex
nonlinear Hamiltonian that includes a dispersive $\chi^{(2)}$ coefficient [60,
61, 62, 63], provided that the Hamiltonian remains quadratic in the electric
field operators. We also emphasize, that in our model, we assume an undepleted
pump approximation, which is a well-justified approximation to almost all
cases with a practically reachable value for the nonlinear gain. In this
approximation, the nonlinear Hamiltonian keeps its quadratic form, allowing us
to formulate an input-output relation for the electric field operator.
Finally, it should be noted, that we do not consider the effect of nonlinear
noise, predicted to have an effect for very strong pump intensities [60, 61].

As will be the convention in the remainder of the paper, repeated Latin
indices are implicitly summed over and the domains of integration for spatial
and frequency integrals are $(-\infty,\infty)$ and $[0,\infty)$, respectively,
unless otherwise noted in the text. We carry out all of our calculations in
the Heisenberg picture and assume that the nonlinear interaction is
adiabatically ”turned on” at $t\rightarrow-\infty$ and adiabatically ”turned
off” as $t\rightarrow\infty$.

To allow for the treatment of open and dispersive systems, the electric field
operator is quantised using the GF quantisation method [42, 43, 44]:

$\displaystyle\mathbf{\hat{E}}^{(\pm)}(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\int\differential{\omega}\mathbf{\underline{\hat{E}}}^{(\pm)}(\mathbf{r},\omega,t),$
$\displaystyle\underline{\hat{E}}^{(+)}_{i}(\mathbf{r},\omega,t)$
$\displaystyle=$ $\displaystyle
i\mathcal{K}\,\frac{\omega^{2}}{c^{2}}\int\differential{\mathbf{r}^{\prime}}\sqrt{\varepsilon^{\prime\prime}(\mathbf{r}^{\prime},\omega)}$
(2)
$\displaystyle\times\,G_{ij}(\mathbf{r},\mathbf{r}^{\prime},\omega)\,\hat{f}_{j}(\mathbf{r}^{\prime},\omega,t),$
$\displaystyle\underline{\hat{E}}^{(-)}_{i}(\mathbf{r},\omega,t)$
$\displaystyle=$
$\displaystyle\left(\underline{\hat{E}}^{(+)}_{i}(\mathbf{r},\omega,t)\right)^{\dagger},$

where $\mathcal{K}=\sqrt{\frac{\hbar}{\pi\varepsilon_{0}}}$,
$\mathbf{\underline{\hat{E}}}^{(\pm)}(\mathbf{r},\omega,t)$ is the amplitude
operator of frequency $\omega$,
$G_{ij}(\mathbf{r},\mathbf{r}^{\prime},\omega)$ are the matrix elements of the
dyadic GF of the dielectric system
$\overleftrightarrow{\mathbf{G}}(\mathbf{r},\mathbf{r}^{\prime},\omega)$ and
$\,\hat{f}_{j}(\mathbf{r}^{\prime},\omega,t)$ are the Cartesian components of
the bosonic annihilation operator
$\mathbf{\,\hat{f}}(\mathbf{r}^{\prime},\omega,t)$.

The GF in (2) is the solution to the classical Helmholtz equation:

$\left(\mathbf{\nabla}\times\mathbf{\nabla}\times-\frac{\omega^{2}}{c^{2}}\varepsilon(\mathbf{r},\omega)\right)\overleftrightarrow{\mathbf{G}}(\mathbf{r},\mathbf{r}^{\prime},\omega)=\overleftrightarrow{I}\delta(\mathbf{r}-\mathbf{r}^{\prime}),$
(3)

together with the condition that it vanishes at infinity. In the above
equation $\overleftrightarrow{I}$ is the identity matrix.

The operator $\mathbf{\,\hat{f}}(\mathbf{r}^{\prime},\omega,t)$ annihilates a
local field-matter excitation of frequency $\omega$, located at
$\mathbf{r}^{\prime}$. Along with its adjoint
$\mathbf{\,\hat{f}}^{\dagger}(\mathbf{r}^{\prime},\omega,t)$, it makes up the
fundamental operator algebra of the GF quantisation method and obeys the
canonical commutation relation:

$\commutator{\,\hat{f}_{i}(\mathbf{r},\omega,t)}{\,\hat{f}\,^{\dagger}_{j}(\mathbf{r}^{\prime},\omega^{\prime},t)}=\delta_{ij}\delta(\mathbf{r}-\mathbf{r}^{\prime})\delta(\omega-\omega^{\prime}).$
(4)

Using Eq. (4), it can be shown that the amplitude operators
$\mathbf{\underline{\hat{E}}}^{(\pm)}(\mathbf{r},\omega,t)$ obey the following
commutation relation:

$\displaystyle\commutator{\underline{\hat{E}}_{i}^{(+)}(\mathbf{r},\omega,t)}{\underline{\hat{E}}_{j}^{(-)}(\mathbf{r}^{\prime},\omega^{\prime},t)}$
$\displaystyle=$
$\displaystyle\mathcal{K}^{2}\frac{\omega^{2}}{c^{2}}\delta(\omega-\omega^{\prime})$
(5)
$\displaystyle\times\text{Im}\left[G_{ij}(\mathbf{r},\mathbf{r}^{\prime},\omega)\right],\quad$

where $\text{Im}\left[G_{ij}(\mathbf{r},\mathbf{r}^{\prime},\omega)\right]$ is
the imaginary part of the Green’s function. The full derivation of the above
relation can be found in Appendix A.

Before proceeding, we note that the explicit results of our non-perturbative
approach, presented in this paper, are strictly valid for fields in non-
magnetic systems. However, the main ideas and steps of our formalism could
potentially be adapted for magnetic systems by using an appropriate field
quantisation, such as one used in Ref. [64].

###  II.1 Equations of motion

The Cartesian components of
$\mathbf{\underline{\hat{E}}}^{(\pm)}(\mathbf{r},\omega,t)$ evolve in time
according to the Heisenberg equation of motion:

$\,\partial_{t}\underline{\hat{E}}_{i}^{(\pm)}(\mathbf{r},\omega,t)=\frac{1}{i\hbar}\commutator{\underline{\hat{E}}_{i}^{(\pm)}(\mathbf{r},\omega,t)}{\hat{H}(t)},$
(6)

where $\hat{H}(t)=\hat{H}_{0}(t)+\hat{H}_{\mathrm{SPDC}}(t)$ and
$\hat{H}_{0}(t)=\hbar\int\differential{\mathbf{r}}\int\differential{\omega}\omega\,\hat{f}\,^{\dagger}_{i}(\mathbf{r},\omega,t)\,\hat{f}_{i}(\mathbf{r},\omega,t)$
is the free-field Hamiltonian in the GF quantisation scheme [43]. The
resulting differential equations thus contain terms describing both the
nonlinear and free-field evolution. The latter can be eliminated and the
equations of motion made simpler by introducing the slowly-varying, rotating-
frame operators, $\mathbf{\,\bar{E}}^{(\pm)}(\mathbf{r},\omega,t)$, where:

$\mathbf{\underline{\hat{E}}}^{(\pm)}(\mathbf{r},\omega,t)&=\mathbf{\,\bar{E}}^{(\pm)}(\mathbf{r},\omega,t)\,e^{\mp
i\omega t},\\\
\mathbf{\hat{E}}^{(\pm)}(\mathbf{r},t)=\int\differential{\omega}\mathbf{\,\bar{E}}^{(\pm)}(\mathbf{r},\omega,t)\,e^{\mp
i\omega t}.$ (7)

The exponential factors account for the free-field component of the evolution
and the rotating-frame operators
$\mathbf{\,\bar{E}}^{(\pm)}(\mathbf{r},\omega,t)$ evolve according to:

$\,\partial_{t}\mathbf{\,\bar{E}}^{(\pm)}(\mathbf{r},\omega,t)=\frac{1}{i\hbar}\commutator{\mathbf{\,\bar{E}}^{(\pm)}(\mathbf{r},\omega,t)}{\hat{H}_{\mathrm{SPDC}}(t)}.$
(8)

A proof of the above equation, as well as a more formal definition of the
rotating-frame operators is given in Appendix B.

Rotating-frame creation and annihilation operators
$\bar{\mathbf{f}}^{(\dagger)}(\mathbf{r},\omega,t)$ are defined equivalently
to Eq. (7) and are related to the rotating frame amplitude operators
$\mathbf{\,\bar{E}}^{(\pm)}(\mathbf{r},\omega,t)$ by the relation:

$\displaystyle\,\bar{E}_{i}^{(+)}(\mathbf{r},\omega,t)$ $\displaystyle=$
$\displaystyle
i\mathcal{K}\,\frac{\omega^{2}}{c^{2}}\int\differential{\mathbf{r}^{\prime}}\sqrt{\varepsilon^{\prime\prime}(\mathbf{r}^{\prime},\omega)}$
(9) $\displaystyle\times
G_{ij}(\mathbf{r},\mathbf{r}^{\prime},\omega)\bar{f}_{j}(\mathbf{r}^{\prime},\omega,t),\quad$

and its adjoint for $\,\bar{E}_{i}^{(-)}(\mathbf{r},\omega,t)$ and
$\bar{f}^{\dagger}_{j}(\mathbf{r}^{\prime},\omega,t)$. The rotating-frame
amplitude operators $\mathbf{\,\bar{E}}^{(\pm)}(\mathbf{r},\omega,t)$, as well
as the rotating-frame creation/annihilation operators
$\mathbf{\bar{f}}^{(\dagger)}(\mathbf{r}^{\prime},\omega,t)$, can be shown to
obey equivalent commutation relations as their Heisenberg picture
counterparts:

$\left[\bar{f}_{i}(\mathbf{r},\omega,t),\bar{f}^{\dagger}_{j}(\mathbf{r}^{\prime},\omega^{\prime},t)\right]=\delta_{ij}&\delta(\mathbf{r}-\mathbf{r}^{\prime})\delta(\omega-\omega^{\prime}),\qquad\\\
\commutator{\,\bar{E}_{i}^{(+)}(\mathbf{r},\omega,t)}{\,\bar{E}_{j}^{(-)}(\mathbf{r}^{\prime},\omega^{\prime},t)}=\mathcal{K}^{2}\frac{\omega^{2}}{c^{2}}\delta(\omega-\omega^{\prime})\\\
\times\,\text{Im}\left[G_{ij}(\mathbf{r},\mathbf{r}^{\prime},\omega)\right].$
(10)

The above relations are further discussed in Appendix B and, from this point
onward, operators will be considered exclusively in the rotating frame without
explicitly noting the fact, e.g., ”rotating-frame amplitude operators” will
simply be referred to as ”amplitude operators” and so on.

We find the equations of motion governing the evolution of the Cartesian
components of $\mathbf{\,\bar{E}}^{(+)}(\mathbf{r},\omega,t)$ from Eq. (8)
(the full derivation is given in Appendix C):

$\displaystyle\,\partial_{t}\,\bar{E}_{i}^{(+)}(\mathbf{r},\omega,t)$
$\displaystyle=$
$\displaystyle\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime}}$
(11) $\displaystyle\times
F_{ij}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t)\,\bar{E}_{j}^{(-)}(\mathbf{r}^{\prime},\omega^{\prime},t),\qquad$

with the corresponding equation for
$\,\partial_{t}\,\bar{E}_{i}^{(-)}(\mathbf{r},\omega,t)$ obtained by simply
taking the adjoint of Eq. (11). Here we also defined:

$\displaystyle
F_{ij}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t)$
$\displaystyle=$
$\displaystyle\frac{2i}{\pi}\frac{\omega^{2}}{c^{2}}\text{Im}\left[G_{il}(\mathbf{r},\mathbf{r}^{\prime},\omega)\right]$
(12) $\displaystyle\times
E^{(+)}_{P,k}(\mathbf{r}^{\prime},t)\chi^{(2)}_{klj}(\mathbf{r}^{\prime})\,e^{i(\omega+\omega^{\prime})t}.\qquad$

Before proceeding with our derivation, we note that the equations of motion in
Eq.(11) are derived under the condition of Kleinman’s symmetry (assuming a
non-dispersive nonlinear susceptibility), where the $\chi^{(2)}_{ijk}$ tensor
is symmetric under the permutation of its indices. Although this assumption is
reasonable in many systems [3], it is not necessary for the development of our
formalism, and only results in simpler expressions. The case of Kleinman’s
symmetry not being applicable is discussed in Appendix C.

###  II.2 Input-output relations

The obtained equations of motion Eqs. (11) are linear in terms of the
amplitude operators $\,\bar{E}_{i}^{(\pm)}(\mathbf{r},\omega,t)$. Along with
the fact that the $\,\bar{E}_{i}^{(\pm)}(\mathbf{r},\omega,t)$ are themselves
canonical, i.e., their commutator Eq. (10) is scalar, this enables us to claim
that, at an arbitrary time $t$, $\,\bar{E}_{i}^{(\pm)}(\mathbf{r},\omega,t)$
can be expressed as a linear combination of
$\,\bar{E}_{i}^{(\pm)}(\mathbf{r},\omega,t\rightarrow-\infty)$, i.e., the
free-field amplitude operators [65, 32, 31]. Thus, we can postulate an ansatz
for the general solution of Eqs. (11) in which the time-dependent (output)
operators are expressed in terms of the free-field (input) operators via an
input-output (IO-) relation. Additionally, due to the linear relation between
the amplitude operators $\,\bar{E}_{i}^{(\pm)}(\mathbf{r},\omega,t)$ and
creation/annihilation operators
$\bar{\mathbf{f}}^{(\dagger)}(\mathbf{r},\omega,t)$, as can be seen from Eq.
(9), either the free-field amplitude operators or the free-field
creation/annihilation operators can be used to formulate the IO-relation. For
the purpose of this work, we use the creation/annihilation operators
$\bar{\mathbf{f}}^{(\dagger)}(\mathbf{r}^{\prime},\omega,t\rightarrow-\infty)$,
which we label $\bar{\mathbf{f}}^{(\dagger)}(\mathbf{r}^{\prime},\omega)$ for
compactness. We note that a completely equivalent formalism is obtained if the
amplitude operators
$\mathbf{\,\bar{E}}^{(\pm)}(\mathbf{r},\omega,t\rightarrow-\infty)$ are used
instead and the main reason we opted for
$\bar{\mathbf{f}}^{(\dagger)}(\mathbf{r}^{\prime},\omega)$ was to have more
streamlined final expressions due to their simpler commutation relation in Eq.
(10). With the above discussion in mind, we assume the following ansatz for
the IO-relation:

$\,\bar{E}_{i}^{(+)}(\mathbf{r},\omega,t)=\int\differential{\bm{\xi}}\int\differential{\nu}\left(\mathcal{B}_{ij}(\mathbf{r},\omega;\bm{\xi},\nu;t)\bar{f}_{j}(\bm{\xi},\nu)+\mathcal{A}^{*}_{ij}(\mathbf{r},\omega;\bm{\xi},\nu;t)\bar{f^{\dagger}}_{j}(\bm{\xi},\nu)\right).$
(13)

In the above expression, $\mathcal{A}_{ij}(\mathbf{r},\omega;\bm{\xi},\nu;t)$
and $\mathcal{B}_{ij}(\mathbf{r},\omega;\bm{\xi},\nu;t)$ are the input-output
(IO)-coefficients. In Eq. (13), $\bm{\xi}\text{ and }\nu$ are spatial and
frequency variables, which enumerate the free-field field-matter local
excitations.

The distribution of the contributing excitations is determined by the IO-
coefficients, which are, in turn, determined by the properties of the system
under consideration. A similar IO-relation is often derived in existing works
on high-gain SPDC in lossless systems, where the output normal-mode photon
creation/annihilation operators are expressed as linear combinations of
normal-mode photon operators at the input, and is associated with frequency
mixing caused by the nonlinear interaction [31, 32]. In the present context,
Eq. (13) does indeed represent the output as a linear combination of
excitation of different frequencies, manifested by the frequency integral over
$\nu$. However, it also shows the output to be a linear combination of
excitations at different spatial positions, manifested by the spatial integral
over $\bm{\xi}$. This is a property inherent to the GF quantisation method and
is a consequence of the local nature of the fundamental bosonic modes, which
itself enables the treatment of arbitrary open optical systems. In the
remainder of the paper, to make expressions more compact, we will combine the
variables $\bm{\xi}\text{ and }\nu$ into the vector
$\bm{\Xi}=\left(\nu,\bm{\xi}\right)$ whenever possible and replace
$\int\differential{\bm{\xi}}\int\differential{\nu}$ with
$\int\differential{\bm{\Xi}}$.

Here, we again emphasise that the linearity of Eqs. (11) in terms of the field
operators, necessary for establishing the IO-relation Eq. (13), is a direct
consequence of the undepleted pump approximation and treating the pump field
as a classical pulse (as shown in Appendix C).

The input-output relation (13) allows us to reconstruct the field operators at
any time $t$, if the IO-coefficients
$\mathcal{A}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ and
$\mathcal{B}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ have been obtained. Moreover,
Eq. (13) can be used to express field-dependent quantities directly in terms
of the IO-coefficients. To illustrate this, we consider the first-order field
correlation function
$g_{ij}^{(1)}(\mathbf{r},\mathbf{r}^{\prime};t,t^{\prime})=\expectationvalue{\hat{E}^{(-)}_{i}(\mathbf{r},t)\hat{E}^{(+)}_{j}(\mathbf{r}^{\prime},t^{\prime})}$
as an example. Here, the expectation value is to be taken in the initial state
of the system, since we are working in the Heisenberg picture. To write
$g_{ij}^{(1)}(\mathbf{r},\mathbf{r}^{\prime};t,t^{\prime})$ in terms of the
IO-coefficients, we expand the field operators using Eq. (7) and then replace
the amplitude operators using the IO-relation Eq. (13). If we assume the
initial state to be the vacuum, we obtain (see Appendix D for details):

$\displaystyle g_{ij}^{(1)}(\mathbf{r},\mathbf{r}^{\prime};t,t^{\prime})$
$\displaystyle=$
$\displaystyle\iint\differential{\omega}\differential{\omega^{\prime}}\,e^{i\omega
t-i\omega^{\prime}t^{\prime}}\int\differential{\bm{\Xi}}$ (14)
$\displaystyle\times\mathcal{A}_{ik}(\mathbf{r},\omega,\bm{\Xi};t)\mathcal{A}^{*}_{jk}(\mathbf{r}^{\prime},\omega^{\prime},\bm{\Xi};t^{\prime}).$

If the system is initially not in the vacuum state, e.g., in the case of a
seeded nonlinear process,
$g_{ij}^{(1)}(\mathbf{r},\mathbf{r}^{\prime};t,t^{\prime})$ (and any higher-
order correlation function) can still be expressed in terms of the IO-
coefficients, however, a non-vacuum initial state will also result in a more
complicated expression, more details on this are given in Appendix D. To
obtain numerical values of
$g_{ij}^{(1)}(\mathbf{r},\mathbf{r}^{\prime};t,t^{\prime})$, one need only
find the IO-coefficients at times $t$ and $t^{\prime}$ and insert them into
Eq. (14).

###  II.3 Comparison with existing methods

To make the benefits of the formalism discussed in this paper more apparent,
before continuing, we will briefly review the differences and similarities
between our approach and existing formalisms for treating high-gain SPDC. The
IO-relation Eq. (13) is conceptually identical to Bogoliubov transformations
used in Refs. [31, 32], as well as the ansatz used in Ref. [33] – all of the
aforementioned relations, including Eq. (13), relate the output field
operators to the creation and annihilation operators for the modes of the
input field in a linear fashion. However, the key difference between these
previous works and ours is the type of quantisation used for the field
operators, which results in different definition of the modes of the quantised
field, which in our case allows for direct inclusion of arbitrary amounts of
linear losses. In Refs. [31, 32, 33], where loss/dissipation is completely
neglected, one can use a quantisation scheme where the creation and
annihilation operators correspond to the normal modes of the photonic system,
which are power-orthogonal and spatially non-local. This results in output
fields being expressed as linear combinations of normal-mode photon operators
(distinguished either by frequency or wavevector). In the present context,
with the GF quantisation method, the field operator is expressed in terms of
local bosonic excitations of the field and matter, sometimes also called
polaritonic modes [45]. Such a description directly includes all the effects
coming from the linear light-matter interactions in the photonic system
(embedded in the GF), with potential loss channels. Some forms of linear
losses can also be treated in the case of normal photonic modes, as has been
done in Refs. [9, 39, 59], through coupling the normal modes to weakly lossy
channels. However, such treatment can only account for weak scattering and
absorption losses, as it relies on assuming lossless/power-orthogonal normal
modes for the underlying photonic system that are only weakly perturbed by the
presence of the loss channels. The GF quantisation method can in turn also
account for cases where loss is very strong and/or is inherently shaping the
modal properties of the system, such as in strongly leaky systems (e.g. Mie
nanoresonators [66]) or strongly lossy systems (e.g. plasmonic systems [46]).
This generality of our method comes at the cost of increased computational
efforts, coming from the local nature of the creation and annihilation
operators in the GF method, which requires carrying around the spatial degree
of freedom in all the calculations. Another notable difference between our
formalism and existing ones is the absence of a weak-dispersion approximation,
which is commonly used to simplify the equations of motion for the field
operators by assuming a linear dispersion relation (between the wave-vector
and frequency of the modes) and using the assumption to replace the time
degree of freedom by the spatial one [9, 31, 32, 33, 39, 59]. In certain
cases, higher-order dispersion can also be included through additional terms
in the equations of motion, corresponding to terms in a power-series expansion
of the wave-vector of the photonic modes of the system [9, 39, 59]. In order
to make our formalism applicable to general nanostructured and open systems,
with potentially highly complex dispersion relations for the involved photonic
modes, we do not make any assumptions regarding the dispersion relation and
solve our equations of motion temporally. This does come at the cost of
increased computational effort when implementing the formalism numerically,
but we note that the prescriptions for simplifying the resulting equations in
cases where the weak-dispersion approximation is valid could also be applied
to our formalism to obtain similar computational benefits. Finally, we
emphasize that generalizing the existing photon-based formalisms for inclusion
of weak losses is usually done on a per-system basis, with focus on a specific
system geometry (e.g. a waveguide or a ring resonator). In contrast, the
formalism presented in our work offers a unified approach, applicable to a
wide variety of open photonic systems, and is formulated in a way that is
compatible with existing numerical methods to facilitate its usefulness.

###  II.4 Coupled-mode equations for the IO-coefficients

In order to calculate the IO-coefficients at a certain time, we find the
coupled differential equations governing their time evolution by inserting Eq.
(13) into the Heisenberg equations of motion Eq. (11) as an ansatz to obtain:

$\displaystyle\,\partial_{t}\mathcal{A}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$
$\displaystyle=$
$\displaystyle\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime}}F^{*}_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t)$
(15a)
$\displaystyle\times\mathcal{B}_{kj}(\mathbf{r}^{\prime},\omega^{\prime},\bm{\Xi};t),\qquad$
$\displaystyle\,\partial_{t}\mathcal{B}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$
$\displaystyle=$
$\displaystyle\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime}}F_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t)$
(15b)
$\displaystyle\times\mathcal{A}_{kj}(\mathbf{r}^{\prime},\omega^{\prime},\bm{\Xi};t).\qquad$

The boundary conditions for the IO-coefficients are derived by equating Eq.
(13) with Eq. (9) and setting $t\rightarrow-\infty$ on both sides:

$\displaystyle\mathcal{A}_{ij}(\mathbf{r},\omega;\bm{\Xi};t\rightarrow-\infty)$
$\displaystyle=$ $\displaystyle 0,$ (16a)
$\displaystyle\mathcal{B}_{ij}(\mathbf{r},\omega;\bm{\Xi};t\rightarrow-\infty)$
$\displaystyle=$ $\displaystyle
i\mathcal{K}\,\frac{\nu^{2}}{c^{2}}\sqrt{\varepsilon^{\prime\prime}(\bm{\xi},\nu)}$
(16b) $\displaystyle\times
G_{ij}(\mathbf{r},\bm{\xi},\nu)\delta(\omega-\nu).\quad$

In most practical cases, Eqs. (15) do not have analytical solutions and have
to be solved numerically. To that end, we can rewrite them in a form more
appropriate for implementation, which has the added benefit of more clearly
distinguishing the different contributions to the output field in the IO-
relation of Eq. (13). We begin by decomposing the coefficient
$\mathcal{B}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ in the following way:

$\mathcal{B}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)=b^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi})+B_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$
(17)

where $b^{(0)}_{ij}$ and $B_{ij}$ satisfy the following conditions:

$\displaystyle b^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi})$ $\displaystyle\equiv$
$\displaystyle\mathcal{B}_{ij}(\mathbf{r},\omega,\bm{\Xi};t\rightarrow-\infty)$
$\displaystyle\,\partial_{t}B_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$
$\displaystyle=$
$\displaystyle\,\partial_{t}\mathcal{B}_{ij}(\mathbf{r},\omega,\bm{\Xi};t),$
$\displaystyle B_{ij}(\mathbf{r},\omega,\bm{\Xi};t\rightarrow-\infty)$
$\displaystyle=$ $\displaystyle 0.$

After inserting the decomposition Eq. (17) into Eq. (13) and relabeling
$\mathcal{A}_{ij}\rightarrow A_{ij}$ for notational convenience, the IO-
relation takes the form:

$\,\bar{E}_{i}^{(+)}(\mathbf{r},\omega,t)=\,\bar{E}^{(0,+)}_{i}(\mathbf{r},\omega)+\int\differential{\bm{\Xi}}\left(B_{ij}(\mathbf{r},\omega,\bm{\Xi};t)\bar{f}_{j}(\bm{\Xi})+A^{*}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)\bar{f}^{\dagger}_{j}(\bm{\Xi})\right),$
(18)

where
$\,\bar{E}^{(0,+)}_{i}(\mathbf{r},\omega)=\int\differential{\bm{\Xi}}b^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi})\bar{f}_{j}(\bm{\Xi})$
is exactly the free-field (input) amplitude operator of frequency $\omega$,
since
$\mathbf{\,\bar{E}}^{(+)}(\mathbf{r},\omega,t\rightarrow-\infty)=\mathbf{\,\bar{E}}^{(0,+)}(\mathbf{r},\omega)$.
The remaining two terms in Eq. (18) quantify the changes that the free-field
component of frequency $\omega$ undergoes due to the nonlinear interaction
through the ”new” set of coefficients – $A_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$
and $B_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$.

The differential equations coupling the IO-coefficient
$A_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ and the newly introduced
$B_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ are derived by inserting the
decomposition Eq. (17) into Eqs. (15):

$\displaystyle\,\partial_{t}A_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$
$\displaystyle=$ $\displaystyle S^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$
(19a)
$\displaystyle+\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime}}F^{*}_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t)$
$\displaystyle\times B_{kj}(\mathbf{r}^{\prime},\omega^{\prime},\bm{\Xi};t),$
$\displaystyle\,\partial_{t}B_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$
$\displaystyle=$
$\displaystyle\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime}}F_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t)$
(19b) $\displaystyle\times
A_{kj}(\mathbf{r}^{\prime},\omega^{\prime},\bm{\Xi};t).$

The coupled quantities in the above system have the initial values:
$A_{ij}(\mathbf{r},\omega,\bm{\Xi};t\rightarrow-\infty)=B_{ij}(\mathbf{r},\omega,\bm{\Xi};t\rightarrow-\infty)=0$
and we defined
$S^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)=\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime}}F^{*}_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t)b^{(0)}_{kj}(\mathbf{r}^{\prime},\omega^{\prime},\bm{\Xi})$,
which acts as a source term for equation (19a). After taking the frequency
integral, it will have the following form:

$\displaystyle S^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ $\displaystyle=$
$\displaystyle\frac{2\mathcal{K}}{\pi}\frac{\omega^{2}\nu^{2}}{c^{4}}\sqrt{\varepsilon^{\prime\prime}(\bm{\xi},\nu)}\,e^{-i(\omega+\nu)t}$
(20)
$\displaystyle\times\int\differential{\mathbf{r}^{\prime}}\chi^{(2)}_{klm}(\mathbf{r}^{\prime})E^{(-)}_{P,k}(\mathbf{r}^{\prime},t)$
$\displaystyle\times\text{Im}\left[G_{il}(\mathbf{r},\mathbf{r}^{\prime},\omega)\right]G_{mj}(\mathbf{r}^{\prime},\bm{\xi},\nu).\quad$

The systems of coupled equations (15) and (19), along with their respective
boundary conditions are equivalent formulations of our theoretical formalism
and are the first part of the main results presented in this work. They are
valid in arbitrary dispersive and open nanostructured systems, such as
photonic crystals, nanoresonators, metasurfaces and waveguides, where the
complex spatial and dispersive properties inherent to such devices can all be
accounted for in the electromagnetic GF, as well as arbitrary loss, e.g.
radiation leakage, material absorption, scattering losses etc.

To summarise the main results derived in this section: in order to find the
output electric field operator of a high-gain SPDC system, our proposed
formalism consists of solving the coupled equations (19), which are defined by
the system’s classical Green’s function
$G_{ij}(\mathbf{r},\mathbf{r}^{\prime},\omega)$, linear susceptibility
$\varepsilon(\mathbf{r},\omega)$ and nonlinear susceptibility
$\chi^{(2)}_{ijk}(\mathbf{r})$. These properties, along with the pump field
$\mathbf{E}_{P}^{(\pm)}(\mathbf{r},t)$ are used to define the integration
kernel $F_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t)$ and
source term $S^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$, given by Eqs. (12)
and (20), respectively. Solving the equations provides the IO-coefficients,
which are then inserted into the IO-relation Eq. (18) to obtain the output
amplitude operators $\mathbf{\,\bar{E}}^{(\pm)}(\mathbf{r},\omega,t)$. These,
in turn enable us to reconstruct the full electric field operator as shown in
Eq. (7). The IO-coefficients can also be used to directly calculate scalar
quantities related to the field, such as correlation functions, as was given
in Eq. (14).

###  II.5 Frequency domain formulation

The formalism developed in the previous section enables straightforward
calculation of output fields and field-related quantities and is valid under
very general considerations. However, the fields (and thus, all field-
dependent quantities) are obtained in the time domain and the calculation of
frequency-domain fields and quantities, e.g., the spectra of output photons,
is inefficient to calculate numerically from the obtained time-domain
quantities. This can be seen by examining, for example, the expression for the
single-photon spectrum of the electric field calculated using the IO-
coefficients used in Eq. (18).

The spectrum of a quantised electric field, measured at a time $t$ and
position $\mathbf{r}_{0}$ and summed over all field polarisations is given by
[67]:

$\displaystyle\sigma(\mathbf{r}_{0},\omega_{0},t)$ $\displaystyle=$
$\displaystyle\iint\limits^{\quad
t}_{-\infty}\differential{t^{\prime}}\differential{t^{\prime\prime}}\,e^{-i\omega_{0}(t^{\prime}-t^{\prime\prime})}$
(21)
$\displaystyle\times\expectationvalue{\mathbf{\hat{E}^{(-)}}(\mathbf{r}_{0},t^{\prime})\cdot\mathbf{\hat{E}^{(+)}}(\mathbf{r}_{0},t^{\prime\prime})}.\qquad$

where $\omega_{0}$ is the frequency at which the spectrum is evaluated.

To obtain an expression for $\sigma(\mathbf{r}_{0},\omega_{0},t)$ in terms of
the IO-coefficients $A_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ and
$B_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$, we use the first-order correlation
function Eq. (14) (remembering that $\mathcal{A}_{ij}\equiv A_{ij}$) and
insert it into Eq. (21):

$\displaystyle\sigma(r_{0},\omega_{0},t)$ $\displaystyle=$
$\displaystyle\iint\limits^{\quad
t}_{-\infty}\differential{t^{\prime}}\differential{t^{\prime\prime}}\iint\differential{\omega}\differential{\omega^{\prime}}\int\differential{\bm{\Xi}}$
(22) $\displaystyle\times
A_{ij}(\mathbf{r}_{0},\omega,\bm{\Xi};t^{\prime})A^{*}_{ij}(\mathbf{r}_{0},\omega^{\prime},\bm{\Xi};t^{\prime\prime})$
$\displaystyle\times\,e^{-i(\omega_{0}-\omega)t^{\prime}}\,e^{i(\omega_{0}-\omega^{\prime})t^{\prime\prime}}.$

Although Eq. (22) theoretically allows one to calculate the spectrum, it can
be impractical in a realistic case, when the IO-coefficients have to be
obtained numerically. Namely, in order to perform the
$t^{\prime}/t^{\prime\prime}\text{ and }\omega/\omega^{\prime}$ integrals with
sufficient numerical accuracy, one requires the values of the coefficient
$A_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ at many narrowly spaced times $t$ in
addition to the other variables $\omega$ and $\bm{\Xi}$, which dramatically
increases the computational resources required for such a calculation.

The above limitation is a consequence of the individual frequency amplitudes
$\mathbf{\,\bar{E}}^{(\pm)}(\mathbf{r},\omega,t)$ not being equivalent to the
spectral components of the field that are detected in a spectrally resolved
measurement in the presence of field sources, such as during a nonlinear field
interaction [68, 67].

To overcome these difficulties and further expand the applicability of our
formalism, we develop a frequency-domain formulation, which can be applied in
the calculation of various spectral quantities, while still rigorously
including the effects of arbitrary loss and dispersion. It focuses on the
”spectrally relevant” parts of the field, i.e. special field operators which
directly contribute to the spectrum. These are defined as time-frequency
transforms of the time-domain fields [68] and, as will be shown in the next
sections, evolve in time in a way quite similar (although not identical) to
$\mathbf{\,\bar{E}}^{(\pm)}(\mathbf{r},\omega,t)$. They also have their own
IO-relation, which enables us to use them as fundamental field variables when
determining the time evolution of frequency-domain quantities.

We begin by noting that the field spectrum Eq. (21) can also be written as:

$\sigma(\mathbf{r}_{0},\omega_{0},t)=\expectationvalue{\mathbf{\tilde{E}}^{(-)}(\mathbf{r}_{0},\omega_{0},t)\mathbf{\tilde{E}}^{(+)}(\mathbf{r}_{0},\omega_{0},t)},$
(23)

where the operators $\mathbf{\tilde{E}}^{(\pm)}(\mathbf{r},\omega,t)$ have
been defined as:

$\mathbf{\tilde{E}}^{(\pm)}(\mathbf{r},\omega,t)=\int_{-\infty}^{t}\differential{t^{\prime}}\,e^{\pm
i\omega t^{\prime}}\mathbf{\hat{E}^{(\pm)}}(\mathbf{r},t^{\prime}).$ (24)

We will refer to them as filtered field operators, due to the definition Eq.
(24) being reminiscent of the action of a causal, ideally monochromatic filter
on the time-domain field. In the absence of an interaction, they reduce to the
free-field amplitude operators
$\mathbf{\,\bar{E}}^{(0,\pm)}(\mathbf{r},\omega)$ as $t\rightarrow\infty$. For
finite $t$ and/or with sources (i.e. an interaction) present, their
interpretation is more involved and is discussed in [68, 67].

The expression for the spectrum in terms of the filtered field Eq. (23) is
immediately more appealing than Eq. (21), as it involves the expectation value
of operators at a time $t$, when the spectrum itself is evaluated. Thus, the
need for knowing the fields (or, equivalently, the IO-coefficients) at all
times prior to $t$ is eliminated, provided the filtered fields
$\mathbf{\tilde{E}}^{(\pm)}(\mathbf{r},\omega,t)$ can be obtained without much
additional computational complexity, which, as will be shown in the following,
is indeed the case.

As can be noted from Eq. (24), the filtered field operators are a linear
transform of the ”regular” $\mathbf{\hat{E}}^{(\pm)}(\mathbf{r},t)$ and thus
of $\mathbf{\,\bar{E}}^{(\pm)}(\mathbf{r},\omega,t)$, as well. This enables us
to formulate an analogous input-output relation for the filtered fields, To
that end, we combine Eq. (18) and Eq. (7) and insert the resulting field
decomposition into Eq. (24). Thus we obtain:

$\tilde{E}^{(+)}_{i}(\mathbf{r},\omega,t)=\tilde{E}^{(0,+)}_{i}(\mathbf{r},\omega)+\int\differential{\bm{\Xi}}\left(\tilde{B}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)\bar{f}_{j}(\bm{\Xi})+\tilde{A}^{*}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)\bar{f}^{\dagger}_{j}(\bm{\Xi})\right),$
(25)

where we defined:

$\displaystyle\tilde{E}^{(0,+)}_{i}(\mathbf{r},\omega)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{t}\differential{t^{\prime}}\int\differential{\omega^{\prime}}$
(26a)
$\displaystyle\times\,e^{i(\omega-\omega^{\prime})t^{\prime}}\bar{E}^{(0,+)}_{i}(\mathbf{r},\omega^{\prime}),$
$\displaystyle\tilde{A}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{t}\differential{t^{\prime}}\int\differential{\omega^{\prime}}$
(26b)
$\displaystyle\times\,e^{-i(\omega-\omega^{\prime})t^{\prime}}A_{ij}(\mathbf{r},\omega^{\prime},\bm{\Xi};t^{\prime}),$
$\displaystyle\tilde{B}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{t}\differential{t^{\prime}}\int\differential{\omega^{\prime}}$
(26c)
$\displaystyle\times\,e^{i(\omega-\omega^{\prime})t^{\prime}}B_{ij}(\mathbf{r},\omega^{\prime},\bm{\Xi};t^{\prime}).$

Here, $\tilde{E}^{(0,+)}_{i}(\mathbf{r},\omega)$ is the filtered analogue of
the free-field amplitudes $\bar{E}^{(0,+)}_{i}(\mathbf{r},\omega)$ and
$\tilde{A}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ and
$\tilde{B}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ are the IO-coefficients of the
filtered field. For conciseness, we will refer to them as filtered IO-
coefficients.

An expression for the spectrum in terms of the filtered coefficients can be
found by grouping time and frequency integrals with their corresponding
$A_{ij}$ coefficient in Eq. (22) and recognising that the resulting integral
quantities exactly match the definitions in Eqs. (26). The expression thus
obtained is:

$\displaystyle\sigma(\mathbf{r}_{0},\omega_{0},t)=\int\differential{\bm{\Xi}}\tilde{A}_{ij}(\mathbf{r}_{0},\omega_{0},\bm{\Xi};t)\tilde{A}^{*}_{ij}(\mathbf{r}_{0},\omega_{0},\bm{\Xi};t).\qquad$
(27)

As expected, we no longer require knowledge of the IO-coefficients at narrowly
spaced times, unlike Eq. (22), - we only need the filtered coefficients at
time $t$. To find their values at a given time, we obtain a new set of coupled
equations by combining the definitions Eqs. (26) with the system Eq. (19) (the
full derivation can be found in Appendix E):

$\displaystyle\,\partial_{t}\tilde{A}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$
$\displaystyle=$
$\displaystyle\tilde{S}^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ (28a)
$\displaystyle+\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime}}\tilde{F}^{*}_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t)$
$\displaystyle\times\tilde{B}_{kj}(\mathbf{r}^{\prime},\omega^{\prime},\bm{\Xi};t),$
$\displaystyle\,\partial_{t}\tilde{B}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$
$\displaystyle=$
$\displaystyle\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime}}\tilde{F}_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t)$
(28b)
$\displaystyle\times\tilde{A}_{kj}(\mathbf{r}^{\prime},\omega^{\prime},\bm{\Xi};t),\qquad$

where $\tilde{S}^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ is derived from the
source term $S^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$, present in Eq. (19a)
and is defined as:

$\displaystyle\tilde{S}^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$
$\displaystyle=$ $\displaystyle
2i\mathcal{K}\,\frac{\nu^{2}}{c^{2}}\sqrt{\varepsilon^{\prime\prime}(\bm{\xi},\nu)}\,e^{-i(\omega+\nu)t}$
(29)
$\displaystyle\times\int\differential{\mathbf{r}^{\prime}}\chi^{(2)}_{klm}(\mathbf{r}^{\prime})G_{mj}(\mathbf{r}^{\prime},\bm{\xi},\nu)$
$\displaystyle\times\int\differential{\omega_{p}}\frac{(\omega_{p}-\nu)^{2}}{c^{2}}\mathcal{E}^{(-)}_{P,k}(\mathbf{r}^{\prime},\omega_{p})$
$\displaystyle\times
G^{*}_{il}(\mathbf{r},\mathbf{r}^{\prime},\omega_{p}-\nu)\,e^{i\omega_{p}t}.$

Here, $\mathcal{E}_{P,k}^{(\pm)}(\mathbf{r},\omega_{p})$ is the Fourier
transform of the pump amplitude, defined by
$\mathbf{E}_{P}^{(\pm)}(\mathbf{r},t)=\int\differential{\omega_{p}}\mathbf{\mathcal{E}}_{P}^{(\pm)}(\mathbf{r},\omega_{p})\,e^{\mp
i\,\omega_{p}t}$. The integral kernel
$\tilde{F}_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t)$
which couples the two coefficients in Eq. (28) is given by:

$\displaystyle\tilde{F}_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t)$
$\displaystyle=$
$\displaystyle\frac{2i}{\pi}\int\differential{\bar{\omega}}\frac{\bar{\omega}^{2}}{c^{2}}\chi^{(2)}_{lmk}(\mathbf{r}^{\prime})\mathcal{E}^{(+)}_{P,l}(\mathbf{r}^{\prime},\omega^{\prime}+\bar{\omega})$
(30)
$\displaystyle\times\text{Im}\left[G_{im}(\mathbf{r},\mathbf{r}^{\prime},\bar{\omega})\right]\,e^{i(\omega-\bar{\omega})t},$

and the initial conditions for the filtered coefficients are easily inferred
from the definitions (26):

$\tilde{A}_{ij}(\mathbf{r},\omega,\bm{\Xi};t\rightarrow-\infty)=\tilde{B}_{ij}(\mathbf{r},\omega,\bm{\Xi};t\rightarrow-\infty)\equiv
0.$

The coupled equations Eqs. (28) are the second part of our main results and
can be summarised as follows: to find the output field in the frequency
domain, one must solve the coupled equations (28) to obtain the filtered IO-
coefficients. As was the case in the time-domain formulation, the properties
of the optical system and the pump field are embedded in the source term
$\tilde{S}^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$, defined in Eq. (29), and
the integration kernel
$\tilde{F}_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t)$,
given by Eq. (30). The filtered IO-coefficients can then be used to
reconstruct the frequency domain field operators using the IO-relation Eq.
(25). Alternatively, the filtered IO-coefficients can also be used to directly
evaluate frequency-domain expectation values, such as the field spectrum Eq.
(27), discussed in the beginning of this section, or field moments of the form
$\expectationvalue{\tilde{E}^{(-)}_{i}(\mathbf{r},\omega,t)\tilde{E}^{(+)}_{j}(\mathbf{r},\omega^{\prime},t)}$
and
$\expectationvalue{\tilde{E}^{(+)}_{i}(\mathbf{r},\omega,t)\tilde{E}^{(+)}_{j}(\mathbf{r},\omega^{\prime},t)}$.
The field moments can be used to reconstruct the joint spectral amplitude
(JSA) of the output state of the system [9], which is extremely useful when
investigating high-gain effects, e.g. squeezing [56].

To verify the results of our frequency-domain approach, we find the low-gain
analytical solutions to Eqs. (28) and insert them into Eq. (27) to obtain the
low-gain signal photon spectrum. These calculations are shown in Appendix F
and the obtained expression for the spectrum, shown in Eq. (F), matches the
results obtained in previous works on low-gain SPDC that used the GF
quantisation formalism [46, 50].

##  III Integrated quantum spectroscopy with undetected photons in the high-
gain regime

Quantum spectroscopy with undetected photons (QSUP) is a spectroscopic
technique which relies on frequency-entangled, non-degenerate photon pairs,
where only one of the entangled photons (e.g., the idler) interacts with the
measured sample, but the effects of the sample’s dispersion and/or absorption
can be studied by detecting the other photon (e.g., the signal). This method
can be used as a way to overcome the limited availability of high-sensitivity
detectors in certain frequency regions, e.g., infrared or terahertz [69, 70,
71].

In Ref. [38], a QSUP scheme was proposed, in which a material to be sensed is
placed in the near-field of a waveguide SPDC source. In such a configuration,
when the material has a spectrally localised absorption line around the
frequency range of the idler photons, it also affects the spectral properties
of the generated signal photons [38, 50]. In the aforementioned work, the
proposed scheme was investigated perturbatively, in the low-gain regime of
SPDC. Using our formalism, we can now extend the investigation into the high-
gain regime through numerical simulations of integrated QSUP at varying
amounts of gain and show that the spectral sensitivity of the scheme improves
at higher parametric gains.

We note here that a similar improvement of sensitivity at high gain has
already been observed in other schemes involving measurements with undetected
photons, such as quantum imaging and optical coherence tomography systems
based on nonlinear interferometers and induced coherence [72, 73, 74]. In
these types of schemes, a sample is placed between two coherently pumped high-
gain SPDC sources, which interacts with the output idler arm of the first
source after generation. In contrast, the integrated QSUP scheme proposed and
investigated in the low-gain regime in [38], involves an analyte sample
interacting with the nonlinear waveguide during the pair-production process,
resulting in the spectroscopic information about the sample being imparted
onto the spectrum of the produced photon pairs. The frequency domain
formulation of our formalism is ideally suited for investigating this kind of
scheme, as it enables us to study the spectral properties of the output
photons at high parametric gains in the presence of significant, spectrally-
varying loss, that is directly affecting the SPDC process at the generation
stage.

Figure 1:  Schematic representation of the integrated quantum spectroscopy
process. A nonlinear waveguide SPDC source of length $L$ produces signal
(blue) and idler (red) photons when excited by a pump pulse (purple). The
idler mode, with a larger wavelength and consequently broader mode profile,
interacts with the analyte (pale red) which has an absorption line in the
frequency range of the idler photons. At the output of the waveguide, the
signal and idler photons are separated by a beamsplitter which directs the
signal photons into a spectrally-resolving detector while the idler photons
remain undetected. The periodic light-blue-green region represent the
periodically poled region of the nonlinear waveguide, allowing efficient
phase-matching within the length $L$.

We consider the system shown schematically in Fig. 1: a periodically poled
nonlinear waveguide, oriented along the $z$-axis, is excited by a pump pulse
with a Gaussian temporal envelope. For simplicity, we neglect any polarisation
or transverse dependence and assume that the waveguide is homogeneous,
dispersive and lossy, characterised by the complex, frequency-dependent
permittivity $\varepsilon(\omega)$, which includes both the effect of the
waveguide plus an analyte that is interacting with the near-field of the
waveguide modes. Under the assumptions noted here, the GF of the waveguide has
the analytical form [43]:

$G(z,z^{\prime},\omega)=\frac{i}{2\frac{\omega}{c}\sqrt{\varepsilon(\omega)}}\,e^{i\frac{\omega}{c}\sqrt{\varepsilon(\omega)}\absolutevalue{z-z^{\prime}}}.$
(31)

Since many relevant analytical properties of the GF are direct consequences of
the Kramers-Kronig relations between the real and imaginary parts of the
dielectric permittivity, absorption and dispersion of the model medium must be
consistent with the aforementioned relations to ensure physically meaningful
results [43]. Accordingly, we use a Lorentz model for the permittivities of
the waveguide dielectric and the analyte [75]. The permittivity of the
combined waveguide and analyte system has the form:

$\varepsilon(\omega)=1+\frac{\Omega_{Pl}^{2}}{(\Omega_{0}^{2}-\omega^{2}-i\Gamma\omega)}+\alpha\,\varepsilon_{\mathrm{loss}}(\omega)$
(32)

where $\Omega_{Pl}$ is the plasma frequency of the waveguide dielectric, while
$\Omega_{0}$ and $\Gamma$ are the frequency and width of the dielectric
resonance, respectively. The term
$\varepsilon_{\mathrm{loss}}(\omega)=\frac{\Omega_{Pl}^{2}}{(\omega_{\mathrm{loss}}^{2}-\omega^{2}-i\gamma\omega)}$,
models the spectrally localised loss of the analyte, whose magnitude is
modulated through the unitless factor $\alpha$. Here, $\omega_{\mathrm{loss}}$
and $\gamma$ are the frequency and width of the analyte resonance,
respectively. The plasma frequency of the analyte is assumed to be the same as
the waveguide for simplicity.

The factor $\alpha$ represents effective overlap of the evanescent tails of
the guided mode at the idler frequency with the surrounding analyte. More
specifically, we are assuming a scenario similar to Ref. [38], where the mode
at the longer-wavelength idler has a much stronger overlap with the analyte,
compared to the shorter-wavelength signal and pump. For a given waveguide
geometry, the value of $\alpha$ can be determined by numerically finding the
eigenmodes of the waveguide and calculating the effective permittivities of
the involved modes [38]. In our simulations, the value was chosen to be
$\alpha=3\times 10^{-8}$, corresponding to only a small portion of the idler
mode interacting with the analyte, which was nevertheless sufficient for our
investigations. We note that in our model, to focus on the main physics, we
are not considering the waveguide-geometry-dependent modal dispersions, which
again, can be calculated through a rigorous eigenmode solver, and can be
directly inserted into our model by substituting $\varepsilon(\omega)$ with a
numerically calculated effective quantity.

For the incoming pump pulse, we use the frequency-domain form:
$\mathcal{E}_{P}^{(+)}(z,\omega)=E_{0}\sqrt{\tau_{p}}\,e^{-2\tau_{p}^{2}(\omega_{p_{0}}-\omega)^{2}}\,e^{ik(\omega)z}$.
Here, $\omega_{p0}$ is the pump central frequency, $\tau_{p}$ is the temporal
width of the pulse, $k(\omega)=\frac{\omega}{c}\sqrt{\varepsilon(\omega)}$ and
$E_{0}$ is a normalisation constant determined by the total energy of the pump
pulse $U_{0}$, where $U_{0}\propto\absolutevalue{E_{0}}^{2}$.

###  III.1 Simulation parameters

All of the simulation parameters given here were normalised in the following
manner: frequency quantities are expressed in units of $\omega_{p0}$, temporal
quantities in units of $\frac{1}{\omega_{p0}}$ and spatial quantities in units
of $\frac{2\pi c}{\omega_{p0}}$ (the central wavelength of the pump in
vacuum). The nonlinear region of the waveguide is centred at the coordinate
origin with the length $L=3\times 10^{4}$ and is periodically poled with a
spatial dependence $\chi^{(2)}(z)=\chi^{(2)}_{m}\cos{\Lambda z}$, where
$\chi^{(2)}_{m}$ is the maximum absolute value of the nonlinear permittivity
and $\Lambda$ is the poling period, chosen such that the central phase-matched
frequencies of the signal and idler photons are $\omega_{s0}=0.7$ and
$\omega_{i0}=0.3$, respectively. The resonance frequency of the analyte was
set to be identical to the idler central phase-matching frequency,
$\omega_{\mathrm{loss}}=\omega_{i0}=0.3$, to make the effects of the loss more
prominent. The width of the loss spectrum was chosen to be $\gamma=2.5\times
10^{-3}$. The waveguide parameters in $\varepsilon(\omega)$ were set to be
$\Omega_{0}=2.1$, $\Omega_{Pl}=0.25$ and $\Gamma=10^{-7}$, in order to satisfy
the following conditions:

  1. i)

the resonance of the dielectric is far above the frequency region of interest
and is narrow enough so the dielectric is effectively lossless for frequencies
around and below $\omega_{p0}$;

  2. ii)

the length of the nonlinear region, in combination with the refractive index
of the dielectric, ensures a phase-matching bandwidth significantly wider than
the width of the loss spectrum;

  3. iii)

the pump temporal width was chosen to be $\tau_{p}=2400$, resulting in a
bandwidth sufficiently narrower than the loss spectrum to allow for observing
the spectral correlation between the signal photon spectrum and the spectrum
of the loss.

The strength of the nonlinear interaction, which contributes to determining
the parametric gain, is characterised by the product
$\chi^{(2)}_{m}E_{0}\propto\chi^{(2)}_{m}\sqrt{U_{0}}$, which is varied to
simulate SPDC at different values of parametric gain. In our case, this is
done by adjusting the pump pulse energy $U_{0}$, while keeping the maximum
nonlinear susceptibility $\chi^{(2)}_{m}$ constant.

Finally, all of our simulations were performed for a sufficiently long time
interval, in which the pump has had enough time to completely pass through the
structure.

###  III.2 Lossless medium

We begin our investigation by studying the single-photon spectrum of SPDC at
different values of parametric gain in the ”lossless” case, i.e. without an
analyte (the loss associated with the waveguide dielectric, although
negligibly small in the frequency region of interest, is nevertheless
accounted for in $\varepsilon(\omega)$). For different values of pump pulse
energy $U_{0}$, we calculate the IO-coefficients by numerically solving Eqs.
(28), which are then used to evaluate the spectrum, as per Eq. (27).

To extract the values of parametric gain associated with different pump
intensities, we follow the prescription detailed in Ref. [33]: we plot the
dependence of the maximal intensity of the single-photon spectrum
$\mathcal{I}_{m}$ (at $\omega_{s}=\omega_{s0}=0.7$) as a function of the pump
energy $U_{0}$ and fit it with the well-known dependence of the single-mode
intensity in the case of a two-mode squeezer:
$\mathcal{I}_{m}=a\sinh^{2}{(b\sqrt{U_{0}})}$, where $a$ and $b$ are fitting
parameters. The parametric gain $\mathcal{G}$ for a particular value of pump
energy is then defined as $\mathcal{G}=b\sqrt{U_{0}}$. In Fig. 2a, we show the
dependence of the spectrum maximum on pump pulse energy for pump pulses of two
different bandwidths, defined by their temporal widths $\tau_{p}=2400$ and
$\tau_{p}=600$, and observe that they both indeed follow the $\sinh^{2}$ law.
The wider-bandwidth pulse, corresponding to $\tau_{p}=600$, was included in
our simulations to make sure our formalism correctly predicts that a wider
pump spectrum results in higher parametric gain for the same total pulse
energy [56, 33]. The inset of Fig. 2a shows that the maximal intensity is
approximately independent of the pump bandwidth at low pump energies, again,
in accordance with previous observations [56, 33]. The dependence of the
parametric gain on the pulse energy is shown in Fig. 2b, where the higher
amount of parametric gain obtained for wider pump bandwidths is evident.

Figure 2: (a) The dependence of the spectral maximum at $\omega=0.7$ on the
total pump energy for the narrow-band pump (red squares) with $\tau_{p}=2400$
and wide-band pump (blue circles) with $\tau_{p}=600$. The solid lines
indicate the fitted $\sinh^{2}{(b\sqrt{U_{0}})}$ and the inset shows the
intensity dependencies at very low pump pulse energies. In this regime, the
behaviour of the spectrum becomes independent of the pump bandwidth. (b) The
dependence of the observed parametric gain on the square root of the pump
energy for the narrow- (red squares) and wide-band pump pulses (blue circles).
In accordance with previous works, the wide-band pump results in higher
parametric gain for a given pump pulse energy.

Aside from the behaviour of the spectrum amplitudes, our formalism also
correctly predicts the broadening of the single- and two-photon spectra,
occurring as gain is increased [56, 33]. The broadening of the two-photon
spectrum can be observed by studying the second-order moment of the output
field
$N(\omega,\omega^{\prime})=\expectationvalue{\mathbf{\tilde{E}}^{(+)}(\mathbf{r}_{0},\omega,t)\cdot\mathbf{\tilde{E}}^{(+)}(\mathbf{r}_{0},\omega^{\prime},t)}$,
which is related to the JSA of the output photons [59]. We can calculate the
moment directly in terms of the IO-coefficients using the expression:

$N(\omega,\omega^{\prime})=\int\differential{\bm{\Xi}}\tilde{B}_{ij}(\mathbf{r}_{0},\omega,\bm{\Xi};t)\tilde{A}^{*}_{ij}(\mathbf{r}_{0},\omega^{\prime},\bm{\Xi};t),$
(33)

which is obtained by expanding the filtered field operators using Eq. (25) and
taking the expectation value. $N(\omega_{i},\omega_{s})$, in the case of wide-
band pump ($\tau_{p}=600$), is shown in Fig. 3a for two values of parametric
gain and we observe that higher gain indeed results in the moment becoming
broader in frequency and, correspondingly, in output photons to exhibit less
frequency correlations. Additionally, in Fig. 3b, we show the normalised
single photon spectrum at different values of gain for the wide-band pump. As
gain is increased, the spectrum broadens, in accordance with previous results
on high-gain lossless SPDC [56, 33].

Figure 3:  (a) The normalised second-order field moment
$N(\omega_{i},\omega_{s})$ for the wide-band pump ($\tau_{p}=600$) at low
parametric gain - $\mathcal{G}=0.0915$ and high parametric gain -
$\mathcal{G}=3.54$. The output state at higher parametric gain exhibits less
frequency correlations between the photons due to gain-induced broadening. (b)
Normalised single-photon spectrum of the signal photons at the output (each
spectrum normalised to its maximum) for the wide-band pump with $\tau_{p}=600$
at varying gain.

###  III.3 Lossy medium

To investigate the effects of increasing gain in QSUP, we introduce the
analyte into the material permittivity $\varepsilon(\omega)$ and calculate
signal photon spectra at the output. In Fig. 4a we show the output spectra
without the analyte and in Fig. 4b we show the spectra with the analyte
present. The effect of the idler loss is immediately evident, as the spectra
show a dip in the signal intensity, centred around the frequency
$\omega^{(c)}=0.7$, which corresponds to the resonance frequency of the
analyte $\omega_{\mathrm{loss}}$ through the conservation of energy
$\omega_{p0}=\omega^{(c)}+\omega_{\mathrm{loss}}$. Both sets of spectra were
obtained for the narrow-band pump with $\tau_{p}=2400$.

Figure 4: Normalised single-photon spectrum of the signal photons for the
narrow-band pump with $\tau_{p}=2400$ at varying amounts of parametric gain,
(a) without the analyte, and (b) with the analyte. The analyte loss is centred
around $\omega_{\mathrm{loss}}=0.3$. In the lossy case, each spectrum is
normalised to its lossless counterpart with the same gain, in order to better
showcase the changes in the spectral shape of the dip.

The appearance of the spectral dip in the presence of idler loss is in
accordance with the perturbative results of Ref. [38], where the perturbative
calculations also suggested the spectral shape of the dip (i.e. its relative
depth and width) to be independent of gain. However, using our non-
perturbative formalism, we observe that the depth and spectral width of the
dip does change as gain is increased. To quantify and investigate this
behaviour we calculate the signal photon extinction spectrum (defined as the
difference between the lossless and lossy spectra) for different values of
gain. The extinction spectra and their properties at different values of gain
are shown in Fig. 5. We observe two main tendencies: the maximal extinction
(equivalent to the relative depth of the dip in the signal spectrum at
$\omega_{s}=\omega^{(c)}=0.7$) increases with gain, while the width of the dip
decreases.

In Fig. 5a, we show extinction spectra obtained at various amounts of gain,
which are normalised to $1$ to better showcase the change in the spectral
width as gain is increased. The lineshape of the analyte absorption spectrum
(i.e. $\text{Im}\left[\varepsilon_{\mathrm{loss}}\right]$) is also shown for
comparison. The dependence of the extinction maximum and its full-width-at-
half-maximum (FWHM) on parametric gain are shown in Fig. 5b. We observe that
the maximum monotonically increases with gain, but with a progressively slower
rate as gain rises and seems to show signs of saturation at very high gain
values. On the other hand, the extinction FWHM decreases with gain (also shown
in Fig. 5b) and shows signs of saturation at higher gain values as well.

Both the monotonic increase of the extinction maximum and the decrease of its
FWHM can be explained as a consequence of self-seeding being impeded by the
analyte loss: idler photons with frequencies within the loss spectrum of the
analyte are removed from the idler field before they have a chance to seed the
production of further photons at those frequencies. This results in signal
photon intensity at frequencies affected by the loss to increase at a lower
rate (with increasing gain) as compared to the lossless case. On the other
hand, the intensity at frequencies unaffected by the loss increases at the
”lossless” rate, resulting in a net increase of the maximum extinction.
Analogously, the decreasing extinction FWHM with gain can also be seen as a
consequence of self-seeding being impeded in the presence of loss, but here,
we also need to take into account the shape of the loss spectrum. The
lineshape of the analyte absorption spectrum dictates that idler photons
experience less loss, the further their frequency is from the analyte
resonance $\omega_{\mathrm{loss}}$. The idler photons further from resonance
thus have a lower chance to be absorbed before seeding further pairs, leading
to the signal photon intensities increasing faster (with gain) for frequencies
further away from $\omega^{(c)}$, resulting in a net narrowing of the
extinction spectrum.

Figure 5:  (a) The extinction spectra of the signal photon at different values
of parametric gain (solid lines) and the lineshape of the idler loss, namely
$\text{Im}\left[\varepsilon_{\mathrm{loss}}\right]$ (dashed line). The plots
are all normalised to their own maximum amplitude to showcase the extinction
spectrum becoming narrower as gain is increased. (b) The dependence of the
extinction maximum (blue circles) and the FWHM of the extinction spectrum on
parametric gain (red squares). The solid lines are interpolation curves to
help illustrate the observed tendencies of the extinction spectrum.

As mentioned at the beginning of this section, sensitivity enhancement in the
high-gain regime has already been observed for undetected photon measurement
schemes based on induced coherence and nonlinear interferometers. The nature
of the enhancement in the case of induced coherence is discussed in detail in
Ref. [72] and the authors conclude that it also occurs due to seeding effects,
where idler photons from the first source, after interacting with the sample,
seed further pairs in the second source and thus impart the effects of the
sample onto those newly generated photons.

The observed ”saturation” of both the extinction maximum and FWHM at very high
gain has, to our knowledge, not been reported before, but it too can be
understood to be another consequence of self-seeding, more specifically, high-
order self-seeding. In such cases, idler (signal) photons of frequency
$\omega$ seed the generation of further signal (idler) photons, not only at
the frequency $\omega_{p0}-\omega$, but also at frequencies around it. At
sufficiently high values of gain, the idler (signal) photons thus generated,
can seed further signal (idler) photons and so on. This enables signal and
idler photons unaffected by the loss to potentially seed the production of
additional photons within the extinction spectrum. Eventually, these high-
order effects compensate the gain reduction experienced by signal photons
within the extinction region, resulting in the spectral shape of the
extinction saturating at very high parametric gain. As these effects
correspond to higher-order processes, they only appear at sufficiently high
gain values ($\mathcal{G}>2$ in our case), but become more prominent as gain
is further increased.

###  III.4 Spatial evolution of the spectrum

The spectra given and discussed in Sec. III.2 and III.3 were obtained assuming
a detector located immediately at the output of the nonlinear waveguide,
however, our formalism also allows the study of the spectrum of the produced
photons throughout the entire length of the nonlinear medium - all obtained by
solving Eqs. (28) once, for a given value of parametric gain.

Although we do not use the evolution of the spectrum within the nonlinear
waveguide to obtain new insight in the context of QSUP, it can be used to
showcase that our non-perturbative formalism can predict the spatial
functionality of the field and field-related quantities. In Fig. 6a, we show
the evolution of the signal photon spectrum as the photons propagate through
the nonlinear waveguide without the analyte present. At the start of the
waveguide, we see the spectrum having a very wide bandwidth and negligible
amplitude; as we go forward in length, we see phase matching taking hold -
reducing the bandwidth and increasing the amplitude. Additionally, in Fig. 6b
we see the effect of the idler loss accumulating during propagation, causing
the dip in spectral intensity, which increases with propagation, according to
the discussion in Sec. III.3.

Figure 6:  The normalised spectral intensity of the signal photons as a
function of their frequency and position within the nonlinear waveguide in the
(a) lossless and (b) lossy cases. Both plots were obtained for the narrow-band
pump ($\tau_{p}=2400$) and parametric gain of $\mathcal{G}=1.33$. (c)
Normalised maximum spectral intensities of the signal photons as a function of
position within the nonlinear waveguide, obtained for the narrow-band pump
($\tau_{p}=2400$) in the lossless case, for different values of gain.

Finally, in Fig. 6c, we show the spatial behaviour of the spectral maximum of
the signal photons for increasing values of gain, in the lossless case. We
observe that, in the low-gain regime, the spectral maximum obeys the well-
established quadratic law, where the intensity is $\propto l^{2}$, $l$ being
the length of the nonlinear region. As gain is increased, the exponential
nature of the length dependence of the maximum spectral intensity becomes more
apparent.

##  IV Conclusion

In summary, we presented a non-perturbative formalism for the description of
high-gain SPDC, applicable to a wide variety of nanostructured and/or open
optical systems with arbitrarily high complex spatial and spectral properties.
Our formalism enables the calculation of the field operators and field-
dependent quantities (e.g., correlation functions and spectra) at the output
and within the nonlinear optical system.

In contrast to previous methods, our formalism is capable of describing
systems with arbitrary amounts of dispersion and loss, through the use of the
Green’s function quantisation method, which intrinsically takes into account
these effects.

We presented both a time-domain, as well as a frequency-domain formulation,
which allows the formalism to be applicable to many types of temporal and
spectral analyses in the field of quantum technologies. As an example, we used
this formalism to investigate quantum spectroscopy with undetected photons in
a nonlinear waveguide surrounded by an analyte with a loss spectrum
corresponding to the frequencies of one of the output photons. We have thus
expanded upon the results of previous low-gain investigations [38] and
discovered that the spectral sensitivity of the QSUP scheme is heavily
dependent on the amount of the parametric gain and can be improved by
operating the nonlinear waveguide in the high-gain regime.

Although derived here for the case of SPDC, the formalism can be generalised
to treat high-gain SFWM in lossy and dispersive systems, since the SFWM
interaction Hamiltonian has an identical operator structure as Eq. (1), when
considered in the undepleted pump approximation [76].

We believe the formalism presented here will advance the design and
implementation of nonlinear sources of quantum light in many emerging quantum
technologies, especially ones implemented on nanostructured platforms. In
addition to integrated squeezed light and high-order Fock state sources in
nanostructured systems with arbitrary loss and dispersion, other examples
include: hybrid systems [5], quantum sensing applications [15] and quantum
frequency conversion in the optical domain [77] or microwave to optical domain
[78].

###### Acknowledgements.

This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under the project identifier 398816777-SFB 1375 (NOA),
the German Federal Ministry of Education and Research (BMBF) under the project
identifiers 13N14877 (QuantIm4Life), 13N16108 (PhoQuant), 13N15956 (LIVE2QMIC)
and the Thuringian Ministry for Economy, Science, and Digital Society and the
European Social Funds (2021 FGI 0043 -Quantum Hub Thuringia). Sina Saravi
acknowledges the funding by the LIGHT profile (FSU Jena).

##  Appendix A Commutation relation for the frequency amplitude operators

We begin by expanding the
$\underline{\hat{E}}^{(\pm)}_{i,j}(\mathbf{r},\omega,t)$ operators in the
commutator according to Eq. (2):

$\displaystyle\commutator{\underline{\hat{E}}_{i}^{(+)}(\mathbf{r},\omega,t)}{\underline{\hat{E}}_{j}^{(-)}(\mathbf{r}^{\prime},\omega^{\prime},t)}$
$\displaystyle=$
$\displaystyle\mathcal{K}^{2}\,\frac{\omega^{2}}{c^{2}}\frac{\omega^{\prime
2}}{c^{2}}\iint\differential{\mathbf{r}_{1}}$
$\displaystyle\differential{\mathbf{r}_{2}}\sqrt{\varepsilon^{\prime\prime}(\mathbf{r}_{1},\omega)}\sqrt{\varepsilon^{\prime\prime}(\mathbf{r}_{2},\omega^{\prime})}$
$\displaystyle\times$ $\displaystyle
G_{ik}(\mathbf{r},\mathbf{r}_{1},\omega)G^{*}_{jl}(\mathbf{r}^{\prime},\mathbf{r}_{2},\omega^{\prime})$
$\displaystyle\times$
$\displaystyle\commutator{\,\hat{f}_{k}(\mathbf{r}_{1},\omega,t)}{\,\hat{f}\,^{\dagger}_{l}(\mathbf{r}_{2},\omega^{\prime},t)}.$

Now, we use the commutation relation Eq. (4) and the properties of the
$\delta$-functions to obtain on the right-hand side:

$\mathcal{K}^{2}\,\frac{\omega^{4}}{c^{4}}\delta(\omega-\omega^{\prime})\\\
\times\int\differential{\mathbf{r}_{1}}\varepsilon^{\prime\prime}(\mathbf{r}_{1},\omega)G_{ik}(\mathbf{r},\mathbf{r}_{1},\omega)G^{*}_{jk}(\mathbf{r}^{\prime},\mathbf{r}_{1},\omega).$

The above expression can be simplified by using the following GF identity
[42]:

$\displaystyle\frac{\omega^{2}}{c^{2}}\int\differential{\bar{\mathbf{r}}}\varepsilon^{\prime\prime}(\bar{\mathbf{r}},\omega)G_{ik}(\mathbf{r},\bar{\mathbf{r}},\omega)$
$\displaystyle G^{*}_{jk}(\mathbf{r}^{\prime},\bar{\mathbf{r}},\omega)=$ (34)
$\displaystyle\text{Im}\left[G_{ij}(\mathbf{r},\mathbf{r}^{\prime},\omega)\right],$

to finally arrive at Eq. (5).

##  Appendix B The rotating frame

The rotating frame operators introduced in Sec. II can be more formally
defined through a decomposition of the unitary evolution operator that
separates the free-field evolution and the evolution induced by the nonlinear
interaction. This procedure is very similar to the definition of the
interaction picture of quantum mechanics.

In general, the Heisenberg equation of motion (6) is equivalent to the
relation:

$\underline{\hat{E}}_{i}^{(+)}(\mathbf{r},\omega,t)=\hat{U}^{\dagger}(t)\underline{\hat{E}}_{i}^{(+)}(\mathbf{r},\omega)\hat{U}(t),$
(35)

where
$\hat{U}(t)=\mathcal{T}\,e^{\frac{1}{i\hbar}\int_{-\infty}^{t}\differential{t^{\prime}}\left(\hat{H}_{0}+\hat{H}_{\mathrm{SPDC}}(t^{\prime})\right)}$,
$\hat{H}_{0}$ and $\hat{H}_{\mathrm{SPDC}}(t^{\prime})$ are the free-field and
nonlinear interaction Hamiltonians, respectively, both considered in the
Schrödinger picture and $\mathcal{T}$ indicates the time-ordering
superoperator. Lastly,
$\underline{\hat{E}}_{i}^{(+)}(\mathbf{r},\omega)\equiv\underline{\hat{E}}_{i}^{(+)}(\mathbf{r},\omega,t\rightarrow\
-\infty)$ represents the field amplitude operators before the nonlinear
interaction, where we omitted the time $t\rightarrow-\infty$ for compactness.
By using Feynman’s disentanglement theorem [79], we can decompose $\hat{U}(t)$
as $\hat{U}(t)=\hat{U}_{0}(t)\hat{U}_{\mathrm{SPDC}}(t)$, where:

$\hat{U}_{0}(t)=&\,e^{\frac{1}{i\hbar}\hat{H}_{0}t}\\\
\hat{U}_{\mathrm{SPDC}}(t)=\mathcal{T}\,e^{\frac{1}{i\hbar}\int_{-\infty}^{t}\differential{t^{\prime}}\tilde{H}_{\mathrm{SPDC}}(t^{\prime})}\\\
\tilde{H}_{\mathrm{SPDC}}(t)=\hat{U}_{0}^{\dagger}(t)\hat{H}_{\mathrm{SPDC}}(t)\hat{U}_{0}(t).$
(36)

We can use the decomposition of $\hat{U}(t)$ to rewrite Eq. (35) as:

$\underline{\hat{E}}_{i}^{(+)}(\mathbf{r},\omega,t)=\hat{U}^{\dagger}_{\mathrm{SPDC}}(t)\hat{U}_{0}^{\dagger}(t)\underline{\hat{E}}_{i}^{(+)}(\mathbf{r},\omega)\hat{U}_{0}(t)\hat{U}_{\mathrm{SPDC}}(t),$

using
$\hat{H}_{0}=\hbar\int\differential{\mathbf{r}}\int\differential{\omega}\omega\,\hat{f}\,^{\dagger}_{i}(\mathbf{r},\omega)\,\hat{f}_{i}(\mathbf{r},\omega)$
and the definition (2) we then obtain:

$\underline{\hat{E}}_{i}^{(+)}(\mathbf{r},\omega,t)=\hat{U}^{\dagger}_{\mathrm{SPDC}}(t)\underline{\hat{E}}_{i}^{(+)}(\mathbf{r},\omega)\hat{U}_{\mathrm{SPDC}}(t)\,e^{-i\omega
t},$

If we now define:

$\,\bar{E}_{i}^{(+)}(\mathbf{r},\omega,t)=\hat{U}^{\dagger}_{\mathrm{SPDC}}(t)\underline{\hat{E}}_{i}^{(+)}(\mathbf{r},\omega)\hat{U}_{\mathrm{SPDC}}(t),$
(37)

we arrive at Eq. (7). The rotating-frame creation/annihilation operators can
be defined in a similar manner to be:

$\bar{f}^{(\dagger)}_{i}(\mathbf{r},\omega,t)=\hat{U}^{\dagger}_{\mathrm{SPDC}}(t)\hat{f}_{i}^{(\dagger)}(\mathbf{r},\omega)\hat{U}_{\mathrm{SPDC}}(t),$
(38)

where again we have
$\hat{f}_{i}^{(\dagger)}(\mathbf{r},\omega)\equiv\hat{f}_{i}^{(\dagger)}(\mathbf{r},\omega,t\rightarrow-\infty)$
for compactness. As with the amplitude operators, the rotating-frame
annihilation operators are related to their Heisenberg picture counterparts
as:

$\hat{f}_{j}(\mathbf{r},\omega,t)=\bar{f}_{j}(\mathbf{r},\omega,t)\,e^{-i\omega
t},$ (39)

with the relation for the creation operators obtained by taking the adjoint of
Eq. (39).

To find the general equation of motion
for$\,\bar{E}_{i}^{(+)}(\mathbf{r},\omega,t)$, we begin by finding the time
derivative of Eq. (37):

$\displaystyle\,\partial_{t}\,\bar{E}_{i}^{(+)}(\mathbf{r},\omega,t)$
$\displaystyle=$
$\displaystyle\left(\,\partial_{t}\hat{U}^{\dagger}_{\mathrm{SPDC}}(t)\right)\underline{\hat{E}}_{i}^{(+)}(\mathbf{r},\omega)\hat{U}_{\mathrm{SPDC}}(t)$
$\displaystyle+\hat{U}^{\dagger}_{\mathrm{SPDC}}(t)\underline{\hat{E}}_{i}^{(+)}(\mathbf{r},\omega)\left(\,\partial_{t}\hat{U}_{\mathrm{SPDC}}(t)\right).$

We then recall Eq. (36) and obtain, after some straightforward calculation:

$\displaystyle\,\partial_{t}\,\bar{E}_{i}^{(+)}(\mathbf{r},\omega,t)=$
$\displaystyle\frac{1}{i\hbar}\hat{U}^{\dagger}_{\mathrm{SPDC}}(t)$
$\displaystyle\commutator{\underline{\hat{E}}_{i}^{(+)}(\mathbf{r},\omega)}{\tilde{H}_{\mathrm{SPDC}}(t)}\hat{U}_{\mathrm{SPDC}}(t).$

Due to the properties of the commutator,
$\hat{U}^{(\dagger)}_{\mathrm{SPDC}}(t)$ can act on the operators inside of it
independently and give, on the right-hand side:

$\frac{1}{i\hbar}\commutator{\,\bar{E}_{i}^{(+)}(\mathbf{r},\omega,t)}{\hat{U}^{\dagger}_{\mathrm{SPDC}}(t)\tilde{H}_{\mathrm{SPDC}}(t)\hat{U}_{\mathrm{SPDC}}(t)},$

where we used Eq. (37). Also recalling Eq. (36) and the decomposition
$\hat{U}(t)=\hat{U}_{0}(t)\hat{U}_{\mathrm{SPDC}}(t)$, we can identify
$\hat{U}^{\dagger}_{\mathrm{SPDC}}(t)\tilde{H}_{\mathrm{SPDC}}(t)\hat{U}_{\mathrm{SPDC}}(t)$
to be the full, Heisenberg-picture SPDC Hamiltonian given in Eq. (1). Thus we
have proven Eq. (8).

The commutation relation for the rotating-frame creation/annihilation
operators is easily found by replacing the definition Eq. (39) in the
commutation relation for the Heisenberg picture operators Eq. (4):

$\left[\bar{f}_{i}(\mathbf{r},\omega,t),\bar{f}^{\dagger}_{j}(\mathbf{r}^{\prime},\omega^{\prime},t)\right]=\delta_{ij}\delta(\mathbf{r}-\mathbf{r}^{\prime})\delta(\omega-\omega^{\prime})\,e^{i(\omega-\omega^{\prime})t}.$

Since the $\delta$-function in the above expression results in
$\omega\equiv\omega^{\prime}$, the exponential factor can be ignored. The same
is valid for the commutation relation of the rotating-frame amplitude
operators, which is found by replacing Eq. (7) in Eq. (10):

$\commutator{\,\bar{E}_{i}^{(+)}(\mathbf{r},\omega,t)}{\,\bar{E}_{j}^{(-)}(\mathbf{r}^{\prime},\omega^{\prime},t)}=\\\
\mathcal{K}^{2}\frac{\omega^{2}}{c^{2}}\text{Im}\left[G_{ij}(\mathbf{r},\mathbf{r}^{\prime},\omega)\right]\delta(\omega-\omega^{\prime})\,e^{i(\omega-\omega^{\prime})t}.$

##  Appendix C Equations of motion for the amplitude operators

To find the commutator
$\commutator{\,\bar{E}_{i}^{(+)}(\mathbf{r},\omega,t)}{\hat{H}_{\mathrm{SPDC}}(t)}$,
we first expand the SPDC Hamiltonian from Eq. (1) in terms of the rotating-
frame operators $\mathbf{\,\bar{E}}^{(\pm)}(\mathbf{r},\omega,t)$:

$\hat{H}_{\mathrm{SPDC}}(t)=-\varepsilon_{0}\int\differential{\mathbf{r}^{\prime}}\iint\differential{\omega^{\prime}}\differential{\omega^{\prime\prime}}\chi^{(2)}_{ljk}(\mathbf{r}^{\prime})E^{(-)}_{P,l}(\mathbf{r}^{\prime},t)\,\bar{E}_{j}^{(+)}(\mathbf{r}^{\prime},\omega^{\prime},t)\,\bar{E}_{k}^{(+)}(\mathbf{r}^{\prime},\omega^{\prime\prime},t)\,e^{-i\omega^{\prime}t}\,e^{-i\omega^{\prime\prime}t}+H.c.$

Due to the linear properties of the commutator and the fact that amplitude
operators with the same sign in the superscript commute, it is sufficient to
calculate:

$\displaystyle\commutator{\,\bar{E}_{i}^{(+)}(\mathbf{r},\omega,t)}{\,\bar{E}_{j}^{(-)}(\mathbf{r}^{\prime},\omega^{\prime},t)\,\bar{E}_{k}^{(-)}(\mathbf{r}^{\prime},\omega^{\prime\prime},t)}=\mathcal{K}^{2}\frac{\omega^{2}}{c^{2}}$
$\displaystyle\Big{(}\delta(\omega-\omega^{\prime})\text{Im}\left[G_{ij}(\mathbf{r},\mathbf{r}^{\prime},\omega)\right]\,\bar{E}_{k}^{(-)}(\mathbf{r}^{\prime},\omega^{\prime\prime},t)$
(41)
$\displaystyle+\delta(\omega-\omega^{\prime\prime})\text{Im}\left[G_{ik}(\mathbf{r},\mathbf{r}^{\prime},\omega)\right]\,\bar{E}_{j}^{(-)}(\mathbf{r}^{\prime},\omega^{\prime},t)\Big{)},$

where we used Eq. (10). As pointed out in the main text, to simplify the
relations, we assume that the Kleinman symmetry condition is valid. This is
applicable when the frequency dependence of the nonlinear susceptibility is
negligible, as is commonly assumed in many practical scenarios of interest
[3], as well as in the present case. In this case, we have the permutation
symmetry of the nonlinear susceptibility $\chi^{(2)}_{ljk}=\chi^{(2)}_{lkj}$,
hence the two terms on the right-hand-side (r.h.s.) of Eq. (41) will result in
identical contributions to the final equation of motion. If Kleinman symmetry
is not valid, we can simply keep both terms in Eq. (41), which adds no
complexity to the calculation, and only makes the expressions longer. Thus we
can write:

$\,\partial_{t}\,\bar{E}_{i}^{(+)}(\mathbf{r},\omega,t)=2i\frac{\varepsilon_{0}\mathcal{K}^{2}\omega^{2}}{\pi
c^{2}}\int\differential{\mathbf{r}^{\prime}}\chi^{(2)}_{ljk}(\mathbf{r}^{\prime})E^{(+)}_{P,l}(\mathbf{r}^{\prime},t)\text{Im}\left[G_{ik}(\mathbf{r},\mathbf{r}^{\prime},\omega)\right]\iint\differential{\omega^{\prime}}\differential{\omega^{\prime\prime}}\delta(\omega-\omega^{\prime\prime})\,\bar{E}_{j}^{(-)}(\mathbf{r}^{\prime},\omega^{\prime},t)\,e^{i\omega^{\prime}t}\,e^{i\omega^{\prime\prime}t}.$

After performing the $\omega^{\prime\prime}$-integral and replacing the value
of $\mathcal{K}$, we obtain the final form of the equation of motion given in
Eq. (11).

##  Appendix D Derivation of Eq. (14)

To find the correlation function
$g_{ij}^{(1)}(\mathbf{r},\mathbf{r}^{\prime};t,t^{\prime})=\expectationvalue{\hat{E}^{(-)}_{i}(\mathbf{r},t)\hat{E}^{(+)}_{j}(\mathbf{r}^{\prime},t^{\prime})}$
in terms of the IO-coefficients, we begin by expanding the electric field
operators in the expectation value according to Eq. (7):

$\expectationvalue{\hat{E}^{(-)}_{i}(\mathbf{r},t)\hat{E}^{(+)}_{j}(\mathbf{r}^{\prime},t^{\prime})}=\iint\differential{\omega}\differential{\omega^{\prime}}\expectationvalue{\,\bar{E}_{i}^{(-)}(\mathbf{r},\omega,t)\,\bar{E}_{j}^{(+)}(\mathbf{r}^{\prime},\omega^{\prime},t)}\,e^{i\omega
t}\,e^{-i\omega^{\prime}t^{\prime}}.$ (42)

we then expand each of the amplitude operators using the IO-relation Eq. (13).
To avoid writing out all of the terms that thus appear, we will only remark
that they are each proportional to one of the products:
$\bar{f}_{\alpha}^{(\dagger)}\bar{f}_{\beta}^{(\dagger)}$, where
$\alpha,\beta$ are arbitrary indices. The expectation value of each product is
then taken in the initial state of the system (as we are working in the
Heisenberg picture). In the present case, we assume the initial state to be
the vacuum, which causes the expectation values to evaluate to zero for all
terms except ones proportional to $\bar{f}_{\alpha}\bar{f}_{\beta}^{\dagger}$.
Thus we are left with:

$\expectationvalue{\hat{E}^{(-)}_{i}(\mathbf{r},t)\hat{E}^{(+)}_{j}(\mathbf{r}^{\prime},t^{\prime})}=\iint\differential{\omega}\differential{\omega^{\prime}}\int\differential{\bm{\Xi}}\differential{\bm{\Xi}^{\prime}}\mathcal{A}_{ik}(\mathbf{r},\omega,\bm{\Xi};t)\mathcal{A}^{*}_{jl}(\mathbf{r}^{\prime},\omega^{\prime},\bm{\Xi}^{\prime};t^{\prime})\expectationvalue{\bar{f}_{k}(\bm{\Xi})\bar{f}_{l}^{\dagger}(\bm{\Xi}^{\prime})}\,e^{i\omega
t}\,e^{-i\omega^{\prime}t^{\prime}},$ (43)

where the remaining expectation value can be evaluated using the commutation
relation Eq. (10) to result in
$\expectationvalue{\bar{f}_{k}(\bm{\Xi})\bar{f}_{l}^{\dagger}(\bm{\Xi}^{\prime})}=\delta_{kl}\delta(\bm{\Xi}-\bm{\Xi}^{\prime})\equiv\delta_{kl}\delta(\bm{\xi}-\bm{\xi}^{\prime})\delta(\nu-\nu^{\prime})$.
The $\bm{\Xi}^{\prime}$ integral can then be evaluated to obtain Eq. (14).

If the system was initially in a state other than vacuum, the expectation
values of the various
$\bar{f}_{\alpha}^{(\dagger)}\bar{f}_{\beta}^{(\dagger)}$ products may result
in multiple nonzero terms, which can also be analytically evaluated using Eq.
(10) by first expressing the initial state using the free-field
creation/annihilation operators $\bar{f}_{\alpha}^{(\dagger)}$.

##  Appendix E Derivation of the coupled equations for the filtered IO-
coefficients

To obtain the differential equation governing the evolution of
$\tilde{A}_{ij}(\dots)$ we first take the time derivative of both sides of Eq.
(26b):

$\,\partial_{t}\tilde{A}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)=\int\differential{\omega^{\prime}}\,e^{-i(\omega-\omega^{\prime})t}A_{ij}(\mathbf{r},\omega^{\prime},\bm{\Xi};t).$
(44)

Then we write out $A_{ij}(\mathbf{r},\omega^{\prime},\bm{\Xi};t)$ as the
formal solution of Eq. (15a). This is found by integrating both sides of the
equation over time:

$A_{ij}(\mathbf{r},\omega,\bm{\Xi};t)=\int_{-\infty}^{t}\differential{t^{\prime}}S^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t^{\prime})+\int_{-\infty}^{t}\differential{t^{\prime}}\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime}}F^{*}_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t^{\prime})B_{kj}(\mathbf{r}^{\prime},\omega^{\prime},\bm{\Xi};t^{\prime}),$

where we used the fact that
$A_{ij}(\mathbf{r},\omega,\bm{\Xi};t\rightarrow-\infty)=0$. When we introduce
this form for $A_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ into Eq. (44), we obtain
two main terms which contribute to the equation of motion. The first term we
label as $\tilde{S}^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ and it has the
form:

$\tilde{S}^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)=\int\differential{\omega^{\prime}}\,e^{-i(\omega-\omega^{\prime})t}\int_{-\infty}^{t}\differential{t^{\prime}}S^{(0)}_{ij}(\mathbf{r},\omega^{\prime},\bm{\Xi};t^{\prime}).$

We can immediately rewrite it using the full expression for
$S^{(0)}_{ij}(\mathbf{r},\omega^{\prime},\bm{\Xi};t^{\prime})$, obtained from
Eq. (20):

$\displaystyle\tilde{S}^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)=\frac{2\mathcal{K}}{\pi}\frac{\nu^{2}}{c^{2}}\sqrt{\varepsilon^{\prime\prime}(\bm{\xi},\nu)}\int\differential{\omega^{\prime}}$
$\displaystyle\frac{\omega^{\prime
2}}{c^{2}}\,e^{-i(\omega-\omega^{\prime})t}\int_{-\infty}^{t}\differential{t^{\prime}}\,e^{-i(\omega^{\prime}+\nu)t^{\prime}}$
$\displaystyle\times\int\differential{\mathbf{r}^{\prime}}\chi^{(2)}_{klm}(\mathbf{r}^{\prime})E^{(-)}_{P,k}(\mathbf{r}^{\prime},t^{\prime})\text{Im}\left[G_{il}(\mathbf{r},\mathbf{r}^{\prime},\omega^{\prime})\right]G_{mj}(\mathbf{r}^{\prime},\bm{\xi},\nu).$

Our next aim is to be able to evaluate the time integral in the above
expression. To do that, we expand the pump field into a Fourier integral:
$E_{P,k}^{(-)}(\mathbf{r},t)=\int\differential{\omega_{p}}\mathcal{E}_{P,k}^{(-)}(\mathbf{r},\omega_{p})\,e^{i\,\omega_{p}t}$,
where $\mathcal{E}_{P}^{(-)}(\mathbf{r},\omega_{p})$ is the pump amplitude in
the frequency domain. The full form of
$\tilde{S}^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ is now:

$\displaystyle\tilde{S}^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)=\frac{2\mathcal{K}}{\pi}\frac{\nu^{2}}{c^{2}}$
$\displaystyle\sqrt{\varepsilon^{\prime\prime}(\bm{\xi},\nu)}\int\differential{\mathbf{r}^{\prime}}\chi^{(2)}_{klm}(\mathbf{r}^{\prime})G_{mj}(\mathbf{r}^{\prime},\bm{\xi},\nu)$
(45)
$\displaystyle\times\int\differential{\omega^{\prime}}\frac{\omega^{\prime
2}}{c^{2}}\text{Im}\left[G_{il}(\mathbf{r},\mathbf{r}^{\prime},\omega^{\prime})\right]\int\differential{\omega_{p}}\mathcal{E}^{(-)}_{P,k}(\mathbf{r}^{\prime},\omega_{p})\,e^{-i(\omega-\omega^{\prime})t}\int_{-\infty}^{t}\differential{t^{\prime}}\,e^{i(\omega_{p}-\omega^{\prime}-\nu)t^{\prime}},$

Now we examine the time integral
$\int_{-\infty}^{t}\differential{t^{\prime}}\,e^{i(\omega_{p}-\omega^{\prime}-\nu)t^{\prime}}$
and introduce a substitution $t^{\prime}=t-\tau$ which allows us to rewrite it
as
$\,e^{i(\omega_{p}-\omega^{\prime}-\nu)t}\int_{0}^{\infty}\differential{\tau}\,e^{i(\omega^{\prime}+\nu-\omega_{p})\tau}\equiv\,e^{i(\omega_{p}-\omega^{\prime}-\nu)t}\zeta(\omega^{\prime}+\nu-\omega_{p})$,
where $\zeta(\omega)$ is a generalised function proportional to the Fourier
transform of the Heaviside step function and is closely related to the
analytical properties of the GF [80, 81]. One notable analytical property of
the GF is [81, 46]:

$\int\differential{\omega}\omega^{2}\text{Im}\left[G_{ij}(\mathbf{r},\mathbf{r}^{\prime},\omega)\right]\zeta(\omega-\omega_{0})=i\pi\omega_{0}^{2}G^{*}_{ij}(\mathbf{r},\mathbf{r}^{\prime},\omega_{0})$
(46)

which we can use to further simplify Eq. (45) by performing the integration
over $\omega^{\prime}$. Thus we arrive at the final form for the filtered
source term (29).

The second term appearing on the r.h.s. of Eq. (44) is:

$\int\differential{\omega^{\prime}}\,e^{-i(\omega-\omega^{\prime})t}\int_{-\infty}^{t}\differential{t^{\prime}}\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime\prime}}F^{*}_{ik}(\mathbf{r},\omega^{\prime};\mathbf{r}^{\prime},\omega^{\prime\prime};t^{\prime})B_{kj}(\mathbf{r}^{\prime},\omega^{\prime\prime},\bm{\Xi};t^{\prime}).$

If we write out the full form of
$F^{*}_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t^{\prime})$
using Eq. (12), we have:

$\displaystyle-\frac{2i}{\pi}\int\differential{\omega^{\prime}}\,e^{-i(\omega-\omega^{\prime})t}\int_{-\infty}^{t}\differential{t^{\prime}}\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime\prime}}\frac{\omega^{\prime
2}}{c^{2}}\chi^{(2)}_{lmk}(\mathbf{r}^{\prime})$
$\displaystyle\text{Im}\left[G_{im}(\mathbf{r},\mathbf{r}^{\prime},\omega^{\prime})\right]$
$\displaystyle\times\left(\int\differential{\omega_{p}}\mathcal{E}_{P,l}^{(-)}(\mathbf{r}^{\prime},\omega_{p})\,e^{i\,\omega_{p}t^{\prime}}\right)\,e^{-i(\omega^{\prime}+\omega^{\prime\prime})t^{\prime}}B_{kj}(\mathbf{r}^{\prime},\omega^{\prime\prime},\bm{\Xi};t^{\prime}),$

where we immediately expanded the pump field into a Fourier integral. After
some reordering of the integrals in the above expression, we can identify:

$\int_{-\infty}^{t}\differential{t^{\prime}}\int\differential{\omega^{\prime\prime}}\,e^{i(\omega_{p}-\omega^{\prime}-\omega^{\prime\prime})t^{\prime}}B_{kj}(\mathbf{r}^{\prime},\omega^{\prime\prime},\bm{\Xi};t^{\prime})=\tilde{B}_{kj}(\mathbf{r}^{\prime},\omega_{p}-\omega^{\prime},\bm{\Xi};t),$

where the r.h.s. was established according to Eq. (26c). Thus, we now write
the term in question as:

$-\frac{2i}{\pi}\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime}}\frac{\omega^{\prime
2}}{c^{2}}\chi^{(2)}_{lmk}(\mathbf{r}^{\prime})\text{Im}\left[G_{im}(\mathbf{r},\mathbf{r}^{\prime},\omega^{\prime})\right]\,e^{-i(\omega-\omega^{\prime})t}\int\differential{\omega_{p}}\mathcal{E}^{(-)}_{P,l}(\mathbf{r}^{\prime},\omega_{p})\tilde{B}_{kj}(\mathbf{r}^{\prime},\omega_{p}-\omega^{\prime},\bm{\Xi};t).$

To arrive at the form of the term given by Eqs. (28a) and (30), we first
introduce a variable substitution $\omega_{p}=\omega^{\prime}+\bar{\omega}$
into the above expression:

$-\frac{2i}{\pi}\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime}}\frac{\omega^{\prime
2}}{c^{2}}\chi^{(2)}_{lmk}(\mathbf{r}^{\prime})\text{Im}\left[G_{im}(\mathbf{r},\mathbf{r}^{\prime},\omega^{\prime})\right]\,e^{-i(\omega-\omega^{\prime})t}\int_{-\omega^{\prime}}^{\infty}\differential{\bar{\omega}}\mathcal{E}^{(-)}_{P,l}(\mathbf{r}^{\prime},\bar{\omega}+\omega^{\prime})\tilde{B}_{kj}(\mathbf{r}^{\prime},\bar{\omega},\bm{\Xi};t).$

The newly introduced $\bar{\omega}$ has the limits
$[-\omega^{\prime},\infty)$, as shown above, which indicates that the filtered
coefficients have to be evaluated at certain negative frequency values during
calculation. Although this is also formally allowed by the definitions in Eqs.
(26), these negative-frequency contributions can be safely neglected and the
lower integration limit set to $0$ without any loss of accuracy. The reason
for this can be seen by examining the source term in (29), where negative
values of the variable $\omega$ would cause the term to oscillate at optical
frequencies, thus averaging to $0$ on the time scales of the SPDC process and
resulting in these negative-frequency values of the IO-coefficients not
contributing in the final equations of motion. Finally, to make the resulting
expression consistent with the notation used in the main text, we will
exchange the variables $\omega^{\prime}$ and $\bar{\omega}$ and write the
final form of the second term of Eq. (44):

$-\frac{2i}{\pi}\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime}}\left(\int\differential{\bar{\omega}}\frac{\bar{\omega}^{2}}{c^{2}}\chi^{(2)}_{lmk}(\mathbf{r}^{\prime})\text{Im}\left[G_{im}(\mathbf{r},\mathbf{r}^{\prime},\bar{\omega})\right]\,e^{-i(\omega-\bar{\omega})t}\mathcal{E}^{(-)}_{P,l}(\mathbf{r}^{\prime},\omega^{\prime}+\bar{\omega})\right)\tilde{B}_{kj}(\mathbf{r}^{\prime},\omega^{\prime},\bm{\Xi};t),$

where the expression in the parentheses exactly corresponds to the definition
Eq. (30) and, when the above expression is combined with Eq. (29), we get
exactly the equation of motion Eq. (28a). An analogous procedure can be
performed to find the equation of motion for the coefficient
$\tilde{B}_{kj}(\mathbf{r},\omega,\bm{\Xi};t)$, after which we obtain the
coupled equations (28).

##  Appendix F Low-gain solution

To test the validity of the coupled equations Eqs. (28), we find the
perturbative solutions for the IO-coefficients in the the low-gain regime,
when the pump amplitude is sufficiently weak, and use them to calculate the
low-gain single-photon spectrum.

The single-photon output spectrum in terms of the IO-coefficients is given in
Eq. (27), where we see that it is entirely determined by
$\tilde{A}_{kj}(\mathbf{r},\omega,\bm{\Xi};t)$. Thus, an analytical expression
for the spectrum in the low-gain regime can be obtained by finding a first-
order perturbative solution for the IO-coefficient
$\tilde{A}_{kj}(\mathbf{r},\omega,\bm{\Xi};t)$ and replacing it in Eq. (27).

To obtain a perturbative expansion of
$\tilde{A}_{kj}(\mathbf{r},\omega,\bm{\Xi};t)$, we begin by integrating both
sides of Eq. (28a) over time, to find the formal form of the solution. This
yields:

$\tilde{A}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)=\int_{-\infty}^{t}\differential{t^{\prime}}\tilde{S}^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t^{\prime})+\int_{-\infty}^{t}\differential{t^{\prime}}\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime}}\tilde{F}^{*}_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t^{\prime})\tilde{B}_{kj}(\mathbf{r}^{\prime},\omega^{\prime},\bm{\Xi};t^{\prime}).$
(47)

Next, we find the formal solution of Eq. (28b) in the same manner:

$\tilde{B}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)=\int_{-\infty}^{t}\differential{t^{\prime}}\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime}}\tilde{F}_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t^{\prime})\tilde{A}_{kj}(\mathbf{r}^{\prime},\omega^{\prime},\bm{\Xi};t^{\prime}),$

and insert it into (47):

$\displaystyle\tilde{A}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)=\int_{-\infty}^{t}\differential{t^{\prime}}\tilde{S}^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t^{\prime})+$
(48)
$\displaystyle\qquad\int_{-\infty}^{t}\differential{t^{\prime}}\int\differential{\mathbf{r}^{\prime}}\int\differential{\omega^{\prime}}\tilde{F}^{*}_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t^{\prime})\left(\int_{-\infty}^{t^{\prime}}\differential{t^{\prime\prime}}\int\differential{\mathbf{r}^{\prime\prime}}\int\differential{\omega^{\prime\prime}}\tilde{F}_{kl}(\mathbf{r}^{\prime},\omega^{\prime};\mathbf{r}^{\prime\prime},\omega^{\prime\prime};t^{\prime\prime})\tilde{A}_{lj}(\mathbf{r}^{\prime\prime},\omega^{\prime\prime},\bm{\Xi};t^{\prime\prime})\right).$
(49)

The full perturbative expansion is obtained by iteratively expanding the above
expression using the formal solutions of Eqs. (28). However, as we are
interested in the first-order solution for
$\tilde{A}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$, Eq. (48) is already
sufficient. This can be confirmed by noting that both
$\tilde{S}^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ and
$\tilde{F}_{ik}(\mathbf{r},\omega;\mathbf{r}^{\prime},\omega^{\prime};t)$ are
proportional to the pump amplitude, according to Eqs. (29) and (30),
respectively, and any further expansion will result in terms of third- and
higher-order in the pump amplitude.

Thus we conclude that the first-order analytical expression for
$\tilde{A}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$, is just the source term
$\tilde{S}^{(0)}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ integrated over time.
Using Eq. (29) we find it has the form:

$\tilde{A}_{ij}(\mathbf{r},\omega,\bm{\Xi};t\rightarrow\infty)=4i\pi\mathcal{K}\,\frac{\nu^{2}\omega^{2}}{c^{2}}\sqrt{\varepsilon^{\prime\prime}(\bm{\xi},\nu)}\int\differential{\mathbf{r}^{\prime}}\chi^{(2)}_{klm}(\mathbf{r}^{\prime})\mathcal{E}^{(-)}_{P,k}(\mathbf{r}^{\prime},\omega+\nu)G^{*}_{il}(\mathbf{r},\mathbf{r}^{\prime},\omega)G_{mj}(\mathbf{r}^{\prime},\bm{\xi},\nu).$
(50)

Now we insert the above expression into (27) to obtain the spectrum in the
low-gain regime. The resulting expression is rather long but, through some
straightforward algebra and an application of Eq. (34), it can be reduced to
obtain:

$\displaystyle\sigma(\omega_{0},t\rightarrow\infty)\propto\frac{\omega_{0}^{4}}{c^{4}}\iint\differential{\mathbf{r}^{\prime}}\differential{\mathbf{r}^{\prime\prime}}$
$\displaystyle\chi^{(2)}_{klm}(\mathbf{r}^{\prime})\chi^{(2)}_{qrs}(\mathbf{r}^{\prime\prime})G_{il}^{*}(\mathbf{r}_{0},\mathbf{r}^{\prime},\omega_{0})G_{ir}(\mathbf{r}_{0},\mathbf{r}^{\prime\prime},\omega_{0})$
$\displaystyle\times$
$\displaystyle\int\differential{\nu}\frac{\nu^{2}}{c^{2}}\mathcal{E}^{(-)}_{P,k}(\mathbf{r}^{\prime},\omega_{0}+\nu)\mathcal{E}^{(+)}_{P,q}(\mathbf{r}^{\prime\prime},\omega_{0}+\nu)\text{Im}\left[G_{ms}(\mathbf{r}^{\prime},\mathbf{r}^{\prime\prime},\nu)\right],$
(51)

##  Appendix G Notes on numerical implementation

In both the time- and frequency-domain formulations, the relatively high
dimensionality of the IO-coefficients
$\tilde{A}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ and
$\tilde{B}_{ij}(\mathbf{r},\omega,\bm{\Xi};t)$ requires a careful approach to
discretising the relevant ranges for each of the variables present in order to
minimise the computational resources required for solving the coupled
equations Eq. (28). In the case of the QSUP scheme investigated in this work,
the one-dimensional nature of the waveguide and the homogeneity of the
waveguide dielectric result in a phase matching function that is well-
approximated by $\text{sinc}^{2}\left(\Delta
k(\omega_{p},\omega_{s},\omega_{i})\frac{L}{2}\right)$, where $\omega_{p}$,
$\omega_{s}$ and $\omega_{i}$ are the pump, signal and idler frequencies,
respectively and $\Delta k(\omega_{p},\omega_{s},\omega_{i})$ is the phase
mismatch. With this in mind, the width of the frequency range of interest
$\Delta\omega$ was chosen to include frequencies around the central and first-
order peaks of the sinc-function for both the signal and idler photons. The
spatial discretisation step $\delta z$ was then chosen to satisfy $\delta
z\ll\frac{L}{2\pi}$ to ensure that phase-matching is correctly reproduced.

The ranges and discretisation steps for the variables $\bm{\xi}$ and $\nu$ are
determined in a more involved manner, but one that is entirely defined by the
properties of the nonlinear waveguide and the behaviour of its GF. As noted in
Sec. II.1, $\bm{\xi}$ and $\nu$ serve to index the initial-time field-matter
excitations taking part in the IO-relations (18) and (25) with the
contributions from each excitation being weighted by the IO-coefficients. In
theory, all possible values for these variables should be taken into account
in order to obtain the exact output fields, however, as will now be discussed,
the relevant contributions come from specific limited ranges determined by the
spatial and spectral properties of the structure in question, as well as the
nonlinear interaction.

To determine which range of the spatial variable $\bm{\xi}$ gives a dominant
contribution to the output we examine the low-gain analytical expression for
the coefficient $\tilde{A}(z,\omega,\bm{\Xi};t\rightarrow\infty)$ (since it is
the only one contributing to the spectrum) given by simplifying Eq. (50)
according to the assumptions made about the waveguide under consideration in
Section III:

$\displaystyle\tilde{A}(z,\omega,\bm{\Xi};t\rightarrow\infty)$
$\displaystyle=$ $\displaystyle
2i\mathcal{K}\,\frac{\nu^{2}\omega^{2}}{c^{2}}\sqrt{\varepsilon^{\prime\prime}(\nu)}$
(52)
$\displaystyle\times\int\differential{z^{\prime}}\chi^{(2)}(z^{\prime})\mathcal{E}^{(-)}_{P,}(z^{\prime},\omega+\nu)\qquad$
$\displaystyle\times G^{*}(z,z^{\prime},\omega)G(z^{\prime},\xi,\nu).\qquad$

The only $\xi$-dependent factor present - $G(z^{\prime},\xi,\nu)$, in the case
of the GF given by (31), is proportional to
$\,e^{i\frac{\nu}{c}\sqrt{\varepsilon(\nu)}\absolutevalue{z^{\prime}-\xi}}$.
Due to $\varepsilon(\nu)$ being complex, this factor represents an oscillating
function with an exponentially decreasing amplitude in
$\absolutevalue{z^{\prime}-\xi}$ with the decay length
$l_{d}(\nu)=\frac{1}{\frac{\nu}{c}\text{Im}\left[\sqrt{\varepsilon(\nu)}\right]}$.
Since the variable $z^{\prime}$ is always confined to the nonlinear region of
the waveguide, the effective range of the variable $\xi$ consists of the
nonlinear region itself and a multiple of the length $l_{d}$ of the linear
regions surrounding the nonlinear waveguide. Physically this means that field-
matter excitations sufficiently far away from the nonlinear region are unable
to influence the nonlinear interaction due to their effects being damped away
by the absorption of the dielectric.

In practice, it is sufficient to find the largest possible decay length for
the frequency range of interest $l_{d,max}$ and fix that as the boundary value
for $\xi$. For materials with small amounts of loss or where the loss is very
narrow-band, these boundary values can be extremely large and could present a
problem in a numerical implementation, where memory limitations are a factor.

The potentially large range of $\xi$ outside of the nonlinear region can be
made more manageable via a variable transformation. We used the double-
exponential transformation of the form [82]:
$\xi=\sinh(\frac{\pi}{2}\sinh{\theta})$ where $\theta$ is the substitution
variable whose range is chosen to correspond to values of $\xi$ outside the
nonlinear region. With such a transformation, we compress the potentially long
range of integration for the variable $\xi$ in the surrounding linear regions
into a much shorter range for the variable $\theta$. This transformation
resembles the coordinate transformation used in implementing perfectly matched
layers (PMLs) in computational photonics [83].

Even though the above conclusions were made using the low-gain result, they
also hold in the high-gain regime as well, as was confirmed by our
simulations. This is due to the spatial/absorbing properties of the structure
being independent of the pump intensity and thus, the spatial distribution of
the initial-time field-matter excitations contributing to the output stays the
same, regardless of gain.

On the other hand, the relevant frequency range for the contributing initial-
time field-matter excitations, i.e. the relevant range of $\nu$, cannot simply
be inferred from the low-gain, due to the seeding effect that couples
neighbouring frequencies in the high-gain regime. Looking again at Eq. (F), we
observe the following in the low-gain regime: the range of frequencies $\nu$
contributing to the output at frequency $\omega_{0}$ is determined by the
spectrum of the pump, i.e. output photons of frequency $\omega_{0}$ are
influenced by field-matter excitations of frequencies within a pump bandwidth
of $\nu_{0}=\omega_{p0}-\omega_{0}$. However, as gain is increasing, the
generated field at frequencies neighbouring $\omega_{0}$ will start
contributing to the intensity at $\omega_{0}$ through higher order seeding
effects, where those neighbouring frequencies themselves are affected by their
own respective neighbouring frequencies further from $\omega_{0}$. Hence, to
have an accurate representation of the intensity at $\omega_{0}$, one has to
include the whole range of frequencies with substantial intensities around
$\omega_{0}$, which are mainly determined by the phase-matching condition.
Hence, in our calculation, for both $\omega$ and $\nu$ that appear in Eqs.
(28), we consider the whole range of frequencies around the main and the first
surrounding peaks of the phase-matching function.

## References

  * O’brien  _et al._ [2009] J. L. O’brien, A. Furusawa, and J. Vučković, Photonic quantum technologies, Nature Photonics 3, 687 (2009).
  * Wang  _et al._ [2020] J. Wang, F. Sciarrino, A. Laing, and M. G. Thompson, Integrated photonic quantum technologies, Nature Photonics 14, 273 (2020).
  * Boyd [2020] R. W. Boyd,  _Nonlinear optics_ (Academic press, 2020).
  * Caspani  _et al._ [2017] L. Caspani, C. Xiong, B. J. Eggleton, D. Bajoni, M. Liscidini, M. Galli, R. Morandotti, and D. J. Moss, Integrated sources of photon quantum states based on nonlinear optics, Light: Science & Applications 6, e17100 (2017).
  * Elshaari  _et al._ [2020] A. W. Elshaari, W. Pernice, K. Srinivasan, O. Benson, and V. Zwiller, Hybrid integrated quantum photonic circuits, Nature Photonics 14, 285 (2020).
  * Anwar  _et al._ [2021] A. Anwar, C. Perumangatt, F. Steinlechner, T. Jennewein, and A. Ling, Entangled photon-pair sources based on three-wave mixing in bulk crystals, Review of Scientific Instruments 92, 041101 (2021).
  * Eisaman  _et al._ [2011] M. D. Eisaman, J. Fan, A. Migdall, and S. V. Polyakov, Invited review article: Single-photon sources and detectors, Review of Scientific Instruments 82, 071101 (2011).
  * Meyer-Scott  _et al._ [2020] E. Meyer-Scott, C. Silberhorn, and A. Migdall, Single-photon sources: Approaching the ideal through multiplexing, Review of Scientific Instruments 91, 041101 (2020).
  * Quesada  _et al._ [2022] N. Quesada, L. G. Helt, M. Menotti, M. Liscidini, and J. E. Sipe, Beyond photon pairs—nonlinear quantum photonics in the high-gain regime: a tutorial, Advances in Optics and Photonics 14, 291 (2022).
  * Menicucci  _et al._ [2006] N. C. Menicucci, P. van Loock, M. Gu, C. Weedbrook, T. C. Ralph, and M. A. Nielsen, Universal quantum computation with continuous-variable cluster states, Physical Review Letters 97, 110501 (2006).
  * Larsen  _et al._ [2021] M. V. Larsen, C. Chamberland, K. Noh, J. S. Neergaard-Nielsen, and U. L. Andersen, Fault-tolerant continuous-variable measurement-based quantum computation architecture, PRX Quantum 2, 030325 (2021).
  * Zhong  _et al._ [2020] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, _et al._ , Quantum computational advantage using photons, Science 370, 1460 (2020).
  * 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).
  * van Loock [2002] P. van Loock, Quantum communication with continuous variables, Fortschritte der Physik: Progress of Physics 50, 1177 (2002).
  * Cutipa and Chekhova [2022] P. Cutipa and M. V. Chekhova, Bright squeezed vacuum for two-photon spectroscopy: simultaneously high resolution in time and frequency, space and wavevector, Optics Letters 47, 465 (2022).
  * Panahiyan  _et al._ [2022] S. Panahiyan, C. S. Muñoz, M. V. Chekhova, and F. Schlawin, Two-photon-absorption measurements in the presence of single-photon losses, Physical Review A 106, 043706 (2022).
  * Lawrie  _et al._ [2019] B. J. Lawrie, P. D. Lett, A. M. Marino, and R. C. Pooser, Quantum sensing with squeezed light, ACS Photonics 6, 1307 (2019).
  * Engelkemeier  _et al._ [2020] M. Engelkemeier, L. Lorz, S. De, B. Brecht, I. Dhand, M. B. Plenio, C. Silberhorn, and J. Sperling, Quantum photonics with active feedback loops, Physical Review A 102, 023712 (2020).
  * Engelkemeier  _et al._ [2021] M. Engelkemeier, J. Sperling, J. Tiedau, S. Barkhofen, I. Dhand, M. Plenio, B. Brecht, and C. Silberhorn, Climbing the Fock ladder: Advancing multiphoton state generation, arXiv preprint arXiv:2105.03720 (2021).
  * Dell’Anno  _et al._ [2006] F. Dell’Anno, S. De Siena, and F. Illuminati, Multiphoton quantum optics and quantum state engineering, Physics Reports 428, 53 (2006).
  * Holland and Burnett [1993] M. J. Holland and K. Burnett, Interferometric detection of optical phase shifts at the Heisenberg limit, Physical Review Letters 71, 1355 (1993).
  * Sperling  _et al._ [2019] J. Sperling, A. Perez-Leija, K. Busch, and C. Silberhorn, Mode-independent quantum entanglement for light, Physical Review A 100, 062129 (2019).
  * 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).
  * 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).
  * Nehra  _et al._ [2022] R. Nehra, R. Sekine, L. Ledezma, Q. Guo, R. M. Gray, A. Roy, and A. Marandi, Few-cycle vacuum squeezing in nanophotonics, Science 377, 1333 (2022).
  * Chen  _et al._ [2022] P.-K. Chen, I. Briggs, S. Hou, and L. Fan, Ultra-broadband quadrature squeezing with thin-film lithium niobate nanophotonics, Optics Letters 47, 1506 (2022).
  * Peace  _et al._ [2022] D. Peace, A. Zappacosta, R. Cernansky, B. Haylock, A. Boes, A. Mitchell, and M. Lobino, Picosecond pulsed squeezing in thin-film lithium niobate strip-loaded waveguides at telecommunication wavelengths, Journal of Physics: Photonics 4, 035002 (2022).
  * Stokowski  _et al._ [2022] H. S. Stokowski, T. P. McKenna, T. Park, A. Y. Hwang, D. J. Dean, O. T. Celik, V. Ansari, M. M. Fejer, and A. H. Safavi-Naeini, Integrated quantum optical phase sensor, arXiv preprint arXiv:2212.09717 (2022).
  * Busch  _et al._ [2007] K. Busch, G. Von Freymann, S. Linden, S. Mingaleev, L. Tkeshelashvili, and M. Wegener, Periodic nanostructures for photonics, Physics Reports 444, 101 (2007).
  * Lalanne  _et al._ [2018] P. Lalanne, W. Yan, K. Vynck, C. Sauvan, and J.-P. Hugonin, Light interaction with photonic and plasmonic resonances, Laser & Photonics Reviews 12, 1700113 (2018).
  * Wasilewski  _et al._ [2006] W. Wasilewski, A. I. Lvovsky, K. Banaszek, and C. Radzewicz, Pulsed squeezed light: Simultaneous squeezing of multiple modes, Physical Review A 73, 063819 (2006).
  * Christ  _et al._ [2013] A. Christ, B. Brecht, W. Mauerer, and C. Silberhorn, Theory of quantum frequency conversion and type-II parametric down-conversion in the high-gain regime, New Journal of Physics 15, 053038 (2013).
  * Sharapova  _et al._ [2020] P. R. Sharapova, G. Frascella, M. Riabinin, A. M. Pérez, O. V. Tikhonova, S. Lemieux, R. W. Boyd, G. Leuchs, and M. V. Chekhova, Properties of bright squeezed vacuum at increasing brightness, Physical Review Research 2, 013371 (2020).
  * Helt and Quesada [2020] L. Helt and N. Quesada, Degenerate squeezing in waveguides: a unified theoretical approach, Journal of Physics: Photonics 2, 035001 (2020).
  * Gigli  _et al._ [2020] C. Gigli, T. Wu, G. Marino, A. Borne, G. Leo, and P. Lalanne, Quasinormal-mode non-hermitian modeling and design in nonlinear nano-optics, ACS photonics 7, 1197 (2020).
  * Lecamp  _et al._ [2007] G. Lecamp, J.-P. Hugonin, and P. Lalanne, Theoretical and computational concepts for periodic optical waveguides, Optics Express 15, 11042 (2007).
  * Saravi  _et al._ [2015] S. Saravi, S. Diziain, M. Zilk, F. Setzpfandt, and T. Pertsch, Phase-matched second-harmonic generation in slow-light photonic crystal waveguides, Physical Review A 92, 063821 (2015).
  * Kumar  _et al._ [2020] P. Kumar, S. Saravi, T. Pertsch, and F. Setzpfandt, Integrated induced-coherence spectroscopy in a single nonlinear waveguide, Physical Review A 101, 053860 (2020).
  * Triginer  _et al._ [2020] G. Triginer, M. D. Vidrighin, N. Quesada, A. Eckstein, M. Moore, W. S. Kolthammer, J. E. Sipe, and I. A. Walmsley, Understanding high-gain twin-beam sources using cascaded stimulated emission, Physical Review X 10, 031063 (2020).
  * Grgić  _et al._ [2012] J. Grgić, J. R. Ott, F. Wang, O. Sigmund, A.-P. Jauho, J. Mørk, and N. A. Mortensen, Fundamental limitations to gain enhancement in periodic media and waveguides, Physical Review Letters 108, 183903 (2012).
  * Saravi  _et al._ [2016] S. Saravi, R. Quintero-Bermudez, F. Setzpfandt, N. A. Mortensen, and T. Pertsch, Effect of loss on slow-light-enhanced second-harmonic generation in periodic nanostructures, Optics Letters 41, 3110 (2016).
  * Vogel and Welsch [2006] W. Vogel and D.-G. Welsch,  _Quantum optics_ (John Wiley & Sons, 2006).
  * Gruner and Welsch [1996] T. Gruner and D.-G. Welsch, Green-function approach to the radiation-field quantization for homogeneous and inhomogeneous Kramers-Kronig dielectrics, Physical Review A 53, 1818 (1996).
  * Dung  _et al._ [1998] H. T. Dung, L. Knöll, and D.-G. Welsch, Three-dimensional quantization of the electromagnetic field in dispersive and absorbing inhomogeneous dielectrics, Physical Review A 57, 3931 (1998).
  * Knöll  _et al._ [2003] L. Knöll, S. Scheel, and D.-G. Welsch, QED in dispersing and absorbing media (2003), arXiv:quant-ph/0006121 [quant-ph] .
  * Poddubny  _et al._ [2016] A. N. Poddubny, I. V. Iorsh, and A. A. Sukhorukov, Generation of photon-plasmon quantum states in nonlinear hyperbolic metamaterials, Physical Review Letters 117, 123901 (2016).
  * Olekhno  _et al._ [2018] N. Olekhno, M. Petrov, and I. Iorsh, Spontaneous parametric downconversion of light by a dielectric nanoparticle, in  _Journal of Physics: Conference Series_ , Vol. 993 (IOP Publishing, 2018) p. 012022.
  * Nikolaeva  _et al._ [2021] A. Nikolaeva, K. Frizyuk, N. Olekhno, A. Solntsev, and M. Petrov, Directional emission of down-converted photons from a dielectric nanoresonator, Physical Review A 103, 043703 (2021).
  * Poddubny and Smirnova [2018] A. N. Poddubny and D. A. Smirnova, Nonlinear generation of quantum-entangled photons from high-Q states in dielectric nanoparticles, arXiv preprint arXiv:1808.04811 (2018).
  * Saravi  _et al._ [2017] S. Saravi, A. N. Poddubny, T. Pertsch, F. Setzpfandt, and A. A. Sukhorukov, Atom-mediated spontaneous parametric down-conversion in periodic waveguides, Optics Letters 42, 4724 (2017).
  * Santos  _et al._ [2022] E. A. Santos, T. Pertsch, F. Setzpfandt, and S. Saravi, Subdiffraction quantum imaging with undetected photons, Physical Review Letters 128, 173601 (2022).
  * Santiago-Cruz  _et al._ [2022] T. Santiago-Cruz, S. D. Gennaro, O. Mitrofanov, S. Addamane, J. Reno, I. Brener, and M. V. Chekhova, Resonant metasurfaces for generating complex quantum states, Science 377, 991 (2022).
  * Santiago-Cruz  _et al._ [2021] T. Santiago-Cruz, A. Fedotova, V. Sultanov, M. A. Weissflog, D. Arslan, M. Younesi, T. Pertsch, I. Staude, F. Setzpfandt, and M. Chekhova, Photon pairs from resonant metasurfaces, Nano Letters 21, 4423 (2021).
  * Zhang  _et al._ [2022] J. Zhang, J. Ma, M. Parry, M. Cai, R. Camacho-Morales, L. Xu, D. N. Neshev, and A. A. Sukhorukov, Spatially entangled photon pairs from lithium niobate nonlocal metasurfaces, Science Advances 8, eabq4240 (2022).
  * Suttorp [2007] L. Suttorp, Field quantization in inhomogeneous anisotropic dielectrics with spatio-temporal dispersion, Journal of Physics A: Mathematical and Theoretical 40, 3697 (2007).
  * Sharapova  _et al._ [2015] P. Sharapova, A. M. Pérez, O. V. Tikhonova, and M. V. Chekhova, Schmidt modes in the angular spectrum of bright squeezed vacuum, Physical Review A 91, 043816 (2015).
  * Yanagimoto  _et al._ [2022] R. Yanagimoto, E. Ng, A. Yamamura, T. Onodera, L. G. Wright, M. Jankowski, M. Fejer, P. L. McMahon, and H. Mabuchi, Onset of non-Gaussian quantum physics in pulsed squeezing with mesoscopic fields, Optica 9, 379 (2022).
  * Xing and Ralph [2023] W. Xing and T. C. Ralph, Pump depletion in optical parametric amplification, Physical Review A 107, 023712 (2023).
  * Quesada  _et al._ [2020] N. Quesada, G. Triginer, M. D. Vidrighin, and J. E. Sipe, Theory of high-gain twin-beam generation in waveguides: From Maxwell’s equations to efficient simulation, Physical Review A 102, 033519 (2020).
  * Scheel and Welsch [2006a] S. Scheel and D.-G. Welsch, Quantum theory of light and noise polarization in nonlinear optics, Physical Review Letters 96, 073601 (2006a).
  * Scheel and Welsch [2006b] S. Scheel and D.-G. Welsch, Causal nonlinear quantum optics, Journal of Physics B: Atomic, Molecular and Optical Physics 39, S711 (2006b).
  * Lindel  _et al._ [2020] F. Lindel, R. Bennett, and S. Y. Buhmann, Theory of polaritonic quantum-vacuum detection, Physical Review A 102, 041701(R) (2020).
  * Lindel  _et al._ [2021] F. Lindel, R. Bennett, and S. Y. Buhmann, Macroscopic quantum electrodynamics approach to nonlinear optics and application to polaritonic quantum-vacuum detection, Physical Review A 103, 033705 (2021).
  * Dung  _et al._ [2003] H. T. Dung, S. Y. Buhmann, L. Knöll, D.-G. Welsch, S. Scheel, and J. Kästel, Electromagnetic-field quantization and spontaneous decay in left-handed media, Physical Review A 68, 043816 (2003).
  * Serafini [2017] A. Serafini,  _Quantum continuous variables: a primer of theoretical methods_ (CRC press, 2017).
  * Weissflog  _et al._ [2023] M. A. Weissflog, R. Dezert, V. Vinel, C. Gigli, G. Leo, T. Pertsch, F. Setzpfandt, A. Borne, and S. Saravi, Nonlinear Nanoresonators are Natural Sources of Bell States, arXiv preprint arXiv:2305.19362 (2023).
  * Cresser [1983] J. D. Cresser, Theory of the spectrum of the quantised light field, Physics Reports 94, 47 (1983).
  * Eberly and Wodkiewicz [1977] J. H. Eberly and K. Wodkiewicz, The time-dependent physical spectrum of light, JOSA 67, 1252 (1977).
  * Kalashnikov  _et al._ [2016] D. A. Kalashnikov, A. V. Paterova, S. P. Kulik, and L. A. Krivitsky, Infrared spectroscopy with visible light, Nature Photonics 10, 98 (2016).
  * Lee  _et al._ [2019] S. K. Lee, T. H. Yoon, and M. Cho, Interferometric quantum spectroscopy with undetected photons via distinguishability modulation, Optics Express 27, 14853 (2019).
  * Lemos  _et al._ [2022] G. B. Lemos, M. Lahiri, S. Ramelow, R. Lapkiewicz, and W. N. Plick, Quantum imaging and metrology with undetected photons: tutorial, JOSA B 39, 2200 (2022).
  * Kolobov  _et al._ [2017] M. I. Kolobov, E. Giese, S. Lemieux, R. Fickler, and R. W. Boyd, Controlling induced coherence for quantum imaging, Journal of Optics 19, 054003 (2017).
  * Giese  _et al._ [2017] E. Giese, S. Lemieux, M. Manceau, R. Fickler, and R. W. Boyd, Phase sensitivity of gain-unbalanced nonlinear interferometers, Physical Review A 96, 053863 (2017).
  * Machado  _et al._ [2020] G. J. Machado, G. Frascella, J. P. Torres, and M. V. Chekhova, Optical coherence tomography with a nonlinear interferometer in the high parametric gain regime, Applied Physics Letters 117, 094002 (2020).
  * Saleh and Teich [2019] B. E. Saleh and M. C. Teich,  _Fundamentals of photonics_ (John Wiley & Sons, 2019).
  * Chen  _et al._ [2005] J. Chen, X. Li, and P. Kumar, Two-photon-state generation via four-wave mixing in optical fibers, Physical Review A 72, 033801 (2005).
  * Zaske  _et al._ [2012] S. Zaske, A. Lenhard, C. A. Keßler, J. Kettler, C. Hepp, C. Arend, R. Albrecht, W.-M. Schulz, M. Jetter, P. Michler, _et al._ , Visible-to-telecom quantum frequency conversion of light from a single quantum emitter, Physical Review Letters 109, 147404 (2012).
  * Han  _et al._ [2021] X. Han, W. Fu, C.-L. Zou, L. Jiang, and H. X. Tang, Microwave-optical quantum frequency conversion, Optica 8, 1050 (2021).
  * Feynman [1951] R. P. Feynman, An operator calculus having applications in quantum electrodynamics, Phys. Rev. 84, 108 (1951).
  * Singh [2020] P. Singh, Novel generalized Fourier representations and phase transforms, Digital signal processing 106, 102830 (2020).
  * Saravi [2018] S. Saravi,  _Photon-pair generation in photonic crystal waveguides_ , Ph.D. thesis, Jena (2018), dissertation, Friedrich-Schiller-Universität Jena, 2018\.
  * Mori and Sugihara [2001] M. Mori and M. Sugihara, The double-exponential transformation in numerical analysis, Journal of Computational and Applied Mathematics 127, 287 (2001).
  * Sauvan  _et al._ [2022] C. Sauvan, T. Wu, R. Zarouf, E. A. Muljarov, and P. Lalanne, Normalization, orthogonality, and completeness of quasinormal modes of open systems: the case of electromagnetism, Optics Express 30, 6846 (2022).