Backscattering in Nonlinear Microring Resonators Via A Gaussian Treatment of Coupled Cavity Modes

Systems of coupled cavity modes have the potential to provide bright quantum optical states of light in a highly versatile manner. Microring resonators for instance are highly scalable candidates for photon sources thanks to CMOS fabrication techniques, their small footprint and the relative ease of coupling many such microrings together, however, surface roughness of the wave-guides, and defects in the coupler geometry routinely induce splitting of the cavity modes due to backscattering and backcoupling. The parasitic back-propagating mode in the microring leads to hybridisation of the modes, altering the linear and nonlinear properties of this system of coupled cavity modes, and ultimately constraining the fidelity of quantum light sources that can be produced. In this paper, we derive a comprehensive general model for Gaussian nonlinear processes in systems of coupled cavity modes, based on an effective field Hamiltonian and a dispersive input-output model. The resulting dynamics of the equations of motion are evaluated in a Gaussian process formalism via the symplectic transformations on the optical modes. We then use this framework to numerically model and explore the problem of backscattering in microring resonators in physically relevant parameter regimes, involving the splitting of various resonances, we calculate the consequent impurity and heralding efficiency of various heralded photon schemes, we explore a perturbative explanation of the observations and assess the correspondence between spontaneous and stimulated processes in these systems.

Systems of coupled cavity modes have the potential to provide bright quantum optical states of light in a highly versatile manner. Microring resonators for instance are highly scalable candidates for photon sources thanks to CMOS fabrication techniques, their small footprint and the relative ease of coupling many such microrings together, however, surface roughness of the wave-guides, and defects in the coupler geometry routinely induce splitting of the cavity modes due to backscattering and backcoupling. The parasitic back-propagating mode in the microring leads to hybridisation of the modes, altering the linear and nonlinear properties of this system of coupled cavity modes, and ultimately constraining the fidelity of quantum light sources that can be produced. In this paper, we derive a comprehensive general model for Gaussian nonlinear processes in systems of coupled cavity modes, based on an effective field Hamiltonian and a dispersive input-output model. The resulting dynamics of the equations of motion are evaluated in a Gaussian process formalism via the symplectic transformations on the optical modes. We then use this framework to numerically model and explore the problem of backscattering in microring resonators in physically relevant parameter regimes, involving the splitting of various resonances, we calculate the consequent impurity and heralding efficiency of various heralded photon schemes, we explore a perturbative explanation of the observations and assess the correspondence between spontaneous and stimulated processes in these systems.
Microring resonators are a promising platform for delivering on-chip quantum light sources 1-3 , offering high brightness [4][5][6] , and increasingly high purities for interference 7,8 . However, even slight fabrication defects of waveguides leads to loss and backscattering 9,10 , which in microrings causes coupling of the forward and backward propagating modes. In the linear optical regime these effects have long been well understood 11,12 , are becoming well characterised experimentally 13,14 , including, asymmetric Fano splitting of the resonances 15,16 and methods to overcome these limitations have been demonstrated [17][18][19] . Microring cavities have great potential for uses in sensing applications [20][21][22] , where backscattering often degrades performance, though counterpropagating modes can be exploited for some sensing applications [23][24][25] . Nonlinear optics remains predominantly constrained to single cavities [26][27][28][29] , with multiple degenerate cavity modes just considered in the continuous-wave (CW) pumping regime 6 , alongside preliminary investigations into interrogation of the back-propagating modes 30 . Similarly, backscatter in cavities generating Kerr frequency combs often plays a significant role in degrading the quality of the resonances and subsequent phase stability [31][32][33] , though these classical models tend to focus on soliton formation in the steady states of the strongly nonlinear regime [34][35][36] , with counter-propagating solitons only recently being proposed 37 and demonstrated 38 .
These systems however differ drastically from the Gaussian regime in consideration here, in which we a) Electronic mail: w.mccutcheon@hw.ac.uk pursue properties of the spectral correlations between quantum fields of few cavity resonances in the presence of weak nonlinearities that can be well approximated by Hamiltonians which are just quadratic in the quantum mode operators, ie. Gaussian nonlinearities. This precludes application to Kerr Hamiltonians, which are quartic, were all the fields to be treated quantum mechanically, however in many instances some bright fields can be treated classically, resulting in quantum fields that evolve under Gaussian (quadratic) Hamiltonians. Beyond backscattering in microrings, a general quantum mechanical treatment of the broadband spectral behaviors of a number of coupled cavity modes, coupled to a number of waveguides has wide application to the increasingly complex photonic devices being fabricated, from the linear dynamics of photonic atoms and molecules [39][40][41] through to next generation quantum light sources comprised of multiple cavities 42,43 . Whilst temporal coupled mode theory (TCMT) has been been demonstrated in the linear regime 44,45 , and for limited nonlinear systems 42,46 , a general Hamiltonian based model for evaluating the Gaussian nonlinear quantum optical mode transformations induced by a several cavities, coupled to several waveguides, with arbitrary pumping and input states, remains to be demonstrated. Moreover, the Gaussian treatment here demonstrated exemplifies the nature of the linear symplectic inversion necessary to arrive at the full Gaussian solutions to such systems (applicable for arbitrary pumping and input fields), it invites the demonstration of a novel perturbative solution, as well as its polar decomposition, which provides insights into spontaneous and stimulated behaviors of general sys-tems.
The result presented herein provides a system of equations that can be directly applied to a wide host of systems, mitigating the need to derive such equations on an ad hock basis. Furthermore, in the pursued treatment different waveguide dispersion can be captured, and the solution to the system is the full Gaussian mode transformation 47,48 easily facilitating the investigation of features beyond the often used first order approximations, such as applications in nonlinear interferometry 49 . In addition, a first order perturbative solution is presented, in the nonlinear interaction, that allows the linear dynamics of the system to be maintained to arbitrary order thereby helping to explain phenomena relating to the stimulated and spontaneous processes. The subsequent application to backscattering in photon pair generation, provides clear understanding of how the various effects manifest and how to accommodate them, which will directly aid in the design of a range of microring based devices.
To pursue a complete quantum optical model of such systems of coupled cavity modes, we draw closely upon the effective field methods developed to treat nonlinear optical processes in dispersive media in a canonical formalism [50][51][52][53] . In particular, its applications to microring resonators 28,29 which closely resembles input-output formalisms for TCMTs often used for such systems. This framework however is also applicable to modes of differing dispersion providing strict notions of mode orthogonality and normalisation 52,54 .

I. The Model
We consider a system of J wave-guide 'channels' with effective field annihilation operatorsŝ j (z, t), with j = 1, ..., J 50,53 . The displacement field in these wave-guides is,D with transverse coordinates x and y and longitudinal coordinate z being dependent on the field j. The channel field j is associated to the central frequency Ω j and the dispersion relation k j (Ω) for which we expand just to first order having group velocity v j . The normalisation prefactors ensure that the effective field operators obey the standard equal-time commutation relations ([ŝ j (z, t),ŝ † j (z , t)] = δ jj δ(z − z ) ) providing the transverse field profiles are normalised accordingly (see Appendix F).
There will exist multiple effective fields indexed by different j, occupying the same spatial mode over disjoint frequency intervals, resulting in effective fields which are narrowband and thus slowly varying in space owing to the fast varying factor e ik(Ωj )z , which leads to additional terms arising in the EM Hamiltonian 53 . The different effective fields are associated to independent first order dispersion relations, meaning that higher-order dispersive properties of the waveguides may be present in the model, however the individual narrowband effective fields are subject to only linear dispersion. As such, we require that each effective field is sufficiently narrowband to neglect higher order dispersion.
These are coupled at z = 0 to N cavities modes with annihilation operatorsâ n (t) with n = 1, ..., N , and energy ω n . The coupling between channel mode j and cavity mode n is given by the channel-cavity coupling elements γ nj . The cavity mode n is coupled to cavity mode m with cavity-cavity coupling g nm . Finally channel mode j is coupled to channel k (at z = 0) by channel-channel couplings C jk . We note that such a point coupling model is only possible when group velocities are matched so insist that C jk = 0 ⇒ v j = v k . The complete linear Hamiltonian is The first term is the free Hamiltonian for the cavity modes and the third, four and fifth terms describe channel-cavity, cavity-cavity and channel-channel coupling respectively. The second term describes the free Hamiltonian of the channels, exhibiting terms arising from the series (linear) expansion of the dispersion relation and a subsequent partial integration 50,55 . Let V , Ω and ω be diagonal matrices of group velocities (v i ), car-rier frequencies of effective fields (Ω i ) and cavity modes (ω i ) respectively. And we denote the matrices C, g and γ having elements C jl , g nm and γ nj respectively so that g = g † and C = C † .
We focus on the nonlinear Hamiltonian resulting from χ (3) response of the material which will be assumed to only occur in the cavities, in which the field amplitudes are large and may be described by, with λ to be determined by the underlying cavity geometry and material properties (see Appendix E). Similarly, χ (2) processes could be tackled in a more straightforward manner. We will neglect nonlinear losses, such as twophoton-and free-carrier-absorption 56 .

Linear Equations of Motion
We construct vectors containing effective field operators in the channel and cavity and denote them by bold operators, s(z, t) = ŝ 1 (z, t),ŝ 2 (z, t), ...,ŝ J (z, t) and The channel mode operators obey the Heisenberg equations of motion (EOMs), (3) The fields entering and exiting the coupling region are given by the one-sided limits s ± (0, t) := lim z→0 ± s(z, t) and at the discontinuity we impose s(0, t) = 1 2 (s + (0, t) + s − (0, t)). We may then compactly express the inputoutput relations (see Appendix A), Note that T :=C −1C † = exp{−2i tan −1 (V −1 C/2)} is explicitly unitary and can asymptotically approach arbitrary passive transformations. We first consider the linear equations of motion (H (N L) = 0) for the cavity modes, which, in terms of the slowly varying envelope operators withΓ Linear transformations admit frequency domain solutions in terms of the Fourier transformed operators, To facilitate moving to a frequency domain model we introduce the phenomenological group velocities, V, for the cavity modes. Then the Fourier amplitudes for the cavity fields are defined Then by induction we find (see Appendix H), where we have defined, This linear system can be solved exactly in terms of the absolute frequency as shown in Appendix D. We often consider however, the case of equal tunings and group velocities everywhereΓ andγ have support, ie. γ nj = 0 ⇒ V jj = V nn and g nm = 0 ⇒ ω n = ω m . In this case both integrands become proportional to δ(k − k ) and we have,ã For which the transmitted fields become, These solutions to the linear equations of motion constitute a system of equations equivalent to those encountered in existing literature on the temporal coupled mode theory of multiple optical cavities. This approach however, derived from an underlying Hamiltonian which is explicitly Hermitian, avoids redundancies in parameterisation 44,58 . The precise relationship between TCMT models and the more general quasi-normal modes of such systems we leave for future work 45,59 .

Nonlinear Equations of Motion
With the inclusion of the nonlinear Hamiltonian (Eqs. 2), we obtain additional terms in the EOMs for the cavity modes, Eqs. 5. In particular, to theâ n (t) element on the l.h.s. of Eqs. 5, we must introduce the terms i n mm λ nn mm â † n â mâm . When all fields are treated quantum mechanically closed form solutions to the EOMs can rarely be found, and for bright quantum fields numeric approaches do not yield well to perturbative methods. At this stage, we move to the Gaussian regime and declare which fields are signal fields which will be assumed relatively weak and treated quantum mechanically. These include both signal and idler fields from conventional treatments of four-wave mixing (FWM), though owing to the multiplicity or degeneracy of these fields they are all considered the signal fields. In contrast the pump fields are those bright fields which may be approximated by their classical mean fields which will obey the same equations of motion, with the operators replaced with the mean fields â pn (t) . These bright pump fields obey cubic equations of motion when including all nonlinear effects, though they can be tractably numerically solved independently providing we neglect the negligible effects arising from back-action from the weak signal fields, ie. we exclude nonlinear coupling to the quantum signal modes in the pump modes' EOMs. This approximation effectively neglects any annihilation of pump photons through coupling to signal fields and is called the undepleted pump approximation. In addition, nonlinear terms involving only signal fields can be safely neglected from their EOMs. This regime prevents the (quartic) Kerr Hamiltonian from acting on the quantum fields directly, and instead is capable only of dealing with instances where two (of the four) fields it acts on need to be treated quantum mechanically, resulting in, at most, quadratic terms in the quantum EOMs.
We have then two systems of equations: 1) the classical EOMs for the pump fields, 2) the quantum mechanical quadratic EOMs for the signal fields driven by timedependent contributions from the pump fields. These signal fields then exhibit time-dependent quadratic nonlinear interactions, hence their Gaussian nature, and these come in two forms. Phase modulation -which contributes terms of the form i pn p λ npn p â † p â p â n which resembles a time varying linear process and includes frequency conversion when n = n. And Squeezing -in which the evolution of the creation and annihilation operators becomes coupled together via terms of the form i n pp λ nn pp â p â p â † n . Note that even in the case of perturbative methods, such as photon pair generation considered below, it is implicit that the system is in the Gaussian regime in order to obtain the quadratic Hamiltonian for declared quantum fields.
To accommodate these quadratic terms in the EOMs, it is convenient to adopt the formalisms of Gaussian quantum optics 47,48 and introduce the capitalised vectors of operators containing both creation and annihilation operators,Ã(t) = (ã (t),ã † (t)) andS ± (0, t) = (s ±T (0, t),s ± † (0, t)) . By doing so we may express the full nonlinear equations of motion for the signal fields in the cavity, where the elements of the nonlinear parts are Solutions to these EOMs are in fact linear symplectic transformations -the central feature of Gaussian quantum optics. We will from hereon use bold uppercase symbols to denote these kinds of vectors of operators and the symplectic operations which act on them. Whilst one could proceed by propagating the system under an appropriate Greens' function, in the absence of phase modulation, it is more straightforward to work again in the frequency domain.
In particular, we pursue a frequency domain Gaussian input-output relation of the form, with the linear symplectic operator M(k, k ) to be defined below. It is convenient to define, so that the frequency domain Fourier transforms arẽ Inline with Eqs.7, by induction we derive (see Appendix H), whereΓ For some cavity modes indexed p r and p r , containing classical pump fields ã pr (t) = dk/ √ 2π ã pr (k) e −ivp r kt (similarly for p r ) and approximately energy matched signal and idler mode s r and i r so that λ prp r srir is non-negligible, the corresponding component can be evaluated where ã pr * ã p r (·) is the convolution. All that remains to be done, is inversion of the scalar valued linear transformation on the l.h.s. of Eqs. 12, resulting in an explicit solution forÃ(k), followed by a substitution into the input-output relatioñ This function, M(k, k ) entirely characterises the linear and nonlinear response of system. It is, by construction, explicitly symplectic and serves as the focus of the forthcoming discussion. On the existence of closed form analytic solutions to this system, we observe that closed form inverses to the function valued matrix in Eqs. 12 limit the practicality of such expressions in all but some particularly convenient examples. Insight can be gleaned by considering the Neumann series expansion of this inverse about Γ (Sq) , which displays the characteristic relationship between Gaussian transformations and their underlying Hamiltonian, in particular, elements of M(k, k ) contain hyperbolic trigonometric functions of Γ (Sq) (along with a linear transformation). For small nonlinearities, perturbative methods can be obtained by truncating this series to first order (see appendix B), which amounts to the approximation sinh(Γ (Sq) ) ≈ Γ (Sq) . We add, that in contrast to many perturbative methods, this approximation can be applied whilst maintaining the full linear response of the system, which is integral for applications involving non-trivial input states to the transformation(to be discussed in Section IV).
We provide a summary of how this framework is applied.
Application of the Framework 1. Identify the pump fields' central frequencies, group velocities and couplings Ω j , V j , ω j , V j , γ ij , g ij , C ij and the input pump fieldss − j (k) for i, j indexing pump fields.
2. Construct Eqs. 7 and solve for the pump cavity fieldsã j (k) with j indexing pump fields, (possibly via Eqs. 8 or Eqs. D.6) 3. Identify the signal fields' parameters: We see that parts 1) and 3) trivially involve gathering the relevant information describing the system being modelled. Part 2) involves solving the dynamics of the pump fields in the cavity modes. Since the pump fields are scalar, this could be achieved by numerically solving Eqs. 10 with the inclusion of phase modulation. Part 4) consists of identifying the strength of the nonlinear interactions between the various pump and signal fields. Parts 5) and 6) then result in the signal fields' Gaussian transformation due to the cavity modes and the nonlinear interaction with the pump fields.
A wide variety of systems, defined by the parameters outlined in parts 1),3) and 4), can be readily solved in this framework. One such example from Chuprina et al 42,43 is given in Appendix G. We turn for now to another application, backscattering in photon pair sources, which despite having widespread ramifications on real world devices, has yet to be presented.

II. Backscattering in microring Resonators
Surface roughness of wave-guides in integrated platforms is routinely a cause of propagation losses. Whilst some of this scattering is lost to bulk/environment modes, some of the field is scattered into backwardpropagating modes in the wave-guides. In microring resonators, this counter-propagating mode constitutes an additional cavity mode, and the scattering causes cou-pling between them. Coupling of cavities induces hybridisation and splitting of their resonance, and this can cause degradation of the quantum light sources they aim to generate. We consider the effects of backscattering in a microring resonator designed to achieve degenerately pumped four-wave mixing. We first solve the linear equations of motion for the pump field in the frequency domain, and use this to drive the nonlinear effects on the signal and idler fields.
The bus wave-guide supports forward and backward propagating effective fields at the pump frequency Ω p with group velocities v p , denotedŝ p1f (z, t) andŝ p1b (z, t). An additional 'phantom' channel is introduced to induce loss in the ring and has effective pump fieldsŝ p2f (z, t) andŝ p2b (z, t), chosen (for convenience) to have equal carrier frequency and group velocity. The microring supports an anti-clockwise (forwards) and clockwise (backwards) mode at pump frequency (ω p = Ω p ) denoted bŷ a pf (t) andâ pb (t) respectively. In total the mode operators we need to consider are a p (t) = â pf (t),â pb (t) and We consider a ring designed to be critically coupled with coupling parameter γ p , perturbed by backscattering in the ring characterised by g p , backcouplings (channel-cavity coupling) δ f b and δ bf (taken proportionally to γ p ), and channel backreflection c p , so that our Hamiltonian is determined by the matrices, From Eqs. 8 and Eqs. 9 we can promptly obtain the pump field in the ring and the resulting input-output relations. We consider Silicon-on-insulator wave-guides with cross section 220nm-by-500nm pumped close to 1.550µm for which the effective index can be found by simulation to be n ef f (λ) = 2.44 − 1.13(λ − 1.55) − 0.04(λ − 1.55) 2 . The pump will have central wavelength 1.546µm (channel 39 of the C-Band ITU grid) so that Ω p = 1.2183 × 10 15 rad·s −1 and v p = 7.1532 × 10 7 ms −1 . The microring will have circumference 202.56µm, and an amplitude round trip transmission (self-coupling) σ r = 0.985, which corresponds to γ p = 8.7277 × 10 8 rad·s −1 Lr . We will take as input pump field a tophat in frequency across the simulation range, which would in practice be well approximated by a sufficiently broadband pulse filtered through a dense wavelength division multiplexer, and permits us to focus on the properties of the system, though we note that more elaborate pumping schemes could be used to tailor the generated states further 60 . Since the power coupled to the ring is dominated only by the overlap with the ring resonance, whereas the extent of the tophat remains a simulation feature, the pump power is not a meaningful quantity in this case. In practice some finite bandwidth pump pulse would be chosen, so, for reference, an example with narrowband pumping is included in Appendix F.
In Fig. 1, we plot the transmitted fields, intra-cavity fields, and their convolutions. The transmitted fields to the bus modes indicate features one could access by interrogation of the ring under linear excitation, whereas the intra-cavity powers, which would in general require fitting of the under-laying linear parameters, indicate toward the resonant enhancement and escape probability of generated photons. The convolution of the intracavity power of the pump field informs us of the nonlinear contribution to Eqs. 14, and together with the forthcoming perturbative treatment, these features dictate the generated states. In part a) we choose a set of parameters which reproduce the routinely observed asymmetric splitting of the transmitted field resonance, as well as slight coupler backreflection, g p = 10 10 rad·s −1 , δ f b = 0.2, δ bf = 0 and c = 0.2 rad·m −1 . We note that these fields are dependent on the phase of the coupling parameters, for instance, when arg(δ f b ) = π/2, the field becomes symmetric about k = 0. With a solution for the pump dynamics, we are free to explore how these generate nonlinear effects on the signal and idler fields. Whilst full knowledge of a systems linear response would allow the nonlinear properties to be evaluated, we will focus on few parameter models, for which parts i) through iv) indicate the various fields for g p = {1/5, 8/5, 22/5, 49/5} * 10 10 rad·s −1 (with other parameters 0).
For the purpose of the present manuscript, exploring photon pair generation, where low pump powers can be used at the expense of low generation rates, we will neglect self-phase modulation of the pump, as well as crossphase modulation, and focus solely on how backscattering affects the generated photons through four-wave mixing. This nonlinear effect only become appreciable close to energy-and phase-matched regions, so we will consider just instances in which some signal and idler cavity fields obey 2Ω p ≈ Ω s + Ω i , and we will include both clockwise and counter-clockwise fields as before. For each of the disjoint spectral intervals, effective fields in the bus and 'phantom' channels are introduced, both forward and backward propagating so our signal fields consist of , so the indexes s (i) relate to signal (idler) spectral intervals, 1 and 2 to the bus channel and phantom channel, and f and b to the forward and backward propagating fields respectively. We consider signal and idler resonance to be next-but-one neighbours of the pump across the frequency comb of the microring, so that Ω s = 1.  troduce a complete set of parameters of the form Eqs. 16 for both signal and idler. We will take self-coupling of the ring to be constant so that γ s = γ i = γ p , ie. the length of the coupling region is sufficiently long that the coupling is achromatic across the region of interest. This achromatic coupling is know to cause impurity in the generated photons even in the absence of backscattering, and whilst methods to overcome these limitations have been proposed 8 , this coupling well approximates the most simple design; a microring evanescently coupled to a waveguide. The complete nonlinear contributions are, We have then a complete model to explore the consequences of backscatter and backcoupling in a physically relevant setting. We will numerically solve Eqs. 15, and discuss the features in the symplectic tensor M(k, k ).
To understand what these transformations imply for photon pair generation in the spontaneous regime, we first perform a polar decomposition, so that M(k, k ) = M h (k, k )M u (k, k ), where M h (k, k ) is Hermitian and M u (k, k ) is unitary. Then in the low squeezing regime (low photon pair probability), in the instance a single pair of photons is generated localised in specific modes, the joint-spectral-amplitude (JSA) describing spectral correlations of this bi-photon state is well approximated by the corresponding component of M h (k, k ). In particular, these complex sympletic transformations have the structure, where the nonlinear Hamiltonian-like operator H N L , is block off-diagonal since M h (k, k ) is Hermitian and the final approximation holds for low squeezing power. In this case, the elements of β h (k, k ) approximate H N L which describes the generation of the bi-photon states, and we may select terms β h ij (k, k ) to associate to the JSA of bi-photons generated across modes i and j. Note that having performed the polar decomposition, by con- We will make extensive use of the Schmidt decomposition which can be achieved by integro-singular value decomposition of the relevant functions β h ij (k, k ) 61 . The singular values, λ i (equivalently Schmidt coefficients) allow the heralded photon purity to be calculated by To numerically solve Eqs. 15 we discretise the wave-number k, sample 201 points on the interval k ∈ (−2515.01, 2515.01)m −1 (0.6 times the free-spectral range) which well resolves the features over the parameter regimes considered. Finally, we can also interrogate the temporal correlations of a bi-photon by performing the appropriate Fourier transform (FT) of the JSA to obtain the Joint-temporal amplitude JTA. We pad this Fourier transform to obtain a temporal resolution 3 times greater than the raw FT.
We focus on how the purity and the heralding efficiency of a bi-photon postselected in various output bus channels is affected by backscattering. We focus on three cases of interest. For each case, the A) Pump Splitting -When only the pump field exhibits splitting, the forward and backward propagating modes of the signal and idler are independent and depend only on the corresponding (forward and back) propagating pump fields. Furthermore, the bus and loss modes are equivalent, in that, all squeezing between loss and bus modes exhibits the same spectral properties, a feature that survives even outside of the critically coupling regime up to a constant factor. There are then, just two relevant two-mode squeezing terms to consider in M h (k, k ), those corresponding to photon pairs that are produced either both propagating forward (f-f ) or backward (b-b).
In Fig.2, a) and b), we observe an increase in purity for the f-f pairs as the backscattering parameter for the pump, g p , is increased, attaining a maximum of ≈ 0.97, followed by a decrease to ≈ 0.89 and a gradual return to the raw purity. The increase at small g p can be attributed to an effective broadening of the pump spectrum, which loosens the energy matching constraint which in standard microrings leads to correlations in the signal-idler JSA, and thus heralded photon impurity. A gentle decrease in efficiency can be seen through the drop in pair rate as the peak pump power in the ring decreases. As splitting gets very large the convolution of the pump's spectrum becomes dominated by the central of the three peaks (see Fig. 1) and the system approaches the non-degenerately pumped case, in which two independent Lorentzian cavity modes drive the process, which is analogous to the non-split case (providing a broad enough injected pump) all but for a factor of 1/4 in the generation rate arising from a combinatorial factor of the process. In the intermediate region, iii), whilst the absolute value of the JSA closely resembles the non-split case (corroborated by the purity of the absolute value of the JSA), the correlations in the phase reduce the purity below the non-split case, which can be readily observed in the JTAs, whereby the secondary peak along t s = t i indicates a suppression and subsequent revival of the forward propagating pump power as the forward and backward cavity modes beat together. This observation helps clarify additional phase-dependency of the system, and demonstrates that increased pump broadening by splitting is not sufficient alone to obtain photon purities beyond 0.97.
Despite the apparent symmetry in the parameters of the set-up, the b-b propagating pairs, depicted in part c), behave very differently due to pump being injected solely in the forward propagating bus mode, and they are consequently nonexistent at g p = 0 (due to the absence of the back propagating field). For small g p , column i), the b-b pairs exhibit strong correlations with minimum purity 0.77 owing to the back-propagating pump having a narrower spectrum, in part due to the filtering it must undergo to couple through the forward propagating cavity mode. Nonetheless, broadening of the pump as g p increases still leads to purities in excess of 0.95, though these reduce to the non-split ring purities with the onset of the non-degenerately pumped regime.  B) Idler Splitting -With only idler mode splitting, all photons are generated by the forward propagating pump field, so the signal field is non-zero only in the forward propagating direction, though the idler scattering allows for pairs to be produced with the idler exiting in the backward mode (f-b), as well as the forward mode (f-f ). In Fig. 3 we see the purity and pair probabilities for the case where only the idler mode is split by parameter g i . The purity of the heralded photon suffers as a result of splitting, and these effects can be readily seen observing the changes to the f-f JSAs. In the case the forward propagating idler mode is used as herald, the heralding efficiency is unchanged by the splitting since every forward propagating herald is accompanied by a heralded signal photon (up to the 0.5 probability caused by the critically coupled phantom/loss channel). If however the signal is used as the herald, some partner photons are lost into the back-propagating modes so the heralding efficiency would be reduced by a factor given by the ratio between the f-f and f-b pair probabilities. These observed single photons (forward propagating signal photons without accompanied idler photons) can be understood as a thermal source when tracing over the back-propagating idler mode, but unlike in the case of loss applied to a two-mode squeezed source, this thermal contribution is characterised by an altogether different spectrum. This spectrum, given by the marginal of the f-b JSA (where the idler marginal is obtained by dk s F (k s , k i )F * (k s , k i )), could in principle, be filtered or mode sorted to some degree, from the f-f JSA of the bi-photon term.
In the case of the f-b bi-photons the pair probability quickly approaches appreciable levels for even relatively slight splitting. The purity of such f-b bi-photons begins at 0.968 at very low rates, owing to the narrower filtered spectrum of the idler photons coupling out backwards, though reduces as the splitting increases and the states become brighter. At large splitting the f-b JSA differs significantly from the f-f case, in that it's non-vanishing at k s = k i = 0, which can be seen to result from the non-vanishing field enhancement of the idler at k i = 0, analogous to the back propagating pump mode in Fig. 1, iii) and iv).
C) All Fields Splitting -When all fields are equally split (g p = g i = g s = g) we can observe photon pairs in all configurations of forward and back propagating and we consider the cases, f-f , f-b, and b-b. The purity of the f-f bi-photons decreases less quickly than in the pure idler splitting case, owing to the pump, and thus energy-matching constraint, broadening. For large g the JSA tends to a four peaked function, whereby the loose energy matching constraints, owing to the broad pump, allow photons to be generated in all four combinations of energies of the bi-peaked signal and idler resonances and the JTA exhibits symmetric beating. In the f-b case, the asymmetry is clear in the JSA and JTA for small g though for larger g the JSA appears not dissimilar to the f-f case, though the JTA still features this asymmetry through a tendency towards larger delays in the idler. The b-b case retains the symmetry of the f-f case, though the JTAs show the tendency towards larger delays due to the slower excitation of the back-propagating pump mode due to coupling via the forward mode.
To better understand the origins of these features we turn to the perturbative solution given in Appendix B, whereby the contributions to the JSA can be understood as arising from the nonlinear part ofΓ Sq (k, k ) being resonantly filtered to the output ports. In this case, the JSA is given by the coherent sum of the contributions from the forward and backward propagating modes in the ring. In Fig. 5 we display the f-f JSA contour plot with theΓ Sq (k, k ) density plot overlaid and the signal and idler resonant filtering functions on the axes, for a (wider) range of all equal splitting parameter, g. The sum  The f-f photons exhibit reduction in both purity and pair probability with increased gi. The double peaked idler spectrum becomes quickly apparent in the JSA, accompanied by beating in the JTA which increases in frequency as the coupling becomes stronger. The f-b pairs however, attain purity of 0.968 for small gi owing to the effective tighter filtering of the idler through the backpropagating mode evidenced in the narrow JSA, c)i). For larger splitting, since the backward out-coupled idler field enhancement doesn't vanish at k = 0, unlike the forward propagating case (see Fig. 1, ii), for an analogous case) the generated idler state exhibits a non-vanishing intensity in the spectrum between the peaks of the split resonance.
of the forward (top) and backward (bottom) terms constitute a perturbative explanation of the observed JSAs. Pair Probability

III. Stochastic Methods
Whilst the trends for particular forms of backscattering are insightful and could inform the design of systems of coupled cavity modes, in practice, fabrication defects are out of the experimenters control and behave randomly from sample to sample, and resonance to resonance 62 . The often observed slight asymmetric splitting example of Fig. 1 a), invites an analysis of the statistics of resonators with qualitatively similar transmission properties. We assess the properties of an ensemble of 1000 devices having each resonance's (pump, signal, and idler) parameters sampled from uniform distributions over the intervals |g| ∈ (0, 10 10 ) rad·s −1 , |δ f | ∈ (0, 0.2), |δ b | ∈ (0, 0.2) and |c| ∈ (0, 0.2) rad·m −1 , with the phase of each random. In Fig. 6 a) a histogram of the observed linewidths is presented and part b) shows the first 80 transmission line shapes (solid) of such cavity modes, along with the ideal case (dashed), to provide a qualitative understanding of the resultant line broadening. Whilst these parameters are not significant enough to often show splitting of the resonance, some line broadening is clearly present in most cases as well as slight wandering of the mean position of the resonances as shown in c). The mean purity of the f-f bi-photons observed from this ensemble was 0.915, only slightly reduced over the ideal (non-split) case at 0.921, as presented in the histogram Fig 6 e). The mean pair probability however drops to 0.00795 from the non-split case 0.01265 (Fig 6 d) ).
Despite this relatively minor reduction in the source purity and probability, the quality of such heralded photons is often only useful when many such sources show high quality interference between one another, which invites consideration of the ensemble purity of these sources, ie. the purity of the mixed state comprised of an equally weighted mixture of the heralded single photons from each. We find ensemble purity of 0.8867 indicating that the visibility of interference from separate sources suffers more significantly.
A discussion of the best and worst cases is in order. The perturbative terms contributing to the best and worst case JSAs are depicted in Figure. 7, a) and b) respectively. The best case, achieves a purity of 0.964, slightly short of the optimum purity, 0.968, achieved by the pure pump splitting case (Section II A)). Figure.  pump resonance lead to higher purities as with the pump split case. Whilst of course this is circumstantial, it provides evidence that purities in excess of this may not be possible for systems of this form. Though it is worth noting that the form of the coupling terms we have imposed to model backscatter, need not hold in more general systems of coupled cavity modes and it would certainly be of interest to consider different coupling models 8 and explore the design of optimal sources under these kinds of defects. The worst case, on the other hand, has a purity of just 0.837, despite the relatively minor splitting parameters. In Figure. 7, b) we see relatively broad, and clearly split, signal and idler resonances along with a narrow pump resonance. We remark that these perturbative methods obtain fidelities 0.999917 and 0.999531 to the full solutions, for the best and worst cases respectively. This is to be expected since the full β h takes the form of a hyperbolic sine in Γ (Sq) , and the truncated series omits terms beyond O(Γ (Sq)3 ), whilst the pair rate satisfies |Γ

IV. Challenges for Stimulated Emission Tomography
Stimulated emission tomography (SET), a common method for performing tomography of such systems, consists of inserting a CW seed laser to drive the signal (idler) mode and recording the spectrum of the resultant conjugate beam in the idler (signal) [63][64][65] . Having pursued a Gaussian treatment, this behaviour is directly described by the transformation Eqs. 15 66 providing the undepleted pump approximation still holds, i.e. the pump fields are not significantly perturbed by displacements of the same order of the conjugate beam displacements.
The standard prescription then consists of the following reasoning: measure elements of the full transformation, Eqs. 15, and use these to infer elements of the Hermitian part of the polar decomposition of the transformation, Eqs. 18, which describes the spontaneous (no seed laser) behaviour of the source. In the instance the entire transformation (Eqs. 15) is measured, the polar decomposition is straightforward. Similarly, if sufficient structure is assumed in Eqs. 15, (such as, for instance, the symmetry between the bus and loss channels in the case of no backscatter) it may be possible to infer properties of the entire transformation, and subsequently its polar decomposition, from incomplete data, ie. in the absence of any characterisation involving seeding or measuring the loss channels. Assumptions of this kind however, are not available in the case of backscattering, ie. SET measurements involving only seeding and measurement of the bus channels, do not uniquely determine the spontaneous behaviour of the source. This becomes apparent when considering that the spontaneous behaviour of the source depends on squeezing of the vacuua that enter from various ports of the device, without sufficient knowledge or assumptions, the squeezing applied to the vacuum entering inaccessible inputs to the device (whilst exiting in the forward bus mode) cannot be determined.
Were we to assume no backscattering was present, and estimate the spontaneous JSA using the conventional method (outlined in Appendix C), our tomography missestimates the outgoing state generated in the spontaneous regime, and a histogram of the fidelity of this estimation procedure (over the ensemble, Section III) is provided in Figure 8 a). We find a mean fidelity of 0.9872 is obtained, which puts strict limits the applicability of such tomography schemes when moving towards source purities near or exceeding this. Finally, we consider whether such an approach to tomography biases our inferred purity towards over-or under-estimating the quality, and present the purity gap (difference in the in SET inferred purity and the actual purity) in Figure. 8d). The mean purity gap is just 0.00049, so whilst we suspect any bias is insubstantial over this particular ensemble, it is by no means apparent how other source defects, or particular distributions over them, may affect tendencies toward infidelity of this estimation process. In light of these challenges, it is clear that to achieve high fidelity assumption free tomography of spontaneous Gaussian processes, including single photon sources, one requires additional data or assumptions, including synthesis of linear properties of the system or seeding and detection of the additional modes at play.

V. Discussion
The resultant transformations demonstrate an important feature of squeezing in the presence of multiple spectral and spatial modes -the states of light exiting the device in some modes of interest (the forward propagating signal and idler modes in the bus for instance) are no longer fully characterised by the single JSA 67 , nor do a single pair of joint-spectra from stimulated emission uniquely define the transformation. Instead, squeezing between different spatial modes has consequences for the modes of interest. Furthermore, whilst the heralded biphoton picture of events applies readily in the low squeezing limit, or where all modes (including loss modes) could be monitored with perfect detection, in general multi-pair events muddy any attempt to isolate a single bi-photon and require more involved analysis. Finally, an unfortunate disparity has occurred in this process, whereby the Schmidt decomposition of the bi-photons we have considered, no longer corresponds to a Bloch-Messiah decomposition of the underlying Gaussian process 47,68,69 . This lack-of correspondence inhibits to a large extent, how multi-photon behaviour can be efficiently modelled in such systems, and whilst some progress in this direction has been made 70,71 , we leave the complete analysis to future work, where we expect the nature of multiphoton terms takes on a rich structure.

VI. Conclusion
We have derived a comprehensive Hamiltonian based input-output model for Gaussian optical processes in the presence of multiple cavity modes for general coupling. This framework can be readily applied to arbitrary systems of coupled cavity modes, and allows the exploration of a wide variety of future quantum light sources. The application of such a framework involves just the input of the systems parameters and the inversion of linear equations.
This was directly applied to the case of a single microring exhibiting backscatter, backcoupling, and backreflection, with relevant Silicon-on-insulator waveguide properties. This model, more entire than previous treatments, reproduces observed asymmetric linear phenomena, and the joint-spectra of the photon pair generated by four-wave mixing in such rings was evaluated for various simple backscattering regimes. The limitations of closed form solutions to such systems have been expressed, and a perturbative solution has been presented shedding light on the behavior of such systems for small nonlinearities.
Additionally, the properties of sources sampled from a representative distribution of parameters was considered, an observed slight decrease in single source purity, whilst a significant drop in ensemble purity was observed. The consequences for stimulated emission tomographic methods were highlighted. This all aids in our understanding of the consequences of backscattering in photon sources, helps inform the development of applicable methods to tomograph these systems, and ultimately how to mitigate the problem of fabrication defects in cavity based quantum light sources.
One further example of the application of this frameworks was considered involving a photonic molecule, demonstrating the straightforward nature of its generalization. Finally, the structure of the obtained transformations invites the consideration of more in-depth analysis of the effects of multi-photon noise in quantum photonic systems via a Gaussian process framework. M.
Borghi, D.P.S.McCutcheon, O.Thomas for fruitful discussion. This work was supported by the UK's Engineering and Physical Sciences Research Council (EP/P510269/1 and EP/L024020/1). The author has no competing interests.

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

A. Channel Equations of Motion
The Heisenberg equations of motion for the channel fields are, or in vectorised notation, which can be solved by a coordinate transform ξ j = x − v j t admitting the general solution 29 , with S j (z, t) any solution to the homogeneous problem. The one-sided limits read, At the coupling region we will impose that the discontinuity obeys, and derive the input-output relations,

B. Perturbative Solutions
SplittingΓ(k, k ) into linear and nonlinear contributions, whereΓ (N L) (k, k ) is small and can be treated as a perturbation, we can expand the inverse in the solution to Eqs. 12 using a Neumann series to obtain the perturbative solution, where we drop the explicit k dependencies and integration for brevity. In the most trivial of circumstances, this first order solution can be seen to be equivalent to solving the linear equations of motion for one of the fields (signal or idler) and using this as a source term for the nonlinear contribution in the conjugate field (idler or signal respectively). The resultant input-output transformation becomes just, This solution maintains the full linear response of a system, whilst including all nonlinear contributions to first order. After imposing the polar decomposition the nonlinear part of this transformation is, In this case, the nonlinear contribution can be understood as the pump field's nonlinear contribution in the cavity modes,Γ (Sq) , resonantly filtered through the linear response of the system coupling the cavities to the output modes, L a→s + .
In the case of the backscattering setting considered here, we find only two relevant nonzero terms contributing to any β ij , those corresponding to pair generation in the forward and backward modes of the ring. If for instance, we take the indexes s1f and i1f , (forward propagating, bus channel, signal and idler modes) The functions L a→s + ij (k, k ) are bijective since the linear system is separable in absolute frequency. We can understand them as a resonant filtering over k, and a transformation of axes. The functions | dk L a→s + (s1n)(in) (k, k )| and | dk L a→s + (i1n)(sn) (k, k )|, with n = f (b) for the forward, f, (backward, b) directions, are plotted on the x and y axes of Figures 5.
A straightforward method for obtaining the spontaneous JSA of a photon source in the absence of additional spatial modes, proceeds by noting that Gaussian mode transformations act on displaced states with the very same transformation as the mode operators themselves. So for some transformation, and a seed coherent state (ie. displacement) is injected into mode j with spectrum d − j (k), the output state contains displacement in mode i, For two-mode squeezing β has the structure, which can be obtained by stimulated emission tomography. The polar decomposition (on the support of β, which is sufficient for these purposes) can be achieved by singular value decomposition of these terms by, Essentially, the left singular vectors of β 12 and β 21 comprise the basis of the symmetric part, β h , of β = β h V (with V some unitary), and thus the Hermitian part, M h of M = M h M u in the polar decomposition. The 'standard' prescription for stimulated emission tomography, then consists of following this method even though β comprises of additional terms beyond those in Eqs. C.4 which are not measured.

D. Linear Solution in Absolute Frequency
The linear dynamics of systems are separable in absolute frequency, which can be readily seen by transforming to, with ξ the absolute frequency. The full expression, becomes separable for each frequency ξ as, The input-output expression can be cast in a similar light, so that the linear system at frequency ξ is described by, (D.6)

E. Nonlinear Coupling Coefficients
The displacement field in a microring cavity supporting modes {â n } n is, where d ⊥ n (x, y) is the transverse field profile, the circumferential coordinate is c ∈ (0, L], and the wave-vector must obey k(ω n ) = 2πL/n. Nonlinear Coupling is then ijkl defined as per reference 50 . For rings we consider, where the signal fields are closely separated from the pump, and not too narrowband we may take ∆k = 2πL( 1 ni + 1 nj − 1 n k − 1 n l ) ≈ 0. In the instance of very high quality factors, or significant group-velocity dispersion over the signal fields spacing, this term can become significant.
Whilst this could couple all cavity fields, the non-negligible terms we consider are the Four-wave mixing (degenerately pumped), where p r , p r , s r , and i r are pump, signal and idler resonances of cavity r, which are close to obeying the energymatching constraint, ω pr + ω p r = ω sr + ω ir . We also consider the Self-Phase Modulation of the pump modes, and the cross-phase modulation terms,

F. Narrowband Pumping: Powers and Rates
The displacement fields' transverse modes must obey the normalisation, with v p (x, y) the local phase velocity, v G (x, y) the local group velocity, and the respective dialectric parameters are to be evaluated at the central frequency for the given effective field. We will take n 0 the linear refractive index in the Silicon waveguide n 0 = 3.48, and in the Silica cladding n 0 = 1.46.
Considering an incident coherent pulse of duration σ p = 20ps (for which a significant fraction of the input field couples to the microring, though slight additional spectral correlations result) of the form, having norm α p , so that it contains mean |α p | 2 photons. One finds that pulses of α p = 0.4 × 10 4 , or 2.1pJ (equivalent to 103µW average power at 50MHz repetition rate) lead to f-f pair probability of 0.0085 per pulse, and the generation rate at low powers is proportional to |α p | 4 .

G. Application: A Photonic Molecule
We present one example of where the general framework could be applied, originating from the work of Chuprina et al 42,43 .
In this example, the pump channel fields s = {a p , f xp , f yp } and cavity modes a = {x p , y p } having central resonance frequencies ω = {ω 0x,p , ω 0y,p } = {ω p , ω p } (where we will take ω p from the main text) are coupled via, γ = κ p γ xp 0 0 0 γ yp , g = 0 g p g * p 0 , C = Nonlinearities occur at λ 22,24 ,λ 22,42 = λ (and we will take λ as in the main text), ie. between the pump field y p (degenerately) and y s and y i . By appropriate choice of Ω and V (taking again from the main text), and bichromatic input pump field (to be determined below), one arrives directly at the dynamics of the given system via the method here presented.
Going a little further one could consider what happens if a backward propagating mode y bi is introduced to generation ring (ring y) at idler frequency along with a corresponding loss channel, f ybi , by appending it to the signal fields, s = {a s , f ys , f zs , a i , f yi , f zi , f ybi }, cavity modes a = {z s , y s , z i , y i , y bi } and coupling, where we introduced the backscattering parameter g ib . Define the pump field, setting the losses in the rings equal γ zs = γ ys = γ zi = γ yi = γ ybi = γ xp = γ yp = γ = 0.5γ (ie. half the loss, γ, from the main text). We set the coupling parameters to κ s = κ i = γ , and κ p = 30γ , v s g s = v i g i = 50γ 2 , and g p = 250γ 2 /v p and the pump bandwidth σ p = 8 ln(2)/(0.6g s ). Then we must choose,

H. Derivation of Frequency Domain Expressions
The time domain equations of motion are linear and thus admit linear solutions in the Fourier domain. We move to the frequency domain by applying the Fourier transform, and substituting the operators' for their Fourier components.
For the linear equations of motion (Eqs. 6) we arrive at Eqs. 7 via, For the nonlinear equations of motion, some care must be taken to ensure the Fourier transform is applied consistently. Given Eqs. 10 we arrive at Eqs. 12 via, for which we indeed find the linear solution.