Recursive quantum detector tomography

Conventional tomographic techniques are becoming increasingly infeasible for reconstructing the operators of quantum devices of growing sophistication. We describe a novel tomographic procedure using coherent states which begins by reconstructing the diagonals of the operator, and then each successive off-diagonal in a recursive manner. Each recursion is considerably more efficient than reconstructing the operator in its entirety, and each successive recursion involves fewer parameters. We apply our technique to reconstruct the positive-operator-valued measure (POVM) corresponding to a recently developed coherent optical detector with phase sensitivity and number resolution. We discuss the effect of various parameters on the reconstruction accuracy. The results show the efficiency of the method and its robustness to experimental noise.


Introduction
Quantum detectors inform our classical world of the underlying quantum world through a set of operators known as positive-operator-valued measure (POVM). In practice, the success of many quantum applications rely on certain knowledge of measurement apparatuses [1,2,3]. Successful applications of sophisticated detectors rely on a complete and accurate knowledge of the detector, i.e. detector characterization. Detector characterization can be implemented in two different ways. One is synthetic, wherein each constituent of a detector is carefully calibrated before being incorporated into a sophisticated physical model of the measurement process. As quantum technologies evolve into increasingly complicated systems, so do quantum detectors, which makes synthetic characterization progressively less feasible. Additionally, any coupling with external degrees of freedom not incorporated into the theoretical model may make the characterization fail [4,5]. A fundamentally different approach is taken in quantum detector tomography (QDT) [6,7,8,9], where the unknown specifics of a detector are characterized in a largely assumption-free way: here, the POVM of a detector are reconstructed from the outcome statistics in response to a set of tomographically complete certified input states.
To date QDT has been successfully applied to avalanche photodiode (APD) [4], timemultiplexed photon-number-resolving detector (TMD) [9,10,11], transition edge sensors [12] and superconducting nanowire detectors [13]. These detectors are phase-insensitive, i.e. they can only measure the mixture of the photon-number states, not the coherence among them. Accordingly, the POVM of these detectors have only non-zero matrix elements on the main diagonals, and the number of parameters to be decided is proportional to d. Here d is the dimensionality of the Hilbert space, and for optical detectors can be estimated as the number of photons to saturate the detector. Yet a large number of effects characteristic for quantum mechanics, including entanglement, violation of local realism [14], measuring non-classical correlations of radiation fields [15], test of macroscopic realism [16] etc., relies on quantum coherence. The effort to harness and exploit quantum coherence brings the prosperity of quantum information processing and quantum metrology. Moreover, exploration and utilization of the full Hilbert space of a quantum system requires a detector capable of implementing a tomographically complete set of measurements [17]. Therefore optical detectors that can access quantum coherence among photon-number states, i.e. phasesensitive detectors, for example strong-and weak-field homodyne detectors [18], are crucial not only for quantum applications, but also for test of fundamental theories of quantum mechanics. Phase-sensitivity comes with the non-zero off-diagonal matrix elements of the POVM. Thus tomography of a phase-sensitive detector requires the estimation of number of parameters proportional to d 2 . For practical detectors, d can range from 10 2 to 10 5 , and d 2 from 10 4 to 10 10 . For example, a weak-field homodyne TMD with 9 time bins requires 1.8 × 10 6 parameters to completely describe its POVM [5], which is about two orders of higher than the largest tomography that had been performed until then [19]. Such large set of parameters represents a considerable challenge to the characterization of phase-sensitive detectors.
In this work we explore potential solutions to QDT of phase-sensitive quantum detectors. In particular we introduce an algorithm that allows to reconstruct the POVM recursively, with no more than d parameters per recursion. Simulations with the QDT of weak-field detectors demonstrate the robustness of this algorithm.

Definition of the problem
QDT is performed by preparing a set of known probe states {ρ m } incident on a quantum detector and observing the detector outcomes. The probability of registering outcome n is given by the Born rule where {Π n } is the POVM of the detector with n = 0, . . . , N − 1, and N is the number of possible outcomes of the detector. In practice the experiment is repeated for each of many identical copies of the probe states, and the frequency f n|m for each measurement outcome n occurring when probe state ρ m is used is recorded. Then p n|m can be estimated from the relative frequency p n|m = f n|m / n f n|m . One can then invert Eq. (1) to find Π n . For a finite number of repetitions, there are always fluctuations in the estimation of p n|m , therefore the inversion should normally be preformed with convex optimization.
A key requirement is that the set of probe states must be tomographically complete. However, it is also important that the set of probe states are experimentally feasible. That means that the states should themselves be well-characterized, and that a large variety should be available with high precision. There are proposed methods to generate the probe states through quantum correlations [8,20]. Yet with current quantum optical sources it is very hard to generate the probe states strong enough to saturate the detector under test. For photodetector measurements, there is a more straightforward option. The set of coherent state vectors |α of an optical beam are ideal candidates, where α = |α|e iθ is the complex amplitude. They are overcomplete in the sense that two different coherent states are not orthogonal with each other yet any quantum state can be decomposed on the set of coherent states. Therefore coherent states can form a tomographically complete set by transforming their amplitude (by means of optical attenuation) and their phase (with a simple delay line). Importantly, they are generated very easily by a laser.
With coherent states as input, the probabilities are given by where Q n (.) is the Q-function of the detector POVM elements Π n . This is equal to the Husimi representation of the POVM, and is uniquely and invertibly related to the POVM.
To reconstruct Π n one can write both |α and Π n in the photon-number basis and truncate the expansion at d − 1, where d − 1 is the number of photons that saturate the detector.
Then Eq. (2) can be written as We can relabel Eq. (5) where P is an M × N matrix with elements P m,n = p n|αm , F is an M × d 2 matrix with elements andΠ is a d 2 × N matrix with elementsπ s,n = π j(s),k(s) n . In practice where the experimental noise is taken into account, the POVM set can be estimated from Eq. (6) with convex optimization subject to the constraints where I is the identity operator. One common approach is the least square estimation where ||A|| 2 = Tr(A † A) is the Frobenius norm. The reconstruction problem effectively deconvolves a coherent state from the statistics to obtain the POVM set. This is an illconditioned problem, as seen by the large ratio between the largest and smallest singular values of the matrix F . This makes the POVM extremely vulnerable to small fluctuations in the measurement statistics. Such instability can be overcome by adding a regularization function g(Π) to the optimization [9,10], therefore the problem is modified as For a phase-insensitive detector with finite detection efficiency, one would expect the variation of the diagonal matrix elements to be smooth, therefore a regularization function known as Tikhonov regularization [21] is applied This limits the variation between adjacent elements along the diagonal matrix elements. Yet for a phase-sensitive detector a regular function is not easy to find: even as each of the leading diagonals are smooth, the relation among different leading diagonals can be arbitrary. An alternative approach for convex optimization is maximum likelihood estimation, which was also proposed for QDT [7]. Maximum likelihood alleviates the requirement of the regularization function. However, its convergence speed is normally not high. Moreover, both the maximum likelihood estimation and the least square estimation in Eq. (11) requires the reconstruction of the whole POVM matrices at the same time. When the size of the matrices becomes large, the problem becomes infeasible. For example, the estimation of a POVM set with 9 elements each of which is a 50 by 50 matrix is already a hard problem for the capability of current multi-processor desktops (2xQuad Core 3GHz, 8GB RAM).
The engineering of large entangled quantum states and development of sophisticated quantum operations has set a challenge for standard quantum tomography techniques. There has been increased interest in the development of novel algorithm with improved efficiency for special situations. In particular, there are process tomography schemes that allows to selectively reconstruct the state or process matrix partially in each run. Several of them use apriori knowledge about the state such as their symmetry [22,23,24], or simply reconstructing a subset of the entire state or process [25,26]. Using improved techniques from classical signal processing have also become common, such as compressed sensing [27,28]. In the following, we introduce a novel algorithm that reconstructs the POVM elements recursively in multiple runs.

The detector model
The algorithm discussed in this work can be universally applied to the tomography of any quantum detector. To illustrate its working in this work we consider a simple example of phase-sensitive detector, the weak-field homodyne APD. The configuration of such detectors together with a schematic tomography setup is given in Fig. (1). The probe state is prepared by the phase modulation and amplitude modulation of the output of a laser system. The weakfield homodyne detector is shown in the black box, where the input state interferes with a local oscillator (LO) and detected by a photon-number-resolving or photon-counting detector. For a weak-field homodyne APD, there are two possible measurement outcomes, no-click and click events, and the corresponding POVM elements Π 0 and Π 1 are given as [18].
where α L is the complex amplitude of the LO and η APD is the detection efficiency of APD.

A selective algorithm with Glauber-Sudarshan P -function
Before we proceed to the recursive algorithm, we consider a more straightforward selective algorithm. Each matrix element π j,k n of Π n is given by Using the Glauber-Sudarshan decomposition of |k j| we have π j,k n = 2 P j,k (α) α|Π n |α In principle we can estimate each individual matrix element π j,k n separately with either the exact form of P j,k (.) which contains the derivative of Dirac-delta function or the approximated regular form [29]. Similar method has been used for quantum process tomography [30,31,32]. However, due to the singularity of P -function, this scheme is extremely sensitive to the noise in the measured Q-function of the POVM element, rendering it infeasible for practical QDT.
As an example, we consider the reconstruction of no-click POVM of a weak-homodyne APD with the reflectivity of the beam-splitter of 0.5, LO intensity |α LO | 2 = 5 and quantum efficiency of the APD 60% (overall detection efficiency 30%). We choose the probe states that sample the phase space from X, Y = −10 to X, Y = 10 with a step size 0.05, where X and Y are the two quadratures of an optical field. We assume that for each probe state we run the experiment f times, then the expected frequency to get the no-click event is f 0|α = πQ 0 (α)f . In practice there are many experimental imperfections that may induce fluctuations to the measurement results. In this work we only consider the most fundamental fluctuation due to the random nature of the outcome of each measurement process, and simulate it by assuming that f 0|α is a random number with a binomial distribution, and assigning the experimentally measured Q-function as Q exp (α) = f 0|α /f . The results are shown in Fig. (2). Without experimental fluctuations, the reconstructed POVM matches almost perfectly with the theoretical prediction. Yet when there is the presence of experimental noise, the results deviates from the theoretical prediction very quickly: for f = 10 5 as in Fig. (2(c)), the reconstructed POVM element is not even a physical measurement operator. Only when the number of measurement is large enough f = 10 10 and thus the experimental fluctuations is small, the reconstructed POVM element is close to the real one. The results presented here is reconstructed up to the 4-photon component. The P-functions of higher photon-components are more singular, and the reconstruction is more sensitive to experimental fluctuations. Therefore, this method has a serious problem for its scalability. For practical QDT, we need to seek another solution.

Outline of the recursive reconstruction
In this section, we discuss a novel recursive method for the tomographic reconstruction of quantum operators. We begin with Eq. (5). Before relabeling, we first integrate over the probe state phase θ. With we have 1 2π The left side of Eq. (19) is a partial integration of the experimental results, while the right side involves only the main diagonals of the POVM. Eq. (19) can be interpreted as using phaserandomized coherent states as input to the detector. Since the input states are completely mixed, the measurement process only involve the main diagonals of the POVM. In a practical For reconstructing the off-diagonals we first note that POVM elements are Hermitian and it is sufficient to reconstruct just the upper or lower off-diagonals. We multiply Eq. (5) by e −ilθ and integrate over θ. Since we have 1 2π Again for practical experiment the integration should be substituted with a summation with an error Eq. (22) Since this is a convex optimization problem, there is only one minimum value which can be found with YALMIP toolbox for Matlab [33] with the solver SeDuMi [34] utilizing primaldual interior point methods [21]. For the examples discussed in this paper, the calculation converges to its minimum in less than 40 iterations which takes about 1 seconds on a multiprocess desktop (Dual Core 2GHz, 2GB RAM). This allows us to reconstruct the POVM recursively for l = 0, . . . , d. For l = 0, the second condition is that the summation of the main diagonals of all the POVM elements equals to 1, while for l = 0, this condition means that the summation of the lth leading diagonals of all the POVM elements equals 0. The positivity condition should be enforced recursively based on Sylvester's criterion, which states that a matrix is positive if and only if all of its principal minors are positive. For l = 0 this requires all the matrix elements on the main diagonals to be positive. Now we derive the condition for l > 0. We start with l = 1. In Eq. (26) we show the matrix Π n where the diagonal elements (green) have been determined using Eq. (19) and the first row of off-diagonals (red) is to be determined with Eq. (23). Any other entry is unknown and will not be reconstructed at this step.
For an input state vector of the form |Φ = a |j + b |j + 1 the effective submatrix of Π n is given by which needs to be positive. We marked these submatrices with blue, orange, and black lines for the j = 0, j = 1, and j = M − 1 cases in Eq. (26). Thus we imposed the additional constraint that all diagonally centered 2 × 2 submatrices of Π j,1 n need to be positive for the reconstruction of the first off-diagonal. This condition is satisfied if and only if the determinant det(Π j,1 n ) is positive, which implies π j,j+1 n 2 ≤ π j,j n π j+1,j+1 n ∀n, j.
For the following reconstruction of the lth leading diagonal we impose a similar constraint on the (1 + l) × (1 + l) submatrices start with π j,j n , illustrated in Eq. (29) where only π j,j+l n is unknown. It is required that Π j,l n is a positive matrix. Since constraints in the previous steps ensure all its principle minors are positive, this condition is equivalent to that its determinant det(Π j,l n ) is positive. To facilitate numerical calculations we derive the bounds on π j,j+l n , which can be done by noticing that det(Π j,l n ) is a quadratic polynomial of π j,j+l It is easy to see that A is positive, since A is the product of the elements along the anti-diagonal and Π j,l n is Hermitian. Therefore det(Π j,l n ) ≥ 0 implies that The value of A, B and C can be easily estimated from Eq. (30) by substituting π j,j+l n = ±1, 0 into Π j,l n and calculating the determinant numerically.

The number of leading diagonals l
To reconstruct the full POVM matrices, we should run the calculation in Eq. (25) until l = d − 1. As can be seen from Eq. (24), for higher l it requires increased number of phases M p to reduce the numerical error. Yet, in practice the number of leading diagonals can be estimated during the calculation. From the positivity condition, one has |π j,j+l n | 2 ≤ π j,j n π j+l,j+l n .
Therefore, after the reconstruction of the principle diagonals, we can put a bound on the number of leading diagonals to be reconstructed. Moreover, in any practical detector there is always a finite fluctuation of the reference phase (with a fluctuation of 2π for a phaseinsensitive detector), which will further reduce the number of leading diagonals, as shown below. In fact, this phase noise will ensure that the entries of the POVM elements decay exponentially away from their main diagonal. Assume the reference phase has a Gaussian distribution with width of δ > 0. Instead of having a POVM Π n we have where R(ξ) is the rotation operator in phase space with angle ξ. The matrix elements of Π n are given by Intuitively, if the fluctuation of the phase reference is small, i.e., δ π, the last integration in Eq. (34) can be approximated as The intuition of exponentially decaying coefficients can be made rigorous as follows. One has, for w.l.o.g. l even, π −π dξ exp(−ξ 2 /(2δ 2 )) cos(lξ) ≤ δ √ 2π exp(−l 2 δ 2 /2) Thus, the matrix elements of Π n satisfy |π j,j+l n | ≤ 2|π j,j+l n | exp(−l 2 δ 2 /2).
The lth leading diagonal is decreased by a factor of 2 exp(−l 2 δ 2 /2). With the increase of l this factor increases therefore reduces the number of significant leading diagonals in Π n , leading to l d. That is to say, the effort of reconstruction up to a constant error is of order O(d) instead of O(d 2 ). For example, with a phase fluctuation of 10 degrees the 18th leading diagonal is reduced to 1% of that with no LO phase fluctuation.
Another reason for the reduction of the required calculation for the leading diagonals comes from one of the major point of performing detector tomography: to predict the response of the detector with various input quantum states. For situations involving input states with a fixed photon number N , like N 00N states [35] or Holland-Burnett states [2], we only require N leading diagonals of the POVM elements to predict all measurement outcomes. Due to the lack of bright quantum sources, N is usually small (less than 8).

Regularization
The numerical stability of a reconstruction algorithm is one of its vital certificates. Numerical instability has been a common problem in tomography [36,37], particularly so in using phase space data from homodyne tomography to reconstruct operators in the Fock spaces [38]. Tools such as pattern functions [39,40,41] exist that can bridge this gap. They are however, hard to identify and cumbersome to work with [42]. The use of maximum likelihood functions has also been suggested for detector tomography [7,8]. Unfortunately, as mentioned earlier, the speed of the convergence of such algorithms is not generally guaranteed to be high, becoming exponentially slow for certain problems.
We strike a balance by developing a recursive algorithm that is efficient by virtue of being cast as a semi-definite programme, as is evident from the convex function to be minimized, and the linear constraints in Eq. (11). Unfortunately, this still leaves us with an ill-conditioned problem, primarily due to extremely large ratio between the largest and the smallest singular values of the matrix F (l) . This is a consequence of the large range of coherent state amplitudes needed to cover the entire dynamical range of the detector in the Fock space. The most common outcome of this ill-conditioning is to result in reconstructed POVMs that have sharp discontinuities [10]. As shown in Eq. (12) this can be resolved by a smoothing function or Tikhonov regularization [21]. We will next discuss how this mathematical technique is physically enforced in realistic detectors.
Most realistic optical detectors have finite efficiencies which enforces a certain degree of smoothness in their corresponding POVM representations. If a lossy optical detector has a POVM element with non-zero amplitude |m n| it will also have a non-zero amplitude in |m + 1 n + 1|, |m + 2 n + 2|, . . . , |m + K n + K|, decreasing with K. If the detector has a finite efficiency η, it will impose some smoothness on the distribution π j,k n . That is because if G(k) is the probability of registering k photons and H(k ) is the probability that k were present, then the loss process will impose This motivates an immediate generalization of Eq. (12) to that in Eq. (25) as While γ is a free parameter introduced into the problem for numerical smoothness, we show that the outcomes of our reconstruction procedure are fairly insensitive to the actual value of the parameter. Fig. (3) presents the effect of the regularization condition for the reconstruction of the no-click POVM of a weak-homodyne APD with the reflectivity of the beam-splitter of 0.5, LO intensity |α LO | 2 = 5 and quantum efficiency of the APD 60% (overall detection efficieny 30%). We vary the weight of the regularization condition for two orders of magnitude. In addition to the fidelity, we also calculate the relative error of the reconstructed POVM ||Π rec 0 − Π the 0 || 2 /||Π the 0 || 2 . The results are presented in Table (1), which show that the change of the reconstructed POVM elements due to the change of regularization strength is small. This confirms that the main effect of the regularization condition is to suppress the ill-conditioning and noise while leaving the POVM fitting unaffected. Table 1. Sensitivity of the reconstruction procedure to the choice of parameter γ. No-click event of a weak-homodyne APD with the reflectivity of the beam-splitter of 0.5, LO intensity |α LO | 2 = 5 and quantum efficiency of the APD 60% (overall detection efficieny 30%). As a comparison, we also calculated the reconstruction of the no-click POVM of a weakhomodyne APD with the reflectivity of the beam-splitter of 0.5, LO intensity |α LO | 2 = 5 and quantum efficiency of the APD 20% (overall detection efficieny 10%) and that of a weakhomodyne APD with the reflectivity of the beam-splitter of 0.1, LO intensity |α LO | 2 = 5 and quantum efficiency of the APD 90% (overall detection efficieny 81%). The results are shown in Figs. (4) and (5). Calculated fidelities and relative errors are presented in Tables (2) and (3). We can see that regularization works very well for moderate and low detection efficiencies, while its performance decreases if the detection efficiency is very high since the corresponding POVM elements are not smooth any more. On the other hand, one can infer the detection efficiency from the differences between the reconstructed results with different regularization strengths. If such difference is large, one should utilize a reduced regularization strength in the reconstruction. Table 2. Sensitivity of the reconstruction procedure to the choice of parameter γ. No-click event of a weak-homodyne APD with the reflectivity of the beam-splitter of 0.5, LO intensity |α LO | 2 = 5 and quantum efficiency of the APD 20% (overall detection efficieny 10%).  Table 3. Sensitivity of the reconstruction procedure to the choice of parameter γ. No-click event of a weak-homodyne APD with the reflectivity of the beam-splitter of 0.1, LO intensity |α LO | 2 = 5 and quantum efficiency of the APD 90% (overall detection efficieny 81%).

Reconstruction of the POVM of a weak-field homodyne APD
To discuss the performance of the recursive reconstruction method, we numerically simulate the reconstruction of the POVM set of a weak-homodyne APD with the reflectivity of the beam-splitter of 0.5, LO intensity |α LO | 2 = 5 and quantum efficiency of the APD 60%. We choose the intensity of the probe state |α u | 2 from 0 to 100 photons with a step size of 0.5 photon. This is sufficient to saturate the detector response. For each intensity we consider probe phases distributed uniformly between 0 and 2π, i.e. θ u,v = {0, 2π/M p , . . . , 2(M p − 1)π/M p }. Again we only consider the fluctuation induced by the random nature of the measurement process. We assume that for each probe state we run the experiment f times and simulate the experimental noise by assume f 0|α is a random number with a binomial distribution. In Fig. (6) which are 87.04%, 98.19% and 98.32% for M p = 5, 20 and 40 respectively. The change in fidelity can be further elucidated by the red bars on top of the reconstructed POVM element which indicate the distance from the theoretical prediction. From the results we can see that although all three phase settings give almost the same results for the principle diagonal, for higher l it requires more probe phases for an accurate reconstruction. This is due to the numerical error for the calculation of the integral given in Eq. (24). This on the other hand shows a practical advantage of the recursive QDT. The probe phase setting can be decided by the elements in the POVM matrices to be reconstructed. If we are only interested in the low leading diagonals, we can greatly reduce the number of probe phases from that needed for a complete reconstruction of the POVM. In Fig. (7) we show the performance of the recursive QDT under different level of experimental fluctuations. Here for each probe intensity we adjust M p = 40 phases. The method discussed in Sec. 4 requires f = 10 10 to achieve a satisfactory accuracy. As a comparison, the recursive QDT is very robust against the experimental fluctuations: a decent accuracy can already be achieved for f = 10 3 (fidelity with the theoretical prediction 98.27%), with further improvement for f = 10 5 (fidelity 98.32%). Depending on the repetition rate of the detector and the laser system for LO and probe state, this requires only several millisecond to one second for each probe state.

Conclusion
Phase-sensitive quantum-optical detectors are crucial to fully exploit the fundamental features of quantum physics and to optimally utilize optical telecommunications channels [43,44,45]. The success of these applications relies on the accurate knowledge of detectors. Yet as quantum-optical detectors become more sophisticated, normal parameters like detectivity, spectral sensitivity and noise-equivalent power are not sufficient to provide a complete specification of the detector. Moreover the complex structures of detectors and the coupling with external degrees of freedom make the conventional characterization of these detectors less feasible. Quantum detector tomography, a black-box or device-independent approach for the complete characterization of quantum detectors, provides a universal solution to this problem. Full characterization enables more flexible design and use of detectors, be they noisy, nonlinear, inefficient or operating outside their normal range. However, the large number of parameters associated with the tomography of coherent quantum detectors presents a technical challenge. This challenge is becoming increasingly typical as quantum devices grow in sophistication. In this work we present a novel recursive reconstruction algorithm to overcome this problem. Aided by numerical simulations, we have demonstrated successful reconstructions of the POVM of a weak-field homodyne APD. The results show the flexibility of the algorithm and its robustness to experimental noise. The capability to fully characterize coherent quantum-optical detectors paves the way to study genuine quantum features, including wave-particle duality, super-sensitivity etc., of a measurement process. It allows the benchmarking of the performance of quantum-optical detectors for various quantum applications and sheds new light on the assessment and verification of more complex detectors. We also hope that recursive quantum tomography provides an efficient procedure for quantum tomography in other quantum state and process characterization problems.