Characterizing multi-photon quantum interference with practical light sources and threshold single-photon detectors

The experimental characterization of multi-photon quantum interference effects in optical networks is essential in many applications of photonic quantum technologies, which include quantum computing and quantum communication as two prominent examples. However, such characterization often requires technologies which are beyond our current experimental capabilities, and today's methods suffer from errors due to the use of imperfect sources and photodetectors. In this paper, we introduce a simple experimental technique to characterise multi-photon quantum interference by means of practical laser sources and threshold single-photon detectors. Our technique is based on well-known methods in quantum cryptography which use decoy settings to tightly estimate the statistics provided by perfect devices. As an illustration of its practicality, we use this technique to obtain a tight estimation of both the generalized Hong-Ou-Mandel dip in a beamsplitter with six input photons, as well as the three-photon coincidence probability at the output of a tritter.


I. INTRODUCTION
Multi-photon quantum interference is a key concept in quantum optics and quantum mechanics.
It has been extensively studied by many authors over the last decades, going from the seminal two-photon interference experiment performed by Hong, Ou and Mandel [1] to more recent experimental demonstrations which involve a higher number of indistinguishable photons in various scenarios [2][3][4][5][6][7][8]. Moreover, besides its indubitable inherent theoretical interest, multi-photon quantum interference also plays a pivotal role on several subfields and applications of quantum information science that use, for example, optical networks (ONs) to interfere photons. These applications include, among others, quantum computing [9], quantum cryptography [10,11], boson sampling [12][13][14][15][16], quantum clock synchronization [17], and quantum metrology [18]. In any practical realization of these applications it is essential to experimentally confirm that the photons interfere as desired [19].
Unfortunately, however, to experimentally characterize multi-photon quantum interference on general ONs is usually a quite challenging task [20]. This is so because, for this, one would ideally need to use high-quality on-demand nphoton sources which are yet to be realized [21,22], together with high-quality photon-number resolving (PNR) detectors, which, besides of being expensive experimental resources, currently can only distinguish up to a certain number of photons and may also introduce noise [23][24][25]. As a result, we have that current experimental techniques to characterize the quantum interference behaviour of ONs at a few photons level typically suffer from inevitable errors due to the use of imperfect sources and detectors [20].
The main contribution of this paper is a novel technique to experimentally estimate the input-output photon number statistics of ONs when the input signals are tensor products of Fock states. For this, we use simple laser sources to generate the input signals to the ON and practical threshold single-photon detectors to measure the output signals [26].
That is, our method is implementable with current technology and allows the estimation of the conditional probability distribution P (x 1 , ..., x M |n 1 , ..., n N ) that describes the behaviour of the ON on the input Fock states |n 1 , ..., n N , where n i (x j ), with i = 1, ..., N (j = 1, ..., M ), denotes the number of photons at the ith (jth) input (output) port of the ON. We emphasize, however, that, in practice, our method is specially suitable to evaluate mainly small-size ONs. This is so because, as we show later, it requires to experimentally estimate the probabilities of certain observable events whose estimation complexity may increase exponentially with the number of input/output ports of the ON [12,27].
The key idea builds on two techniques that are extensively used in the field of quantum cryptography: the decoystate method [28][29][30] and the so-called detector-decoy technique [31,32]. We use the former at the input ports of the ON to estimate the statistics provided by ideal n-photon sources. Besides standard quantum key distribution [28][29][30][33][34][35], the decoy-state method has also been used for example to estimate the yield of two single-photon pulses in measurement-device-independent quantum key distribution [10,36,37], to simulate single-photons sources with imperfect light sources [38], and to perform singlephoton quantum state tomography with practical sources [39]. That is, so far the use of the decoy-state method has been limited to evaluate the behaviour of ONs when they receive single-photon pulses at their input ports. Here we extend its use to estimate the behaviour of ONs in the general case where they receive as input signals multi-photon pulses. Furthermore, we employ the detector-decoy method at the output ports of the ON to estimate the statistics provided by ideal PNR detectors [31,32].
To illustrate the practicality of our technique to study ONs, we evaluate two simple examples of interest. In the first one, we estimate the generalized Hong-Ou-Mandel (HOM) dip [1] in a beamsplitter when the total number of input photons is six for two different conditional probabilities, P (3, 3|3, 3) and P (5, 1|5, 1). The first case has been experimentally studied in [5], where the authors used for this a spontaneous parametric down-conversion source in combination with a measurement setup with six threshold single-photon detectors. The second case, however, (to the best of our knowledge) has not been experimentally implemented yet due to the difficulty of generating fivephoton states to input the beamsplitter. In both scenarios we use our method to estimate the HOM dip by means of just two laser sources and two threshold single-photon detectors. In the second example, we estimate the three coincidence detection probability in a tritter [40] when there is just one single-photon pulse in each of its input ports, i.e., we estimate P (1, 1, 1|1, 1, 1). This example is also used to obtain a high precision estimation of the dependence of that probability with the triad phase, which arises when one considers more than two input photons [40]. While these two examples correspond to evaluating linear ONs, we remark that our method could also be used to study multi-photon quantum interference in non-linear ONs.
The paper is organized as follows. In Sec. II, we present our method in detail. Then, in Sec. III we evaluate the two practical examples described above. Finally, we summarize the content of the paper in Sec. IV. The paper also includes two appendixes with additional information.

II. METHOD
As already mentioned above, we use the decoy-state (detector-decoy) method at the input (output) ports of the ON to estimate the statistics provided by ideal n-photon sources (PNR detectors). Of course, in contrast to the case where one really uses perfect n-photon sources and PNR detectors, the use of decoy settings does not provide single shot resolution about how many photons input and output each port of the ON each given time. However, it permits to estimate the full statistics that such perfect devices could give, which is enough for our purposes.
More precisely, we use as input signals to the ON Fock diagonal states with different photon-number statistics. This type of signals could be generated, for instance, with attenuated laser diodes emitting phase-randomised weak coherent pulses (WCPs), triggered spontaneous parametric down-conversion sources or practical single-photon sources, together with variable attenuators to vary the intensity of the different light pulses. To implement the detector-decoy method, on the other hand, we place variable attenuators also on the output ports of the ON together with threshold single-photon detectors. This general scenario is illustrated in Fig. 1. In so doing, as we show below, we have that the probability of each possible detection pattern observed on the threshold single-photon detectors can be written as a sum of linear terms where the only unknowns are the  [28][29][30] and the detector-decoy technique [31,32]. More precisely, we place at each input port i = 1, ..., N, of the ON a source of phase-randomized WCPs together with a variable attenuator of transmittance γi. At each output port j = 1, ..., M of the ON we place a variable attenuator of transmittance ωj and a threshold single-photon detector Dj. The input (output) spatial modes are denoted in the figure with the letters ai (bj), and the output detection pattern of click and no-click events given by the threshold single-photon detectors is denoted by θ.
probabilities P (x 1 , ..., x M |n 1 , ..., n N ) ≡ P (x|n) (where, for easy of notation, we use x ≡ x 1 , ..., x M and n ≡ n 1 , ..., n M in what follows). As a result, we obtain a set of linear equations which are function of the probabilities P (x|n) and, in principle, one can estimate these quantities accurately. The more decoy-state/detector-decoy settings we use, the higher the number of linear equations that we obtain and, thus, the better the accuracy of the estimation. Indeed, in the asymptotic limit where one uses an infinite number of decoystate/detector-decoy settings then the probabilities P (x|n) could be estimated precisely. Importantly, however, as we show below, already a small number of decoy-state/detectordecoy settings can typically provide a quite tight estimation of P (x|n) for small values of n and x.
Our starting point is the input state to the ON. As shown in Fig. 1, this is the state of the N spatial modes after the input attenuators of transmittance γ i . This state can be written as where ρ µi i = ∞ ni=0 p µi ni |n i n i | is the Fock diagonal state at the ith input spatial mode of the ON, which in the case of phase-randomized WCPs satisfies p µi ni = e −µi µ ni i /n i !. Here, the mean photon number µ i = γ i µ , with µ being the initial intensity of the laser sources. The quantity P µ n = N i=1 p µi ni , on the other hand, represents the conditional probability of having the input state |n ≡ |n 1 , ..., n N given the set of input intensities µ = {µ 1 , ..., µ N }.
Let us now consider the output state ρ µ out = U ρ µ in U † of the ON, where U denotes the evolution unitary operator applied by the network. We can write this state in terms of the probabilities P (x|n). For this, for convenience, we first combine the effect of each output attenuator ω j with the detection efficiency of each threshold single-photon detector D j (see Fig. 1). By doing so, we can conceptually consider that at the jth output port of the ON there is now a threshold single-photon detector with efficiency κ j = ω j η D , for j = 1, ..., M , where η D is the detection efficiency of the threshold single-photon detector D j in the original scenario (note that here, for simplicity, we assume that all detectors D j have the same detection efficiency η D ). This is so because when a detector has some finite detection efficiency η D it can be mathematically described by a beamsplitter of transmittance η D combined with a lossless detector [41]. Importantly, since the positive-operator valued measure (POVM) that characterizes the behaviour of a typical threshold singlephoton detector is diagonal in the Fock bases, it follows that the resulting measurement statistics when measuring ρ µ out remain unchanged if, before the actual measurements, we perform a quantum nondemolition (QND) measurement of the total number of photons at each output mode of the ON. This means, in particular, that for any ρ µ out , there is always a Fockdiagonal state, which we shall denote byρ µ out , of the form that provides exactly the same measurement statistics as ρ µ out . In Eq. (2), |x ≡ |x 1 , ..., x M is the Fock state after the QND measurements on ρ µ out , and P (x|n) = | x|U |n | 2 denotes the conditional probability of having such state |x given that the input state to the ON is |n . Note that here, for simplicity, we consider passive networks that do not create photons, and therefore we assume that That is, the total number of output photons cannot be greater than the total number of input photons to the ON. This last condition is expressed in Eq. (2) with the symbol x ≤ n. However, we remark that our method could be applied as well to evaluate active ONs.
Finally, to estimate the unknown probabilities P (x|n) we need to relate them with some observable quantities. For this, we use the fact that the probability P µ,κ θ of observing the detection pattern θ ≡ (θ 1 ...θ M ), where θ j is equal to zero (one) for a no-click (click) event in the threshold singlephoton detector D j , given the stateρ µ out and the detectors' efficiencies κ = κ 1 , ..., κ M , is given by with the POVM elements Π κj θj given by where p dark denotes the dark-count probability of the detector D j , which for simplicity we assume is equal for all j = 1, ..., M . That is, the operator Π κj 0 (Π κj 1 ) is associated to a no-click (click) event at the detector D j . After substituting Eqs. (2) and (4) in Eq. (3), we finally obtain Here, P κ (θ|x) = x| ⊗ M j=1 Π κj θj |x denotes the probability of observing the detection pattern θ given the output state |x , the detection efficiencies κ and the dark count probability p dark . If the detectors D j are well-characterized, this quantity is known. Importantly, Eq. (5) relates the observed probabilities P µ,κ θ , which can be directly measured in the actual experiment, to the unknown probabilities P (x|n) via the statistics P µ n and P κ (θ|x), which are both known a priori given the experimental parameters µ , η D and p dark together with the attenuator settings γ = {γ 1 , ..., γ N } and ω = {ω 1 , ..., ω M }. Indeed, as already mentioned previously, each decoy/detector-decoy setting provides a new linear equation which has the same unknowns P (x|n) but different coefficients P µ n and P κ (θ|x) and constant terms P µ,κ θ . Then, by solving the set of linear equations given by Eq. (5) one can, in principle, estimate any conditional probability P (x|n).
In what follows, we illustrate this method with two simple examples of practical interest.

A. First example: beamsplitter
In this case, we have that the creation operators,â † 1 andâ † 2 , for the input modes of a beamsplitter and those,b † 1 andb † 2 , for its output modes satisfy the relationsb † where the parameters r, t, r and t fulfill |t| 2 + |r| 2 = 1, |t| = |t |, |r| = |r | and t r + r t = 0 [42]. That is, if the state at the input spatial modes a 1 and a 2 is say |n 1 , n 2 a1,a2 (i.e., it consist in n 1 and n 2 indistinguishable photons respectively), the state at the output modes b 1 and b 2 is given by the following coherent superposition of Fock states where, for simplicity, we have considered the particular case in which r = − √ 1 − η, r = −r and t = t = √ η, with η being the transmittance of the beamsplitter. From Eq. (6) one could directly theoretically calculate the probability distribution P (x 1 , x 2 |n 1 , n 2 ) = | ψ out |x 1 , x 2 | 2 of finding, respectively, x 1 and x 2 photons at the output ports b 1 and b 2 of the beamsplitter given that there are n 1 and n 2 photons at its input ports a 1 and a 2 . Importantly, according to quantum mechanics the value of this probability strongly differs from that of a classical scenario, where the photons are considered distinguishable particles which do not interfere. The HOM dip [1] is a well-known example of this fact. Indeed, when two photons input a 50 : 50 beamsplitter through a different input port, classical mechanics predicts a probability equal to 1/2 of finding the two photons at different output ports of the beamsplitter, while quantum mechanics predicts (for indistinguishable photons) that this probability is equal to zero. In general, this difference between the predictions of quantum and classical mechanics can be quantified by means of the visibility, which is defined as where the subindex c denotes the classical case, i.e., when the photons are perfectly distinguishable. Eq. (7) has been experimentally evaluated in many different experiments over the last years. For instance, in [4] and [5], the authors obtain visibilities V 2,2|2,2 equal to 88% for a four-photon interference scheme within an asymmetric beamsplitter and V 3,3|3,3 equal to 92% for a six-photon interference scheme, respectively. For this, they use type-II parametric down-conversion sources to generate pairs of entangled photons and a measurement setup with four ans six threshold single-photon detectors, respectively, in combination with beamsplitters. Also, in the experiment reported in [8], the authors interfere two bosonic atoms (instead of photons) and they observe a visibility equal to about 65%.
We now apply our method based on two sources of phaserandomized WCPs and two threshold single-photon detectors to evaluate the visibility V x1,x2|n1,n2 . Like in the general case considered in the previous section, it is straightforward to show that by varying the intensity µ i of the input signals at the i-th input port of the beamsplitter as well as the attenuator's transmittance ω j (and thus the effective detector's efficiency κ j ) at its j-th output port, with j = 1, 2, one can generate an arbitrary number of inequalities that involve the unknown probabilities P (x 1 , x 2 |n 1 , n 2 ). The final system of linear equations, particularized from Eq. (5), is given by for each one of the four possible detection patterns θ ≡ (θ 1 θ 2 ) ∈ {00, 01, 10, 11}. Again, in a real experiment the probabilities P µ n1,n2 = e −(µ1+µ2) µ n1 1 µ n2 2 /(n 1 !n 2 !) and P κ (θ|x 1 , x 2 ) = x 1 , x 2 | ⊗ 2 j=1 Π κj θj |x 1 , x 2 , with Π κj θj given by Eq. (4), are known given the experimental sets µ and κ, as well as the value of the dark count probability of the detectors, while the probabilities P µ,κ θ can be directly observed in the experiment, once performed. For our simulations we use as observed values P µ,κ θ those predicted by quantum mechanics (see A for more details).
To solve the set of linear equations given by Eq. (8) one can use analytical or numerical tools. For simplicity, here we solve Eq. (8) numerically. For this, we first transform the set of equalities given by Eq. (8), which contains an infinite number of unknowns P (x 1 , x 2 |n 1 , n 2 ), into a set of inequalities with a finite number of unknowns, as shown in B. Also, we use the linear programming solver Gurobi [43] and the Matlab interface Yalmip [44].
Just as an example, Fig. 2 shows our results for the conditional probabilities P(3, 3|3, 3) and P(5, 1|5, 1) in a beamsplitter with transmittance η = 1/2 and η = 5/6 respectively, as a function of the relative delay dT /∆T between the arrival times of the phase-randomized WCPs at the two input ports of the beamsplitter. Here dT denotes the absolute delay between the arrival times of the optical pulses at each input port of the beamsplitter and ∆T is the full-widthhalf-maximum (FWHM) of the pulses, which for simplicity we assume is equal for all of them. In these simulations, the efficiency of the threshold single-photon detectors is set equal to 80% [45], and the dark count probability is p dark = 10 −6 . We have chosen these particular examples because quantum mechanics predicts that these probabilities are equal to zero (i.e., complete destructive interference) when dT /∆T = 0. As we can see from Fig. 2, our estimations approximate very well the theoretical value, and the simulated lower bounds for the visibilities V 3,3|3,3 and V 5,1|5,1 are very close to one. To be precise, we obtain V 3,3|3,3 ≥ 0.99994 and V 5,1|5,1 ≥ 0.99996. The reasons for the slightly noisy behaviour of the estimated values as well as for the small discrepancy between these and the theoretical values predicted by quantum mechanics (especially when dT /∆T = 0) are mainly twofold. First, as we have already mentioned above, in our simulations we use a relatively small number of decoy-state/detector-decoy settings. In particular, for each value of dT /∆T , we choose an optimized set of six possible values for the input parameters µ 1 and µ 2 and five possible values for the output parameters κ 1 and κ 2 . By using a larger number of settings one could in principle approximate the theoretical value as much as desired. The second reason is the limited numerical precision of the linear solver as well as the fact that, as explained in B, to solve Eq. (8) numerically we reduce the number of unknowns P (x 1 , x 2 |n 1 , n 2 ) to a final set. Also, we emphasise that the upper and lower bounds illustrated in Fig. 2 depend on the absolute value of dT /∆T . This is because the experimental data P µ,κ θ that we use in our simulations depend on |dT /∆T | (see A for further details).
Finally, let us remark that when we try to estimate the conditional probabilities P (x 1 , x 2 |n 1 , n 2 ) for higher total input photon numbers, the accuracy of the estimation decreases. This is so because the value of the coefficients P µ n P κ (θ|x) decreases very rapidly when n increases, which bound for these probabilities obtained with our method based on the use of two laser sources emitting phase-randomized WCPs and two threshold single-photon detectors. In our simulations we consider that the efficiency ηD of the detectors is 80% [45], the dark count probability is p dark = 10 −6 , and the transmittances γi (ωj) take six (five) different values.
renders the estimation problem difficult to solve numerically even with strong scaling methods. Moreover, increasing the value of the intensity setting µ is not of much help here, since it entails an increase of the leftover term (see B). Possible solutions might be to try to solve the set of linear equations analytically by means of say Gaussian elimination, or to develop more efficient numerical estimation methods. It would be definitively interesting to further investigate these two options.
The first scenario that we consider is shown in Fig 3(a). In this case, the input pulses to the tritter have the same polarization state, but their arrival times to the different input ports of the tritter vary. The result predicted by quantum theory for P (1, 1, 1|1, 1, 1) in this situation is shown with a solid line in Fig 3(a), while our estimations are shown with dots. Again, we can see that the estimated upper and lower bounds for P (1, 1, 1|1, 1, 1) fit tightly the theoretical probability. In the second scenario, now the polarization states of the input light pulses are chosen to compensate the temporal distinguishability between the arriving photons and they might be different for the signals at each input port. The motivation for this scenario is to observe the dependence that the three photon coincidence probability P (1, 1, 1|1, 1, 1) has on the triad phase φ by keeping constant all those terms r jk that affect such probability but arise from two-photon distinguishability [40]. The results are shown in Fig 3(b), where once again we can see that our method provides a tight estimation of the theoretical values, thus showing its practicality. Also, we remark that, as in the case of Fig. 2, the upper and lower bounds illustrated in Fig. 3 depend on |dT /∆T | because in our simulations the experimental data P µ,κ θ depend on |dT /∆T |. Moreover, like in the previous example of the beamsplitter, in our simulations we consider that the efficiency η D of the threshold single-photon detectors D j is 80% [45] and their dark count probability p dark = 10 −6 . Also, for the observables P µ,κ θ we use the expected values predicted by quantum mechanics. Furthermore, for each value of dT /∆T in Fig. 3(a), and for each value of φ in Fig. 3(b), we choose three different values for the intensities of the phaserandomized WCPs as well as two possible values for the output attenuators.

IV. CONCLUSION
In this paper we have proposed a simple method to experimentally characterize the behaviour of small-size optical networks (ONs) on input signals that are tensor products of Fock states. More precisely, our method could be used to obtain a tight estimation of the input-output photon number statistics of a ON. Importantly, our technique could be easily implemented with current technology like, for instance, phase-randomized weak coherent pulses together with threshold single-photon detectors. The main idea of the method is rather simple: it estimates the statistics provided by ideal n-photon sources at the input ports of a ON by means of decoy-state techniques and it estimates the statistics provided by ideal photon-number resolving detectors at its output ports by means of detector-decoy techniques. To illustrate the practicality of the method we have evaluated two simple examples.
In the first one, we have estimated the generalized Hong-Ou-Mandel dip in a beamsplitter for a total number of six input photons, while in the second example we have estimated the three coincidence detection probability in a tritter when it receives one singlephoton pulse in each of its input ports. In both cases we have obtained tight estimations that approximate very well the theoretical values.

V. ACKNOWLEDGMENTS
The authors wish to thank Daniel J. Gauthier, Hoi-Kwong Lo and Norbert Lütkenhaus for very useful discussions on the topic of this paper. This work was supported by the Galician Regional Government (consolidation of Research Units: AtlantTIC), the Spanish Ministry of Economy and Competitiveness (MINECO), the Fondo Europeo de Desarrollo Regional (FEDER) through grant TEC2014-54898-R, and the European Commission (Project QCALL). A.N. gratefully acknowledges support from a FPU scholarship from the Spanish Ministry of Education. F. X. acknowledges support from the 1000 Young Talents Program of China.
Appendix A: Toy model for the experimental data In order to evaluate the performance of our technique we need to generate the experimental data P µ,κ θ which is required to run the simulations. For this, and in the absence of a real experiment, we use a simple mathematical model that we detail below. In particular, let A † (B † ) be the creation operators for the input (output) spatial modes of the ON. That is, A † = [â † 1 , ...,â † N ] T and B † is defined similarly. These vectors satisfy where U is the unitary transformation that describes the behaviour of the ON.
In the case of WCPs, the input state to the ON can be written as |Ψ in = N k=1 |ψ in,k , where is the coherent state at the kth input mode [47]. Here, the parameters α k (ω) are defined as That is, for simplicity we assume that each |α k (ω)| 2 follows a Gaussian distribution of mean zero and standard deviation σ which is multiplied by the intensity µ k to guarantee that the condition |α k (ω)| 2 dω = µ k holds. The temporal parameter t k represents the arrival time of the optical pulse that enters the ON through its kth input port. We remark that in the definition of the states |ψ in,k we have not included yet the fact that their phases φ k are randomized. We will return to this point later.
Let {u jk } be the elements of the unitary matrix U . By applying Eq. (A1), and due to the linearity of the integral, we have that the state at the output ports of the ON can be written as Ψ out = M k=1 |ψ out,k , where and β k (ω) = N j=1 α j (ω)u jk . This means that the state |Ψ out at the output ports of the attenuators of efficiency κ is given by For convenience, note that here, like in the main text, we have included the effect of the efficiencies η D of the threshold single-photon detectors into the efficiency of the attenuators. The probability of having vacuum in a specific output mode k is related to the mean photon numbern k = |β k (ω)| 2 dω of the coherent state in that mode by P 0 = e −n k . In order to calculaten k , let ϕ jk be the phase of the element u jk of U , i.e., u jk = |u jk |e iϕ jk . Then, it is straightforward to show that This means, in particular, that where τ ij = t j − t i represents the delay between the arrival times of the pulses that enter the ON through its input ports ith and jth, and ∆T is their FWHM. Finally, we have that the joint probability of detecting a certain pattern θ on the threshold single-photon detectors is given by where φ = {φ 1 , ..., φ N } represents the dependence of that probability on the phase of each input coherent pulse. This is so because the probability of having no click at the output port k (that is, θ k = 0) is given by P µ,κ k ,φ 0 = (1 − p dark )e −κ knk , and thus the probability of having a click (θ k = 1) has the form P µ,κ k ,φ 1 = 1 − (1 − p dark )e −κ knk . If we consider now the fact that the input coherent states are phase-randomized, we obtain that the probability of detecting the pattern θ on the threshold single-photon detectors D j is given by which can be calculated numerically or even analytically for the simplest cases.

Appendix B: Numerical estimation with linear programming
For small values of the intensities µ = {µ 1 , ..., µ N } we have that the coefficients P µ n P κ (θ|x) of the set of linear equations given by Eq. (5) drop quickly to zero when the number of photons n ≡ n 1 , ..., n M increases. Therefore, one can neglect some of the terms in Eq. (5) to decrease the number of unknowns P (x|n) to a finite set. For instance, one can discard all the summation terms that satisfy N i n i > M cut , for a certain prefixed parameter M cut . In this way, we obtain that P µ,κ θ ≥ n∈Scut x≤n P µ n P (x|n)P κ (θ|x), where S cut is the subset that contains all possible n such that N i n i ≤ M cut . Similarly, one could also obtain an upper bound on P µ,κ θ that depends on the same finite number of unknowns P (x|n). For this, note that P µ,κ θ = n∈Scut x≤n P µ n P (x|n)P κ (θ|x) + n / ∈Scut x≤n P µ n P (x|n)P κ (θ|x) ≤ n∈Scut x≤n P µ n P (x|n)P κ (θ|x) + n / ∈Scut P µ n x≤n P (x|n) = n∈Scut x≤n P µ n P (x|n)P κ (θ|x) + 1 − n∈Scut P µ n = n∈Scut x≤n P µ n P (x|n)P κ (θ|x) + Λ µ Scut , where the first inequality is due to the fact that P κ (θ|x) ≤ 1 and the second equality comes from x≤n P (x|n) = 1 and n P µ n = 1, ∀n. Obviously, the leftover term Λ µ Scut = 1 − n∈Scut P µ n should be as small as possible. By using this result, one can numerically obtain an upper bound for the probability P (x|n) by solving the following linear program max P (x|n) s.t. P µ,κ θ ≤ n∈Scut x≤n P µ n P (x|n)P κ (θ|x) +Λ µ Scut , ∀µ, κ, θ P µ,κ θ ≥ n∈Scut x≤n P µ n P (x|n)P κ (θ|x), ∀µ, κ, θ 0 ≤ P (x|n) ≤ 1, ∀x ≤ n, n ∈ S cut x≤n P (x|n) = 1, ∀n ∈ S cut .
The lower bound can be estimated by simply replacing the max with a min.