Proper error bars for self-calibrating quantum tomography

TL;DR

This paper develops a method for accurately estimating quantum states and measurement device parameters simultaneously, providing reliable error regions, with applications to quantum cryptography and satellite-based quantum key distribution.

Contribution

It introduces a joint estimation framework for self-calibrating quantum tomography using optimal error regions, improving uncertainty quantification.

Findings

Successfully applied to BB84 measurement scenario

Reconstructed quantum state and device parameters simultaneously

Demonstrated relevance for satellite-based quantum key distribution

Abstract

Self-calibrating quantum state tomography aims at reconstructing the unknown quantum state and certain properties of the measurement devices from the same data. Since the estimates of the state and device parameters come from the same data, one should employ a joint estimation scheme, including the construction and reporting of joint state-device error regions to quantify uncertainty. We explain how to do this naturally within the framework of optimal error regions. As an illustrative example, we apply our procedure to the double-crosshair measurement of the BB84 scenario in quantum cryptography and so reconstruct the state and estimate the detection efficiencies simultaneously and reliably. We also discuss the practical situation of a satellite-based quantum key distribution scheme, for which self-calibration and proper treatment of the data are necessities.

Tables

Table 1. Table 1: Mock-true and ML-estimated values of the 10 parameters used for the simulated data of Fig.  4 . Since there are only 66 detection events in this simulation, the estimated values do not agree well with the mock-true values.

Parameter	Mock-true value	ML-estimated value
$⟨1\otimes {\sigma }_x⟩$	-0.1201	-0.2658
$⟨1\otimes {\sigma }_z⟩$	-0.0803	-0.0578
$⟨{\sigma }_x\otimes 1⟩$	-0.0592	0.2200
$⟨{\sigma }_z\otimes 1⟩$	0.3783	0.1643
$⟨{\sigma }_x\otimes {\sigma }_x⟩$	-0.0182	-0.0736
$⟨{\sigma }_x\otimes {\sigma }_z⟩$	0.4009	0.5693
$⟨{\sigma }_z\otimes {\sigma }_x⟩$	-0.0434	0.0488
$⟨{\sigma }_z\otimes {\sigma }_z⟩$	0.1359	-0.1060
${\eta }_{\text{left}}$	0.6755	0.5831
${\eta }_{\text{right}}$	0.7746	0.6565

Table 2. Table 2: Mock-true and ML-estimated values for the 11 parameters used for the simulated data of Fig.  6 . There are 300 898 detection events in this simulation and, owing to their large count, some of the ML estimates of the parameters are comparable with their mock-true values.

Parameter	Mock-true value	ML-estimated value
$⟨1\otimes {\sigma }_x⟩$	-0.4062	-0.4095
$⟨1\otimes {\sigma }_z⟩$	-0.0478	-0.0421
$⟨{\sigma }_x\otimes 1⟩$	-0.1985	-0.3878
$⟨{\sigma }_z\otimes 1⟩$	0.0085	-0.0909
$⟨{\sigma }_x\otimes {\sigma }_x⟩$	0.3595	0.3190
$⟨{\sigma }_x\otimes {\sigma }_z⟩$	-0.0607	-0.1802
$⟨{\sigma }_z\otimes {\sigma }_x⟩$	-0.0664	-0.1180
$⟨{\sigma }_z\otimes {\sigma }_z⟩$	0.4192	-0.1524
$T{\eta }_{\text{left}}$	$7.3771\times {10}^{-5}$	$9.7099\times {10}^{-5}$
${\eta }_{\text{right}}$	0.7238	0.7435
$\nu$	500 000	486 868

Equations

$$p_k=\mathrm{tr}\left\{{\Pi }_k(\mathbf{\alpha })\rho \right\}=⟨{\Pi }_k(\mathbf{\alpha })⟩.$$

$$L(D|\rho ,\mathbf{\alpha })=\sum _{n_0=0}^{\infty }L(D,n_0|\rho ,\mathbf{\alpha }),$$

$${\mathcal{R}}_{\lambda }(D)=\{(\rho ,\mathbf{\alpha })|L(D|\rho ,\mathbf{\alpha })\ge \lambda L_{\text{max}}(D)\}$$

$$s_{\lambda }(D)=\underset{{\mathcal{R}}_{\lambda }(D)}{\int }(\mathrm{d}\rho )(\mathrm{d}\mathbf{\alpha })$$

$$(\mathrm{d}\rho )$$

$$(\mathrm{d}\mathbf{\alpha })$$

$$c_{\lambda }(D)=\frac{1}{L(D)}\underset{{\mathcal{R}}_{\lambda }(D)}{\int }(\mathrm{d}\rho )(\mathrm{d}\mathbf{\alpha })L(D|\rho ,\mathbf{\alpha }),$$

$$L(D)={\int }_{{\mathcal{R}}_0}(\mathrm{d}\rho )(\mathrm{d}\mathbf{\alpha })L(D|\rho ,\mathbf{\alpha })$$

$$c_{\lambda }=\frac{\lambda s_{\lambda }+{\int }_{\lambda }^1\mathrm{d}{\lambda }^{\prime }{s}_{{\lambda }^{\prime }}}{{\int }_0^1\mathrm{d}{\lambda }^{\prime }{s}_{{\lambda }^{\prime }}}$$

$${\lambda }_{\text{crit}}(D)\equiv \frac{L(D)}{L_{\text{max}}(D)}={\int }_0^1\mathrm{d}\lambda s_{\lambda }.$$

$$\displaystyle\left.\begin{array}[]{l}p_{1}\\ p_{2}\end{array}\right\}=\frac{1}{4}(1\pm z)\,,\quad\left.\begin{array}[]{l}p_{3}\\ p_{4}\end{array}\right\}=\frac{1}{4}(1\pm x)\,,$$

$${\eta }_2{\eta }_3^{\prime }⟨\frac{1-{\sigma }_z}{2}\otimes \frac{1+{\sigma }_x}{2}⟩$$

$$L(D|\rho,\bm{\eta})=\mathrm{e}^{\mbox{\footnotesize$\nu p_{0}$}}\prod\limits_{k=1}^{24}{p_{k}}^{n_{k}}\,,$$

$$p_0=1-\sum _{k=1}^{24}{p}_k$$

$$L(D|\lambda {\mathbf{s}}_1+(1-\lambda ){\mathbf{s}}_2,\mathbf{k})\ge$$

$$\text{mbox}+(1-\lambda )L(D|{\mathbf{s}}_2,\mathbf{k}),$$

$${\eta }_1^{\prime }:{\eta }_2^{\prime }:{\eta }_3^{\prime }:{\eta }_4^{\prime }=0.4172:0.5510:1:0.6777,$$

$${\eta }_1:{\eta }_2:{\eta }_3:{\eta }_4=0.6595:1:0.6287:0.7619,$$

$$L(D|\rho,\bm{\eta},\nu)=\prod_{k=1}^{24}(\nu p_{k})^{n_{k}}\mathrm{e}^{\mbox{\footnotesize$-\nu p_{k}$}}\,,$$

$$w_0({\eta }_{\text{right}})$$

$$w_0(T{\eta }_{\text{left}})$$

$$w_{0}(\nu)\propto\nu^{99}\mathrm{e}^{\mbox{\footnotesize$-\frac{\nu}{5000}$}}\,.$$

$${\eta }_1^{\prime }:{\eta }_2^{\prime }:{\eta }_3^{\prime }:{\eta }_4^{\prime }=0.7064:0.5242:1:0.3419,$$

$${\eta }_1:{\eta }_2:{\eta }_3:{\eta }_4=0.7518:0.7520:0.6969:1.$$

$$L(S|\rho ,\mathbf{\eta },n_0)=\frac{(N+n_0)!}{N!n_0!}\prod _{k=0}^{24}{p}_k^{n_k},$$

$$L(S|\rho ,\mathbf{\eta },\nu )$$

$$L(D|\rho ,\mathbf{\eta },\nu )$$

Full text

Keywords:

self-calibration, optimal error regions, plausible region, space-based quantum key distribution

I Introduction

Quantum tomography LNP649 ; Teo:16 is a basic yet crucial element in most quantum information processing tasks. In the typical quantum state estimation scenario—the focus of our paper—a finite number of quantum systems, all prepared in the same unknown state, are probed by a measurement of one’s choosing. From the gathered data, comprising the observed sequence of measurement detector clicks, one tries to arrive at a best guess of the unknown state. The notion of “best” here depends on one’s choice of figure of merit. A popular choice is to maximize the likelihood of the data, giving the maximum-likelihood (ML) estimator as the best guess of the state MLEreview .

Standard quantum tomography strategies rely on measurement devices that are fully characterized and well calibrated in advance, so that the unknown quantities are solely those of the state to be estimated. In practice, there are situations where this in-advance calibration is not possible for all aspects of the measurement devices. For example, the measurement device may have some properties that drift over time, and hence require frequent recalibration. In this case, it may be useful to do the calibration of those properties of the measurement device at the same time as the actual tomography of the state. Self-calibrating quantum state tomography deals precisely with such a situation, where the measurement devices are not fully characterized. The properties of the measurement devices that are not calibrated in advance and the state parameters can be reconstructed from the same gathered data.

The first initiative towards self-calibrating quantum state tomography was taken by Mogilevtsev et al. PRA79.020101 , with a scheme for reconstructing the quantum state and for quantifying the mismatch between the signal and reference states. Later, in Ref. PRA82.021807 , Mogilevtsev presented another scheme for reconstructing the state and calibrating the single-photon detectors simultaneously by exploiting some partial knowledge about the state. By having squeezed noisy signal states in an on/off detection scheme, this work enables absolute calibration of single-photon detectors in a very simple way. The first experimental realization was achieved by Brańczyk et al. NJP14.085003 , who estimated the unknown rotation angle of the measurement basis and the state parameters simultaneously. There, the states of polarization-encoded photonic qubits were reconstructed using wave plates with unknown retardance. In Ref. PRA87.062118 , Quesada et al. showed how to circumvent the requirement for well-characterized unitary operations in quantum state estimation by treating unknown parameters in the state and the unitary operations on an equal footing. Recently, Williams and Lougovski njp19.043003 used a Bayesian mean estimation-based method for the simultaneous reconstruction of the unknown state and determination of the overall detection efficiencies.

In reporting the estimates of the various parameters, one has to include error bars that quantify the uncertainty, and thereby make a statement about the quality and quantity of the data, together with available prior information. A common strategy is to deduce error bars for the measurement device parameters and error bars (more precisely, error regions in the state space) for the estimated state separately. For precalibrated schemes, this is a reasonable approach as the device parameters and the estimated state are inferred from different data; in self-calibrating quantum tomography, however, the device and state parameters are estimated from the same data. A proper reporting of uncertainty should hence involve error regions that are regions in the combined device-state parameter space, rather than separate error bars.

In this paper, we discuss how to deduce proper error bars for self-calibrating quantum tomography schemes. The notion of optimal error regions (OERs) OER13 permits consistent treatment of error regions for device and state parameters within a single, rigorous framework; since the OERs are regions around the multiparameter ML estimator—this is an implied property, not a matter of definition—they naturally justify the use of ML estimators as best guesses. The Bayesian foundation of OERs further provides the smooth incorporation of any prior information about the measurement device and the source of the state. We work out examples on this matter, applying our methods, in particular, to the scenario of satellite-based quantum key distribution (QKD), for which a self-calibrating approach is a necessity. As a point of caution, we also note the additional structure in the likelihood function, in this situation of self-calibration, that can present difficulties for standard ML approaches.

Below, we first begin with an overview on the self-calibration procedure, and explain the basic notion of constructing joint ML estimators and OERs. We then illustrate our approach in Sec. III using the example of the $x$-$z$ measurement scheme of the BB84 protocol for QKD BB84 . We explore in greater detail the specific scenario of space-based QKD experiments in Sec. III.4, and offer concluding remarks in Sec. IV.

II Basic concepts

II.1 Self-calibration scheme

In the typical scenario, as shown in Fig. 1, a finite number of identically prepared quantum systems, the relevant properties of which are described by the unknown quantum state $\rho$, are measured by the measurement apparatus, correspondingly described by a probability-operator measurement (POM), also known as a positive operator-valued measure. The POM consists of $K$ non-negative probability operators $\Pi_{1}(\bm{\alpha})$, $\Pi_{2}(\bm{\alpha})$, …, $\Pi_{K}(\bm{\alpha})$, one for each of the $K$ measured outcomes, and another probability operator $\Pi_{0}(\bm{\alpha})$ for the null event. The elements of the POM resolve the identity, ${\sum_{k=0}^{K}\Pi_{k}(\bm{\alpha})=1}$. The unknown parameters $\bm{\alpha}$ represent one’s lack of knowledge about the measurement devices. Then, the probability for the $k$th outcome to occur is given by the Born rule,

[TABLE]

The data $D$ are made up of the sequence of counts of the detection events, ${D=\{n_{1},n_{2},\dots,n_{K}\}}$, with ${N=\sum_{k=1}^{K}n_{k}}$ total counts of detection events. In scenarios with imperfect detectors, there are also missing counts denoted by $n_{0}$. If we know the total number of copies sent into the apparatus, then $n_{0}$ is known once we know $N$; more typically, one does not have precise information about the total number of copies, and $n_{0}$ is an unknown parameter that enters the model for the experiment.

In quantum state estimation, $\bm{\alpha}$ is precalibrated and the task is to reconstruct $\rho$ from the data. In quantum parameter estimation, $\rho$ is known and the task is to reconstruct $\bm{\alpha}$ from the data. In a self-calibration scheme, one aims at reconstructing the state $\rho$ and the device parameters $\bm{\alpha}$ from the same data.

II.2 Maximum-likelihood estimation

The ML estimators for the state and device parameters—denoted by $\hat{\rho}_{\text{ML}}$ and $\hat{\bm{\alpha}}_{\text{ML}}$, respectively—can be obtained by maximizing the likelihood function over the joint state and device parameter space using an iterative algorithm MLEreview . For the scenarios with missing counts, the likelihood of obtaining the observed data $D$, given the state $\rho$ and the device parameters $\bm{\alpha}$, is

[TABLE]

where $L(D,n_{0}|\rho,\bm{\alpha})$ is the corresponding likelihood of obtaining the data $D$ and also having $n_{0}$ null events.

For our present discussion of general concepts and methods, we are content with the scenario sketched in Fig. 1. We note, however, that there can also be unknown parameters of the source that enter the likelihood function but are not accounted for by the quantum state $\rho$, and source parameters of this kind can be handled analogously to the apparatus parameters $\bm{\alpha}$. An example is provided by the situation discussed in Sec. III.4.

II.3 Optimal error regions

The ML estimators are point estimators which represent our best guess for the unknown quantum state and device parameters. The point estimator calculated from a finite amount of data will not coincide exactly with the true parameters. To be statistically meaningful, a point estimator should be endowed with error bars—error regions in higher dimensions—and these regions should be optimally chosen by appropriate criteria. In the Bayesian approach used here, the OERs can be equivalently specified as having either largest credibility ($\equiv$ posterior content) for the prechosen size ($\equiv$ prior content) or smallest size for the prechosen credibility; see Refs. Evans+2:06 ; OER13 . The OERs can be characterized very easily, as they are bounded-likelihood regions—regions containing all the points with their likelihood larger than or equal to some threshold value,

[TABLE]

with $0\leq\lambda\leq 1$, where $L_{\text{max}}(D)=L(D|\hat{\rho}_{\text{ML}},\hat{\bm{\alpha}}_{\text{ML}})$ is the maximum value of the likelihood function. $\mathcal{R}_{0}$ is the same for all data $D$; it contains all thinkable quantum states and device parameter values. Owing to the simple inequality in Eq. (3), it is easy to check whether a particular $(\rho,\bm{\alpha})$ of interest is inside the $\mathcal{R}_{\lambda}$ in question, even if the high-dimensional OERs are difficult to visualize.

For any OER $\mathcal{R}_{\lambda}$, there are two important quantities, its size $s_{\lambda}$ and its credibility $c_{\lambda}$. The size of a region is the assigned probability of finding the true state and device parameters in the region, prior to acquiring the data. For the OERs, then,

[TABLE]

with

[TABLE]

where $q_{1},\dots,q_{J}$ denote the state parameters; the factor $w_{\text{cstr}}(q_{1},\dots,q_{J})$ accounts for all the constraints that the state parameters have to satisfy; $w_{0}(q_{1},\dots,q_{J})$ and $w_{0}(\alpha_{1},\dots,\alpha_{M})$ are the prior densities which represent our knowledge on the state and device parameters.

After the data have been obtained, we update our belief about the state and device parameters by multiplying the prior density with the likelihood function to obtain the posterior density. The credibility of a region is the probability of finding the true state and device parameters in the region, conditioned on the observed data. For the OERs, then,

[TABLE]

where

[TABLE]

is the likelihood of the data $D$. We have ${c_{\lambda}\geq s_{\lambda}}$, with the equal sign usually only holding for ${\lambda=0}$ and $1$, and there is the link

[TABLE]

between $c_{\lambda}$ and $s_{\lambda}$, so that sometimes only the first of the high-dimensional integrals in Eqs. (4), (6), and (7) needs to be computed numerically.

A related concept is the plausible region Evans:15 ; OEI16 which is based on the principle of evidence. The principle of evidence states that if the posterior probability of a $(\rho,\bm{\alpha})$ pair is larger than its prior probability then the data give evidence in favor of the $(\rho,\bm{\alpha})$ pair, and the data give evidence against the pair if the posterior probability is less than the prior probability. The plausible region is composed of all the $(\rho,\bm{\alpha})$ pairs, in favor of which the data give evidence—all the pairs with $L(D|\rho,\bm{\alpha})>L(D)$. Accordingly, the plausible region is the OER for the critical $\lambda$ value

[TABLE]

Once a prior has been chosen, the plausible region is solely determined by the data. As more and more data are obtained, the size of the plausible region decreases while its credibility increases. If the size $s_{\lambda_{\text{crit}}}$ of the plausible region is small and its credibility $c_{\lambda_{\text{crit}}}$ is large, “then this suggests an accurate estimate has been obtained” Evans-essay-2019 . One can quantify the strength of the evidence associated with an OER in terms of a suitably chosen function of the ratio $c_{\lambda}/s_{\lambda}$; see Ref. Evans-Guo-2019 for a discussion of various aspects thereof.

To report the error region for the observed data $D$, following the scheme in Ref. OER13 , $s_{\lambda}$ and $c_{\lambda}$ are plotted as functions of $\lambda$; see Figs. 4 and 6 below. For any desired level of credibility, the corresponding $\lambda$ and the size of the OER can be determined from the graph. Alternatively, if the experimenter chooses to report the plausible region, the critical $\lambda$ value can be calculated and the size and credibility of the plausible region can be read off from the graph. The difference ${c_{\lambda}-s_{\lambda}}$ is largest for ${\lambda=\lambda_{\text{crit}}(D)}$.

To calculate the size and credibility of a region for the observed data $D$, one has to resort to Monte Carlo integration (see Refs. Shang+1:15 ; Seah+1:15 ; YS19 for the application of Monte Carlo integration to this context of OERs). Random samples in the joint state and device parameters space have to be generated to perform the high-dimensional integrals for the size and credibility. We do this by employing the Hamiltonian Monte Carlo strategy described in Ref. Seah+1:15 .

III Self-calibration for a BB84 experiment

To explain our approach, and illustrate the significance of proper error regions in self-calibration schemes, we consider a concrete example: a BB84 QKD experiment BB84 . We are envisioning an experimenter wishing to set up a BB84 QKD protocol, who needs to calibrate the entangled photon source, as well as the detectors to be used in the protocol. Such knowledge is needed for protocols like BB84 that do not have the device-independence properties of more recent QKD schemes; yet, even experimenters setting up device-independent schemes would potentially require such self-calibration experiments to understand their own setup, even if the information is not used in the security analysis. We first set up the problem, and discuss the general approach to constructing proper error regions; it should be emphasized that, while we discuss our approach for this specific BB84 example, our methods apply to other situations. We discuss, as a final subsection, a practical application to space-based QKD, where such self-calibration is a necessity.

III.1 Setup

In one version of the BB84 scheme, a photon source emits entangled photon pairs that follow a Poissonian distribution with mean number $\nu$; see Fig. 2. One of the photons in each pair is sent to a measurement apparatus on the left; the other photon is sent to the apparatus on the right. On each side, the photon is measured by a four-outcome “crosshair” POM with the ideal detection probabilities

[TABLE]

where $x={\left\langle{\sigma_{x}}\right\rangle}$ and $z={\left\langle{\sigma_{z}}\right\rangle}$ are the expectation values of two components of the Pauli vector operator $\bm{\sigma}$. These expressions would apply to the photon detection on the right in Fig. 2 if the detectors had unit efficiency. When accounting for the finite detection efficiencies, we have expressions such as

[TABLE]

for the probability of detecting the photons of a pair by detectors $2$ on the right and $3^{\prime}$ on the left in coincidence. Here, $\eta_{k}$ is the probability that detector $k$ functions correctly, i.e., detects a photon that falls on it. For simplicity, we are here assuming that there is a negligible chance of losing the photon on the way from the source to the detector; this is markedly different in the situation of Sec. III.4.

Owing to the imperfection of the detectors, there are actually five possible outcomes on each side including a null event, thus 25 different joint outcomes in total. However, the double-null events where both photons escape detection are not recorded and, since the actual number of entangled photon pairs is not known, we also do not know how many double-null events have occurred. Therefore, the data $D$ are made up of the sequence of 24 counts of detection events, $D=\{n_{1},n_{2},\dots,n_{24}\}$; see the table in Fig. 2. In this scenario, the likelihood for obtaining the data $D$, given the state $\rho$ and detector efficiencies $\bm{\eta}$, is

[TABLE]

up to an overall factor of no consequence, where $p_{k}$ is the probability of detecting an event of the $k$th kind and

[TABLE]

is the probability of getting a double-null event; see the Appendix for a derivation.

In this example, the only unknown parameters of the measurement devices are the detector efficiencies; other aspects of the measurement devices are assumed to have been precalibrated. So, our task is to infer the eight state parameters—${\left\langle{1\otimes\sigma_{x}}\right\rangle}$, ${\left\langle{1\otimes\sigma_{z}}\right\rangle}$, ${\left\langle{\sigma_{x}\otimes 1}\right\rangle}$, ${\left\langle{\sigma_{z}\otimes 1}\right\rangle}$, ${\left\langle{\sigma_{x}\otimes\sigma_{x}}\right\rangle}$, ${\left\langle{\sigma_{x}\otimes\sigma_{z}}\right\rangle}$, ${\left\langle{\sigma_{z}\otimes\sigma_{x}}\right\rangle}$, and ${\left\langle{\sigma_{z}\otimes\sigma_{z}}\right\rangle}$—and the eight detector efficiencies—$\eta_{1}$, $\eta_{2}$, $\eta_{3}$, $\eta_{4}$, $\eta^{\prime}_{1}$, $\eta^{\prime}_{2}$, $\eta^{\prime}_{3}$, and $\eta^{\prime}_{4}$—from the data, and not only report our best guess for these 16 numbers but also quantify the accuracy of the inferred values. As it turns out, estimating all eight detector efficiencies using the ML approach cannot be done in a straightforward manner: The resulting likelihood function can have multiple maxima, which renders ML estimation ambiguous. We elaborate on this matter in the next section.

To focus our discussion on the issue of proper error regions, rather than resolving this ambiguity in the ML estimation scheme (worthy of further investigation elsewhere), we instead assume a simplification in the form of prior information about the detectors. Specifically, we assume that the detector efficiencies on each side of the setup are in stable, and precalibrated, ratios with one another. In Sec. III.4, we describe a physical scenario for which such prior information about the detectors is natural. Then, what is unknown and to be estimated in the self-calibration scheme, are the maximum detector efficiencies, one for each side, $\eta_{\textrm{right}}=\max\{\eta_{1},\eta_{2},\eta_{3},\eta_{4}\}$ and $\eta_{\textrm{left}}=\max\{\eta^{\prime}_{1},\eta^{\prime}_{2},\eta^{\prime}_{3},\eta^{\prime}_{4}\}$. Once these efficiencies are known, we know the individual detector efficiencies from our precalibrated ratios.

III.2 Multiple local maxima of the likelihood function

When we do not take the prior information on the ratios of detector efficiencies into account and try to reconstruct the eight efficiencies and the eight state parameters with ML estimation, we observe multiple maxima in the likelihood function in some cases, which are neither typical nor rare. These maxima have approximately the same height, but refer to very different situations as they are usually far from each other in the joint parameter space. We illustrate this feature with an example from a simulated experiment, see Fig. 3. The data obtained in the simulation are $D=\{$597, 515, 658, 345, 2012, 1039, 804, 1081, 506, 3308, 1091, 795, 990, 600, 3446, 740, 558, 826, 326, 2334, 4228, 3123, 4295, 2143$\}$. As can be seen, the three maxima are far away from each other. In particular, if we choose two of them, the fidelity between the two states is not high, ${F(\rho_{1},\rho_{3})=\mathrm{tr}{\left\{{\left|\sqrt{\rho_{1}}\sqrt{\rho_{2}}\right|}\right\}}=0.84}$, and the efficiencies are rather distinct, $|\bm{\eta}_{1}-\bm{\eta}_{3}|=0.97$.

The presence of multiple maxima indicates that the log-likelihood function is not concave in the joint space of quantum states and efficiencies. A unique “best guess” based on maximizing the likelihood is then not available, even if there is—mathematically speaking—a single global maximum because local maxima may have parameter values that are equally plausible for a physicist. This is the situation illustrated by the example in Fig. 3, where the self-calibrating scheme for the double-crosshair measurement does not yield an unambiguous answer, and the OERs $\mathcal{R}_{\lambda}$ consist of disjoint pieces when $\lambda\lesssim 1$; the plausible region is also of this kind. While it is true that the global maximum tends to dominate when sufficiently many events are observed, and then the resulting plausible region is convex, the situation can easily be inconclusive for a large, but not that large, number of observations.

Yet, an ambiguity of this kind can often be resolved by taking additional information into account. Rather than using it for choosing among the candidate estimators identified by maximizing the likelihood, it is preferable, if possible, to reduce the number of parameters. As explained earlier, in the example studied here, this reduction is achieved by prior knowledge about the ratios of the four detection efficiencies on each side of the scheme in Fig. 2.

Ratios of the detector efficiencies can be determined much more easily than the absolute values. For example, by counting photons from the same source for the same period of time by each of the detectors, the ratios of their efficiencies are simply the ratios of the photon counts. Therefore, we modify our scheme such that ratios of the detector efficiencies on each side of the setup are measured prior to the BB84 experiment.

When the prior information of knowing ratios of the detector efficiencies is included, we are left with the eight state parameters, plus two efficiencies, i.e., $\eta_{\text{right}}=\max\{\eta_{1},\eta_{2},\eta_{3},\eta_{4}\}$ and $\eta_{\text{left}}=\max\{\eta^{\prime}_{1},\eta^{\prime}_{2},\eta^{\prime}_{3},\eta^{\prime}_{4}\}$ to be estimated from the data. All the other detector efficiencies can then be deduced from the estimates for $\eta_{\text{right}}$ and $\eta_{\text{left}}$.

Simulations provide strong numerical evidence that the likelihood function for the new scheme has only one maximum; it is log-concave indeed. For the example in Fig. 3, a unique estimator $(\hat{\rho}_{\text{ML}}$, $\hat{\bm{\eta}}_{\text{ML}})$ is found after the prior knowledge of ratios of the detector efficiencies is taken into account. Moreover, this estimator is very close to the mock-true state $\rho_{\text{true}}$ and the efficiencies $\bm{\eta}_{\text{true}}$ used for generating the simulated data. Specifically, we have the fidelity $F(\hat{\rho}_{\text{ML}},\rho_{\text{true}})=0.9999$, and $|\hat{\bm{\eta}}_{\text{ML}}-\bm{\eta}_{\text{true}}|=0.0019$.

For the examples in this paper, we assume the specific values of ratios of detector efficiencies can be obtained from precalibration measurements. In the more general situation, rather than prescribing specific values for the ratios of detector efficiencies, one can also remove the multiple local maxima from the problem by prescribing a distribution for the ratios. We first reparametrize $(\rho,\bm{\eta})$ as $(\bm{s},\bm{k})$, where $s=(\rho,\eta_{\text{left}},\eta_{\text{right}})$ as before, and $\bm{k}$ are the ratios of detector efficiencies. We then regard $\bm{k}$ as nuisance parameters to be integrated out, according to a prescribed distribution deduced from precalibration measurements, so that we are again left with only $\bm{s}$ to estimate from our data. Our numerical investigations suggest that, for chosen fixed values of $\bm{k}$, the likelihood function $L(D|\bm{s},\bm{k})$ is concave in the space of $\bm{s}$,

[TABLE]

for $\lambda\in[0,1]$ and $\bm{s}_{1}$ and $\bm{s}_{2}$ are two values of $\bm{s}$. Integrating over the ratios $\bm{k}$ according to a prior distribution $(\mathrm{d}\bm{k})w(\bm{k})$ gives the marginal likelihood, $L(D|\bm{s})\equiv\int(\mathrm{d}\bm{k})w(\bm{k})L(D|\bm{s},\bm{k})$, which, by linearity, inherits this concavity property and has a single maximal point.

III.3 ML estimation and proper error regions

In one simulated experiment for the double-crosshair measurement of the BB84 scheme, we obtained the data $D=\{$1, 2, 1, 1, 7, 1, 1, 0, 0, 2, 3, 3, 1, 1, 7, 1, 0, 1, 3, 1, 5, 12, 3, 9$\}$. In this example, we assume that we do not know anything about the mock-true state, $\eta_{\text{left}}$, and $\eta_{\text{right}}$ before the data were taken. Thus, the prior we choose is uniform in the eight state parameters, and also uniform in $\eta_{\text{left}}$ and $\eta_{\text{right}}$. Figure 4 shows the plot for the size and credibility of the ten-dimensional OERs $\mathcal{R}_{\lambda}$ as a function of $\log_{10}\lambda$, calculated by a Monte Carlo integration that uses a random sample with $500\,000$ points. To ensure the physicality of the states in the Hamiltonian Monte Carlo sampling, we use the parametrization described in Sec. 4.3 of Ref. Seah+1:15 . The size and credibility of the plausible region are ${s=0.0378}$ and ${c=0.9826}$, respectively. The mock-true state and efficiencies used for the simulation are contained in the OERs with ${\lambda<8.27\times 10^{-2}}$ and ${c>0.249}$. Thus, they are in the plausible region.

More specifically, this simulation used detector efficiencies with the following ratios:

[TABLE]

and

[TABLE]

and the parameters of the mock-true state are reported in the middle column of Table 1. This table also shows the parameter values of the ML estimators $\hat{\rho}_{\text{ML}}$ and $\hat{\eta}_{\text{left}}$, $\hat{\eta}_{\text{right}}$. We note that, although the estimators are in the plausible region, the estimated parameter values are rather different from the ones used for the simulation. This is not unexpected for so few data, namely, only 66 detection events.

As an illustration that the error regions are regions in the joint device-state parameter space, Fig. 5 shows the two-dimensional slice of the ten-dimensional plausible region with the other parameters set to their mock-true values. Observe that there is no natural way to break up the error region into a device-parameter-only region, and a state-parameter-only region. A naive approach might be to report the maximum width of a OER of a chosen credibility or plausible region along the $\eta_{\textrm{left}}$ axis as the error bar for that device quantity, and then construct the state-only error region as usual. However, this is not representative of the actual situation where both state and device parameters are inferred from the same data.

The correct procedure for determining an optimal error range for a single parameter is based on the likelihood of observing the data, conditioned solely on the parameter of interest. This inferred likelihood is obtained from $L(D|\rho,\bm{\alpha})$ by marginalizing all other parameters; see Ref. OEI16 . We shall not elaborate on this matter here, except for noting that, usually, the best guess for a singled-out parameter is not equal to the best guess for this parameter when it is estimated together with all others.

III.4 Application to space-based QKD experiments

We turn to a practical scenario where self-calibration is crucial for the QKD experiment to function properly: space-based QKD. One of the main challenges in QKD is to extend the coverage towards a global scale. Optical fibers and free-space links between ground stations have distance limits due to losses in fibers or the need of line-of-sight locations. To establish a global quantum communication network, the usage of satellites as transmitters and receivers has been proposed space , and there have been successful experiments demonstrating the feasibility China-Micius ; China-Micius2 ; China-Austria .

We consider a variation of the previous setup in Fig. 2 applicable to satellite-based QKD experiments. We have the set of detectors on the right-hand side and the photon source located on a satellite, whereas the set of detectors on the left-hand side are located on Earth. Within the satellite, the environment can be stabilized so that the efficiencies of the detectors in the satellite relative to one another are stable over time; the same is true for those located in the Earth laboratory. These efficiency ratios can be precalibrated and recorded ahead of time. What cannot be controlled and can vary to a large degree over time is the relative efficiencies between the detectors on Earth and those in space. In particular, the photons have to pass through the atmosphere before reaching the Earth-bound detectors. The atmosphere acts as an absorber which only allows a rather small fraction of the photons to pass through while leaving the polarization of the photons essentially unchanged, and we can consider this absorption as part of the loss in efficiency of the Earth-based detectors. The fraction $T$ of the photons which passes through the atmosphere depends on a lot of factors such as temperature, humidity, atmospheric turbulence, etc. space2 ; space3 ; space4 . Thus, it is hard to determine the true value of $T$; it also keeps changing over time. When applying the self-calibration procedure, however, the value of $T$ does not need to be determined in advance. Instead, it is treated as an unknown parameter which can be determined from the data. In this scenario, the overall efficiencies of the set of detectors on Earth depend both on their imperfection and on losses in the atmosphere, specified by two parameters $\eta_{\textrm{right}}$ and $T\eta_{\textrm{left}}$, to be estimated in a self-calibration experiment.

In the previous sections, we have assumed that the mean number of photons $\nu$ is known. This assumption is unrealistic in the current scenario as the power of the photon source might deteriorate over time and it might be hard to calibrate it since it is located on the satellite. Thus, we get rid of this assumption in this scenario and treat $\nu$ as another unknown to be determined. The likelihood for this scenario is

[TABLE]

up to an overall factor of no consequence (see the Appendix). The unknowns to be determined are the eight state parameters, the mean number of photons $\nu$, and the two efficiencies, $\eta_{\text{right}}$ and $T\eta_{\text{left}}$. As $T$ is a small number, the number of detections on the satellite side is much much greater than that on the Earth side, which makes this problem highly asymmetric.

We obtained the following data ${D=\{}$0, 1, 0, 1, 2, 1, 2, 0, 0, 0, 1, 0, 0, 4, 5, 2, 2, 0, 3, 3, 65188, 70928, 37230, 127525$\}$ from a simulated experiment. In this example, we assume that we do not know anything about the mock-true state, but we know $\eta_{\text{right}}$, $T\eta_{\text{left}}$, and $\nu$ to a certain extent before the data are taken. Thus, we choose a uniform prior for the eight state parameters and a beta prior for both $\eta_{\text{right}}$ and $T\eta_{\text{left}}$,

[TABLE]

and a gamma prior for $\nu$ ,

[TABLE]

For these choices, the expected values of $\eta_{\text{right}}$, $T\eta_{\text{left}}$, and $\nu$ are $0.778$, $1.87\times 10^{-4}$, and $500\,000$ with standard deviations of $0.049$, $1.53\times 10^{-4}$, and $50\,000$, respectively. The shortest intervals where these values lie with 0.95 prior probabilities are $[0.681,0.870]$, $[1.97\times 10^{-7},4.88\times 10^{-4}]$, and $[403716,599105]$ respectively. This quantifies our prior guesses about these parameters and our trust in the guessed values. The true values of $\eta_{\text{right}}$, $T\eta_{\text{left}}$, and $\nu$ used in the simulation are $0.724$, $7.38\times 10^{-5}$, and $500\,000$, respectively.

Figure 6 shows the plot for the size and credibility of the OERs as a function of $\log_{10}\lambda$, from a Monte Carlo integration that uses a random sample with $500\,000$ points. The size and credibility of the plausible region are $s\simeq 4.43\times 10^{-6}$ and $c\lesssim 1$, respectively; there are so few sample points in the very small plausible region that its size cannot be determined accurately here (the sampling error is indicated in the figure noteSamplingError ). The mock-true state and efficiencies are contained in the OERs with $\lambda<7.73\times 10^{-3}$ and $c>0.7477$. Thus, they are inside the plausible region. The sizes of the OERs shown in Fig. 6 decrease much more rapidly than those in Fig. 4 because the likelihood function is very sharply peaked in certain directions. This is due to the large number of detections on the satellite side. Therefore, the plausible region is a very small region around the maximum likelihood estimator. Figure 7 shows the two-dimensional slice of the 11-dimensional plausible region with the other parameters set to their mock-true values. Again, one observes that such joint device-state error regions are much more informative and representative of the joint error than reporting separate error bars.

The ratios of detector efficiencies in this simulation are

[TABLE]

and

[TABLE]

The parameters of this simulation are reported in Table 2. Owing to the large count of detection events on the satellite side, some of the estimated parameter values are quite close to their mock-true values, while others are not. More importantly, however, the ML estimates are inside the plausible region, which is of very small size.

IV Conclusions

Self-calibration is needed whenever precalibration of the experimental devices is not feasible. The data from the self-calibration experiment yield information about both the state as well as the measurement apparatus. One estimates state and device parameters from the same data. In this paper, we explain how to do this jointly, treating the state and device parameters on equal footing. In particular, we construct state-device optimal error regions, properly quantifying their joint uncertainty. We applied our methods to the example of the crosshair measurement in BB84 QKD experiments, and discussed the case of satellite-based QKD, where self-calibration is a necessity due to the physical situation.

Our methods can, of course, be applied to many other situations that involve parameter estimation. Examples include the absolute calibration of photon detectors in Klyshko’s scheme Klyshko80 ; RarityRidleyTapster87 (with or without the simultaneous calibration of the photon-pair source), the determination of an interferometer phase QuPhase1 ; QuPhase2 , and the precise estimation of the down-conversion efficiency in experiments with entangled photons arXiv:1808.06863 .

Acknowledgements.

This work is funded by the Singapore Ministry of Education (partly through Tier-2 Grant No. MOE2016-T2-1-130) and the National Research Foundation of Singapore. J.S. acknowledges support by the Beijing Institute of Technology Research Fund Program for Young Scholars, the National Natural Science Foundation of China through Grant No. 11805010, the European Research Council (Consolidator Grant No. 683107/TempoQ), and the Deutsche Forschungsgemeinschaft.

Appendix: Derivation of the likelihoood functions in Eqs. (16) and  (21)

In the lapse of time, during which we detect $N=n_{1}+n_{2}+\cdots+n_{24}$ events, there are also $n_{0}$ double-null events that are not recorded, and we do not know when they occur. For a given sequence of $N$ detected events, interspersed by $n_{0}$ unrecorded events, there are $\frac{(N+n_{0})!}{N!\,n_{0}!}$ possible sequences of ${N+n_{0}}$ events. Then, the likelihood for observing the actual sequence $S$ of detected events, given the quantum state $\rho$, the detection efficiencies $\bm{\eta}$, and the count $n_{0}$ of unrecorded double-null events, is

[TABLE]

with the $p_{k}$s related to $\rho$ and $\bm{\eta}$ through Born’s rule, as in Eq. (15).

The number ${N+n_{0}}$ of down-converted photon pairs follows a Poissonian statistic with an average number of $\nu$ pairs during the period of data taking. Accordingly, the likelihood for $S$, conditioned on $\rho$, $\bm{\eta}$, and $\nu$, is

[TABLE]

The actual sequence $S$ does not matter here, as only the event counts $D=(n_{1},n_{2},\dots,n_{24})$ enter. The likelihood $L(D|\rho,\bm{\eta},\nu)$ for this minimal statistic differs from $L(S|\rho,\bm{\eta},\nu)$ by the combinatorial factor $\frac{N!}{n_{1}!\,n_{2}!\,\cdots\,n_{24}!}$, so that

[TABLE]

This is the product of 24 independent Poisson distributions, one for each kind of detection event, with the average count $\nu p_{k}$ for the event of the $k$th kind. After removing the $p_{k}$-independent factors, which cancel in Eq. (6) and therefore have no bearing on the error regions of Sec. III, we arrive at Eq. (16). For the scenario in Sec. III.4, where $\nu$ is also treated as unknown, we have Eq. (21) after removing only the factorial factors.