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

We develop a non-perturbative formulation based on the Green-function quantization 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 eﬀects in the performance of the system.


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 highgain 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 [9,34], 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 [9,32,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 [38,50,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 nonperturbative 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 [15,16,38,51].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 ε(r, ω) = ε (r, ω) + iε (r, ω), where ε and ε 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 χ (2) ijk (r).In the Heisenberg picture, the SPDC interaction Hamiltonian has the following form [56]: where Ê(+) j,k (r, t) are Cartesian components of the positive-frequency part of the generated quantum fields, E (−) P,i (r, t) are the Cartesian components of the negativefrequency part of the pump field and ε 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 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, results in a more complex dynamics for generation of non-Gaussian quantum states [57,58].Advancing the GF formalism to go beyond the undepleted pump approximation could be an interesting direction for future works.
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 (−∞, ∞) and [0, ∞), 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 → −∞ and adiabatically "turned off" as t → ∞.
To allow for the treatment of open and dispersive systems, the electric field operator is quantised using the GF quantisation method [42][43][44]: is the amplitude operator of frequency ω, G ij (r, r , ω) are the matrix elements of the dyadic GF of the dielectric system ← → G (r, r , ω) and fj (r , ω, t) are the Cartesian components of the bosonic annihilation operator f (r , ω, t).
The GF in 2 is the solution to the classical Helmholtz equation: together with the condition that it vanishes at infinity.In the above equation ← → I is the identity matrix.The operator f (r , ω, t) annihilates a local field-matter excitation of frequency ω, located at r .Along with its adjoint f † (r , ω, t), it makes up the fundamental operator algebra of the GF quantisation method and obeys the canonical commutation relation: Using Eq. 4, it can be shown that the amplitude operators Ê(±) (r, ω, t) obey the following commutation relation: where Im [G ij (r, r , ω)] 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. [59].

A. Equations of motion
The Cartesian components of Ê(±) (r, ω, t) evolve in time according to the Heisenberg equation of motion: where Ĥ0 (t) = dr dω ω f † i (r, ω, t) fi (r, ω, 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, Ē(±) (r, ω, t), where: The exponential factors account for the free-field component of the evolution and the rotating-frame operators Ē(±) (r, ω, t) evolve according to: 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 f ( †) (r, ω, t) are defined equivalently to Eq. 7a and are related to the rotating frame amplitude operators Ē(±) (r, ω, t) by the relation: and its adjoint for Ē(−) i (r, ω, t) and f † j (r , ω, t).The rotating-frame amplitude operators Ē(±) (r, ω, t), as well as the rotating-frame creation/annihilation operators f ( †) (r , ω, t), can be shown to obey equivalent commutation relations as their Heisenberg picture counterparts: 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 Ē(+) (r, ω, t) from Eq. 8 (the full derivation is given in Appendix C): with the corresponding equation for ∂ t Ē(−) i (r, ω, t) obtained by simply taking the adjoint of Eq. 11.Here we also defined:

B. Input-output relations
The obtained equations of motion Eqs.11 are linear in terms of the amplitude operators Ē(±) i (r, ω, t).
Along with the fact that the Ē(±) i (r, ω, t) are themselves canonical, i.e., their commutator Eq. 10b is scalar, this enables us to claim that, at an arbitrary time t, Ē(±) i (r, ω, t) can be expressed as a linear combination of Ē(±) i (r, ω, t → −∞), i.e., the free-field amplitude operators [31,32,60].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 freefield (input) operators via an input-output (IO-) relation.Additionally, due to the linear relation between the amplitude operators Ē(±) i (r, ω, t) and creation/annihilation operators f ( †) (r, ω, 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 f ( †) (r , ω, t → −∞), which we label f ( †) (r , ω) for compactness.We note that a completely equivalent formalism is obtained if the amplitude operators Ē(±) (r, ω, t → −∞) are used instead and the main reason we opted for f ( †) (r , ω) was to have more streamlined final expressions due to their simpler commutation relation in Eq. 10a.With the above discussion in mind, we assume the following ansatz for the IO-relation: In the above expression, A ij (r, ω; ξ, ν; t) and B ij (r, ω; ξ, ν; t) are the input-output (IO)-coefficients.In Eq. 13, ξ and ν 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 normalmode 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 ν.However, it also shows the output to be a linear combination of excitations at different spatial positions, manifested by the spatial integral over ξ.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 ξ and ν into the vector Ξ = (ν, ξ) whenever possible and replace dξ dν with dΞ.
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 A ij (r, ω, Ξ; t) and B ij (r, ω, Ξ; t) have been obtained.Moreover, Eq. 13 can be used to express field-dependent quantities directly in terms of the IO-coefficients.To il-lustrate this, we consider the first-order field correlation function g ) 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 (1) ij (r, r ; t, t ) in terms of the IO-coefficients, we expand the field operators using Eq.7b 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): If the system is initially not in the vacuum state, e.g., in the case of a seeded nonlinear process, g ij (r, r ; t, t ) (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 (1) ij (r, r ; t, t ), one need only find the IO-coefficients at times t and t and insert them into Eq.14.

C. 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: The boundary conditions for the IO-coefficients are derived by equating Eq. 13 with Eq. 9 and setting t → −∞ on both sides: 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 B ij (r, ω, Ξ; t) in the following way: where b (0) ij and B ij satisfy the following conditions: After inserting the decomposition Eq. 17 into Eq.13 and relabeling A ij → A ij for notational convenience, the IOrelation takes the form: where Ē(0,+) The remaining two terms in Eq. 18 quantify the changes that the freefield component of frequency ω undergoes due to the nonlinear interaction through the "new" set of coefficients -A ij (r, ω, Ξ; t) and B ij (r, ω, Ξ; t).

The differential equations coupling the IO-coefficient
A ij (r, ω, Ξ; t) and the newly introduced B ij (r, ω, Ξ; t) are derived by inserting the decomposition Eq. 17 into Eqs.15: The coupled quantities in the above system have the initial values: kj (r , ω , Ξ), which acts as a source term for equation 19a.After taking the frequency integral, it will have the following form: 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 (r, r , ω), linear susceptibility ε(r, ω) and nonlinear susceptibility χ (2) ijk (r).These properties, along with the pump field E (±) P (r, t) are used to define the integration kernel F ik (r, ω; r , ω ; t) and source term S (0) ij (r, ω, Ξ; 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 ampli-tude operators Ē(±) (r, ω, t).These, in turn enable us to reconstruct the full electric field operator as shown in Eq. 7b.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.

D. 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 fielddependent 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 r 0 and summed over all field polarisations is given by [61]: where ω 0 is the frequency at which the spectrum is evaluated.
To obtain an expression for σ(r 0 , ω 0 , t) in terms of the IO-coefficients A ij (r, ω, Ξ; t) and B ij (r, ω, Ξ; t), we use the first-order correlation function Eq. 14 (remembering that A ij ≡ A ij ) and insert it into Eq.21: 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 /t and ω/ω integrals with sufficient numerical accuracy, one requires the values of the coefficient A ij (r, ω, Ξ; t) at many narrowly spaced times t in addition to the other variables ω and Ξ, which dramatically increases the computational resources required for such a calculation.The above limitation is a consequence of the individual frequency amplitudes Ē(±) (r, ω, 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 [61,62].
To overcome these difficulties and further expand the applicability of our formalism, we develop a frequencydomain 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 [62] and, as will be shown in the next sections, evolve in time in a way quite similar (although not identical) to Ē(±) (r, ω, 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: where the operators Ẽ(±) (r, ω, t) have been defined as: 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 Ē(0,±) (r, ω) as t → ∞.For finite t and/or with sources (i.e. an interaction) present, their interpretation is more involved and is discussed in [61,62].
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 IOcoefficients) at all times prior to t is eliminated, provided the filtered fields Ẽ(±) (r, ω, 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" Ê(±) (r, t) and thus of Ē(±) (r, ω, 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.7b and insert the resulting field decomposition into Eq.24.Thus we obtain: where we defined: Here, Ẽ(0,+) i (r, ω) is the filtered analogue of the freefield amplitudes Ē(0,+) i (r, ω) and Ãij (r, ω, Ξ; t) and Bij (r, ω, Ξ; t) are the IO-coefficients of the filtered field.For conciseness, we will refer to them as filtered IOcoefficients.
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: As expected, we no longer require knowledge of the IOcoefficients 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): where S(0) ij (r, ω, Ξ; t) is derived from the source term S (0) ij (r, ω, Ξ; t), present in Eq. 19a and is defined as: × dr χ (2) Here, E P,k (r, ω p ) is the Fourier transform of the pump amplitude, defined by E (±) The integral kernel Fik (r, ω; r , ω ; t) which couples the two coefficients in Eq. 28 is given by: lmk (r ), ( 30) and the initial conditions for the filtered coefficients are easily inferred from the definitions 26: 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 timedomain formulation, the properties of the optical system and the pump field are embedded in the source term S(0) ij (r, ω, Ξ; t), defined in Eq. 29, and the integration kernel Fik (r, ω; r , ω ; 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 Ẽ(−) i (r, ω, t) Ẽ(+) j (r, ω , t) and Ẽ(+) i (r, ω, t) Ẽ(+) j (r, ω , 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.F5, 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 [63][64][65].
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 [66][67][68].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.
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 ε(ω), 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]: 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 [69].The permittivity of the combined waveguide and analyte system has the form: where Ω P l is the plasma frequency of the waveguide dielectric, while Ω 0 and Γ are the frequency and width of the dielectric resonance, respectively.The term ε loss (ω) = , models the spectrally localised loss of the analyte, whose magnitude is modulated through the unitless factor α. Here, ω loss and γ 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.
For the incoming pump pulse, we use the frequencydomain form: e ik(ω)z .Here, ω p0 is the pump central frequency, τ p is the temporal width of the pulse, k(ω) = ω c ε(ω) and E 0 is a normalisation constant determined by the total energy of the pump pulse U 0 , where U 0 ∝ |E 0 | 2 .

A. Simulation parameters
All of the simulation parameters given here were normalised in the following manner: frequency quantities are expressed in units of ω p0 , temporal quantities in units of 1 ωp0 and spatial quantities in units of 2πc ω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 × 10 4 and is periodically poled with a spatial dependence χ (2) m cos Λz, where χ (2) m is the maximum absolute value of the nonlinear permittivity and Λ is the poling period, chosen such that the central phase-matched frequencies of the signal and idler photons are ω s0 = 0.7 and ω i0 = 0.3, respectively.The resonance frequency of the analyte was set to be identical to the idler central phase-matching frequency, ω loss = ω i0 = 0.3, to make the effects of the loss more prominent.The width of the loss spectrum was chosen to be γ = 2.5 × 10 −3 .The waveguide parameters in ε(ω) were set to be Ω 0 = 2.1, Ω P l = 0.25 and Γ = 10 −7 , in order to satisfy the following conditions: 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 ω p0 ; 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; iii) the pump temporal width was chosen to be τ 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 χ (2) √ 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 χ (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.

B. Lossless medium
We begin our investigation by studying the singlephoton 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 ε(ω)).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 I m (at ω s = ω 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: , where a and b are fitting parameters.The parametric gain G for a particular value of pump energy is then defined as G = b √ 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 τ p = 2400 and τ p = 600, and observe that they both indeed follow the sinh 2 law.The wider-bandwidth pulse, corresponding to τ 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 [33,56].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 [33,56].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.
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 [33,56].The broadening of the two-photon spectrum can be observed by studying the second-order moment of the output field N (ω, ω ) = Ẽ(+) (r 0 , ω, t) • Ẽ(+) (r 0 , ω , t) , which is related to the JSA of the output photons [70].We can calculate the moment directly in terms of the IO-coefficients using the expression: which is obtained by expanding the filtered field operators using Eq. 25 and taking the expectation value.N (ω i , ω s ), in the case of wide-band pump (τ 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 highgain lossless SPDC [33,56]. ) 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.

C. Lossy medium
To investigate the effects of increasing gain in QSUP, we introduce the analyte into the material permittivity ε(ω) 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 ω (c) = 0.7, which corresponds to the resonance frequency of the analyte ω loss through the conservation of energy ω p0 = ω (c) + ω loss .Both sets of spectra were obtained for the narrow-band pump with τ p = 2400.
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 ω s = ω (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.Im [ε loss ]) is also shown for comparison.The dependence of the extinction maximum and its fullwidth-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 (a) The normalised second-order field moment N (ωi, ωs) for the wide-band pump (τp = 600) at low parametric gain -G = 0.0915 and high parametric gain -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 τp = 600 at varying gain.FIG. 4. Normalised single-photon spectrum of the signal photons for the narrow-band pump with τp = 2400 at varying amounts of parametric gain, (a) without the analyte, and (b) with the analyte.The analyte loss is centred around ω 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.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 ω 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 ω (c) , resulting in a net narrowing of the extinction spectrum.
As mentioned at the beginning of this section, sensitiv- ity 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. [66] 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 ω seed the generation of further signal (idler) photons, not only at the frequency ω p0 − ω, 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 (G > 2 in our case), but become more prominent as gain is further increased.

D. Spatial evolution of the spectrum
The spectra given and discussed in Sec.III B and III C 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 C.
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 wellestablished quadratic law, where the intensity is ∝ 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 pho-tons.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 [71].
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 highorder 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 [72] or microwave to optical domain [73].

ACKNOWLEDGMENTS
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).We begin by expanding the Ê(±) i,j (r, ω, t) operators in the commutator according to Eq. 2: Now, we use the commutation relation Eq. 4 and the properties of the δ-functions to obtain: The above expression can be simplified by using the following GF identity [42]: 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.
The rotating-frame creation/annihilation operators can be defined in a similar manner to be: where again we have f ( †) i (r, ω) ≡ f ( †) i (r, ω, t → −∞) for compactness.As with the amplitude operators, the rotating-frame annihilation operators are related to their Heisenberg picture counterparts as: fi (r, ω, t) = fj (r, ω, t) e −iωt , (B5) with the relation for the creation operators obtained by taking the adjoint of Eq.B5.
To find the general equation of motion for Ē(+) i (r, ω, t), we begin by finding the time derivative of Eq.B3: We then recall Eq.B2b and obtain, after some straightforward calculation: Due to the properties of the commutator, Û ( †) SPDC (t) can act on the operators inside of it independently and give: where we used Eq.B3.Also recalling Eq.B2c and the decomposition Û (t) = Û0 (t) ÛSPDC (t), we can identify Û † SPDC (t) HSPDC (t) ÛSPDC (t) to be the full, Heisenbergpicture 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.B5 in the commutation relation for the Heisenberg picture operators Eq. 4: Since the δ-function in the above expression results in ω ≡ ω , the exponential factor can be ignored.The same is valid for the commutation relation of the rotatingframe amplitude operators, which is found by replacing Eq. 7a in Eq. 10b: Appendix C: Equations of motion for the amplitude operators To find the commutator Ē(+) i (r, ω, t), ĤSPDC (t) , we first expand the SPDC Hamiltonian from Eq. 1 in terms of the rotating-frame operators Ē(±) (r, ω, t): 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: where we used Eq.10b.To simplify the relations, we assume that the Kleinman symmetry condition is valid (there is negligible dispersion in the nonlinear response), which is a commonly used assumption in many practical scenarios of interest [3].In this case, we have the permutation symmetry of the nonlinear susceptibility χ (2) lkj , hence the two terms on the right-hand-side (r.h.s.) of Eq.C1 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.C1, which adds no complexity to the calculation, and only makes the expressions longer.Thus we can write: After performing the ω -integral and replacing the value of K, we obtain the final form of the equation of motion given in Eq. 11.
Now we examine the time integral t −∞ dt e i(ωp−ω −ν)t and introduce a substitution t = t − τ which allows us to rewrite it as e i(ωp−ω −ν)t ∞ 0 dτ e i(ω +ν−ωp)τ ≡ e i(ωp−ω −ν)t ζ(ω + ν − ω p ), where ζ(ω) 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 [75,76].One notable analytical property of the GF is [46,76]: which we can use to further simplify Eq.E2 by performing the integration over ω .Thus we arrive at the final form for the filtered source term 29.
The second term appearing on the r.h.s. of Eq.E1 is: If we write out the full form of F * ik (r, ω; r , ω ; t ) using Eq. 12, we have: where we immediately expanded the pump field into a Fourier integral.After some reordering of the integrals in the above expression, we can identify: where the r.h.s. was established according to Eq. 26c.Thus, we now write the term in question as: To arrive at the form of the term given by Eqs.28a and 30, we first introduce a variable substitution ω p = ω + ω into the above expression: P,l (r , ω + ω ) Bkj (r , ω, Ξ; t).
The newly introduced ω has the limits [−ω , ∞), 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 ω 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 ω and ω and write the final form of the second term of Eq.E1: lmk (r )Im [G im (r, r , ω)] e −i(ω−ω)t E (−) P,l (r , ω + ω) Bkj (r , ω , Ξ; 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 Bkj (r, ω, Ξ; t), after which we obtain the coupled equations 28.
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 ξ gives a dominant contribution to the output we examine the low-gain analytical expression for the coefficient Ã(z, ω, Ξ; t → ∞) (since it is the only one contributing to the spectrum) given by simplifying Eq.F4 according to the assumptions made about the waveguide under consideration in Section III: × dz χ (2) (z )E .Since the variable z is always confined to the nonlinear region of the waveguide, the effective range of the variable ξ 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 ξ.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 ξ outside of the nonlinear region can be made more manageable via a variable transformation.We used the double-exponential trans-formation of the form [77]: ξ = sinh π 2 sinh θ where θ is the substitution variable whose range is chosen to correspond to values of ξ outside the nonlinear region.With such a transformation, we compress the potentially long range of integration for the variable ξ in the surrounding linear regions into a much shorter range for the variable θ.This transformation resembles the coordinate transformation used in implementing perfectly matched layers (PMLs) in computational photonics [78].
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 ν, 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. F5, we observe the following in the low-gain regime: the range of frequencies ν contributing to the output at frequency ω 0 is determined by the spectrum of the pump, i.e. output photons of frequency ω 0 are influenced by field-matter excitations of frequencies within a pump bandwidth of ν 0 = ω p0 −ω 0 .However, as gain is increasing, the generated field at frequencies neighbouring ω 0 will start contributing to the intensity at ω 0 through higher order seeding effects, where those neighbouring frequencies themselves are affected by their own respective neighbouring frequencies further from ω 0 .Hence, to have an accurate representation of the intensity at ω 0 , one has to include the whole range of frequencies with substantial intensities around ω 0 , which are mainly determined by the phase-matching condition.Hence, in our calculation, for both ω and ν 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.

FIG. 1 .
FIG.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.

FIG. 2 .
FIG. 2. (a) The dependence of the spectral maximum at ω = 0.7 on the total pump energy for the narrow-band pump (red squares) with τp = 2400 and wide-band pump (blue circles) with τp = 600.The solid lines indicate the fitted sinh 2 (b √ U0) 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.

FIG. 3 .
FIG. 3.(a) The normalised second-order field moment N (ωi, ωs) for the wide-band pump (τp = 600) at low parametric gain -G = 0.0915 and high parametric gain -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 τp = 600 at varying gain.

FIG. 5 .
FIG. 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 Im [ε loss ] (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.

FIG. 6 .
FIG. 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 (τp = 2400) and parametric gain of 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 (τp = 2400) in the lossless case, for different values of gain.