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

Almost all tasks in optical implementations of quantum communication,^1–3 quantum cryptography,^4–6 quantum sensing, and quantum information processing^7–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, 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¹⁴ and can in many cases be more readily incorporated into integrated devices with 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 toward making a deterministic source is to instead use single quantum emitters such as semiconductor quantum dots,²¹ defect centers in crystals and 2D materials,¹² 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–26 However, there are often efficiency–purity trade-offs using such sources,²⁷ 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²¹ leave significant questions regarding their utility in information processing applications.

For these reasons, the advantages offered by parametric sources motivate a detailed study of their operation beyond the weak excitation limit. While 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,²⁸ 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 toward 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²⁹ and number resolving detectors,³⁰ 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 and 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 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,35 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 they provide 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.

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 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

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ωNf)⊤$ 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_sN_f. The mode operators satisfy the usual commutation relations $[âiω,âjω′†]=δijδωω′$⁠.

A variety of candidates exist for parametric photon sources, such as free-space systems,³⁶ cavity-based systems,³⁷ in-line or wave-guided systems,³⁸ and inter-modal phase-matched systems.¹⁴ 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 non-linearity 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 can in all generality be written as

where $H$ is a matrix of scalar coefficients, which may be a discrete approximation to a continuous function under appropriate regularity conditions.³⁸ Although we are interested here in optical systems, this form of Hamiltonian and the analysis below is rather general and may also be applied to any quadratic bosonic Hamiltonian, such as those arising in condensed matter systems³⁹ as well as quadratically coupled oscillator models in studies of quantum Brownian motion and decoherence.^40,41 The Hamiltonian has a corresponding unitary time evolution operator $Û=exp[−iĤ]$⁠, and using $Ĥ=Ĥ†$ and the basic commutation relations for the creation and annihilation operators, one can show⁴² 

where

with $1N$ being 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⁻¹, where M is defined as a linear symplectic matrix.⁶⁵ Equation (3) provides us with a way to use the underlying Hamiltonian parameters to propagate the collection of mode operators in the Heisenberg picture.

Rather than working directly with the quantum state of the optical modes $ρ̂$⁠, we instead consider the corresponding characteristic functions.⁴² We define the s-ordered characteristic function as

where $D̂(Λ)=exp†KΛ$ and Λ is a 2N-dimensional vector of complex variables. Equation (6) describes a mapping from the quantum state of the field $ρ̂$ in Hilbert space to a multivariate function in an associated phase space. In this phase space representation, the nature and properties of the field (e.g., its energy and coherence) are reflected in the specific functional form of the characteristic function. 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 $12N$⁠, 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 $vacvacS$⁠. 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³¹ 

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=12|S|$⁠, we find P_off(S) = 1 as expected.

We now consider the probabilities associated with photon detection using threshold detectors, extending the existing results that include only multiple spectral modes³¹ or multiple spatial modes.²⁹ 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 in which the detection events take place, while $Θ={ω1,…,ωNf}$ 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 $S={1,3}×{ω1,…,ωNf}={1ω1,…,1ωNf,3ω1,…,3ωNf}$⁠.

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 $2S$ 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 involving the Torontonian function in Ref. 29, which gives the probability for threshold detector click patterns in the absence of unresolved spectral degrees of freedom (i.e., for N_f = 1). We note that the number of spectral degrees of freedom increases the size of the reduced covariance matrices in Eq. (8), while the number of explicit terms in the sum in Eq. (10) depends exponentially only on the number of spatial modes in which photon(s) are detected. In this way, the computational difficulty of these threshold detector probabilities asymptotically scales exponentially and is dominated by the total number of detection events. This exponential scaling is the origin of the intractability of Gaussian boson sampling with threshold detectors.²⁹ 

We now present expressions for probabilities associated with number resolving detection in spatial modes. Our method takes inspiration from the way in which pseudo-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 one 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 that 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. 44.

1. 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 as $σS×S=σS⊕12|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 the initial amplitude in σ_S and has the matrix elements $(UQFT)nn′=exp[2πinn′/m]/m$⁠. The covariance matrix following the action of the QFT is given by

   $σS×S′=UQFT⊗12|S|σ̃S+12|S|⋯0⋮⋱⋮0⋯0+12|S|UQFT†⊗12|S|=1mσ̃S⋯σ̃S⋮⋱⋮σ̃S⋯σ̃S+12|S|⋯0⋮⋱⋮0⋯12|S|=1mσ̃S⃗(12|S|⃗)⊤+12m|S|,$

   (13)

   where we have defined $σ̃S=σS−12|S|$ and $σ̃S⃗=(σ̃S,…,σ̃S)⊤$ as the block-constant vector with each of its m elements equal to $σ̃S$ and $12|S|⃗$ is similarly a block-constant vector with all entries being the identity matrix $12|S|$⁠. We see that this procedure results in a block rank 1 form for the total covariance matrix $σS×S′$⁠.
2. We next determine the probability to detect zero clicks (the vacuum) in some subset of the supermodes $B$⁠. The quantity of interest is

   $PoffB×S=det(12|S‖B|+σB×S)/2−1/2,$

   (14)

   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 $12|S|⃗$ instead having length $|B|=b$⁠. Using the matrix determinant lemma and commutative subring properties, we can write⁶⁶ 

   $Poff(B×S)=poff(b)=det12|S|+b2mσ̃S−1/2,$

   (15)

   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.
3. 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,

   $Pon,off(X,Y)=∑B∈2X(−1)|B|Poff((B∪Y)×S).$

   (16)

   Equation (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

   $Pon,off(X,Y)=∑l=0k(−1)lklpoff(m−(k−l)),$

   (17)

   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

   $∑|X|=kPon,off(X,Y)=mk∑l=0k(−1)lklpoff(m−(k−l)).$

   (18)
4. 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 recognize 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,

   $PPNR(n)=limm→∞∑|X|=nPon,off(X,Y)=(−1)nn!∂tndet12|S|+t2σ̃S−1/2t=1.$

   (19)
5. 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=(n1,…,n|S|)$ in spatial modes $S$ is given by

   $PPNR(S;n)=∏i=1|S|(−1)nini!∂tini×det12|S|+TSσ̃STS2−1/2T=12|S|,$

   (20)

   where $TS=12⊗(t⊗1|S|)$ with $t=diag(t1,…,t|S|)1/2$⁠.

Equation (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 unresolved spectral degrees of freedom may be present. It should be noted that although our expression in Eq. (20) is valid for any finite detector number pattern n and 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.²⁹ 

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 affecting the scaling of the computation with respect to photon number. To do so using the expressions from 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.

The derivatives in Eq. (20) can in principle be evaluated analytically using repeated applications of the Jacobi identity for derivatives of determinants. This may be a promising approach to pursue in order to draw conclusions regarding the computational complexity of photon number probabilities in the presence of multiple spectral modes. Being primarily interested in photon number probabilities for relatively small total photon numbers, in this work, we use a finite difference approximation to numerically evaluate the derivatives by taking linear combinations of the vacuum probabilities, in which the relevant function must be evaluated at $∏i=1|S|(ni+1)$ points for different values of T_S.

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$ is 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 italicized 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, $PThres(S;n)≠PPNR(S;n)$⁠, so this common notation only indicates a correspondence and not equality.

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 ν. We will refer to the function F(ν₁, ν₂) 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,^45,46 namely,

where λ_l are positive coefficients known as the Schmidt coefficients and the functions {ψ_l(ν)} and {ϕ_l(ν)} are sets of orthonormal functions (though not necessarily equal or mutually orthonormal). In the low squeezing limit ‖F(ν₁, ν₂)‖ ≪ 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(ν₁, ν₂) means that the Schmidt basis defined by the functions {ψ_l(ν)} and {ϕ_l(ν)} must be found by solving integral eigenvalue equations. In practice, it is often easier to instead discretize the integral in Eq. (22), which, again, is valid under reasonable regularity conditions.³⁸ 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 discretization 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 the 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).

Equation (26) allows the Hamiltonian coefficients to be written as $H=UFDU†$⁠, 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 as $M=UMDU†$⁠, where

We see that a multimode two-mode squeezer is simply a set of independent two-mode squeezers acting on the appropriate Schmidt modes.

The formalism presented here, which represents optical elements as symplectic matrices, requires that the frequency degree of freedom is discrete. In order that frequency features in the (continuous) JSA are accurately captured, we must choose our discretization step small enough that these features are resolved. We note, however, that the discretization need not necessarily be performed uniformly, and in fact, the Schmidt decompositions above provide one with a formal way in which to identify a small number of dominant frequency modes even in broadband cases where a large span of the frequency range may be important. Although we consider non-frequency resolving detection in this work, we also note that a discretization of frequencies can quite naturally and accurately reflect the timing or spectral resolution of detectors used.

In addition to two-mode squeezers that 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 α^† = α⁻¹. For a dispersionless (frequency independent) beam splitter, we have

while a dispersionless phase shifter acting on say spatial mode 1 of 2 is described by $αPS(ϕ)=diag(exp[iϕ],1)⊗1Nf$⁠. 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(ατ,1Nf)$ with spectral matrix elements $(ατ)ω1ω2=δω1ω2⁡exp[iω1δντ]$⁠. 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 that 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 $USSUSS†+USLUSL†=12|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 $USS=diag(USS,USS*)$⁠, and we have omitted modes that we will shortly trace out. The top left block pertains to our modes of interest and can be written using $σ̃S≔σS−1|2S|$ as $USS(σ̃S+12|S|)USS†+USL12USL†$⁠. Then, using the unitary condition, we can rearrange to find

This map, described uniquely by the sub-matrix $USS$ 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 modeling frequency-independent loss on all spatial modes, we can take $USS=1−ϵ1|S|$ with 0 ≤ ϵ ≤ 1 being the loss parameter. For loss acting on, for example, spatial mode 1 of 2, we have $USS=diag(1−ϵ,1)⊗1Nf$⁠. Spectral filtering can be included as a frequency-dependent loss with an associated filter function f(ν) such that $(USS)ω1ω2=δω1ω2f(ω1δν)$⁠. For example, a bandpass filter with central frequency ν₀ and bandwidth Δν_f is described by

This provides us with all the optical components necessary to model a HOM interference experiment.

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¹⁴ 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₃ have a second order non-linearity and are able to use the polarization degree of freedom to create a time delay.⁴⁷ 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 the interference between heralded photons from two sources, which we term heralded HOM interference. The optical circuits are shown in Figs. 3(a) and 3(d), and are 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 fourfold coincidence terms,

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,

We envisage simultaneously pumping both sources with, for example, a common optical pump field. In this case, the probabilities in both Eqs. (36) and (37) correspond to measured counts per excitation pulse, assuming the specified level of losses.

In Figs. 3(b), 3(c), 3(e), and 3(f), we plot the fourfold 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 [(b) and (e)] correspond to coincidences measured with threshold detectors, while the purple curves [(c) and (f)] correspond to the same quantity using number resolving detectors. For the varying beam splitter case, referred to as the Mach-Zehnder HOM interferometer, the fourfold coincidence probability has a maximum for a beam splitter angle of θ = 0 (acting trivially on modes 2 and 3), which we label $(PD4)max$⁠. This is compared to the corresponding value for a balanced beam splitter with θ = π/4, which we label $(PD4)min$⁠. The visibility is typically defined as⁴⁸ 

On the other hand, for a fixed beam splitter angle of θ = π/4 and variable time delay, the coincidences are typically normalized with respect to a large (effectively infinite) delay, which we label $(PD4)τ=∞$⁠. This is then compared to the zero time delay coincidence probability $(PD4)τ=0$⁠, which is equal to $(PD4)min$⁠. We then have the visibility metric,⁴⁹ 

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 $(PD4)min$ term in the denominator of Eq. (38) results in different expressions, as the maximum coincidence probability for variable beam splitter $(PD4)max$ is not equal to the maximum coincidence probability for variable time delay $(PD4)τ=∞$⁠, leading to a different normalization in each case. While in the low power limit $(PD4)max=2(PD4)τ=∞$⁠, which reflects the relative number of single photon pathways leading to a coincidence and gives a simple 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=(PD4)max/(PD4)τ=∞$⁠, which allows us to write

We see that with access only to coincidence probabilities measured with variable beam splitter, the 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 $VPNRHOM$ 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 Eq. (39) as our visibility metric for the remainder of this paper, but note that careful consideration of how this figure of merit affects a specific application will be required.

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 to heralded photons that are spectrally pure. To do so, we take as an example an idealized JSA that takes the functional form of the product of two Gaussian functions,¹⁴ 

where $Δνi=νi−ν̄i$ with $ν̄i$ being 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 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 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 the efficiency of one of the sources per excitation pulse. It is defined as

and

for number resolving detectors (yellow curve) and threshold detectors (blue dashed curve), respectively. We see that for pure sources, the visibility with number resolving detectors 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 exhibit perfect HOM interference 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. As can be seen in Eq. (36), with threshold detectors, all events involving four or more photons are registered, and this decrease in visibility is due to the higher order multiphoton terms contaminating the “single photon” state; threshold detectors register events in which one or both sources produce more than one photon in an idler mode, for which perfect HOM interference does not occur. The heralding rate monotonically increases as it does not distinguish between the single 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 the parameter $L̃$ captures the fields’ propagation over the effective length of the two-mode squeezing interaction. Such a JSA is achieved for spontaneous parametric down-conversion using a χ⁽²⁾ non-linearity where ζ represents the pump bandwidth or using degenerately pumped spontaneous four-wave mixing under a χ⁽³⁾ non-linearity with ζ being the width of 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 $k̄s+k̄i−(n−1)k̄p+2π/P=0$⁠, for possible poling period P, and the fields’ group velocities obey $vi−1=2vp−1−vs−1$ so that $L̃=L(vs−1−vp−1)$ with L being 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 normalization factor in Eq. (43) ensures the JSA has Schmidt coefficients of the form λ_l = ξα_l with $∑lαl2=1$⁠, and the multiplying factor ξ is the (multi-mode) squeezing parameter, which is related to the strength of the pump laser. While in the high-power regime, other non-linear effects such as phase-modulation can become appreciable,³⁸ in this work, we neglect these effects and vary the pump strength solely through the parameter ξ. In this way, we can separate out the impact of the photon number statistics independently of the more application or material specific nonlinear effects that may also occur.

In Fig. 6, we again show the JSA of each source on the left and the heralded HOM visibility and joint 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.²⁸ 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 as $Ĥ=∑lλlĈl†D̂l†+h.c.$⁠, where $Ĉl†=∫dνψl(ν)â1†(ν)$ and $D̂l†=∫dνϕl*(ν)â2†(ν)$ are the generalized Schmidt modes as before. Since these generalized 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 $Ŝl(2)(z)=exp[zĈl†D̂l†−z*ĈlD̂l]$ is a two-mode squeezer acting on spectral mode l in spatial modes 1 and 2, and which in normally ordered form is written as⁵⁰ 

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=(nl1,nl2,…)$⁠, while $|n〉i=⨂l|nl〉il$ with $|nl〉1l=Ĉl†nl|vac〉1l/nl!$ (and similarly for $|nl〉2l$⁠) representing the corresponding Fock state in spatial mode i, and the coefficients are $c(n)=∏lsechλl(−i⁡tanhλl)nl$⁠.

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=∑l111l$⁠. The probability for such a detection is $PPNR(1;1)=Tr[ΨΨΠ1]$⁠, giving

while the post measurement state is $ρ̂=Tr1[ΨΨΠ1]/PPNR(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. Equation (49) is plotted with a black curve in Fig. 6 and perfectly matches the green circular data points as expected.