Benchmarking single-photon sources from an auto-correlation measurement

Here we argue that the probability that a given source produces exactly a single photon is a natural quantity to benchmark single-photon sources as it certifies the absence of multi-photon components and quantifies the efficiency simultaneously. Moreover, this probability can be bounded simply from an auto-correlation measurement -- a balanced beam splitter and two photon detectors. Such a bound gives access to various non-classicality witnesses that can be used to certify and quantify Wigner-negativity, in addition to non-Gaussianity and P-negativity of the state produced by the source. We provide tools that can be used in practice to account for an imperfect beam splitter, non-identical and non-unit detection efficiencies, dark counts and other imperfections, to take finite statistical effects into account without assuming that identical states are produced in all rounds, and optionally allow one to remove the detector inefficiencies from the analysis. We demonstrate the use of the proposed benchmark, non-classicality witness and measure using a heralded single-photon source based on spontaneous parametric down-conversion. We report on an average probability that a single photon is produced $\geq 55\%$ and an average measure of the Wigner negativity $\geq 0.004$ with a confidence level of $1-10^{-10}$.


Introduction
Single-photon sources [1,2] are key resources for quantum communication [3], photonic quantum computation [4] or radiometry [5,6]. Not all single-photon sources are alike and to be scaled up, most applications require efficient sources of true single photons (single-photon Fock/number states). The quality of single-photon sources is usually quantified from an auto-correlation measurement [7], that is, by sending the photons to a balanced beam splitter and checking that the ratio between the detected twofold coinci- o X 6 6 S J K U 8 R N q e 5 N p 5 Z Z 8 P + 5 n 1 t P c 3 d h v Q P c 6 8 B s R p n x P 6 l + 8 z 8 r 8 7 U o t H D r q 1 B U E 2 J Z U x 1 L H f J b F f M z d 0 v V W l y S I g z u E t x R Z h Z 5 W e f X a t J b e 2 m t 4 G N v 9 p M w 5 o 9 y 3 M z v J l b 0 o D 9 n + M c B Y 2 t i u 9 V / M P t c n U v H / U 0 V r G G D Z r n D q o 4 Q A 1 1 8 r 7 A A x 7 x 5 B w 5 V 8 6 t c / e R 6 h R y z Q q + L e f + H T f Y l I 0 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " A u H h 8 o X 6 6 S J K U 8 R N q e 5 N p 5 Z Z 8 P + 5 n 1 t P c 3 d h v Q P c 6 8 B s R p n x P 6 l + 8 z 8 r 8 7 U o t H D r q 1 B U E 2 J Z U x 1 L H f J b F f M z d 0 v V W l y S I g z u E t x R Z h Z 5 W e f X a t J b e 2 m t 4 G N v 9 p M w 5 o 9 y 3 M z v J l b 0 o D 9 n + M c B Y 2 t i u 9 V / M P t c n U v H / U 0 V r G G D Z r n D q o 4 Q A 1 1 8 r 7 A A x 7 x 5 B w 5 V 8 6 t c / e R 6 h R y z Q q + L e f + H T f Y l I 0 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " A u H h 8 o X 6 6 S J K U 8 R N q e 5 N p 5 Z Z 8 P + 5 n 1 t P c 3 d h v Q P c 6 8 B s R p n x P 6 l + 8 z 8 r 8 7 U o t H D r q 1 B U E 2 J Z U x 1 L H f J b F f M z d 0 v V W l y S I g z u E t x R Z h Z 5 W e f X a t J b e 2 m t 4 G N v 9 p M w 5 o 9 y 3 M z v J l b 0 o D 9 n + M c B Y 2 t i u 9 V / M P t c n U v H / U 0 V r G G D Z r n D q o 4 Q A 1 1 8 r 7 A A x 7 x 5 B w 5 V 8 6 t c / e R 6 h R y z Q q + L e f + H T f Y l I 0 = < / l a t e x i t > ? Figure 1: Schematic representation of the measurement that is considered to characterize an unknown photon source which supposedly produces single photons. It is realized with a beam splitter and two non-photon-number-resolving detectors, as in a standard auto-correlation measurement. At each round, each detector either clicks • or not •. By analyzing the frequency of these events, the probability that the source actually produces exactly a single photon can be lower bounded. Furthermore, if a single radiation mode is detected, various forms of non-classicality can be witnessed and quantified.
dences and the product of singles vanishes, see Fig. 1 and the discussion in Appendix A. This ensures that the source produces no more than one photon. The result is, however, insensitive to loss as the efficiency cancels out of the ratio. These two aspects -the capacity of a source to produce no more than one photon and its efficiency -are thus considered separately. Both aspects are, however, important and are quantified jointly by the probability that the source actually produces exactly a single photon. Characterizing this probability is a direct and more complete way to benchmark single-photon sources.
Interestingly, this probability can be bounded by reconsidering the statistics of detector counts in an auto-correlation measurement. This suggests a systematic way to benchmark the quality and quantify the efficiency of single-photon sources, and to witness and quantify their quantum nature. To motivate benchmarking single-photon sources by the probability that the source actually produces exactly a single photon, we provide a detailed analysis which includes a simple statistical tool to account for finitesize effects without assuming that identical states were produced in all rounds of the experiment. We show how to include imperfections in the measurements apparatus and how to remove the detector ef-ficiencies from the analysis to facilitate the use of the proposed benchmark. An experimental demonstration is presented, illustrating the quality of heralded single-photon sources based on spontaneous parametric down-conversion.
The auto-correlation measurement is also known to be valuable for witnessing various forms of nonclassicality, including the non-positivity of the Pdistribution and quantum non-Gaussianity [8][9][10].
Here we show that for single-mode states or measurements a bound on the single-photon probability can be readily used to quantify the negativity of the Wigner representation [7], the strongest form of non-classicality, with respect to a measure proposed in [11]. We apply this method to verify Wignernegativity in the reported experiments.

Measurement apparatus
The measurement apparatus we consider is similar to the one used for the second order auto-correlation measurement. It is a simple measurement consisting of sending the photonic state to be measured (labeled ρ) on to a beam splitter and recording the photoncount correlations between the two outputs, see Fig. 1. We consider that photon detections are made with typical non-photon-number-resolving detectors. In order to draw conclusions from such photon-counts, we introduce a simple quantum model for such a measurement setup. A non-photon-number-resolving detector of efficiency η can be modeled with a two element positive operator-valued measure (POVM) {E • , E • } corresponding to click (•) and no-click (•) outcomes. When the measurement acts on a single mode characterized by bosonic operators a and a † , the POVM elements take the following form This model describes a non-photon-number-resolving detector, for which every incident photon can trigger the detection event with probability η. The no-click outcome then corresponds to the event where none of the incident photons triggered a click, and occurs with probability p • = tr ρ (1 − η) a † a , see e.g. [12]. From now on we assume that the detectors are accurately described by Eq. (1) for some value of η. In section 3.1.4 we will discuss how to account for small deviation from this model. When two such detectors are placed after a beam splitter with reflectance r (transmittance t), it is straightforward to see that the four possible outcomes are given by the POVM elements where the first (second) label •/• refers to the detector after the reflected (transmitted) output of the beam splitter, see Appendix B for a formal derivation. The events where a fixed detector does not click are modeled by the two POVM elements The corresponding probabilities are labeled p • and p • . Note that the case where the two detectors do not have the same efficiency η R = η T can be accounted for by replacing t with t = tη T tη T +rη R , r with r = rη R tη T +rη R , and setting η = tη T + rη R in Eq. (2).

Benchmarking a single-photon source
For simplicity, through the main part of this section we will be considering the case of single-mode sources. In section 3.3, we show that all the presented tools also apply to multi-mode sources.
With the measurement we just described, any single-mode state ρ incident on the beam splitter can be associated with a probability vector governing the occurrence of clicks. Our goal is to construct an estimatorP 1 (p) that relates this vector p to the photon number statistics of the state and in particular to the weight of the single-photon component P 1 = 1| ρ |1 . Directly lower bounding P 1 is a natural way to benchmark a single-photon source. In particular, it sets a bound on the trace distance between the state ρ prepared by the source and an ideal single photon |1 which can be readily used to bound the errors for various applications of single-photon sources.

Measurement calibration independent benchmark
We proceed step by step, considering first an ideal measurement apparatus, then considering two identical non-unit detector efficiencies and finally focusing on the most imperfect measurement, where we consider an unbalanced beam splitter and two detectors having different efficiencies. At this point we only assume that the measurement is described by the POVM in Eq. (2) for some parameters η, t and r, that do not need to be known. The benchmark that we will derive in this section is thus independent of the calibration of the measurement apparatus. Finally, we show that by adding a small correction term our benchmark can be applied in situations where the POVM in Eq. (2) is only an approximate description of the measurement apparatus.

Ideal measurement apparatus
For clarity, we first assume that the source produces an identical state ρ at each round, the detectors have unit detection efficiency η = 1 and the beam splitter is balanced t = r = 1/2. In this ideal case, the probabilities p •• , p •• are equal, and p is described by two independent real parameters. For convenience, we intro- First, we note that p • ≥ p •• holds by definition. Furthermore, the points (1, 1) and (0, 0) are attained by the vacuum and the state with infinitely many photons, respectively. Thus, the line p • = p •• is also attainable by mixtures of aforementioned states. Then, we look for the maximum value of p • = n P n 1 2 n for a fixed p •• . We have to solve As (1/2) n is decreasing with n, the maximum is attained by saturating the values of P n starting with P 0 . Hence, it equals The set of possible values (p • , p •• ) is thus included in a convex polytope with four vertices Q P = Polytope{(0, 0), 1+P 4 , 0 , 2−P 2 , 1 − P , (1, 1)}, sketched in Fig. 2. The only nontrivial facet of this polytope is the edge connecting 1+P 4 , 0 and 2−P 2 , 1 − P which is associated to the inequality 4p • − 3p •• − 1 ≤ P , and is given by the colored lines in Fig. 2. Thus, without loss of generality, the condition 1| ρ |1 ≤ P implies that the elements of p satisfy the linear constraint Conversely, by measuring the pair (p • , p •• ) and by computing the resulting value ofP T 1 , we can guarantee that for any value of P such thatP T 1 > P , 1| ρ |1 > P holds, that is we get a lower bound on the probability that the source to be benchmarked produces exactly a single photon.

Identical non-unit efficiency detectors
To move away from the ideal case, we still consider a perfectly balanced beam splitter and focus on a situation where non-unit efficiency detectors are used. We consider the case where the detector efficiency η is unknown. In the measurement setup shown in Fig. 1, non-unit efficiency detectors can be modeled by taking ideal detectors and placing a beam splitter with transmission η before the balanced beam splitter. As a consequence, observing a violation of Ineq. (8) proves that the state produced by the single-photon source and undergoing losses satisfies 1| |1 ≥P T 1 (p). This provides a valid benchmark even though the intrinsic quality of the source is estimated with a lossy measurement apparatus. It is interesting to note that for any state ρ with P 1 ≥ 2/3, the probability of the single-photon weight P 1 can only decrease with loss, see Appendix C. Therefore, showing that P 1 ≥ 2/3 with lossy detectors implies that the original state also satisfies P 1 ≥ 2 3 . This is not the case when P 1 < 2/3, i.e. for specific states the single-photon weight P 1 can be increased by loss (an intuitive example is the twophoton Fock state).

Unbalanced beam splitter and different non-unit efficiency detectors
We now relax the assumptions that the beam splitter is balanced and the detector efficiencies are the same, i.e. we consider the case with a measurement performed with a beam splitter having an unknown transmission t and reflection r = 1 − t and two detectors having different efficiencies labeled η R and η T . In this case, the observed statistics would be equivalently obtained with a beam splitter having a transmission coefficient t = tη T /(tη T + rη R ) and two detectors with the same efficiency η = tη T + rη R , as already mentioned below Eq. (2). This means that the measurement can be modeled with a first unbalanced beam splitter with transmission coefficient η corresponding to loss on the state to be characterized, an unbalanced beam splitter with transmission t and two detectors with unit detection efficiency. In this case, the relation between p •• and P 0 is unchanged, are however no longer the same. They are now given by p • = n P n (1 − r ) n and p • = n P n (1−t ) n . Hence, the quantity n P n 1 2 n is no longer directly related to the probability p. Nevertheless, it can be bounded from observable quantities, as n P n Using the definition of p • and p • , we rewrite them in terms of probabilities of disjoint events aŝ With this notation in hand, we conclude that the quantityP is a benchmark for single-photon sources, without assumptions on the detector efficiencies and on the fact that the beam splitter is balanced. This means that from the outcome probabilities p of a usual autocorrelation measurement, we can computeP T 1 (p) and P R 1 (p), deduce their minimumP 1 (p) and guarantee the tested source produces states with the weight of the single-photon component satisfying P 1 = 1| ρ |1 ≥P 1 (p).

Dark counts and other imperfections
Finally let us briefly consider general passive detectors described by a POVM (11) where e • (n)+e • (n) = 1. Consider combining two such detectors in the g (2) setup, with the POVM {Ẽ •/• } performed on the transmitted mode and {Ẽ •/• } on the reflected one. The resulting POVM elements for with P BS (k|n) = n k t k r n−k . Now let us assume these detectors are not too different from the textbook single-photon detector model of Eq. (1). Concretely we assume that for some η, t and r where ρ can be any state susceptible to be prepared in the experiment. In particular, consider the case where both detectors satisfy for some η, η , and all n on which the input state is supported.
It follows that the probabilitiesp With the help of Eq. (9) it is then straightforward to lower bound the value of the benchmark given the coincidence probabilitiesp observed with any detectors abiding to Eq. (13).
Here, it is worth noting that the POVM model of Eq. (1) does not account for dark counts. These can be modeled by setting a nonzero click probability for the vacuum e • (0) = p dc , implying |e • (n) − (1 − η) n | ≤ δ = p dc . Then accordingly to Eqs. (15,17), to account for the effect of dark counts one can use the bound where we only took the leading order correction in p dc , which is typically fairly small.

Measurement calibration dependent benchmark
The benchmark we proposed relies on no assumption on the characteristics of the beam splitter or the two detectors used in the auto-correlation measurement. This prevents a miscalibration of the measurement apparatus that could result in an overestimation of the quality of the single-photon source. Nevertheless, as one would expect, the bound on the quality of the tested source (see Eq. (10)) is reduced when using an unbalanced beam splitter or inefficient detectors. We now discuss a way to characterize the intrinsic quality of the source by making additional assumptions on the measurement apparatus. The basic idea is to exploit estimations of different losses in photonic experiments. In particular, we assume that the detector efficiencies η R(T ) and the beam splitter reflectivity r are bounded, that is η R(T ) ≤η R(T ) and r ∈ [1 −t,r]. We show in Appendix D that the condition 1| ρ |1 ≤ P implies that the two following inequalities hold where the coefficients C 1 and C 2 , defined as have been optimized such thatP T * 1 (p) andP R * 1 (p) give the tightest bound on P . Finally, one can choose the best among the two bounds, giving rise to the benchmark for the single-photon probability, which takes advantage of the additional experimental knowledge.
To account for dark counts one may use ∆ = 2p dc .
is the vector of probabilities observed with the real detectors, while p gives the probabilities that would have been observed with the textbook model and that are used in the benchmark.

Multi-mode sources
We have so far considered sources emitting light in a single mode, or equivalently that the emitted light is filtered in all the auxiliary degrees of freedom so that a single mode of light is detected. We will now briefly consider the situation where the detected state is multi-mode, each mode being associated to an annihilation operator a k satisfying [a k , a † ] = δ k . The no-click and click events for a multi-mode input are associated to the POVM elements does not click only if none of the modes triggers a click. To a multi-mode state ρ one associates the distribution P n of the total photon number operatorn = k a † k a k . Assuming that the beam splitter acts identically on all modes, the measurement apparatus is only sensitive to the total number of photonsn. That is, the POVM ele- albeit with a † a replaced byn, see Appendix B. We thus conclude that the quantitiesP 1 (p) andP * 1 (p) that have been derived, can be readily used to benchmark the probability that a multi-mode source emits a single photon It is worth noting that a high P 1 for a multi-mode source does not guarantee that the single-photon probability is high in any of the individual modes. As a practical example, consider the multi-mode singlephoton state This state has exactly one photon in totaln = 1. Nevertheless, it is not a good single-mode single-photon state. The probability to find exactly one photon in any single mode is only 1/N , therefore two such states would only exhibit a limited two-photon interference (bunching). In other words is not a pure single-photon state |1 1| 1 , the degree of freedom that distinguishes between the different modes is in a highly mixed states. Depending on the application, the mode-purity (or single-mode character) of the source may be either irrelevant (e.g. for some quantum random number generators) if the interference between different sources plays no role, or crucial (e.g. quantum repeaters, boson sampling or photonic quantum computation) if it is at the heart of the task.
In Appendix G, we show that one can guarantee that the source produces a single-mode state with high P 1 if it is reasonable to assume that the multimode state is a product = k ρ k . In general, however, the auto-correlation measurement is intrinsically insensitive to the multi-mode characteristics of the source, as illustrated by the example of Eq. (24). In principle, it is possible to extract a single-mode state from a multi-mode source by filtering the auxiliary degrees of freedom, e.g. using a single-mode fiber for spatial degrees of freedom and a spectral filter for frequency domain. In any case, the mode-purity of the detected state has to be verified with a different set of measurements, e.g. via Hong-Ou-Mandel interference [13] or based on a physical model of the source. For example, in the case of a heralded single-photon source, the spectral purity can be determined by measuring the signal-idler joint spectral intensity [14], as will be explained in Sec. 6.3.
A systematic analysis of the characterization of the mode-purity of a source is beyond the scope of this paper. Nevertheless, it is worth mentioning that with the additional information on the multi-mode characteristics of the source and the knowledge of the total single-photon probability P 1 it is possible to bound the probability to find a single photon in the mode of interest P [1] 1 (denoted mode 1). In particular, given a bound on the probability that all the modes but one are empty tr ( since the maximal contribution to P 1 from the other modes is ε. Here, ρ 1 = tr 2,...,N is the marginal state of the first mode, which can in principle be filtered from the source. A bound on P [1] 1 can then be used to, e.g., quantify the Wigner-negativity of the state ρ 1 , see Sec. 5.

Finite statistics
In this section, we analyze finite size effects for the benchmark.
We sketch an analysis to account for finite statistics in any experiment aiming to evaluateP 1 (p), the benchmark for single-photon sources derived in Sec. 3. For a measurement round described by p we associate a random variable X T that takes different real values depending on the measurement result This random variable satisfies E(X T ) =P T 1 (p) in Eq. (9) (here and further E denotes the expected value of a random variable). Analogously, we define X R by exchanging the role of the two detectors, such that it satisfies E(X R ) =P R 1 (p). In general, the source may prepare a different state ρ (i) at each round, corresponding to different probabilities p (i) of the measurement outcomes. This means that in each round, we sample different random variables X (i) T and X (i) R , which are independent between rounds given the sequence of states ρ (1) , . . . , ρ (n) produced in the experiment. In this case, a reasonable figure of merit is the average quality of the state prepared by the sourceP 1 is the probability of the single-photon component of the state ρ (i) . Because the transmission and reflection coefficients of the beam splitter can be considered to be constant, either E(X holds for all i. This means that the average single-photon weight fulfills . Finally, we use the Hoeffding 1963 theorem [15] to show that is a one-sided confidence interval onP 1 with confidence α (see Appendix F). Precisely, with probability 1 − α the observed value ofq α lower bounds P 1 . It might be convenient to note that the quantity min{X T ,X R } of this confidence interval can be computed using Analogously, one can derive a one-sided confidence intervals associated to the calibration-dependent benchmark derived in Sec. 3.2. With very similar arguments one can show that botĥ where the functions C 1 and C 2 defined in Eq. (20) are confidence intervals forP 1 . Since they are derived from the same data, for a statistically meaningful statement one has to chose parameter α before computing the confidence interval. On the other hand, it is straightforward to see that is also a confidence interval at confidence level 1 − 2α. See App. F for a detailed derivation.

Relation to the non-classicality of the source
The data obtained from auto-correlation type measurements are known to be valuable for witnessing and quantifying various forms of non-classicality, including the non-positivity of the P-function and quantum non-Gaussianity [8][9][10]. We now show that the knowledge of P 1 in a single bosonic mode (as e.g. provided by our benchmarksP 1 (p) andP * 1 (p)) can reveal Wigner-negativity [16], arguably the strongest form of non-classically for a bosonic mode. In particular, Wigner-negativity implies the non-positivity of the P-function [17]. Similarly, it implies that the corresponding state is non-Gaussian, as a Gaussian state has a Gaussian (and thus positive) Wigner function 1 . Thus, demonstrating Wigner-negativity for a light source brings evidence of its strong quantum nature. Note that it has been shown recently that witnesses of Wigner-negativity can be derived systematically using a hierarchy of semidefinite programs [19]. Our contribution is more specific and aims at witnessing Wigner-negativity simply and directly fromP 1 (p) orP * 1 (p) observed on a single mode.

Wigner-negativity witness
The Wigner function is a representation of a singlemode state ρ in terms of the following quasiprobability distribution [20] with dβ 2 W ρ (β) = 1. Here, D β = e a † β−aβ * is the displacement operator with a complex amplitude β. Applying Eq. (32) to a Fock state gives [17] W |n n| (β) = 2(−1) n π e −2|β| 2 L n 4|β| 2 where L n is the Laguerre polynomial. Note that the following bound on the Laguerre polynomials e −x/2 |L n (x)| ≤ 1, see e.g. Eq. (18.14.8) in [21], leads to a bound on the Wigner function of Fock states |W |n n| (β)| ≤ 2 π . Note also that L 1 (x) = 1 − x. With the help of Eq. (33), the upper bound on the Wigner function of Fock states and the definition of the Laguerre polynomial L 1 (x), it is easy to see that the Wigner function of any mixture of Fock states 1 In addition, Hudson's theorem [18] tells us that any pure state with a positive Wigner function is Gaussian.
Focusing on the origin β = 0, we get W ρ (0) ≤ 2 1−2P1 π which is negative if P 1 is larger than 1 2 . Hence, if one concludes from the measurement of p thatP 1 (p) > 1 2 , one can conclude that the measured state is Wignernegative (recall that we assumed that the state ρ is single-mode).

Wigner-negativity measure
A natural way to quantify the negativity of the Wigner representation of a given state ρ is to measure the total quasi-probability for which the function W ρ (β) takes negative values [11], i.e.
which is manifestly zero for states with a positive Wigner function. In the Appendix E we show that N W (ρ) is non-increasing under Gaussian operations, which justifies its use as a measure of Wignernegativity. Note that with the help of Ineq. (34), we show that N W (ρ) satisfies where w 0 is the principal branch of the Lambert W function. The function F (P 1 ) is non-decreasing. Hence, from the measurement of p, we get a lower boundP 1 (p) on P 1 that can be used to lower bound N W (ρ) using F (P 1 (p)). The bound (36) is tight by construction in the ideal case N W (|1 ) = F (1) = 9 4 √ e − 1 ≈ 0.36. By computing F (P 1 ) ≥ 0 we show that the function F (P 1 ) in Eq. (36) is convex. This property will be used in the following section, where we discuss the finite statistics effects.

p-value to witness Wigner-negativity
We now present the extension of finite statistics analysis presented in Sec. 4 to the task of Wigner-negativity detection and quantification.
First, let us now consider the witness of Wignernegativity discussed in Sec. 5.1, that is W ρ (0) ≥ 0 =⇒P 1 (p) ≤ 1/2, and quantify the statistical significance of its contrapositive given the measurement data. This can be done by computing the p-value associated to the hypothesis that the Wigner function of the state is positive. As before, we consider the general case where a different state ρ (i) may be prepared at each run. Nevertheless, at each round the bound holds. Therefore, for the sequence of states prepared in the experiment, the average Wigner function at the origin satisfies and is negative ifP 1 > 1 2 . Given some value of min{X T ,X R } recorded after n measurement rounds, we show in Appendix F that for any collection of n states withW (0) ≥ 0, the probability that the results are equal or exceed the observed value of min{X T ,X R } is given by for min{X T ,X R } > 1 2 . In other words, given the observed value of min{X T ,X R }, the probability that it is coming from states that are Wigner-positive on average is bounded by the right-hand side of Ineq. (38). In App. F one finds a bound on the p-value for the calibration-dependent setting.

Confidence interval on the measure of Wigner-negativity
Finally, the convexity of the function F (P 1 ) in Eq. (36) implies that the average Wigner-negativity is a one-sided confidence interval onN W , that is, with probability 1 − α, the average Wigner-negativity as quantified byN W is lower bounded by nw α = F (q α ). This can be used both in the calibrationindependentnw α = F (q α ) and calibration-dependent nw * α = F (q * α ) settings.

The crucial role of the single-mode hypothesis
It is important to emphasize that the single-mode hypothesis is crucial in order to relate P 1 (or its esti-mated valueP 1 (p),P * 1 (p)) to Wigner-negativity. In particular, the multi-mode single-photon state of Eq. (24) becomes Wigner-positive for N > 2. In section 3.3 we have discussed how the single-mode character of the emitted radiation can be verified in practice. Here we merely recall that a bound of the form P [1] 1 ≥ P 1 − ε, where P [1] 1 is the single photon probability for a given mode, can be readily used to verify Wigner-negativity. One simply has nw [1] α if min{X T ,X R } > 1 2 + ε, for Wigner-negativity of the said mode.
It is worth mentioning that a possibility to ensure that the detected radiation is single-mode is to perform homodyne measurements. In such a measurement the incoming beam is mixed with a strong local oscillator on a beam splitter, the intensity of the output beams are then measured with linear detectors and subtracted. Under the assumption that the local oscillator is single-mode 3 the obtained signal is only sensitive to the single input mode identified by the local oscillator. Photon number statistics, and P 1 in particular, can be reconstructed from the statistics of a phase-averaged quadrature measurement [22], i.e. a homodyne measurement with a phase-randomized local oscillator. This offers the possibility to use our bound on the Wigner-negativity N W (ρ) in Eq. (36) with measurements that are guaranteed to pick up a single mode.

Experiment
To demonstrate the feasibility of our tools, we experimentally benchmark, witness and quantify the nonclassical nature of a heralded single-photon source [23] that is optimized for high efficiency of the heralded photon [24]. A periodically poled potassium titanyl phosphate (PPKTP) crystal is pumped by a Ti:Sapphire laser at λ p = 771.8 nm in the picosecond pulsed regime with a repetition rate of 76 MHz to create nondegenerate photon pairs at λ s = 1541.3 nm (signal) and λ i = 1546.1 nm (idler) via type-II spontaneous parametric down-conversion (SPDC). The pair creation probability per pump pulse is set to P pair ≈ 1.0 × 10 −3 and high-purity heralded signal photons are ensured by spectrally filtering the heralding idler photons using a dense wavelength division multiplexer at ITU channel 39. From a joint spectral intensity measurement [14], we estimate the spectral purity of the heralded photon to be 98.59% ± 0.04%. In this way, we herald signal photons at a rate of 19.1 kcps.
For the heralded auto-correlation measurement, the signal photon is sent to a 50/50 fiber coupler (AFW FOBC). All photons are detected by MoSi superconducting nanowire single-photon detectors [25] and time-correlated single-photon counting in a programmable time-to-digital converter (ID Quantique ID900) is used to register the detection events. Data are acquired for 200 s in order to evaluate p = (p •• , p •• , p •• , p •• ) for the signal photons after the 50/50 beam splitter.
The dark count probabilities of the detectors we used for the auto-correlation measurement are fairly small, p dc ≤ 4 × 10 −7 , and can be completely neglected. That is, the dark count corrections (≈ 4 × 10 −6 ) to the estimated values of the benchmark in Eqs. (18,22) are more than two orders of magnitude lower than the statistical noise, see Tables 1 and 3.

Calibration-independent benchmark
In a first step we apply our benchmark to the experimental results without taking the splitting ratio of the beam splitter and the detector efficiencies into account. The overall efficiency is 25 % for the heralding idler photons and 62 % for the heralded signal photons. In order to simulate a less efficient singlephoton source, we introduce loss by inserting a fiber coupled variable attenuator (JDS Uniphase MV47W) into the heralded photon path before the 50/50 beam splitter and repeat the auto-correlation measurement for eight different transmission efficiencies η att . Each transmission efficiency leads to a value for the pair (min(p • , p • ), p •• ) that is represented by a black cross in Fig. 2. In the same figure, we represent the polytope Q P (defined after Eq. (7)) containing all possible values (min(p • , p • ), p •• ) associated to states with 1| |1 ≤ P . Four polytopes are represented corresponding to the values P ∈ {0.25, 0.5, 0.75, 1}. A measurement result associated to a black cross lying outside a polytope Q P is guaranteed to come from a state with a single-photon component satisfying P 1 > P .  For the measurements with the three highest transmission efficiencies, we give the results of our benchmark in Tab. 1. We conclude for the highest transmission for example, that the measured states have on average a single-photon component with a weight P 1 ≥ 0.554 with a confidence level of 1 − 10 −10 .

Calibration-dependent benchmark
To compute the value of the calibration-dependent benchmark one needs to estimate the detector efficiencies and the reflection/transmission coefficient of the beam splitter. In order to characterize the detectors, we use the standard method, see e.g. [25] for a detailed description. For our setup we find that the beam splitter coefficients are bounded bŷ r ∈ [0.49, 0.50] and the detector efficiencies are upper bounded by (η R ,η T ) = (0.95, 0.88). The upper bounds for the detector efficiencies are obtained from results of the measured detection efficiencies given in Tab. 2 by adding three times the measurement uncertainty of around 0.01, see Supplementary Material of [25].
Under the assumptions that the (r,η R ,η T ) belong to these intervals, the values ofP T * 1 (p) andP R * 1 (p) as measured in our experiment are given in Tab. 3 for the three highest transmission efficiencies. The confidence intervalq * α ≤P 1 is also reported for a confidence level of 1 − α = 1− 10 −10 .   (19). The confidence intervalsq * α ≤P1 in the finite statistics analysis are calculated for a confidence level of 1 − α = 1 − 10 −10 .

Non-classicality of the source
As already mentioned, there is no guarantee that the single-mode assumption is exactly satisfied in our experiment, hence we have to estimate the mode purity of the source. We assume that the spatial mode purity is guaranteed by coupling the photons into a single-mode fiber. The polarization-purity is ensured by the fact that the signal and idler photons are separated with a polarizing beam splitter. To estimate the spectral purity, we apply the standard approach that relies on a physical model of the source, which we believe to properly describe the experiment. Precisely, we assume that the SPDC process responsible for the generation of the photon pairs is of the form where a(w s ) and b(w i ) are the frequency field-modes of the signal/idler photons with [a(w), a † (w )] = [b(w), b † (w )] = δ(w − w ). At low pumping power, we reconstruct f (w s , w i ) with a signal-idler joint spectral intensity measurement. Via a 2D-Gaussian fit we perform a singular value decomposition of f (w s , w i ) to rewrite the interaction in the form where [a k , a † ] = [b k , b † ] = δ k now describe discrete spectral modes. With this procedure the largest Schmidt coefficient λ 1 is computed to be λ 1 = 0.992 92 (18), where the standard deviation σ λ1 = 1.8 × 10 −4 is obtained from a Monte Carlo method assuming Poissonian count statistics in the joint spectral intensity measurement. Therefore, assuming that the efficiency of the trigger detector is the same for all idler modes, we obtain the leading order estimate of the probability that the modes a k≥2 are empty conditional to the detection of an idler photon. This corresponds to the mode-purity of k λ 2 k ≈ 98.59%. With the help of Eq. (25) we can take this into account for the quantification of Wigner-negativity, resulting innw α=10 −10 = 0.0046 for the case of no added loss on the heralded single-photon state and no assumptions on the calibration of the measurement apparatus. The corresponding p-value and the results for the measurement-apparatus-dependent case are given in Tab. 4.

Conclusion
Auto-correlation measurements are commonly used to assess the quality and the quantum nature of singlephoton sources, that is, they are used to check that a given source does not emit more than one photon η s,totn w α=10 −10 p-valuenw * α=10 −10 p-value * 62 % 0.0046 10 −603 0.053 10 −6420 52 % 0 × 0.0072 10 −894 42 % 0 × 0 × Table 4: Wigner-negativity in our experiment for the three highest transmission efficiencies ηs,tot of the heralded singlephoton state. The confidence interval on the measure of Wigner-negativitynwα ≤NW is obtained from Eq. (39) for the confidence level α = 10 −10 , assuming the reduced singlemode P [1] 1 as given in Eq. (25). Further, we give the p-value according to Eq. (38) associated with the hypothesis that the measured states are on average Wigner-positive. The quantities with a * are taking the detector efficiencies into account and are obtained accordingly from Eqs. (19) and (25). and its emission is non-classical in the sense that its P-distribution is non-positive or that its state is non-Gaussian. We have shown that the statistics obtained from these measurements is actually richer. They can be used to lower bound the probability that a given source actually produces a single photon. We argued that this probability is a good benchmark for single-photon sources as it captures both its quality and its efficiency. Moreover, we showed that if the mode purity of the source can be assessed the lower bound on the single-photon emission probability can be used to witness and quantify the negativity of the Wigner function, a stronger form of nonclassicality than the negativity of the P-distribution and the non-Gaussianity. We have proposed practical tools to benchmark single-photon sources and characterize its Wigner-negativity this way. With this material in hand, we hope that the community which is developing single-photon sources could exploit the statistics of their auto-correlation measurements in a more enlightening way. 000910-S), Fundació Cellex, Fundació Mir-Puig, Generalitat de Catalunya (CERCA, AGAUR SGR 1381) and from the ERC AdGCERQUT.

Data availability
The data supporting the experimental results within this paper are available on the Zenodo data repository [26].

A Auto-correlation function and autocorrelation measurement
Here we briefly recall the definitions of the autocorrelation function g (2) . We will use the notation introduced in Section 2, 3 and 3.1.1 of the main text, and assume that the beam splitter is balanced and the two detectors have equal efficiency, i.e. t = r = 1 2 in Eq. (2).
Historically [27], the auto-correlation function was defined as the ratio and the efficiency of the source can be characterized by the average number of photons it emits I = a † a .
Then it is not difficult to see that the bound is a tight benchmark. To see this note that in the a † a , (a † a) 2 plane the quantity w = 2 a † a − (a † a) 2 measures the distance from the line connecting the points (0, 0) and (2, 4) corresponding to Fock states |0 and |2 . With w = 1 for the point (1, 1) corresponding to the single-photon state |1 . In practice, one can not directly measure (a † a) 2 and a † a , but can approximate these values by increasing the loss artificially, as . Such an approach thus requires a precise control of the efficiency η and is statistically inefficient, since additional losses are introduced.
Alternatively [28], the auto-correlation function can by directly defined as The efficiency can also be characterized by probability that a source produces a clickĨ = E Given the two valuesg (2) andĨ it is then possible to reconstruct the full probability distribution so that p is defined by two parameters. As argued above the two functions coincide g (2) =g (2) in the limit η → 0.
In both cases the auto-correlation measurements relies on the setup of the Fig. 1. The measurement data can thus be readily used to estimate p and compute the benchmarks proposed in this paper.
To finish the discussion of the auto-correlation functions we recall that both g (2) andg (2) are witnesses of the non-classicality of the state ρ [28], i.e. g (2) ,g (2) < 1 is only possible for states whose P-function admits negative values. In fact, this is true in a more general context, as given by the following observation.
Observation. For any two binary POVMs {M • , M • } and {M • , M • } measured at the two outputs of a beam splitter (Fig. 1), the inequality is a witness of non-classicality (P-function  For the sake of completeness we prove the above statement here. A coherent state splits into two coherent states on a beam splitter |α → BS | √ rα √ tα . Hence for a coherent state and G (2) = 1. Then for any mixture of coherent states ρ cl = d 2 αP(α) |α α| one finds (49) To shorten the equations let us denote z = |α|, f (z) = p • ( r|α| 2 ), g(z) = p • ( t|α| 2 ), µ(z) = 2z dϕP ze iϕ and dµ(z) = dzµ(z) with dµ(z) = 1, such that Without loss of generality consider z ≥ z, by p • (|n ) ≥ p • (|m ) for n ≥ m it follows that f (z ) ≥ f (z) for z ≥ z (using the fact that a coherent state has a Poissonian photon number distribution). Hence, one can write and (f (z)g(z) + f (z )g(z )) − (f (z)g(z ) + f (z )g(z)) for each z and z . Therefore the inequality is also true for the integrals Showing that G (2) ≥ 1 for any state with a positive P-function. Notably, the same proof works if both n| M • |n and n| M • |n are decreasing functions of n, e.g. if they are replaced with M • and M • . The above observation has been used to propose an experiment where non-classicality of light would be demonstrated by directly using human eyes as detectors [29].
B POVM corresponding to photon detection preceded by a beam splitter Since two single-photon detectors after a beam splitter cannot detect any coherence between Fock states of an incoming state, we can consider a state ρ = n cn n! (a † ) n |0 0| a n without loss of generality. The state ρ arrives at a beam splitter with transmittance t, the resulting state is (ignoring again the coherence between the Fock states) One can then compute the probability of the different events for the state ρ r , for example where and The resulting POVM is thus E Let us consider a two-mode state (generalization to n modes is straightforward) ρ = n1,n2 The POVM element corresponding to no-clicks after the beam splitter on each mode now reads E (2) • = E • ⊗ E • . The state after the beam splitter is now We focus on the event click on the first detector and no click on the second We retrieve that the POVM elements in the multimode case are given by C The effect of loss on P 1 We show here that the set of states ρ with P 1 ≥ 2 3 is closed under losses. Consider a state ρ associated with a single-photon component P 1 = P . Let us apply infinitesimal transmission losses η = 1 − d . After the loss, the photon number distribution P n = n|ρ|n reads In particular, which is negative for P ≥ 2/3. On the other hand, there are states with P 1 < 2/3 for which the single-photon probability can be increased substantially by losses. Consider a channel with transmission efficiency η and apply it to the state The blue line depicts maxη P1(ρη) as the function P1, and the black line is simply P1. One sees that for this example the single-photon component can be increased by losses if P1 < 2/3, as illustrated by the shaded area. Conversely, we show that for any single-mode state ρ with P1 ≥ 2/3, the single-photon component can not be increased by losses. ρ = P 1 |1 1| + (1 − P 1 ) |2 2|. This leads to a state ρ η having a single-photon component This quantity is maximized at the maximum being depicted in Fig. 3 as a function of P 1 . To give a concrete example, for the initial ρ = 1 2 (|1 1| + |2 2|), the weight of the singlephoton component can be increased to P 1 (ρ η=3/4 ) = 1/2 + 1/16 = 0.562 while for the Fock state ρ = |2 2|, it is possible to reach P 1 (ρ η=1/2 ) = 1/2.

D Parameter dependent witness
We consider the case where the two detectors have efficiencies η T and η R and the beam splitter has reflectance r and transmittance t (with t + r = 1). For an incoming Fock state |n the probabilities of clicks are given by For a mixture of Fock states ρ = n P n |n n| one has from which we get In order to set a lower bound on P 1 we thus need to upper bound the term n≥2 P n f n . We also have Therefore, to derive a benchmark for P 1 we are looking for a function f * (p •• ) such that To find f * (p •• ) we define q n = P n g n , so that the maximization can be rewritten as Using n q n fn gn ≤ ( n q n ) max n fn gn we see that the solution of Eq. (67) satisfies We prove right after that the maximum is achieved for n = 2. By plugging n≥2 P n f n ≤ f * (p in Eq. (65), we get the desired inequality The same bound holds with p •• instead of p •• and h n instead of f n . The right-hand side of these two inequalities are a function of the observed probabilities p and a lower bound on P 1 , hence defining two benchmarksP The best option is to consider the larger value of p •• and p •• . The proof of max n fn gn = f2 g2 . To maximize the ratio fn gn express it as Manifestly, maximizing fn gn is equivalent to minimizing In other words we want to show that for n ≥ 2 the fraction with x = η T t and y = 1 − η R r satisfying 0 < x < y < 1, is minimized at n = 2. However, it is enough to show that the expression 1 n −(1−x) n y n −(y−x) n is increasing with n. To do so, let us derive this quantity with respect to n. We have d dn To show that it is positive we can omit the denominator (y n − (x − y) n ) 2 . Labeling a = (1 − x) n , b = y n , c = (y − x) n and noting that log For b > a we set c = a and get To show that the last expression is positive, we divide it by and note that the function (80) Hence, the fraction 1−(1−x) n y n −(y−x) n is increasing with n and attains its minimum at the boundary n = 2 of the interval [2, ∞). Therefore, fn gn is maximized at n = 2, which concludes the proof.

E Wigner-negativity measure
For a single-mode state ρ, the Wigner function W ρ (β) is a quasi-probability distribution satisfying dβ 2 W ρ (β) = 1. The negativity of the Wigner function (W (β) < 0 for some β ∈ C) is an important nonclassical feature of the state, as argued in the main text. A natural way to quantify this negativity is to measure the total quasi-probability where the function W ρ takes negative values, that is to compute This intuitive quantity was introduced in [11]. We now show that N W (ρ) is a good "measure" of Wigner negativity in the sense that it can not be increased by Gaussian operations. Pure Gaussian operations are displacements D γ = e γa † −γ * a , single-mode squeezing SMS g = e g 2 (a †2 −a 2 ) , phase rotations e iϕa † a , or combination thereof. Consider a single-mode state ρ with its Wigner function W ρ (β) and its Wigner negativity measure N W (ρ). The effect of a displacement = D γ ρD † γ on the Wigner function is a mere translation in phase space W (β) = W ρ (β − γ), which does not affect the Wigner negativity measure N W (ρ) = N W ( ). The same goes for a phase rotation, which merely transform W (β) = W ρ (βe iϕ ). For a squeezing operation = SMS g ρSMS † g , the Wigner function is transformed as where β = β + iβ andβ = e g β + e −g iβ . This implies for the Wigner negativity measure that Hence, N W (ρ) is also unchanged by squeezing. Moreover, the quantity is manifestly convex as Hence, N W (ρ) is nonincreasing under mixtures of pure Gaussian operations. We conclude that N W (ρ) is a reasonable measure of Wigner negativity.
Let us now show how the Wigner negativity measure of a given state can be related to the weight of its single-photon component. The Wigner function for an arbitrary Fock state |n reads [17] W n (β) = 2(−1) n π e −2|β| 2 L n 4|β| 2 ) , with |W n (β)| ≤ 2 π since e −x/2 |L n (x)| ≤ 1 (L n are Laguerre polynomials). Hence, for any mixture of Fock states ρ = P n |n n|, we have For P 1 ≥ 1/2, the Wigner function is negative in the phase-space region with where w 0 (x) is the principal branch of the Lambert W function. Eq. (87) defines a disk disk(P 1 ) = β ∈ C |β| 2 ≤ (P 1 ) (88) in phase-space centered at the origin, where the Wigner function is negative. We can now compute the integral over this region that is plotted in Fig. 4. Notably, for P 1 = 1 the bound becomes tight N W (|1 ) = F (1) = 9 4 √ e − 1 ≈ 0.36. To show that the function F (P 1 ) is convex, one computes which is positive since w ≥ 0. In the main text, we defined the function F (P 1 ) continued on the whole interval P 1 ∈ [0, 1] by simply setting F (P 1 ) = 0 for P 1 ≤ 1 2 . Obviously, the continued function remains convex. Furthermore, at P 1 = 1 2 the derivative of the function F is zero F (1/2) = 0, and since F (P 1 ) ≥ 0 we can conclude that F (P 1 ) is non-decreasing on the whole interval.

F Finite statistics
Consider n independent random variables X (i) ∈ [a, b] with its mean E(X) = E( 1 n X (i) ). The Hoeffding theorem [15] gives a simple bound on the deviation of the observed averageX after n trails from the expected value E(X) In our case, the observables X (i) takes values in the interval [−1, 3] so that (b − a) 2 = 16.
Let us now defined X (i) as the minimum of two variables X (i) = min{X We have for the probability We now useP 1 + t ≥ x = min{E(X T ), E(X R )} + t such that P X ≥P 1 + t ≤ P X ≥ min{E(X T ), E(X R )} + t (93) and consider two cases.
If E(X T ) ≤ E(X R ) we use Otherwise, we do the same withX R . For both cases we find that (96) Writing the last expression in the form the latter being a confidence interval forP 1 .
For the calibration-dependent benchmark, one naturally defines the random variable for the quantityP T * 1 (p) of Eq. (19), with positive constants C 1 , C 2 ≥ 0 giving rise to the confidence interval (100) on the average single-photon weightP 1 . Defining Z R similarly (with detector's roles exchanged) gives rise to the confidence interval by exchanging the roles of the detectors. Both are confidence intervals onP 1 , that is It follows that is also a confidence interval onP 1 with confidence level 1 − α.
Finally, let us discuss the witness of Wigner negativity. First, we label by Q the measured value of min{X T ,X R } after n measurement rounds and consider the case Q > 1/2. Given that W ρ (i) (0) ≤ For such a collection we thus have In the calibration-dependent setting with the same argument, the p-value can be obtained analogously. For any value Q and an ensemble of random variables Z (1) Hence, the hypothesisP 1 ≤ 1/2 is associated to the p-value T observed in the experiment. Exchanging the roles of the detectors we also obtain (109) Given the observed data (n •• , n •• , n •• , n •• ) both bounds are valid statements about all collections of states withP 1 ≤ 1 2 . One can then simply choose the most favorable bound, that is p-value ≤ min exp − 2n(Z T − 1 2 ) 2 (C 1 (t,η T ) + C 2 (t,η T ,η R )) 2 , exp − 2n(Z R − 1 2 ) 2 (C 1 (r,η R ) + C 2 (r,η R ,η T )) 2

G Multi-mode product states
We now consider multi-mode product states of the form = k ρ k , • is a lower bound on n P [k] n (λ) 1 2 n . We deduce that max k P [k] 1 ≥ P 1 (p) = min{ P T 1 (p), P R 1 (p)}. (113) Under the assumption that the source produces a multi-mode product state of the form given in Eq. (111), we thus deduce that there is a mode k * , that can in principle be filtered out, such that the corresponding state ρ [k * ] satisfies 1|ρ [k * ] |1 ≥ P 1 (p).

G.1 The proof of Eq. (112)
Consider a multi-mode product state = k ρ k with P [k] 1 ≤ P in each mode. We label (p Our first aim is to analyze the possible set of values Q ∞ P = {(p • , p •• )} for all n from 1 to ∞ and in particular, to show that Q ∞ P ⊂ Q P for P ≥ 1 2 . Naturally, we are interested in the extreme points of this set. Eq. (114) is linear in all points (p • ≤ 1+P 4 . This means that the point (p • , p •• = 0) remains inside the original polytope Q P . We can thus remember that (0, 0) and 1+P 4 , 0 are points of Q ∞ P , but ignore these vertices in the further construction. Analogously, all modes with (p It remains inside Q P if and only if Naturally, if (p • , p •• ) 2 ∈ Q P the next points (p • , p •• ) n are also in Q P . Therefore, Q ∞ P ⊂ Q P for P ≥ 1 2 . This means that for any state of the form = k ρ k , with P In order to extend the analysis to any value of P , we would need to analyze Q ∞ P for a arbitrary P , which is cumbersome. Instead, we analyze the maximal value thatP T 1 (p) takes on Q ∞ P . We know that it takes its maximum value on one of the vertices (p • , p •• ) n , and denote these valueŝ for n ≥ 1. Its maximal value in Q ∞ P is thus given bŷ Let us now look atP T 1 (n) as a function of a continuous parameter n ∈ [0, ∞), and compute its derivative .
It is quite an irregular function, as can be seen in Fig. 5. The boundary valueP T ↑ 1 (0) = lim P →0P T ↑ 1 (P ) = 1 3 can be computed analytically. We show numerically that it is upper-bounded by a simple linear function also can be seen in Fig. 5. By inverting the last inequality we find that Therefore, for any value P such that P T 1 (p) > P , we can conclude that at least one mode satisfies max k 1|ρ [k] |1 > P .

G.2 Finite statistics
Here, in order to obtain a confidence interval we define an independent random variable with E(Y T ) = P T 1 (p). Similarly, we define Y R by exchanging the role of the two detectors and get E(Y R ) = P R 1 (p) with Y T (R) ∈ [−2, 4]. The same exact analysis as above yields a one-sided confidence interval q α (Ȳ T ,Ȳ R ) = min{Ȳ T ,Ȳ R } − 36 log(1/α) 2n (125) on the quantity max k P [k] 1 averaged over all states produced by the source. In the general multi-mode case, the same quantity is a one-sided interval on the quantity λ p(λ) max k P