Estimating the photon-number distribution of photonic channels with realistic devices and applications in photonic quantum information processing

Characterising the input-output photon-number distribution of an unknown optical quantum channel is an important task for many applications in quantum information processing. Ideally, this would require deterministic photon-number sources and photon-number-resolving detectors, but these technologies are still work-in-progress. In this work, we propose a general method to rigorously bound the input-output photon number distribution of an unknown optical channel using standard optical devices such as coherent light sources and non-photon-number-resolving detectors/homodyne detectors. To demonstrate the broad utility of our method, we consider the security analysis of practical quantum key distribution systems based on calibrated single-photon detectors and an experimental proposal to implement time-correlated single photon counting technology using homodyne detectors instead of single-photon detectors.


I. INTRODUCTION
Quantum photonics is the art of using low-light optical signals to exchange and process information in the quantum regime [1,2]. Today, photonic systems represent one of the most promising platforms to implement quantum technology, with already several well-established applications ranging from quantum cryptography [3,4] and communications [5], to sensing and metrology [6][7][8], lithography [9] and imaging [10].
In the most general setting, one considers the preparation, transmission, and detection of optical signals. Here, the photonic channel of interest accepts an N -mode input state and returns an M -mode output state, which is then measured by a series of photon-counting devices. More formally, let q( m| n) be the probability of obtaining m = (m 1 , m 2 , . . . , m M ) photons across the M output modes given n = (n 1 , n 2 , . . . , n N ) photons are injected into the channel (see Figure 1). Here, we note q( m| n) can be characterised independently of the channel's dynamics. That is, the knowledge of the channel is not needed to estimate q( m| n), i.e., we can treat the channel as a black box with N inputs and M outputs and sample accordingly. To keep the characterisation as general as possible, we also allow cases in which the channel is not photonnumber preserving. This can happen when the channel experiences loss and/or suffers from background noise. Correspondingly, if the channel is known to be photonnumber preserving, then we have In practice, the input-output photon-number distribution is central to a multitude of information processing tasks. A good example is boson-sampling-a type of non-universal quantum computation [11,12]. Here, a linear optical interferometer with N inputs and N outputs is considered, where the input is injected with a fixed number of single photons and the output measured with photon-counting devices. In Ref. [11], it was shown that * emilien.lavie@u.nus.edu  Each of the N input modes of the photonic channel is connected to a light source that supplies the input photons and each of the M output modes is connected to a detector that gives an outcome related to the number of photons leaving the photonic channel. Note that the sources and detectors considered here may include some form of active modulation devices, such as intensity modulators.
q( m| n) evaluates directly the permanents of sub-matrices of the interferometer's matrix. On the other hand, solving such matrix permanents with a classical computer is known to be computationally hard. Another example is quantum key distribution (QKD), particularly those using discrete variable encoding [13,14]. For such a communication system (connected by a quantum channel with one input and one output), the estimation of single-photon statistics is essential for protocol security. Take for instance q(1|1), which quantifies how often the untrusted channel behaves as a true singlephoton channel. Having this information strengthens the QKD security analysis by allowing one to assume that (1) the adversary forwards exactly one photon to the receiver and (2) the trusted detector noise and untrusted channel noise are separated. The former is especially powerful as it enables the security analysis of practical QKD under the assumption of a qubit channel (it also applies generally to any qudit channel of interest and hence to high-dimensional QKD as well [15][16][17]).
As a final example, we consider time-correlated singlephoton counting (TCSPC) [18,19], an optical waveform measurement technique that is widely used in fundamental physics research (e.g., ranging [20,21], imaging [22,23], light source characterisation, and life sciences experiments (e.g., fluorescence-lifetime imaging microscopy [24] ). In this setting, the input and output optical modes of the channel are temporal modes corresponding to different time intervals. Here, we have N i=1 n i = 1 and M j=1 m j ≤ 1 since only one photon is deployed in a single trial. In recent years, to improve the efficiency and speed of optical waveform measurement, the idea of using photon-number resolving measurements has been proposed as well [25,26].
The input-output photon-number distribution can be easily estimated if one injects fixed photon-number states |n 1 |n 2 . . . |n N into the channel and measure the output using photon-counting devices (which counts the number of photons in each output mode). However, this would require deterministic photon-number sources and photon-number resolving detectors (PNRDs), which at the moment are still in development [1,2]. The more realistic options are probabilistic photon-number sources (e.g., coherent lasers) and non-photon-number resolving detectors (e.g., threshold detectors and homodyne detectors); these optical devices are not only highly reliable and cost-effective but also widely available.
To this end, it is natural to ask if one can use these standard optical devices to estimate the input-output photon-number distribution of an unknown photonic channel.
It is worthwhile to mention that the problem we are interested in can be seen as a special case of coherentstate quantum process tomography (csQPT), which uses coherent states (probe states) and homodyne measurement to reconstruct the process matrix of an unknown optical channel [27][28][29][30]. Indeed, by looking at only the diagonal components of the process matrix, one can recover the photon-number distribution of the unknown channel. Importantly, since we are only interested in the photon-number distribution, the implementation can be significantly simplified, i.e., the probe states and local oscillators can come from independent laser sources, as we will show later. Note that in the case of csQPT, the relative phase between the probe states and the local oscillator has to be calibrated (needed to fully recover the underlying process matrix), which may be an issue in long-distance quantum communication protocols such as QKD [31,32].
The problem of estimating the photon-number distribution of an unknown optical channel is not new and has been studied before in the field of QKD using threshold detectors. On the input side, decoy-state method has been proposed, which uses phase-randomised light pulses with different intensities to estimate single-photon statistics [33][34][35]. The method essentially entails solving a system of linear equations constrained by the different expected detection rates of the protocol. As such, there are two approaches towards solving the problem, namely one can do it analytically via Gaussian elimination [36,37] or numerically with linear programming [38,39]. The same principle can also be applied to threshold detectors to estimate the output photon-number distribution of the channel [40]. In this approach, called detector-decoy method, one randomly varies the detection efficiency with a variable optical attenuator or intensity modulator to generate a system of linear equations; likewise, these are constrained by the different detection rates effected by the variation of detection efficiency. Given that both decoy-state and detector-decoy methods are based on the same concept, it is thus natural to consider the combination of these two approaches. This direction was recently pursued by the authors of Ref. [41], who used the direct combination of decoy-state and detector-decoy methods to characterise multi-photon quantum interference patterns. Alternatively, the authors of Ref. [42] used a source modulation along with PNRDs to characterise multi-photon quantum interference.
Here, based on the above ideas, we provide a systematic approach to analyse optical communication systems using a linear estimation of the photon-number statistics, extending the decoy-state, detector-decoy and homodyne based linear estimation methods. Our theoretical contributions are three-fold: (1) the extension of Ref. [41] to homodyne detectors, (2) the security analysis of practical QKD systems based on calibrated single-photon detectors, and (3) an experimental proposal to implement TCSPC technology using homodyne detectors instead of single-photon detectors. Concerning the latter, there are two practical advantages in using homodyne detectors: (1) these detectors are typically much more costeffective than single-photon detectors and (2) no active intensity modulation is required to achieve the same effect as detector-decoy. More generally, our extended approach with homodyne detectors provides a simpler and more cost-effective implementation path for applications that requires only the knowledge of q( m| n) instead of single-shot information (see the examples above). We also present two different methods to estimate the desired input-output photon-number probabilities: one based on Gaussian elimination and the other based on linear programming.
Additionally, we highlight that our work is focused on providing interval estimates on the photon-number statistics instead of point estimation. Indeed, our approach is essentially motivated by how parameters are estimated in QKD: there, it is imperative to provide reliable upper and lower bounds on parameters such as bit error rates and detection rates, which characterise the amount of key information leaked to the unknown channel. Therefore, the methods that we describe later include a guarantee on the statistical distance to the true value, unlike other estimation techniques such as maximum-likelihood [43] and least square estimation [44].
The paper is organised as follows. In Section II, we introduce a general channel model and the optical device models used in the estimation. Then in Section III, we present two methods for bounding the desired photonnumber probabilities. Finally, in Section IV we show how our method can be used to analyse the security of practical QKD with calibrated detectors and TCSPC using homodyne detectors.

II. CHANNEL MODELLING
In the following, we keep our analysis to a single-mode photonic channel (i.e. a channel with one input mode and one output mode); the generalisation to multi-mode channels is straightforward. The starting point of our approach is the measurement function, {f (x, y)} x,y , which characterises the observed statistics depending on two designated control parameters x and y owned respectively by Alice (transmitter side) and Bob (receiver side). Here, the parameters are quantities used to test the unknown photonic channel. In the case of active schemes, they are random optical modulations operated by the users. In the case of passive schemes (e.g. passive decoy states [45][46][47][48][49] or homodyne detection as presented later in this paper), they are random variables whose outcomes are correlated to the behaviour of the unknown channel.
In the most general setting, the measurement function is modelled by where p n (x) is the input photon-number distribution (representing correlations between a n-photon state transmission event and Alice's parameter x), q(m|n) is the probability of the channel emitting m photons given it has received n photons, and r m (y) is the measurement response representing the correlation between a mphoton reception event and Bob's parameter y.
As mentioned above, our goal is to estimate certain elements of the unknown channel's input-output photonnumber distribution, q(m|n). To that end, we suppose the input photon-number distribution p n (x), the measurement response r m (y), and the measurement function f (x, y) are fully characterised for any x and y. That is, we assume the user has complete knowledge of the underlying optical devices and has made enough measurements to accurately infer f (x, y).
Similar to standard decoy-state method implementations, we use a phase-randomised coherent-wave laser to generate photon-number states at the channel's input. In this case, the light field entering the channel is described by a Poisson distribution of photon-number states where µ is the mean photon number of the field. The random input x is achieved by modulating the mean photon number with an intensity modulator. As such, the probability model of the source is fully characterised by x and given by For the measurement model, we can use either a threshold detector or a homodyne detector. In the former case, the detector only fires if some photons are detected. As such, there are only two possible outcomes, detection and no detection. Here, we consider only the no detection outcome since the detection outcome is simply the complement event. Following Ref. [40], the response of a practical threshold detector can be modelled using where p dc is the probability of dark count, η det is the single-photon efficiency of the detector, and ν is the transmission efficiency of the intensity modulator (placed in front of the detector) controlled by input y. We note that this model is general and applies to most of today's standard single-photon detection techniques, e.g., single-photon avalanche diodes (SPADs) and superconducting nano-wire single-photon detectors (SNSPDs); see Ref. [50] for an overview of single-photon technology. This model is simple and might be inaccurate under specific operating conditions like fast repetition rate under which other effects like after-pulses might appear. However it is easy to replace the simple model used in Eq. (4) by a more refined model like the one suggested in [51] to account for such effects. For simplicity in this paper, we stick to the simple model to avoid unnecessary complication in the understanding of the underlying method.
In the case of homodyne detection, it does not count the number of photons in the incoming light field but rather gives an outcome whose probability density function is correlated to the number of photons [52]. This relation becomes more apparent when the local oscillator is phase-randomised and the response of the detector when given m photons is given by [43,44]: where y is a real number and {H k (y)} k are Hermite polynomials [53,54].
Here, two observations are in order. Firstly, unlike threshold detectors, one can obtain any number of discrete outcomes by binning y (in practice, an Analogto-Digital Converter (ADC) is used). Secondly, notice that no additional intensity modulation is required here. This is because the response density function of a homodyne detector is sensitive to the range of y and hence one can optimise the binning function (i.e., the ADC) to assign different weights to different input photon-number states. Essentially, this is the same as the detector-decoy method, which assigns different detection probabilities to different photon-number states via the variation of the detection efficiency. Again here, the model we use in Eq. (5) is relatively simple and could be refined to include additional imperfection of realistic detectors like electronic noise [57].

III. METHODS
In most quantum information processing tasks, one is only interested in elements of q( m| n) that are small in the input photon number and output photon number. Additionally, the possible values for the controlled parameters x and y are limited to fixed sets x ∈ X := {x 0 , x 1 , . . . , x n0 } and y ∈ Y := {y 0 , y 1 , . . . , y m0 }. We consider a practically-relevant finite subset of q(m|n) by focusing on n ∈ N 0 := { 0, 1, . . . , n 0 } and m ∈ M 0 := { 0, 1, . . . , m 0 }. Our objective is to derive upper and lower bounds on a specific element or a linear combination of different elements from the set {q(m|n)} n∈N0,m∈M0 . As mentioned, this problem is essentially a linear optimisation problem with constraints given by positivity, normalisation, and the measurement distribution. More specifically, for positivity one has q(m|n) ≥ 0 for any n and m, for sub-normalisation m0 m=0 q(m|n) ≤ 1 for any n, and for measurement distribution for any x and y. In addition, one could also exploit the knowledge of characterised functions p n (x) and r m (y) to construct linear constraints like where h(x, y) is some positive function depending on the optical devices used in the application. We will provide some examples later in Section IV and more technical details in Appendix A. In the following, we present two methods to estimate the desired input-output photon-number statistics.
Linear programming method : Let q(m * |n * ) be the quantity of interest to which an upper bound is desired, then the linear programming (LP) problem is Evidently, the idea behind LP is to use the various constraints on q(m|n) or some linear combination of them to provide bounds for the possible values of q(m|n). Also, since the optimisation is numerical, it is useful to first narrow down to a set of m and n of interest, which can be done via the truncation of both m and n up till some suitable choice of m c ≥ m 0 and n c ≥ n 0 . Notably, LP is performed by first defining a feasible region where the set of input-output distribution {q(m|n)} m≤mc,n≤nc satisfies the constraints. One can then maximise (resp. minimise) the desired probability q(m * |n * ) over the feasible region to obtain the upper (resp. lower) bound on q(m * |n * ).
Analytical method : The basic idea of the second method is to leverage the complete knowledge of the characterised devices to estimate q(m * |n * ) without using any cutoff condition n ≤ n c or m ≤ m c . To that end, we consider a linear combination of the measurement function (see Eq. (1)) over a finite set of evaluation points in X and Y. This gives us a real-valued quantity Λ which is defined as where coefficients {α i } n0 i=0 and {β j } m0 j=0 are real numbers. Also, we write to capture the summation over all the considered evaluation points.
Here, we want to get Λ as close as possible to q(m * |n * ). To do that, we set u n = δ n,n * and v m = δ m,m * for all values of n ∈ N 0 and m ∈ M 0 , where δ a,b is the Kronecker delta function, and solve for α i and β j . In essence, this step requires solving two systems of linear equations, namely one for the source device, and one for the measurement device, Notice that this step does not require the knowledge of f (x i , y j ) and hence can be seen as part of the calibration process prior to characterising the channel. Solving Eq. (11) and Eq. (12) hence gives As one can see, Λ is now expressed in terms of the desired quantity, q(m * |n * ), and some other irrelevant terms that emanate from higher photon number contributions, i.e. those from n ≥ n 0 + 1 and m ≥ m 0 + 1. In Appendix A, we show that these terms can be rigorously bounded by using known information of the optical devices. This in turn provides upper and lower bounds on q(m * |n * ). More concretely, the idea is to establish bounds on u n and v m using the characterised input photon-number distribution and detection response function. Therefore, these bounds are specific to the types of light sources and detectors used in the setup; in Appendix A, we provide standard bounds for common optical devices such as phase-randomised lasers, threshold detectors and homodyne detectors with phase randomised local oscillators. We also note that these bounds can be made arbitrarily tight by selecting large enough n 0 and m 0 values. Indeed, a key condition is to ensure that the derived bounds on the extra terms in Eq. (13) are small when compared to q(m * |n * ); and this can be achieved by using bigger values of n 0 and m 0 . This linear method is conceptually similar to early papers in homodyne tomography using pattern functions to recompute photon-number statistics [55,56]. They considered the use of pattern functions M n (x) such that: Indeed, there is an obvious similarity with our method when using only one input mode: and Λ is a good approximation of p n up to some deviation we can bound. Our computation in Eq. (15) can be seen as a discretized version of Eq. (14). The main benefit is that it can be generalised easily beyond homodyne tomography, for instance using threshold detectors or other light sources, and as such may provide a unified understanding of photon-number probability estimation.

IV. APPLICATIONS
We present here two applications to illustrate the utility of our framework. In the first application, we consider the security of practical prepare-and-measure QKD with realistic photon sources and single photon detectors. More specifically, we show how to rigorously bound the single-photon channel security of the protocol. In the second application, we show how to use our framework to enable TCSPC with homodyne detection instead of single photon detection.
A. Prepare-and-measure QKD with single-photon channel security As a first application of our method, we analyse the asymptotic security analysis of discrete-variable QKD protocols based on practical optical devices such as lasers and threshold detectors. Here, we consider the celebrated Bennett-Brassard 1984 (BB84) protocol [13] and the sixstate protocol [58,59]. These two protocols are formulated as qubit protocols, i.e. the preparation, evolution and measurement can be described in a Hilbert state of dimension 2 [14].
In practice, however, most QKD systems use weak coherent laser sources and threshold detectors [14] to implement qubit states and measurements. This is because, as mentioned above, deterministic single-photon sources and PNRDs are not yet available and weak coherent laser sources and threshold detectors are the closest one can get to achieving qubit states and measurements in practice (at least with regards to cost and practicality). However, one critical drawback is that there are some fundamental differences between the qubit models (assumed in the protocol) and these optical devices. Some of these differences are irreconcilable and hence cannot be applied to qubit protocols, while some may lead to implementation loopholes such as photon-number splitting attacks [60,61].
Fortunately, for most qubit protocols (including BB84), the gap between theory and practice can be mitigated using innovative techniques such as the decoy-state method [33][34][35] and squashing [62][63][64][65]. In the case of the former, assuming that the prepared optical signals are diagonal in the photon-number basis, Eve's attacks can be categorised according to the emitted photon numbers. This allows us to focus on the single-photon component of emitted optical signal and hence view the states prepared by Alice as qubits. On the other hand, squashing models provide an elegant method to map the actual full (infinite dimensional) mode measurement onto a finite dimensional Hilbert space followed by the ideal measurement. In the case of the BB84 protocol, a squashing model exists and hence one could assume qubit models for Bob's measurements in practice.
Therefore, using the decoy-state method together with a squashing model, one can derive statistical bounds on the single-photon error rates of the BB84 protocol and hence compute its secret key rate. However, this approach is quite restrictive and does not apply readily to other qubit protocols. For instance, it has been shown that squashing does not immediately apply to the sixstate protocol [63,66]; on the other hand, it has been shown that by relaxing certain statistical constraints, it is possible to define squashing maps for a wide range of finite-dimensional protocols [65].
By contrast, our method allows us to analyse any qubit (more generally, any higher dimension) protocol [15][16][17] without using a squashing model. Assuming that Alice and Bob prepares and receives a single photon defined across two orthogonal optical modes, the states and measurements can be described by qubit states and qubit measurements, respectively. Our method allows us to bound the probability of Bob receiving a single photon (just before the measurement) when Alice prepares and sends a single photon as well as the corresponding error rate in each basis.
Roughly speaking, our method provides three practical advantages over existing methods. Firstly, since we consider only secret key contributions from events in which Alice prepares a single photon and Bob receives a single photon, we can directly use any security proof technique for qubit models without applying any squashing model. As such, our method can be applied to most practical QKD systems under the condition that these systems randomly vary their detection efficiency as specified by the detector-decoy method. Secondly, since squashing models typically require mapping double-detection events to random outcomes [63,66], applying a squashing model would likely introduce some additional errors from the detector background noise. On top of that, squashing also does not differentiate between clicks due to true single photon detections and empty detections, which may also introduce additional errors. Hence, our method, which can rigorously bound the true channel error rates, could give an enhanced bound on the secret key rate especially in the high loss regime where the dark count rate is not negligible. Finally, as mentioned, our method allows us to analyse the security of the protocol based on single-photon channel security. To appreciate this feature better, we comparatively note that when using the decoy state method combined with a squashing model, one actually evaluates the security of the channel together with the detector noise, which is normally trusted. In this case, the single-photon error rates includes the trusted detector noise. By contrast, our method allows the rigorous separation of channel noise and detector noise and thus provides a concise method to derive lower bounds on the secret key rate in the calibrated detector setting; in fact, the security of QKD with calibrated devices is known to be an open problem [14] .
We now demonstrate how our method can be applied to the security analyses of practical QKD systems. We emphasise that our method can be used for most discretevariable protocols, but as concrete examples, we only apply our method to BB84 and six-state protocol. To that end, we introduce some notations that we use in this subsection. We denote the basis choice of Alice and Bob by x and y respectively. x and y are randomly chosen from the set X and Y respectively where X = Y = {X, Z} for BB84 and X = Y = {X, Y, Z} for the six-state protocol. The symbol value encoded by Alice is denoted by a ∈ A and Bob's detection pattern is denoted by b ∈ B. For both the BB84 and six-state protocol, a ∈ {0, 1} and b = b 0 b 1 is a two-bit string where b i indicates whether the detector in mode i clicks (we have b i = 0 when the detector in mode i does not click and b i = 1 when the detector clicks). Based on b, Bob would then map the observed click pattern into the decoded symbol or he could choose to discard inconclusive events (such as no-click or double-click events). Finally, we denote Alice's intensity setting by µ and Bob's detection efficiency by η (composed of η 0 for detector 0 and η 1 for detector 1) which are chosen randomly in each round. We will now focus our attention on the case where |A| = 2, i.e., the protocol uses binary symbols. The generalisation to a higher dimensional protocol is straightforward.
In a typical discrete-variable protocol, Alice would randomly choose a basis x and bit value a. She would then choose an intensity setting µ and then prepare a phaserandomised coherent state in the corresponding mode. For phase-randomised coherent source, the emitted photon number n would follow a Poisson distribution with mean µ. Similarly, Bob randomly chooses a basis choice y as well as detection efficiency setting η and he obtains the outcome b. In the parameter estimation step, Alice and Bob can estimate the following conditional probabilities (for all possible combination of parameters): y, a, µ, η).
This function can be expanded as where p n (µ) is the probability of the source emitting n photons (defined in Eq. (3)) while q {xya} (kl|n) is the probability of k photons arriving in mode 0 and l photons arriving in mode 1 given that the source emitted n photons, Alice chooses the basis x and bit value a and Bob chooses measurement basis y. r {b0} k (η 0 ) denotes the probability of the detector in mode 0 clicking (b 0 = 1) or not (b 0 = 0), given that k photons arrived in that mode and r {b1} l (η 1 ) is defined similarly. Note that the probability of a threshold detector not clicking is given in Eq. (4).
As such, Eq. (17) represents the conditional probability as a product of different system elements (transmitter, channel and receiver) in the form given in Eq. (1). The parameters x, y, a and b are considered as fixed parameters when estimating the input-output photonnumber distribution. Here, we remark that these parameters serve only as a means for Alice and Bob to organise their measurement data, i.e. they categorise the inputoutput photon-number distributions according to x, y, a and b. Of note, the input-output photon-number distributions have to be independent of these parameters: this is needed to ensure that the single-photon channel behaviour is basis-independent. In other words, Eve's attacks on the quantum channel have to be independent of Alice's and Bob's basis choices [62]. In addition to this, we also need that the input-output photon-number distributions are independent of µ and η. In practice, these conditions can be reasonably enforced by using phaserandomised coherent lasers and photon-counting detectors, which is the case in our consideration.
Notice that the n-photon yield, which is denoted by Y {xyab} n here, is the usual quantity of interest in decoystate QKD [34]. More precisely, in Ref. [34], the authors considered events in which Alice and Bob choose the same basis, i.e. x = y, and the cases in which Bob observes at least one click. They also average the yield over Alice's bit value a. Hence, the n-photon yield is the probability of observing a click given that Alice's laser emits n photons. Clearly this would depend on both the channel and the trusted detectors which are located in Bob's lab. In contrast, our method bounds the probability of k and l photons arriving at Bob's measurement device given that n photons are prepared by Alice. This would depend only on the behaviour of the channel and not on the detectors.
To analyse the security of the BB84 and six-state protocol, the quantity of interest is q {xya} (kl|n) when n = 1 and k + l = 1, i.e. the probability of the channel outputting a single photon given that a single photon enters the channel. In this case, the conventional security proofs [67,68] that rely on the qubit models can be directly applied without the need for squashing model [62,63,65]. Key rate per channel use 6 states with source and detector modulation BB84 with source and detector modulation BB84 with source modulation only Figure 2. We simulate the achievable secure key rate with two avalanche photo-diode detectors, each one featuring a dark count rate of 10 −6 and a fixed channel error rate of 5% in all bases (depolarising channel). The secure key rate in the source and detector modulation case is K ≥ p0(µ)q 00|0 1 − r (0) 0 (η0)r  . We compare the exact channel error rate value used for the simulation (5% here) to the upper bound provided by our proposed analytical method and the one proposed in Ref. [36]. The definition in Ref. [36] is including the noise from the detector dark counts while our is considering only the noise on the channel. As a result, the definition of Ref. [36] is increasing while ours is staying close to the exact value in the high loss regime, hence a slight improvement in distance for the key rate. Note that this improvement is only affecting privacy amplification, since the noise due to dark counts still has to be corrected in the error correction step.
ple, the single-photon error rate given that Alice and Bob have chosen the same basis (i.e. x = y) is given by where the qubit detection probability is Once the single-photon error rates are determined, the secret key rate can be easily computed. Here, we present the bound on the secret key rate while we defer the detailed security analysis of the six-state protocol to Appendix B. In the asymptotic limit and under the assumption of collective attacks, the secret key rate K of BB84 and six-state protocols is given by (20) where H(A|E) is single-photon conditional entropy given Eve's quantum side information; h 2 (·) denotes the binary entropy function. Q(µ, η) and E(µ, η) is the observed gain and quantum bit error rate when Alice and Bob choose intensity µ and detection efficiency η respectively. Hence, the first term is the contribution due to the events in which Alice prepares vacuum state. Since no quantum information is leaked whenever Alice prepares the vacuum state, all the detected events due to the transmission of vacuum states are secure. The second term is the single-photon contribution and the third term is the leakage due to error correction assuming that the error-correcting code saturates the Shannon limit. Using {e X , e Y , e Z } as short-hand for the single-photon error rate in the respective basis, the single-photon conditional entropy H(A|E) for the BB84 protocol and the six-state protocol is given by Here, λ = (λ 0 , λ 1 , λ 2 , λ 3 ) is a real vector containing the unique solution to the following simultaneous equations and H(λ) is the corresponding Shannon entropy. Therefore, by substituting the appropriate H(A|E) to Eq. (20), we obtain the bound on the secret key rate of the corresponding protocol. We present the simulated secret key rates in Fig. 2 assuming standard SPADs parameters. Here, we compare against the asymptotic secret rate based on the standard decoy-state method [36]. As mentioned above, the key difference is that our method directly evaluates the single-photon channel security whereas the standard decoy-state method would include the detectors' background noise (dark counts). Indeed, in Fig. 3, we see that the single-photon error rate of our method does not include the detectors' dark counts in the channel error rate and hence our proposed upper bound on the error rate remains close to the exact value while that based on Ref. [36] is dominated by the dark count noise in the high loss regime.

B. Time-Correlated Single Photon Counting with
Homodyne detection Here we propose TCSPC with homodyne detection instead of single photon detection. There are significant practical benefits in doing so; largely one could reduce the implementation cost and footprint of TCSPC through integrated photonics platforms. The basic idea of TCSPC is to measure the single photon emission time profile in the nanosecond time scale, e.g. fluorescence decays of excited samples [24]. The current method achieves this via the small time resolution (hundreds of picoseconds) of Single Photon Detectors (SPDs; e.g. photo-multiplier tubes, micro-channel plates, SPADs) to measure the time difference between a reference "start" signal and a "stop" signal triggered by a single photon emission event [18,19]. As such, by using a pulsed laser the intensity profile can be sampled repeatedly and a histogram of photon arrivals per time bin over the intended time domain can be constructed; assuming the probability of multi-photon emission is negligible. This concept is depicted in Fig. 4.
However, SPDs typically suffer from finite recovery time (or dead-time); consequently, the detector becomes inactive for a period of time after a first detection event (hundreds of nanoseconds to dozens of microseconds [50,69]). As such, any optical signal arriving during this time window will not be detected and this problem tends to bias the measurement results towards earlier detection events. This is a well known issue called pulse pileup [70]. In practice, to mitigate this problem, a popular approach is to keep the multi-photon emissions low and the counting rate below 2%-5% or lower [70], e.g. by restricting the excitation power. In this case, nothing is detected most of the time and once in a while a unique photon is detected and recorded. This common approach of limiting the excitation power solves the multi-photon emission issue but leads to a longer acquisition time.
Yet, TCSPC is not making full use of the single-shot information of a photon being detected or not at a particular time window; it extracts only the average count rate Original TCSPC with SPD Figure 4. The original TCSPC concept is based on the recording of single-photon detection events and construction of the corresponding histogram as shown above. As a result, there is an implicit post-selection of the conclusive outcomes. That is, the experiment is repeated until the bin with the maximum number of events reaches a certain level. The histogram is then normalised to derive the probability of having a detection event in a particular bin t ∈ T given that there was a conclusive outcome C, i.e. Pr(t|C). Nevertheless, due to the dead time effect, the recorded probability is not exactly the one aforementioned, but rather the probability of observing a detection in a bin and not recording any detection before. When the counting rate is low, the two probabilities are close though [24].
for each time window. This suggests that other forms of detection technology could be used instead, for instance, homodyne detection. Indeed, this possibility has already been discussed in Ref. [71]: the authors therein described a linear method to compute the moments of an unknown probability distribution using the moments of the outcome function obtained experimentally. Our proposal with homodyne detection essentially follows this idea: that we can recover the photon-number statistics with phase-randomised homodyne detection. To help fix ideas, in the following we first briefly describe an ideal version based on PNRDs. We then provide a proof of concept simulation of the homodyne TCSPC technique.

TCSPC with perfect photon-number-resolving detection
It is useful to first consider an intermediate ideal TCSPC protocol to illustrate the main ideas of our homodyne-based protocol. Here, we assume a perfect PNRD is used to measure the photons arrival time. By perfect, we mean that the detector has zero dead-time, perfect detection efficiency, and able to tell how many photons are detected in a given time window. Mathematically, the outcome of the protocol is described by a sequence of time-ordered random variables, Y {t} ∈ N, where each random variable counts the number of pho-TCSPC with perfect Photon Number Resolving Detector Figure 5. We consider here an ideal version of TCSPC where perfect PNRDs are used to record the exact number of incoming photons in each time bin. We label Y {t} ∈ N the corresponding random variable for each time bin. After enough events have been recorded, it is possible to estimate the photon number probability distribution at each time bin. By keeping only the single photon events, it is easy to recover the same TCSPC information as in Fig. 4. tons in the time bin t ∈ T as shown in Fig. 5. Evidently, this ideal TCSPC protocol can recover the original TC-SPC protocol's information by keeping only the events in which Y {t} ≥ 1. We also highlight that there will be the same number of events recorded (regardless of the value of Y {t} ) in each time bin. As a result, the probability of recording an event in a certain time window is Pr(t) = 1/|T |.
In the limit of many repetitions, the data from Y {t} allows computation the probability of detection of n photons in the time bin t: We label C a conclusive event ; it is the set of photon number values leading to a conclusive outcome, here all values n ≥ 1. We further denote the probability of a conclusive event within a time window: From this information only, it is possible to recompute the same probability as in the original version of TCSPC using the uniformity of T and Bayes' rule:

Homodyne TCSPC
We now replace the perfect PNRD in the ideal version described above by a homodyne detector that sequentially measures all the time bins one after the other, up TCSPC with Homodyne Detector Figure 6. We consider here a practical implementation of TC-SPC using a binned homodyne detector. We label X {t} ∈ X the corresponding random variable for each time bin. After enough events have been recorded, it is possible to estimate the measurement PDF at each time bin. Then from this information, we show that it is possible to recompute the same photon number distribution as in Fig. 5.
to the time resolution and speed of the ADC. In this case, the homodyne detector gives a continuous outcome which is also binned depending on the ADC resolution.
We label X the set of bins. Now the sequence of random variables Y {t} ∈ N is replaced by X {t} ∈ X as drawn in Fig. 6.
In the limit of many repetitions and small bins in X , we can recompute the probability density function (PDF) of every time bin which is assumed to have the following structure: with A n (x) representing the homodyne measurement response as we define in Section II Eq. in Eq (26) are defined in Eq (23) and represent the contribution to the outcome due to n photons on the detector. From the value of f {t} (x) at each x ∈ X , the direct application of Section III allows us to recompute a reasonably good estimation of the weights q {t} n in front of the A n (x) in Eq (26) for the low photon number events. Then recomputing the probability distribution as in the original TCSPC can be done with Eq (25).
We consider a simple physical experimental model to highlight the feasibility of our homodyne TCSPC. Here, we consider the intensity profile of a fluorescence decay after a delta excitation happening at t 0 = 50 ns. Following Ref. [72] and assuming an exponential decay for the intensity profile, the photon number distribution is where E(t) is the energy arriving on the detector for t ≥ t 0 and T = 5 ns is the time bin duration, τ = 100 ns is the decay time, α = 0.9 is a coefficient including the excitation power and the detector sensitivity. Then, the measurement PDF is computed according to Eq. (26) and a direct application of Section III allows us to compute upper and lower bounds on the single and two photon emission probabilities. Here, we use 16 bins evenly spaced over the range [0, 5] and apply the analytical method presented in Section III to obtain the results shown in Fig. 7.
The bounds on the single photon emission probability are very close to the exact value for any time considered in our simulation while the bounds for the two photon emission probability are less tight due to its lower value that cannot be estimated well via our method. This is due to the finite value of the extra terms in Eq. (13) as discussed in Section III. One way to tighten this bound would be to consider more bins for the PDF.

V. CONCLUSION
In this paper, we show the possibility of using realistic light sources and detectors to recompute relevant information about the input-output photon number distribution of any unknown channel. We describ a simple linear framework to model characterised sources, detectors and any multi-mode unknown channel. Then, we present two computational methods to derive upper and lower bounds on the input-output photon number distribution. Such information can be used for various applications in quantum optics and quantum information processing. To that end, we highlight two applications: the singlephoton channel security of practical QKD and TCSPC with homodyne detection. For example, the application to QKD shows that this framework is a bridge between practical implementations and theoretical qubit security proofs. Here, it is useful to mention that even though we present only one class of qubit protocols and an associated security proof, other finite-dimensional QKD protocols can also be proven secure in the same fashion, such as the Reference Frame Independent [73], Loss-Tolerant [74] and Tomography-based [75,76] protocols. We can also highlight that the single photon detectors could possibly be replaced by homodyne detectors in certain schemes [77]. The TCSPC example also suggests that this technology could possibly be deployed with homodyne detection for more cost-effective implementations. Therefore, many applications relying on TCSPC as a module could potentially be upgraded to use a homodyne TCSPC instead. For instance, we can think of applications based on time-of-flight measurement [78] such as ranging [20,21] and low light imaging [22,23].
The goal of this section is to show that the estimate given in Eq. (13) approximates well the quantity of interest q(m * |n * ). We give practical bounds on the extra terms. First we describe the general strategy and arguments that we use regardless of the actual hardware under use, and then we give the explicit bounds for a laser source, a threshold detector and a homodyne detector.
We recall a few notations: Let us define the residuals:

General strategy
The general strategy is to compute upper and lower bounds on these individual residuals to obtain bounds on R = R n0 + R m0 + R n0m0 = Λ − q(m * |n * ). We highlight below a few general remarks that are useful for the subsequent derivation.
1. The sequences u n and v m depend respectively on n 0 and m 0 , since it is the number of terms in the summation.
The values of α and β also depend on n 0 and m 0 . Hence one has to be careful when considering n 0 or m 0 → +∞ 2. The series u n has a finite zeroth moment ( n≥0 u n < +∞) and first moment ( n≥0 nu n < +∞). That is because it is related to the probability distribution of a source with finite energy.
3. The series v m does not have a finite zeroth moment in general since it is related to a conditional probability. However, the sequence is always bounded.
4. ∀n, m ≥ 0 : 0 ≤ q(m|n) ≤ 1 5. The series q(m|n) has a finite zeroth moment in m for all n: The best we can do for bounding R n0 is to bound q(m * |n) by +1 or −1 depending on the sign of u n . Let us denote: Then the bound is related to the remainder of the convergent series u ⊕ n and u n : For this one, we need to combine the series q(m|n) and v m together since the former is convergent in m while the latter is not. Similarly, let us denote: Additionally we assume that the sequence v ⊕ m is non-increasing for m ≥ m 0 + 1 and similarly, v m is non-decreasing for m ≥ m 0 + 1. If it is not the case, we assume that we can find an upper (lower) bound on v ⊕ m (v m ) that satisfies this property, and we use it instead. This assumption helps to compute an upper bound and lower bound on v m by simply considering the first element: In that case, we can derive the following bound: We can see that this bound depends on our ability to find a good bound on: As a first approximation, we can useq ≤ 1. It is possible to start from this rough estimate to obtain a valid upper and lower bounds on a few q(m|n * ) in the range m ∈ { 0 . . . m 0 } and then use this information to update the bound onq. After a few iterations of this procedure, the bounds on R m0 usually become of the same order of magnitude as those on R n0 . Alternatively, it is also possible to analyse some geometric-arithmetic sequence to compute the limit after many iterations.
c. Bound on Rn 0 m 0 This bound is easy to compute since we do not have much information anyway. The best we can do is to bound m≥m0+1 q m|n by 1.
The reader can notice that the bound on R n0m0 is roughly the product of the previous bounds on R n0 and R m0 hence it is small in practice, and the main contribution comes from the residuals related to source only R n0 and detector only R m0 . Similar to the previous cases, it is also possible to consider quantities like u n v m ⊕ to refine R n0m0 but the improvement is limited.

Phase randomised laser
We consider here a Poisson distribution for the source p n (x) = e −x x n n! , and n 0 + 1 evaluation points 0 ≤ x 0 < x 1 < · · · < x n0 ≤ x max with x max ≤ n 0 .
We notice that p n (x) is non-increasing in n when n is larger than x max . It is useful to keep the sign in α otherwise the bound quickly becomes loose, therefore it is not satisfactory to simply bound α by its 1-norm or positive or negative part. However, due to the simple expression of p n (x), it is easy to check that u n has a constant sign. Indeed, we find that u n ∼ n→+∞ α n0 p n (x n0 ) so for n large enough, the sign of u n is given by the sign of α n0 . For the first few n before that property is satisfied, it is possible to numerically check a finite number of values to establish a correct upper and lower bound. For the rest, we can take u ⊕ n = u n and u n = 0, if u n0+1 ≥ 0 (A16) u ⊕ n = 0 and u n = u n , if u n0+1 < 0 (A17) In our examples, the sign is always constant and depending on the parity of n * − n 0 , similar to what Ref. [37] reported for decoy states with their bounds X n and Z n alternatively being lower or upper bound.
We first notice that r m (y) is non-increasing in m for any y. We also have the following property: from which we find that v m is decreasing if it is positive and increasing if it is negative. Similar to the laser case, we find v m ∼ m→+∞ β m0 r m (y m0 ) where y m0 is the lowest value for y. Therefore the sign of v m is given by the sign of β m0 for m large enough, then by monotonicity we find that the sign is constant for all m. Therefore similar to the laser case, we take:

Homodyne detectors
We use a homodyne detector at the receiver. The detection function is as follows [43,44]: Since the functions A m (y) are even, we can assume without loss of generality that the outcome y is non-negative. We further assume that y ∈ [0; y max ]. This is motivated by the finite range of the ADC that will restrict the maximum observable value for the outcome.
We recall a few useful properties of a m (y) [53,54,79]: √ 2ma m (y) = ya m−1 (y) − a m−1 (y) (A25) a m (y) = ya m−1 (y) − √ 2ma m (y) (A26) We define two functions for m large enough to ensure 2m − y 2 max > 0: g m (y) = a m (y) 2 + a m (y) 2 2m + 1 − y 2 (A27) h m (y) = a m (y) 2 + a m (y) 2 2m − y 2 (A28) g m (y) was defined by Szego in Ref. [80] to prove Sonin's theorem for Hermite functions. More specifically for our purpose, g m (y) is non-decreasing in y. We can show using Eq. (A25) and (A26) that h m (y) also satisfies [79]: h m (y) = a m−1 (y) 2 + a m−1 (y) 2 2m − y 2 (A29) and then g m+1 (y) ≤ h m+1 (y) ≤ g m (y) ≤ h m (y) hence g m (y) is non-increasing in m. Evidently, we also have a m (y) 2 ≤ g m (y). We define as usual: On the interval of interest, the functions A m (y) are oscillating between 0 and some local maxima. More precisely, the plot of |a m | has m + 1 "bumps" before quickly decreasing to zero. When the detector efficiency is lower than 1, the oscillations are attenuated for the first few local extrema, and the last bump remains predominant. As a result, the series v m is also oscillating somehow, and it is difficult to analyse.
To simplify the analysis, we use an upper bound defined as follows for m large enough to ensure 2m + 1 − y 2 max > 0: Then it is sufficient to take: min(0, β j )G m (y j ) (A33) efficiency, inconclusive outcome b = 00, Bob will disclose all outcome results b. By doing so, they can estimate the statistics f {xyab} (µ, η) and use Section III to compute q {xya} (10|1), q {xya} (01|1) for all relevant x, y, a and then the achievable secure key rate the way we describe below. If the latter is positive, they proceed to step 5, otherwise they abort the protocol.

5.
Post-processing: For the remaining rounds, Alice and Bob apply suitable error correction and privacy amplification procedures to extract the secret key.
The quantities q {xya} (kl|n) when k + l = n = 1 are enough to characterise the proportion of events that can be analysed using conventional single photon security proof. We show below how to relate them to the security of the qubit six-state protocol.
We define the qubit detection probability for a basis match as: In other words, a qubit detection happens when there is only one photon in any arm of the receiver and any conclusive detection pattern occurs.
Let us assume that the detectors are labelled in a way that in the absence of noise and x = y, there could be photons arriving only on detector number a, according to Alice's bit. This allows us to define an ideal detection, hence any other result would be an erroneous detection. We can define the success rate for the match basis: where ⊕ is the addition modulo 2 and the error rate is defined as p x=y err = 1 − p x=y suc . For simplicity and to relate to existing notation in the literature, we denote: We rephrase the results summarised in Appendix A of Ref. [14] that is if we denote λ = (λ 0 , λ 1 , λ 2 , λ 3 ) to be the unique solution of: and H(λ) = − 3 i=0 λ i log 2 (λ i ) the Shannon entropy, then the conditional entropy on Alice's bit value given Eve's side information is as follows: We also define the signal detection rate and error rate for the key basis 0: In the asymptotic analysis, we can always consider an efficient implementation where the maximum intensity µ and efficiency η are used most of the time. Otherwise, the framework allows to extract key from any -key state-