A general framework for multimode Gaussian quantum optics and photo-detection: application to Hong-Ou-Mandel interference with filtered heralded single photon sources

The challenging requirements of large scale quantum information processing using parametric heralded single photon sources involves maximising the interference visibility whilst maintaining an acceptable photon generation rate. By developing a general theoretical framework that allows us to include large numbers of spatial and spectral modes together with linear and non-linear optical elements, we investigate the combined effects of spectral and photon number impurity on the measured Hong--Ou--Mandel interference visibility of parametric photon sources, considering both threshold and number resolving detectors, together with the effects of spectral filtering. We find that for any degree of spectral impurity, increasing the photon generation rate necessarily decreases the interference visibility, even when using number resolving detection. While tight spectral filtering can be used to enforce spectral purity and increased interference visibility at low powers, we find that the induced photon number impurity results in a decreasing interference visibility and heralding efficiency with pump power, while the maximum generation rate is also reduced.


Introduction
Almost all tasks in optical implementations of quantum communication [1,2,3], quantum cryptography [4,5,6], quantum sensing and quantum information processing [7,8,9] require sources of single photons [10,11].Although the requirements of such a source depend on the particular application [12,13], an ideal source would produce single photons with high spectral purity, and either deterministically or with a very high efficiency.To this end, parametric non-linear processes can produce pairs of correlated photons in distinct modes, with the detection of a photon in one mode heralding the presence of a single photon in another.Such sources benefit from potentially very high photon spectral purity [14], and can in many cases be more readily incorporated into integrated devices with a high density and reproducibility [15,16].Such sources, however, are not deterministic, and in order to achieve high efficiencies, it is necessary to multiplex many sources together [17,18].While many of the impressive purity metrics of these sources have been measured in the weak-excitation limit, an analysis of the level of multiplexing necessary to achieve a given efficiency obviously requires operation beyond this limit, and in conjunction with any effects of loss, including filtering [19,20].This problem is further complicated when one also considers the types of detectors employed, as beyond the weak excitation limit photon number resolving detectors project states in a way not possible with threshold (or 'bucket') detectors, which instead only indicate the presence of one or more photons.
Another route towards making a deterministic source is to instead use single quantum emitters such as semiconductor quantum dots [21], defect centres in crystals and 2D materials [12], and single organic molecules [22,23].Using these systems, an excited state can be populated on demand, and radiative relaxation can then occur with very high internal quantum efficiency, with the emitted photon in many cases having a high spectral purity [24,25,26].However, there are often efficiency-purity trade-offs using such sources [27], as it is typically necessary to strongly modify the surrounding photonic environment to boost extraction efficiencies or increase spectral purity (using e.g.cavities), which can simultaneously have detrimental effects on other figures of merit.These considerations, and the challenges involved in incorporating many such sources into a scalable platform [21] leaves significant questions regarding their utility in information processing applications.
For these reasons, the advantages offered by parametric sources motivates a detailed study of their operation beyond the weak excitation limit.Whilst it is possible to write down the quantum state of the field produced by a parametric source in the Fock basis, and from this evaluate the heralding rate and fidelity to a pure single photon state including spectral impurity and multi-photon events [28], it is not a straightforward matter to propagate many such multi-photon multi-modal states through subsequent optical elements and calculate measured detection probabilities.This challenge arises due to the cumbersome nature of dealing directly with large Fock states and points towards the difficulties in the concatenation of multiple systems and tracking of errors.This represents a problem even when seeking to accurately model Hong-Ou-Mandel (HOM) interference of two heralded photons on a subsequent beam-splitter, and particularly so in the presence of loss or if spectral filters are used, since strongly frequency dependent photon number states evolve non-trivially through latter devices in the system.Regarding detection, recent work on Gaussian boson sampling provides expressions for the detection probabilities using threshold detectors [29] and number resolving detectors [30], while multiphoton contributions for threshold detectors have been investigated for single-mode sources [31,32], although in all cases ignoring the photon spectral properties.Spectrally multi-moded sources have been investigated but only for measurements of intensity correlation functions [33,34], which are not in general the appropriate measurements to investigate single photon interference.
In this work we present a general framework of multimode Gaussian optics which includes arbitrary spatial and spectral degrees of freedom, arbitrary photon numbers, and readily accounts for linear optical elements such as beam-splitters and filters, as well as non-linear squeezing operations which represent parametric sources.The Gaussian formalism offers a means to efficiently model an entire system, requiring propagation of a covariance matrix of size only linear in the number of modes, followed by a threshold or photon counting detection model which extends previous results to the multi-spectral moded case.We put this new framework to use by performing a thorough evaluation of single photons heralded from a parametric photon pair source valid for all excitation powers and corresponding photon generation rates, and in which we simultaneously account for spectral impurity and the effects of loss and filtering.We assess the degree of quantum interference from a simulation of the experimentally measured HOM interference statistics between heralded photons from two sources, elucidating the quantitative and qualitative differences found using both detector types.As has been previously established, for threshold detectors we find that the HOM interference visibility decreases with increasing photon generation rate, as at higher powers multiphoton components contaminate the heralded 'single photon' state [28].When spectral impurity is included, this problem is also present when using number resolving detectors and post-selecting only one-photon events; even in the single photon subspace increasing the photon generation rate necessarily and detrimentally affects the distribution of spectral modes in the heralded (truly) single photon state.With the inclusion of spectral filtering, while the interference visibility can be made arbitrarily high in the weak excitation limit, any increases in power to improve heralding rates deteriorate both the interference visibility and heralding efficiency.These results indicate that considerable care should be taken when designing a source with targeted figures or merit, but at the same time provides a general framework to efficiently describe the larger systems for which the sources are intended, allowing for the consequences of multiple and varied source imperfections to be captured.

Gaussian symplectic formalism and photon detection
The general scenario to which our formalism applies is depicted in Fig. 1, and consists of a collection of optical modes, which can be excited, coupled, and the states of which can be detected.As we V n m 7 s r O 7 t 3 9 Q P T x q m S T T j P s s k Y n u h N R w K R T 3 U a D k n V R z G o e S t 8 P x 3 c x v P 3 F t R K I e c Z L y I K Z D J S L B K F r J 7 w 0 S N P 1 q z a 2 7 c 5 B V 4 h W k B g W a / e q X z b E s 5 g q Z p M Z 0 P T f F I K c a B Z N 8 W u l l h q e U j e m Q d y 1 V N O Y m y O f L T s m Z V Q Y k S r R 9 C s l c / Z 3 I a W z M J A 7 t Z E x x Z J a 9 m f i f 1 8 0 w u g l y o d I M u W K L j 6 J M E k z I 7 H I y E J o z l B N L K N P C 7 k r Y i G r K 0 P Z T s S V 4 y y e v k t Z F 3 b u q u w + X t c Z t U U c Z T u A U z s G D a 2 j A P T < l a t e x i t s h a 1 _ b a s e 6 4 = " f G u g O v z m + h q w y p R 8 S t 3 K 6 1 + Z g q s = " > A A A B 7 X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e z 6 Q I 9 B L x 4 j m A c k S 5 i d z C Z j Z m e X m V 4 h L P k H L x 4 U 8 e r / e P N v n C R 7 0 M S C h q K q m + 6 u I J H o e S t 8 P x 3 c x v P 3 F t R K I e c Z L y I K Z D J S L B K F r J 7 w 0 S N P 1 q z a 2 7 c 5 B V 4 h W k B g W a / e q X z b E s 5 g q Z p M Z 0 P T f F I K c a B Z N 8 W u l l h q e U j e m Q d y o e S t 8 P x 3 c x v P 3 F t R K I e c Z L y I K Z D J S L B K F r J 7 w 0 S N P 1 q z a 2 7 c 5 B V 4 h W k B g W a / e q X z b E s 5 g q Z p M Z 0 P T f F I K c a B Z N 8 W u l l h q e U j e m Q d y o e S t 8 P x 3 c x v P 3 F t R K I e c Z L y I K Z D J S L B K F r J 7 w 0 S N P 1 q z a 2 7 c 5 B V 4 h W k B g W a / e q X z b E s 5 g q Z p M Z 0 P T f F I K c a B Z N 8 W u l l h q e U j e m Q d y

E↵ects of photon spectral and number impurity in Hong-Ou-Mandel interference 9
We imagine that Eq. ( 31) describes modes of interest, and ancillary loss modes are initially in the vacuum state.Coupling into the loss modes is a unitary process, and introducing an ancillary mode for each mode in S, the unitary can be written in a general form where unitarity is ensured by , and we are now working in a basis in which ancillary modes are ordered after all modes of interest.The action of this unitary on the total covariance matrix is then where we have omitted modes which we will shortly trace out.The top left block pertains to our modes of interest, and can be written using ˜ Then using the unitary condition we can rearrange to find This map is valid for all transformations that can be constructed as a unitary transformation acting on our system of interest and ancillary modes.In our usual basis the relevant component of the global unitary takes the form U SS = diag(C, C).
For the purposes of modelling frequency-independent loss on all spatial modes we can can take C = p 1 L1 |S| , with 0  L  1 the loss parameter.For loss acting on, for example, spatial mode 1 of 2, we have C = diag( p 1 L, 1) ⌦ 1 N f .Spectral filtering can be included as a frequency-dependent loss with an associated filter function f(!) such that (C) For example a bandpass filter with central frequency ! 0and bandwidth ! is described by This provides us with all the optical components necessary to model a Hong-Ou-Mandel interference experiment.We note that multiple sources can be Schmidt decomposed individually, exponentiated to construct their corresponding symplectic transformations, and then combined straightforwardly as they act on distinct spatial modes.

Hong-Ou-Mandel interference visibilities
modes of i nterest, and can be wri tten usi ng ˜ Then usi ng the uni tary condi ti on we can rearrange to find Thi s map i s val i d f or al l transf ormati ons that can be constructed as a uni tary transf ormati on acti ng onoursystemofi nterestandanci l l arymodes.Inourusualbasi stherel evantcomponentofthegl obal uni tary takes the f orm U SS =di ag(C,C).
For the purposes of model l i ng f requency-i ndependent l oss on al l spati al modes we can can take C = p 1 L1 |S| , wi th0L1 the l oss parameter.For l oss acti ng on, f or exampl e, spati al mode 1 of2,wehaveC =di ag( p 1 L,1)⌦1 N f .Spectralfil teri ngcanbei ncl udedasaf requency-dependent l oss wi th an associ ated fil ter f uncti on f(!) such that (C) For exampl e a bandpass fil ter wi th central f requency ! 0and bandwi dth !i s descri bed by Thi s provi des us wi th al l the opti cal components necessary to model a Hong-Ou-Mandel i nterf erence experi ment.We note that mul ti pl e sources can be Schmi dt decomposed i ndi vi dual l y, exponenti ated to construct thei r correspondi ng sympl ecti c transf ormati ons, and then combi ned strai ghtf orwardl y as Figure 1: The general optical circuit that our presented formalism applies to consists of a collection of modes which each carry one of N s spatial labels i indicating, for example, in which waveguide they are confined, and one of N f spectral labels ω indicating their colour or frequency.A succession of multimode Gaussian transformations labelled M j (which may be linear or non-linear) couple the modes.Threshold or number resolving photon detection then takes place in a subset S (for illustrative purposes here 1 and 3) of the spatial modes.are interested in accurately describing the spectral properties of photons, we take particular care to distinguish between spatial and spectral mode labels.A spatial mode label indicates where (e.g. in which waveguide) a mode is defined, while a spectral label refers to the frequency or wavelength of that mode.To accommodate both of these properties, we label the creation and annihilation operator of a mode with two indices, generically i and ω, referring to spatial and spectral degrees of freedom, respectively.We list all of our mode operators in a vector, which takes the form Â = (â 1 , . . ., âNs , â † 1 , . . ., â † Ns ) , (1) where N s is the number of spatial mode labels, and for each spatial mode label i we have the vector âi = (â iω1 , . . ., âiω N f ) , with âiω being the annihilation operator for a mode with spatial label i and spectral label ω.The number of spectral (equivalently frequency) modes is N f , giving a total number of modes N = N s N f .The mode operators satisfy the usual commutation relations [â iω , â † jω ] = δ ij δ ωω .A variety of candidates exist for parametric photon sources, including free-space systems [35], cavity-based systems [36], in-line or wave-guided systems [37], and inter-modal phase-matched systems [14].Crucially, each shares the common feature that their dynamics are generated by Hamiltonians which are at most quadratic in the quantum field mode operators.While the underlying electromagnetic nonlinearity may be cubic or quartic, in all cases mentioned above only two of the field operators in the Hamiltonian are treated quantum mechanically, and the remaining bright fields contribute as time-dependent scalars.Consequently, once the dynamics are solved, the solutions are Gaussian transformations, which are described by an effective Hamiltonian Ĥ which is quadratic in the field mode operators of interest, and which can in all generality be written where H is a matrix of scalar coefficients which may be a discrete approximation to a continuous function under appropriate regularity conditions [37].The Hamiltonian has corresponding unitary time evolution operator Û = exp[−i Ĥ], and using Ĥ = Ĥ † and the basic commutation relations for the creation and annihilation operators one can show [38] Û and where with 1 N the N -dimensional identity matrix.Since the commutation relation between the mode operators is preserved under a unitary transformation we have or equivalently M † K = KM −1 , and which defines M as a linear symplectic matrix ‡.Eq. ( 3) provides us with a way in which to use the underlying Hamiltonian parameters to propagate the collection of mode operators in the Heisenberg picture.

Gaussian states and transformations
Rather than working directly with the quantum state of the optical modes ρ, we instead consider the corresponding characteristic functions [38].We define the s-ordered characteristic function as where D(Λ) = exp Â † KΛ and Λ is a 2N -dimensional vector, the elements of which are complex variables which describe the quantum state in an associated phase-space.As we detail in the appendix, a Gaussian state is one having a Gaussian characteristic function of the complex variables Λ, and is therefore uniquely determined by a 2N × 2N matrix of scalar coefficients σ which is known as the covariance matrix, and a 2N -dimensional vector d called the displacement vector.In this work we will exclusively consider states for which d = 0.If a Gaussian state undergoes a linear symplectic transformation as in Eq. ( 3), it follows that its covariance matrix transforms as By noting that the vacuum has corresponding covariance matrix 1 2N , we see that Eq. ( 7) allows us to construct a description of a quantum state generated by successive symplectic transformations.
In what follows we show that it is possible to model both threshold ('bucket') and number resolving detectors using only projections of a state onto the vacuum for different subsets of the modes.For this reason, let us therefore consider a subset of the total N modes, which we label S. We are then interested in the probability associated with the projector onto the vacuum for all modes in S, which we label |vac vac| S .As shown in the appendix, the projection of a general Gaussian state described by a covariance matrix σ onto this state can be found to be [31] where |S| is the number of modes in S and σ S is the 2|S| × 2|S| covariance matrix pertaining only to modes in S. Setting σ S = 1 2|S| we find P off (S) = 1 as expected.

Photon detection using threshold detectors
We now consider the probabilities associated with photon detection using threshold detectors, extending the existing results which include only multiple spectral modes [31] or multiple spatial modes [29].Threshold detectors are currently the most widely used experimentally, returning a signal (or 'click') with a probability which depends on the presence of one or more photons in a given spatial mode.These detectors do not typically resolve spectral degrees of freedom, and as such the relevant projection operator corresponding to a detector click in any spectral mode with spatial label i is which is the projection onto all states with spatial mode label i except the vacuum on all N f spectral modes.In the following we use the calligraphic notation S to represent the set of spatial mode labels ‡ The symplectic condition is usually written in which the detection events take place, while Θ = {ω 1 , . . ., ω N f } is the set of all spectral mode labels.As such, a complete set of (spatial and spectral) labels is S = S × Θ.For example, we may be interested in spatial modes S = {1, 3} (see Fig. 1), which including spectral labels gives the set Using this, the probability to detect at least one photon in each of the spatial modes i ∈ S (and with any spectral mode labels) is where 2 S is the power set (the set of all subsets) of S and with B = B × Θ, and for which we can use Eq. ( 8) for Gaussian states.The general form for threshold detectors allows us to also calculate joint probabilities of clicks and vacuum detection events in sets of spatial modes.For detection events in spatial modes X and vacuum in spatial modes Y we have This expression extends the previously known result for threshold detector probabilities involving the Torontonian function in Ref. [29].We note that the number of explicit terms in the sum in Eq. ( 10) depends only on the number of spatial modes in which photon(s) are detected.The number of spectral degrees of freedom increases the size of the reduced covariance matrices in Eq. ( 8).

Photon number resolving detection
We now present expressions for probabilities associated with number resolving detection in spatial modes.Our method takes inspiration from the way in which pesudo-number resolving detectors can be constructed experimentally, namely by 'fanning' out a mode across multiple modes and using threshold detectors.Interestingly, we find that these photon number statistics can be constructed only considering projections onto the vacuum for different combinations of modes.
In our derivation, for each spatial mode we distribute its amplitude across m fictitious ancillary modes using a unitary which has equal amplitudes across its range.Doing so creates an equal superposition of the input state diluted by the vacuum over m modes, and at the end of each we conceptually place a threshold detector.Provided that the number of modes m is much larger than the number of photons present in the initial state described by covariance matrix σ S , we can assume that each of the m diluted modes contains at most 1 photon.In this way the probability for a total number of threshold detector clicks, which we label k, gives an approximation to the probability to detect n photons in the initial state described by σ S .Experimentally this leads to a correspondence between the number of ancillary modes used and the number of photons which can be detected accurately.In our case the ancillary modes are conceptual, allowing us to take the limit m → ∞ and give exact expressions for the probability to detect a fixed number of photons n.
We derive our expressions using the following steps, in part inspired by Ref. [39].
(i) We consider a subsystem of our total Gaussian state, which is also Gaussian and described by a 2|S| × 2|S| covariance matrix σ S , which may contain both spatial and spectral degrees of freedom.
As depicted in Fig. 2, we fan this state out into m supermodes (which may each contain spatial and spectral degrees of freedom), and we label the set of supermodes with boldface calligraphic symbols, i.e. S = {1, . . ., m}.The complete set of mode labels is then given by S ×S ×Θ = S ×S, of which there are in total |S||S| = m|S|.All |S|(m−1) ancillary modes are initially in the vacuum state.The total initial covariance matrix is written σ S×S = σ S ⊕ 1 2|S|(m−1) .The fanning out transformation is achieved using the quantum Fourier transform (QFT) acting on the supermodes, although the specific unitary is unimportant.The QFT evenly distributes initial amplitude in σ S , where we have defined σS = σ S − 1 2|S| and # » σS = (σ S , . . ., σS ) as the block-constant vector with each of its m elements equal to σS , and # » 1 2|S| is similarly a block-constant vector with all entries the identity matrix 1 2|S| .We see that this procedure results in a block rank 1 form for the total covariance matrix σ S×S .(ii) We next determine the probability to detect zero clicks (the vacuum) in some subset of the supermodes B. The quantity of interest is where σ B×S is the covariance matrix pertaining to the subset of modes with supermode labels in B, which takes on precisely the same form as in Eq. ( 13), but with the constant vectors # » σS and # » 1 2|S| instead having length |B| = b.Using the matrix determinant lemma and commutative subring properties we find we can write § which we see depends only on the number of supermodes in subset b.This result also reduces the determinant of the potentially high dimensional covariance matrix in Eq. ( 13) to essentially the determinant of one block σS which is only 2|S|-dimensional, meaning we can evaluate an arbitrary number of fanned out modes without changing the complexity of the determinant.(iii) We now extend Eq. ( 12) to write down the probability to detect at least one photon in each supermode in the subset X and vacuum in all other supermodes Y, namely Eq. ( 15) tells us that the terms in the summation above depend only on the size of the set of supermode labels (B ∪ Y).We then collect terms involving subsets B of the same size to write where k = |X | is the total number of detection events.As we might expect, this probability depends only on the number of detection events, and not the particular pattern.We then sum over all patterns X with fixed length k to give the probability for all detection patterns with k clicks as (iv) In the limit that the number of supermodes m becomes very large, the probability that any supermode contains more than one photon becomes vanishingly small, and a threshold detector click amounts to the detection of exactly one photon.Our expression above for k clicks then becomes equal to the probability to detect exactly n photons.In the limit m → ∞ we have k → n, we can recognise the right-hand side of Eq. ( 18) as a derivative, and we find the number resolving detection probability associated with the detection of n photons in the state described by covariance matrix σ S : (v) Similar steps can be used to find probabilities associated with number resolving detection in distinct spatial modes.The probability for a spatial mode detection pattern described by the vector n = (n 1 , . . ., n |S| ) in spatial modes S is given by where Eq. ( 20) constitutes one of the major results of this work.It allows for the calculation of number resolved detection probabilities across multiple spatial modes, within which multiple spectral degrees of freedom may be present.It should be noted that although our expression in Eq. ( 20) appears relatively compact, the presence of the derivatives means that there is in principle an exponentially large number of terms involved in the limit of large photon numbers.Indeed, this can be seen in Eq. (18).This is to be expected, however, as it is known that the calculation of photon number probabilities from Gaussian states is in the #P computational complexity class [29].
However, the utility of our expression above in fact lies in the way in which the spectral mode degrees of freedom are included.In our expression the size of the matrix entering the determinant scales only with the number of modes, and is fixed with respect to photon number.This should be compared to the other number resolving detection probability involving matrix Hafnians in Ref. [29], for which no distinction between spatial and spectral modes is made.As such, we can here include many spectral modes by simply (linearly) increasing the size of the covariance matrix σ S , without effecting the scaling of the computation with respect to photon number.To do so using the expressions in Ref. [29], one would need to calculate matrix Hafnians for all of the different patterns of spectral modes in which photons could have been detected, thus exponentially increasing the number of terms being calculated.Our expressions permit, for example, one to consider arbitrary photon spectra, purity, and indistinguishability.
Finally we note that making the set of modes explicit in the argument in Eq. ( 20) clarifies detection patterns when other spectator modes are to be traced out.In order that we have a similar notation for threshold detectors, we introduce a list n for threshold detectors which is a click pattern, analogous to the photon number resolving (PNR) click pattern n in Eq. ( 20), but each element can take on values of only 'on' or 'off'.We then define where it is understood that the set X is those modes for which n i = 'on', and Y those modes for which n i = 'off'.This notation allows us to specify, for either detector type, the set of modes we keep from the whole system S and the particular detector pattern n.For example we may wish to consider the first four (of potentially greater than four) spatial modes which we label with italicised numbers, S = {1 , 2 , 3 , 4 }.With number resolving detectors we then specify the photon numbers in each mode, e.g.n = (1, 0, 2, 0), while for threshold detectors we specify the click pattern, e.g.n = (on, off, on, off), and in Eq. ( 21) we should take X = {1, 3} and Y = {2, 4}.Of course, in general P Thres (S; n) = P PNR (S; n), so this common notation only indicates a correspondence and not equality.

Two-mode squeezers
Having introduced Gaussian states and photon detection in general terms, we now explore how to describe specific optical elements within this formalism.We begin with the parametric photon pair sources themselves, which arise from Hamiltonians that take the form of multimode two-mode squeezers.These have the general form where â † i (ν) is the creation operator for a mode with spatial mode i and frequency ν.Note that at this stage the spectral degree of freedom is treated as a continuous variable.We will refer to the function F (ν 1 , ν 2 ) as the joint-spectral-amplitude (JSA), the properties of which have a significant effect on the quality of a single photon source based on a two-mode squeezer.To see this we write the JSA in its Schmidt decomposition [40,41], namely where the λ l are positive coefficients known as the Schmidt coefficients, and the functions {ψ l (ν)} and {φ l (ν)} are each sets of orthonormal functions (though not necessarily equal or mutually orthonormal).
In the low squeezing limit F (ν 1 , ν 2 ) 1 the propagator can be expanded to first order to give the bi-photon state where we have introduced the broadband mode operators, Ĉ † l = ´dν 1 ψ l (ν 1 )â † 1 (ν 1 ) and D † l = ´dν 2 φ * l (ν 2 )â † 2 (ν 2 ), which themselves are mutually orthogonal.We can therefore see that when the Schmidt decomposition has more than one term the biphoton state is entangled, and when detecting one of the photons in an unknown spectral mode, the other will be left in a spectrally mixed state.Conversely, JSA separability guarantees a pure quantum state after detection of one of the photons.For these reasons we will refer to a source described by a separable JSA as a spectrally pure source.
The continuous nature of the JSA function F (ν 1 , ν 2 ) means that the Schmidt basis defined by the functions {ψ ν (ω)} and {φ ν (ω)} must be found by solving integral eigenvalue equations.In practice, it is often easier to instead discretise the integral in Eq. ( 22) which, again, is valid under reasonable regularity conditions [37].In doing so we find the Hamiltonian can be written in precisely the form of Eq. ( 2), where the matrix of coefficients takes the form and the matrix F has elements F ω1ω2 = δνF (δνω 1 , δνω 2 ), with δν being the discretisation step in frequency.The Schmidt decomposition of Eq. ( 23) in the discrete case is the singular value decomposition: where F D is a diagonal matrix of the singular values of F , and U and V are unitary matrices.We note that the factor of δν in the definition of the matrix elements of F means the singular values in F D approximate the true singular values λ l in Eq. ( 23).
Eq. ( 26) allows the Hamiltonian coefficients to be written H = 1 2 UF D U † , where with and we note that U and U are unitary.With H written in this way we can perform the exponentiation in Eq. ( 3) straightforwardly, and find that the symplectic transformation for the multimode two-mode squeezer can be written M = UM D U † where We see that a multimode two-mode squeezer is simply a set of independent two-mode squeezers acting on the appropriate Schmidt modes.

Unitary and passive transformations
In addition to two-mode squeezers which describe parametric sources, we also require unitary mode transformations such as beam-splitters and phase shifters.In terms of symplectic transformations as described in Eq. ( 5), these take the general block-diagonal form M = diag(α, α * ) with α † = α −1 .For a dispersionless (frequency independent) beam-splitter we have while a dispersionless phase-shifter acting on say spatial mode 1 is described by α PS (φ) = diag(exp[iφ], 1) ⊗ 1 N f .A linear frequency dependent phase shift corresponds to a delay in the time domain.Such a transformation (again acting for illustrative purposes here on mode 1 ) is described by α delay = diag(α τ , 1 N f ) with spectral matrix elements (α τ Here ω is a discrete (integer) index, and we note that the delay τ should, in our context, be thought of relative to the bandwidth of photons generated by the two-mode squeezing process, which itself is determined by the typical width of the JSA in Eq. (22).As well as these unitary transformations, we will also be interested in more general non-unitary yet passive transformations, and in particular those which correspond to loss or filtering.To implement such transformations we add ancillary loss mode(s), then use a unitary beam-splitter type transformation as in Eq. ( 22) acting on the mode of interest and ancillary modes, and then trace out the ancillary modes.This can be done analytically and we do not need to explicitly include the ancillary modes in our calculations.From the condition in Eq. ( 5), we know that a general pure covariance matrix will take the form, We imagine that Eq. ( 31) describes modes of interest, and ancillary loss modes are initially in the vacuum state.Coupling into the loss modes is a unitary process, and introducing an ancillary mode for each mode in S, the unitary can be written in a general block form where unitarity is ensured by U SS U † SS + U SL U † SL = 1 2|S| , and we are now working in a basis in which ancillary mode operators are ordered after all mode operators of interest.The action of this unitary on the total covariance matrix is then where U SS = diag(U SS , U * SS ) and we have omitted modes which we will shortly trace out.The top left block pertains to our modes of interest, and can be written using σS : Then using the unitary condition we can rearrange to find This map, described uniquely by the sub-matrix U SS of the unitary, is valid for all transformations that can be constructed as a unitary transformation acting on our system of interest and vacuum ancillary modes, followed by a trace over the ancillary modes.
For the purposes of modelling frequency-independent loss on all spatial modes we can take U SS = √ 1 − 1 |S| , with 0 ≤ ≤ 1 the loss parameter.For loss acting on, for example, spatial mode 1 of 2, we have U SS = diag( √ 1 − , 1) ⊗ 1 N f .Spectral filtering can be included as a frequencydependent loss with an associated filter function f (ν) such that (U SS ) ω1ω2 = δ ω1ω2 f (ω 1 δν).For example a bandpass filter with central frequency ν 0 and bandwidth ∆ν f is described by This provides us with all the optical components necessary to model a HOM interference experiment.

Heralded Hong-Ou-Mandel interference visibilities
To begin our investigation of HOM interference visibilities, let us first discuss how such a measurement may be performed.To measure a HOM visibility in a manner closest to the original experiment two signals are interfered on a balanced beam-splitter, and coincidences at the outputs are compared when the signals are made to be as indistinguishable as possible (by, for example, tuning their arrival times to be equal), and when they are made distinguishable by varying some degree of freedom in one of the signals (typically delaying the arrival time of one of the signals).However, in practice it is not always straightforward to toggle a distinguishability degree of freedom in this way.In particular, typically HOM visibilities measured in integrated silicon platforms are performed by instead scanning through beam-splitter angles using a Mach-Zehnder interferometer [14], as there is no straightforward way to create a temporal delay with the third order non-linearity in silicon.The maximum and minimum coincidences as a function of beam-splitter angle are then used to derive a visibility.Conversely, other platforms such as Ti:LiNbO 3 have a second order nonlinearity and are able to use the polarisation degree of freedom to create a time delay [42].We introduce the term interference parameter to allow us to compare different degrees of freedom in one framework.A measured HOM visibility is then potentially dependent on a) the interference parameter used, b) the detector type used, and c) the visibility function used to combine raw counts into a figure of merit.
We consider interference between heralded photons from two sources, which we term heralded HOM interference.The optical circuits are shown in Figs. 3 a) and d), composed primarily of two sources, a beam-splitter and four detectors.In the ideal case, when the detectors in the two herald equivalently signal modes (1 and 4) register the presence of photons, photons are then necessarily present in the two idler modes (2 and 3), which then interfere and bunch, leading to a detection event in either mode 2 or 3. We are therefore interested in the four-fold coincidence terms P 4 PNR = P PNR (1 , 2 , 3 , 4 ; 1, 1, 1, 1), P 4  Thres = P Thres (1 , 2 , 3 , 4 ; on, on, on, on),  Thres (black), and parts c) and f) using number resolving detectors, P 4 PNR (purple).Parameters correspond to sources with non-separable JSAs and operating beyond the weak-excitation limit.which should vary from large to small as the interference parameter is increased.For later convenience we will also introduce the bunching terms, which indicate successful heralded HOM interference: In Figs. 3 b), c), e), and f) we plot the four-fold coincidences as defined in Eq. ( 36) for identical non-separable sources, as a function of time delay (top row) and beam-splitter angle (bottom row).The black curves (parts b) and e)) correspond to coincidences measured with threshold detectors, while the purple curves (parts c) and f)) correspond the same quantity using number resolving detectors.For the varying beam-splitter case, referred as the Mach-Zehnder HOM interferometer, the four-fold coincidence probability has a maximum for a beam-splitter angle of θ = 0 (acting trivially on modes 2 and 3) which we label (P 4 D ) max .This is compared to the corresponding value for a balanced beam splitter with θ = π/4, which we label (P 4 D ) min .The visibility is typically defined as [43], On the other hand, for a fixed beam-splitter angle of θ = π/4 and variable time delay, the coincidences are typically normalised with respect to a large (effectively infinite) delay, which we label (P 4 D ) τ =∞ .This is then compared to the zero time delay coincidence probability (P 4 D ) τ =0 , which is equal to (P 4 D ) min .We then have the visibility metric [44], In both cases the visibility can depend on the type of detector used, and we see that in general these two visibilities are not equal.Even removing the (P 4 D ) min term in the denominator of Eq. ( 38) results (a)  in different expressions, as the maximum coincidence probability for variable beam-splitter (P 4 D ) max , is not equal to the maximum coincidence probability for variable time delay (P 4 D ) τ =∞ , leading to a different normalisation in each case.While in the low power limit (P 4 D ) max = 2(P 4 D ) τ =∞ , which reflects the relative number of single photon pathways leading to a coincidence, and gives a single way to equate the two visibilities, this does not hold true away from the low power limit.
To explore this further we can introduce the ratio R = (P 4 D ) max /(P 4 D ) τ =∞ , which allows us to write We see that with access only to coincidence probabilities measured with variable beam-splitter, knowledge of R allows a visibility calculated in this way to be equated to that measured using a fixed beam-splitter and a time delay.In Fig. 4 we show how R varies with increasing pump power as captured by the squeezing parameter ξ (defined below), and depending on the type of detector used and level of loss.For threshold detectors we see R decreases with increasing power as coincidence probabilities begin to saturate.However R also decreases for number resolving detectors beyond the low squeezing regime and for non-separable sources.These findings demonstrate that some care must be taken in deducing a HOM interference visibility when using any given experimental setup.Although it would seem from Fig. 4 that either interference parameter can be used provided the measurements are taken in the low power limit, as we will see in the remainder of this paper, interference visibilities are not in general constant with power, even when using number resolving detectors.Moreover, when considering larger scale systems involving more photons, it is unlikely that their successful operation or fidelity with target states will be linear functions of, or even uniquely defined in terms of, these simple HOM interference visibilities.We will see in what follows, however, that in the limiting case of two identical sources and in the absence of loss, it is the visibility V HOM PNR measured with variable time-delay and number resolving detectors which gives a direct measure of the heralded photon purity.As such, for the sake of concreteness, we will use use Eq.(39) as our visibility figure of merit for the remainder of this paper, but note that careful consideration of how this figure or merit affects a specific application will be required.

Effects of spectral and number impurity
We now investigate how the pump power simultaneously affects the heralded HOM interference visibility and heralding rate.We first consider sources which have separable JSAs, which give rise .We show the visibility using number resolving detectors V HOM PNR (green circles) and threshold detectors, V HOM Thres (purple squares), together with the joint heralding rate for number resolving detectors P Herald PNR (yellow solid curve) and threshold detectors P Herald Thres (blue dashed curve).We use the parameters ζ = 0.1 THz.
to heralded photons that are spectrally pure.To do so we take as an example an idealised JSA that takes the functional form of the product of two Gaussian functions [14], where ∆ν i = ν i − ν i with ν i are the central frequencies of the signal and idler photons and ζ represents their bandwidth.The denominator here represents the Frobenius norm of the exponential factor.In this way the JSA features a single Schmidt coefficient given by ξ.
In Fig. 5 we plot the real part of the JSA of both sources on the left, and on the right we show the heralded HOM visibility and joint heralding rate as a function of the squeezing parameter ξ which represents the pump power.Visibilities are calculated using both threshold (purple squares) and number resolving detectors (green circles).The heralding rate is the probability that both heralding (signal) arms (modes 1 and 4) register photons, and can be thought of as the square of efficiency of one of the sources per excitation pulse.It is defined as for number resolving detectors (yellow curve) and threshold detectors (blue dashed curve), respectively.We see that for pure sources the number resolved visibility is 1 for all values of the squeezing parameter.This reflects the fact that regardless of the power, a heralding event selects single photon states in the idler modes, which then perfectly interfere to give no coincidences.The heralding rate, however, begins to decrease for large squeezing parameters as the probability of only single photon events decreases, and the optimal squeezing parameter is approximately ξ = 0.9, which corresponds to 7.8dB of squeezing.Using threshold detectors, we see that the HOM visibility decreases with increasing squeezing parameter, which is due to the higher order multiphoton terms contaminating the 'single photon' state.The heralding rate monotonically increases as it does not distinguish between the single show the real part of the JSAs (here being equal), while c) is the HOM visibility for number resolving detectors, V HOM PNR (green circles) which matches exactly with Eq. ( 49), and threshold detectors V HOM Thres (purple squares), together with the joint heralding rate for number resolving detectors P Herald PNR (yellow curve) and threshold detectors P Herald Thres (blue dashed curve).We use the parameters ζ = 0.1 THz, L = 29.0ps.photon and multiphoton subspaces, and eventually saturates at 1.As we might expect, threshold detectors do reasonably approximate the number resolving heralding rate for squeezing parameter ξ < 0.2 (squeezing of 1.7dB).Since the sources here produce spectrally pure photons, Fig. (5) serves only to highlight the effects of photon number purity, which we see can be mitigated with the use of number resolving detectors.
We now also include the effects of spectral impurity.To do so we consider two typical waveguide sources which gives rise to non-separable JSAs.The functional form we use is where here the parameter L captures the fields' propagation over the effective length of the twomode squeezing interaction.Such a JSA is achieved for spontaneous parametric down-conversion using a χ (2) non-linearity where ζ represents the pump bandwidth, or using degenerately pumped spontaneous four-wave mixing under a χ (3) non-linearity with ζ being the autoconvolution of the pump.The functional form above is a consequence of symmetric phase-matching with the central frequencies obeying νs +ν i −(n−1)ν p = 0, the central wave-vectors satisfy ks + ki −(n−1) kp +2π/P = 0, for possible poling period P , and the fields' group velocities obey ) with L the interaction length.In all cases the subscripts indicate the signal, idler and pump, and n = 2, 3 is the order of the non-linearity.The normalisation factor in Fig. (43) ensures the JSA has Schmidt coefficients of the form λ l = ξα l with l α 2 l = 1, and the multiplying factor ξ is the (multi-mode) squeezing parameter which is directly related to the strength of the pump laser.
In Fig. (6) we again show the real part of the JSA of each source on the left, and the HOM visibility and heralding rate on the right.For threshold detectors, we see similar trends with increasing squeezing as for the pure sources case, though now with the HOM visibility starting below unity, reflecting the fact that even in the single photon subspace interference is imperfect owing to the spectral impurity of the heralded photons.In contrast to the case explored above however, when using number resolving detectors we see that the HOM interference visibility decreases with increasing squeezing parameter [28].This suggests that even when using number resolving detectors to herald photons only in the single photon subspace, the power, as captured by the squeezing parameter, cannot be increased without detrimentally affecting the interference probability of the photons produced.
We can see this decreasing of the interference visibility by following Ref.[28] and using the Schmidt decomposition in Eq. ( 23), as it allows the two-mode squeezing Hamiltonian in Eq. ( 22) to be written Ĥ = l λ l Ĉ † l D † l + h.c., where Ĉ † l = ´dνψ l (ν)â † 1 (ν) and D † l = ´dνφ * l (ν)â † 2 (ν) are the generalised Schmidt modes as before.Since these generalised modes are independent, the unitary time evolution operator corresponding to a multi-mode two-mode squeezing operation is in fact just a product of independent two-mode squeezers, and we have where Ŝ( 2) Ĉl Dl ] is a two-mode squeezer acting on spectral mode l in spatial modes 1 and 2 , and which in normally ordered form is written [45] Ŝ( 2) where we have written z = re iφ .To find the joint signal-idler state produced by a source including all spectral modes and photon numbers, we can act this on the vacuum to give where the sum runs over all possible integer tuples indicating the number of photons in each spectral Schmidt mode, i.e. n = (n l1 , n l2 , . ..), while |n i = l |n l il with |n l 1l = Ĉ †n l l |vac 1l / √ n l !(and similarly for |n l 2l ) represents the corresponding Fock state in spatial mode i, and the coefficients are To find the interference probability of the idler photons, we consider the conditional state obtained when a single photon in any Schmidt mode is detected in spatial mode 1 , i.e. the signal mode.The corresponding measurement operator is the projector Π 1 = l |1 1| 1l .The probability for such a detection is P PNR (1 ; 1) = Tr[|Ψ Ψ|Π 1 ], giving while the post measurement state is ρ = Tr 1 [|Ψ Ψ|Π 1 ]/P PNR (1 ; 1) with the trace only over modes with spatial label 1 , which gives The purity of this single photon state is a measure of its spectral indistinguishability.For this we find which we see depends on the distribution of the Schmidt coefficients.If one Schmidt mode dominates and is much larger than all others, we have simply Tr[ρ 2 ] ≈ 1.However, if this is not the case, then we see that linearly increasing each coefficient causes a decrease in the idler photon purity.Eq. ( 49) is plotted with black curve in Fig. 6, and perfectly matches the green circular data points as expected. (a)

Effects of photon loss and filtering
While the above analysis in the Fock basis does allow the state produced by a source to be scrutinised in this way, it is not a straightforward matter to describe its evolution through subsequent optical elements.As a pertinent example of this, and one of the main advantages of the formalism presented in this work, we now investigate the effect of frequency selective loss, i.e. filtering.The results so far have demonstrated that non-separability of the JSA is a significant factor affecting source figures of merit, and it is natural to ask to what extent separability of a JSA can be imposed by filtering.
In order to gain insight let us first consider the case in which a fixed amount of frequency independent loss is present on all modes, and investigate the measured interference visibility when using number resolving detectors.The results are shown in Fig. (7), with parts a) and b) corresponding to, respectively and as above, pure sources and non-separable sources.For each the different curves correspond to different loss levels as indicated.We see that with loss there is a greater decrease in visibility with increasing power as compared to the cases without loss.In fact, we see from part a) that even in the case for which the sources are pure and number resolving detectors are used, when loss is included there is a decrease in interference visibility with increasing power, which is not the case without loss.In all cases, the detrimental effects of loss can be understood as compromising the ability for a number resolving detector to herald a truly single photon state in the idler modes.
So far the rate or efficiency of the sources investigated has been characterised by the heralding rate, which is the probability that there is a detection event in both of the heralding (signal) modes.In the absence of loss this quantity is precisely the probability that photon(s) are heralded in both of the idler modes since the photon numbers in the signal and idler modes are perfectly correlated.With the inclusion of loss, however, this is not the case, and it is instructive to also consider the heralding efficiency, defined as the conditional probability of detection events in both of the heralded (idler) modes (2 , 3 ) when no beam-splitter is present, given both the herald (signal) modes (1 , 4 ) registered a detection event.In our case this quantity can be calculated as the ratio of fourfold coincidence events measured with no interference to the heralding rate, and this is therefore written where is the sum of the bunching and antibunching terms and is independent of the beam-splitter angle.Writing the heralding efficiency in this way allows it to be consistently defined even in the case where the beam-splitter angle is fixed.In Fig. 8 we plot the heralding rate and heralding efficiency when including loss as a function of the squeezing parameter for both a) number resolving detection and b) threshold detectors, and in both cases consider non-separable sources with the parameters in Fig. 6.We see that the peak in number resolved heralding rate is shifted to higher squeezing values with increasing loss, suggesting that decreases in photon generation rates causes by losses can in principle be easily overcome by increasing the pump power.However, while Fig. 7 demonstrates that this compensation will necessarily decrease the interference probability, Fig. 8 further shows that increasing the power will decrease the heralding efficiency, meaning that photons are less likely to be present in the idler modes when heralding events are registered.Loss therefore detrimentally affects the effective purity, heralding rate, and heralding efficiency of a source, with none of these figures of merit independently compensated by changing the pump power.With threshold detectors we similarly see that loss-induced decreases in the heralding rate can be compensated by increasing the squeezing parameter.Although the heralding efficiency appears to increase with increasing power for threshold detectors, this reflects only the fact that multiphoton terms increasingly tend to saturate these detectors.
Let us now consider the effect of spectral filtering.We investigate the case in which two sources each described by a non-separable JSA are subject to spectral filtering on both spatial modes (signal and idlers).Each mode is filtered with the same ideal bandpass filter through which photons are either fully transmitted or fully removed depending on their frequency.The bandpass filters are centred at the centre of the JSAs (ν 0 = ν = 193.1 THz) and the sources have the same parameters as in Fig. 6.The results are shown in Fig. 9, where the left column corresponds to no filtering, the middle column to a filter width of ∆ν f = 2×10 15 rad/s, and the right column to ∆ν f = 1.2×10 15 rad/s.The widths of the filers are indicated by the white lines in the plots of the real part of the JSA in the top row, and in all cases the different curves in the plots below correspond to the different levels of loss.As one might expect, we see that tight filtering of the sources means that at low powers an interference visibility of up to 99.58% can be achieved for the initially non-separable imperfect source.The trade off is that the filtering process has reduced the maximum heralding rate from 0.07 to 0.063, and the heralding efficiency decreases to from 100% to 88% at the point of maximal heralding rate.
More importantly, however, is that it appears spectral filtering of this sort does not allow unit visibility for all squeezing parameter values, even with number resolving detectors; a tightly filtered source does not behave like a separable source, even if perfect interference visibilities are found in the low power regime.For a truly pure source, as explored in Fig. 5, the use of number resolving detectors allows the heralding rate to be increased by increasing the pump power at no cost to the interference visibility.In the present case, however, this is no longer true, even when there is no Figure 9: Source figures of merit when using number resolving detectors and no filtering (left column), a band-pass filter with width 2×10 15 rad/s (middle column), and similarly with width 1.2×10 15 rad/s (right column).Unfiltered source has a non-separable JSA as in Fig. 6, and the different curves correspond to different loss levels as in Fig. 8.We see that tight filtering (right column) achieves near unit visibility in the low power limit, however this does not persist with increasing power, and becomes very sensitive to loss.
(frequency independent) loss.Once again, the problem here is that filtering has introduced photon number noise, and number resolving detection is no longer able to herald a truly single photon state.

Structured and non-identical sources
As a final demonstration of the utility of our presented formalism, we will now consider heralded HOM interference between structured and non-identical sources.The pairs of sources that we will consider are shown in the first two rows of Fig. 10.The first column shows a pair of identical sources, with the JSA of each consisting of the sum of two Gaussian functions of the form in Eq. ( 41) with   43), and a separable double lobed source.The right column shows a filtered non-separable source and a double lobed source with each lobe of opposite sign.All sources have squeezing parameter of ξ = 0.4.The bottom row shows the associated four-fold coincidences as measured using threshold detectors, P 4  Thres (black, dashed), and number resolving detectors, P 4 PNR (purple).
widths ζ = 0.03 THz.The four-fold coincidences are shown in the bottom plot, calculated for number resolving detectors (solid, purple curve) and threshold detectors (dash, black).The heralded photons here have double peaked temporal (and spectral) envelops, resulting in a HOM dip with multiple features reflecting the Fourier transform of the product of the two heralded photon spectra.In the middle column we keep one source the same, and replace the other with a filtered waveguide source with a non-separable JSA of the form in Eq. ( 43).The corresponding HOM dip shown in the plot at the bottom again has interesting temporal features, and here the dip minima are higher than in the previous case, reflecting both the spectral impurity of the first source and the poor overlap of the two heralded photons.Finally, in the column on the right we consider again the filtered non-separable waveguide source, now interfering with a double peaked Gaussian with opposite phase, i.e. a JSA that is the difference of two terms of the form in Eq. ( 41).Here we see almost no interference at 0 time delay, as the two heralded photons now have approximately orthogonal spectra.However, for time delays of ∆τ ≈ ±4 ps we see the interference is improved, as for these parameters one of the peaks of the second source is centred on that of the first.We stress that in all cases here, while it may be possible to write down the conditional state of the heralded photons in the Fock basis and describe the qualitative features of these HOM dips, calculating the actual four-fold coincidence count rates and deriving a corresponding interference visibility beyond the weak excitation limit is considerably more involved.Within our continuous variable formalism, on the other hand, the complete state of the system to all photon number orders and with all spectral degrees of freedom is explicitly propagated through all optical elements, allowing for photon detection rates to be calculated exactly.

Conclusion
In this work we have introduced a Gaussian quantum optics formalism that explicitly includes a tensor product structure of the spatial and spectral modes.This formalism allows us to model realistic and relevant highly multimoded fields generated by collections of parametric sources, together with non-Gaussian threshold and number resolving detection.This has allowed us to extend the known expressions for non-frequency resolving threshold and number resolved detection probabilities to include spectral degrees of freedom.We applied this formalism to the study of heralded Hong-Ou-Mandel interference, in which we elucidated the inter-dependencies of the source heralding rates, heralding efficiencies and interference visibilities.In particular, we demonstrated that any nonseparability of the source joint-spectral amplitude results in decreases of the interference visibility beyond the weak excitation limit, and further decreases in the interference visibility and heralding efficiency when any loss is present.Furthermore, we found that while spectral filtering can improve a source's interference visibility, the filtering process necessarily introduces photon number noise which detrimentally affects the interference visibility and heralding efficiency at higher powers.
The HOM setup considered here constitutes a fundamental primitive for many more complex systems, and as such understanding these effects will become ever more important as demonstrations (and applications) are increased in scale.Consider for instance, current state-of-the-art 4-photon experiments achieving detection rates of 10 −2 s −1 [46], 1 s −1 [47,48], 1 − 2 s −1 [49], 21 s −1 [50], and 46 s −1 [51], with target state fidelities not exceeding 0.94.Whilst state infidelity is routinely attributed to multi-photon events, with heuristic or idealised error models of such often presented [52], an exact analysis including multiphoton events alongside spectral effects and detection mechanisms has yet to be included.As systems are required to provide higher count rates, source figures of merit including purity and multi-photon contamination are traded in favour of achieving acceptable count rates.Recently a 10-photon experiment was demonstrated at a rate 4 per hour [53] at the expense of a state fidelity of 0.573, and a subsequent 12-photon experiment achieved fidelity of 0.576 at a rate 1 per hour [54].Of course, the more efficient generation of higher photon-number states for computation applications requires multiplexing [55], and one would necessarily have to chose an operation regime in which generation rates are traded off against multi-photon contamination, and a realistic noise model capturing multi-photon and spectral impurity is then required in order to demonstrate that acceptable error thresholds are indeed overcome.In particular, it is clear that leakage errors, associated to the higher dimensional state space of multi-photon states (outside the desired photon number subspaces) will present a key challenge for both the experimental systems, and the theory of their operation [56].
In this work we have primarily focused on the generation of heralded single photon states, which could in turn to be used to generate larger states with fixed photon number.However, our formalism also naturally lends itself to applications where states without fixed photon number are of interest, most notably those of Gaussian boson sampling, which are based on the observation that calculating photon detection probabilities from Gaussian states with many modes is computationally hard [29,57,58].Using the present formalism it would be interesting to investigate the extent to which photon spectral impurity or mutual distinguishability can affect the computational complexity of multi photon detection probabilities.It is in turn worth mentioning that the number resolving detection probabilities derived here make use only of derivatives of vacuum projections and combinatorics, and it would be interesting to explore the extent to which these new expressions could lead to insights into the development of classical algorithms for simulating photonic measurements of this sort.

<
l a t e x i t s h a 1 _ b a s e 6 4 = "Q o / 4 h / t K o 9 J p f M p C r v m k K F 8 Q y x 8 = " > A A A B 7 H i c b V B N S 8 N A F H y p X 7 V + V T 1 6 W S y C p 5 K I o s e i F 4 8 V T F t o Q 9 l s N + 3 S z S b s v g g l 9 D d 4 8 a C I V 3 + Q N / + N 2 z Y H b R 1 Y G G b e s O 9 N m E p h 0 H W / n d L a + s b m U b g r f 4 8 l / S P K l 6 5 1 X 3 7 q x S u 8 r j K M I B H M I x e H A B N b i B O j S A w Q C e 4 A V e H e E 8 O 2 / O + 7 y 1 4 O Q z + / A L z s c 3 z n m N e w = = < / l a t e x i t > M 3

Figure 2 :
Figure 2: Conceptual fan-out circuit used to count the number of photons in an arbitrary Gaussian state σ S using threshold detectors.and has the matrix elements (U QFT ) nn = exp[2πinn /m]/ √ m.The covariance matrix following the action of the QFT is given by fold coincidences Beamsplitter angle θ (radians)

Figure 3 :
Figure 3: Variable time delay (top row) and variable beam-splitter (bottom row) heralded Hong-Ou-Mandel (HOM) measurements.Parts a) and d) show schematics of the optical circuits, while parts b) and e) show the four-fold coincidences as measured using threshold detectors, P 4Thres (black), and parts c) and f) using number resolving detectors, P 4 PNR (purple).Parameters correspond to sources with non-separable JSAs and operating beyond the weak-excitation limit.

Figure 4 :
Figure 4: The ratio of maximum four-fold coincidences, measured with large time delay ∆τ in Fig. 3 a), to that measured with zero beam-splitter angle θ = 0 in Fig. 3 d), i.e. (P 4 D ) τ =∞ /(P 4 D ) max = R(ξ, L), and for (a) threshold detectors and (b) number resolving detectors, for different loss values as indicated.

Figure 5 :
Figure 5: Figures of merit for interference between two sources with separable joint-spectral-amplitudes (JSAs).Parts a) and b) show the real part of the source JSAs (which are here equal), part c) shows the joint heralded HOM visibility and heralding rate as a function of squeezing parameter ξ (i.e.pump power).We show the visibility using number resolving detectors V HOM PNR (green circles) and threshold detectors, V HOM Thres (purple squares), together with the joint heralding rate for number resolving detectors P Herald PNR (yellow solid curve) and threshold detectors P Herald Thres (blue dashed curve).We use the parameters ζ = 0.1 THz.

Figure 6 :
Figure 6: Figures of merit corresponding to two sources with non-separable (JSAs).Parts a) and b)show the real part of the JSAs (here being equal), while c) is the HOM visibility for number resolving detectors, V HOM PNR (green circles) which matches exactly with Eq. (49), and threshold detectors V HOM

Figure 7 :
Figure 7: Interference visibility using number resolving detectors as a function of squeezing parameter (power), with each curve corresponding to different loss values as indicated.Part a) corresponds to two pure sources with separable JSAs, and part b) corresponds to sources with non-separable JSAs.Parameters as in Fig. 5 for part a) and Fig. 6 for part b).

Figure 8 :
Figure 8: Heralding rates and efficiencies when using part a), number resolving detectors and part b), threshold detectors.The curves indicate the heralding efficiencies with different values of loss as indicated, and the plot points indicate the corresponding heralding rates.

Figure 10 :
Figure10: Variable time delay HOM measurements for more exotic sources.The left columns shows the real part of two identical and separable sources, each consisting of two Gaussian lobes with width 0.03 THz.The middle column corresponds to a filtered non-separable waveguide source as in Eq. (43), and a separable double lobed source.The right column shows a filtered non-separable source and a double lobed source with each lobe of opposite sign.All sources have squeezing parameter of ξ = 0.4.The bottom row shows the associated four-fold coincidences as measured using threshold detectors, P4Thres (black, dashed), and number resolving detectors, P 4 PNR (purple).