Quantum tomography of light states by photon-number-resolving detectors

We address state reconstruction by photon-number-resolving detectors, and demonstrate that they may be effectively exploited to perform quantum tomography of states of light. In particular, we find that the pattern function technique, originally developed for optical homodyne tomography, may be also applied to discrete data. Our results open new perspectives for quantum-state reconstruction in the mesoscopic regime, and pave the way to the use of photon-number-resolving-based detection schemes in Quantum Information science.


Introduction
A quantum tomographic technique is a scheme to evaluate the expectation value of any observables as well as to reconstruct the density matrix or the Wigner function of a physical system. To this aim, the technique processes the outcomes of a quorum of observables, measured on repeated preparations of the state under investigation [1][2][3][4]. In the continuous-variable domain, quantum states are usually reconstructed by means of optical homodyne tomography (OHT) [5]. From the optical point of view, OHT is based on an interferometric scheme (see the left panel of figure 1) in which a state ρ (the signal) is mixed at a balanced beam splitter (BS) with a highintensity coherent state b b ñ = ñ f | || |e i , usually referred to as the local oscillator (LO). The two outputs of the interferometer are then detected by two p-i-n photodiodes, whose difference photocurrent is formed, suitably amplified, rescaled by the LO amplitude b | |, and finally recorded as a function of the LO phase f. The values of any observable are then obtained by properly processing the data [6], some of them by averaging pattern functions [7,8] or using maximum-likelihood reconstruction techniques [9][10][11][12].
In order to avoid the use of a strong LO and the subtleties of these reconstruction protocols [13] (e.g. due to binning), in the last decade some efforts have been devoted to the development of hybrid detection schemes, involving photon counting [14][15][16][17][18]. In this paper, we focus our attention on photon-number-resolving detectors [19][20][21]. The implemented scheme bears some resemblance to homodyne detection and leads to similar results. However, at variance with the standard scheme, the p-i-n photodiodes are replaced by photonnumber-resolving detectors and the LO bñ | is a low-intensity (few tens of photons) coherent state. The output of the apparatus is thus the difference between the effective number of detected photons [22][23][24], rather than the difference between two macroscopic photocurrents. Note that the use of photon-number-resolving detectors instead of p-i-n photodiodes gives direct access to the statistics of light at each BS output. This gives the possibility to easily investigate some properties of the state under exam, such as the mean value and the variance of the photon number, and, at the same time, it makes it possible the self-consistent calibration of the employed photon-number-resolving detectors [25]. In addition, the homodyne-like scheme is based on a low-intensity LO, thus demonstrating that, even in the presence of a small imbalance between signal and LO, the state reconstruction is guaranteed. This result allows a deep insight into the physics of quantum tomography by investigating, for instance, the effect of the LO intensity on the state reconstruction. We note that this kind of study is not possible with the standard homodyne scheme, since p-i-n photodiodes are macroscopic detectors and cannot be used to count few photons. At the same time, the exploitation of a detection scheme based on a weak LO could be crucial in the case in which such a field is not directly available from the laser source and must be obtained by means of nonlinear interactions. The only limitation of the homodyne-like scheme is given by the lower quantum efficiency of the photon-number-resolving detectors compared to p-i-n photodiodes. Currently, only superconducting detectors, such as transition-edge sensors and nanowire detectors, can reach a detection efficiency larger than 80%.
Up to now, our scheme has been exploited for state-discrimination [26,27], and proved useful to outperform standard homodyne detection in quantum key distribution with continuous variables [28]. In addition, the same scheme (in the absence of the LO) can be employed to generate negative photon-number correlations, thus demonstrating the analogous of antibunching at the mesoscopic intensity level [29]. We also notice that the scheme adopted in this work is completely different from the so-called weak-homodyne detection ones, in which photon-number-resolving detectors have been employed. Indeed, in those schemes the reconstruction of the states has been achieved by considering only a single output of the interferometer [30][31][32][33][34][35][36][37], whereas here we deal with the difference between the two outputs.
In this paper, we address quantum state reconstruction by the scheme of figure 1 and prove that it may be efficiently employed to retrieve the density matrix of single-mode quantum states, as well as to evaluate the expectation values of the first moments of the quadrature operator. The homodyne-like scheme can thus be considered as a valid alternative to the standard one. Since the reconstruction method is applied to detected quantities, namely to the difference between the number of photons detected at the outputs of the interferometer, the retrieved density matrices correspond to the detected states and not to the ones at the input. We also investigate the role played by quantum efficiency, which is a crucial parameter in the detection of nonclassical features of quantum states [38]. In all cases examined in the paper, robust and reliable results have been obtained with a relatively modest imbalance between signal and LO.
The paper is structured as follows. In section 2 we describe our reconstruction scheme and compare its features with those of standard OHT. Section 3 reports the experimental results, whereas section 4 is devoted to Monte Carlo simulated experiments. Section 5 closes the paper with some final remarks.

State reconstruction by photon-number-resolving detectors
As a matter of fact, information coming from photon-number-resolving-based homodyne-like scheme is of different nature compared to that obtained from homodyne detection. In principle, this would require the design of a new tomographic procedure, in order to reconstruct the quantum state under investigation from experimental data. Nevertheless, the parametric approximation (i.e. treating the LO as c-number, rather than as an operator) used to describe homodyne detection is known to hold also when the intensity of the LO is low [39], as far as its coherence is preserved. In turn, this fact suggests that the pattern functions technique used in OHT The second-harmonic pulses of a ps-laser amplified at 500 Hz are sent to a Mach-Zender interferometer (consisting in two balanced BS, two right-angle prisms, and two high-reflective mirrors (M)). The light pulses in one arm correspond to the signal under investigation, whereas those in the other arm correspond to the LO. Note that in both arms neutral density filter wheels (ND) are used for the fine tuning of the intensities. The length of one arm can be micrometrically changed by means of a piezoelectric movement (Pz). At each output of the interferometer the light is focused by an achromatic doublet (L) into a multi-mode fiber (MF) and delivered to a hybrid photodetector (HPD), whose output is amplified, synchronously integrated and acquired. The photon-number difference is post-processed offline. [1-5, 7, 8] may be exploited also in our scheme. One of the main outcomes of this paper, substantiated by experimental and numerical results, is indeed the proof of this equivalence.
We assume that the state of the input modeâ, with = [ˆˆ] † a a , 1 , is described by the density matrix ρ. In OHT, given the homodyne probability distributions , where ñ f |x are the eigenstates of the quadrature operator = + , the elements of the density matrix ρ in the photon number basis can be reconstructed as ò ò The functions F nm (x, f) are a set of sampling functions and can be written as [7] q where the pattern functions ( ) f x nm can be expressed in terms of regular and irregular wave functions (ψ n (x) and f m (x), respectively) with  m n. Quite simple expressions for the numerical computation of these functions can be found in the original paper [7] and are not reported here explicitly.
Usually, the homodyne probability distributions p(x, f) are retrieved from the difference of the two macroscopic photocurrents exiting the homodyne detector [40][41][42]. In our scheme, we replace p( and Δ is the (phase-dependent) difference between the number of photons detected at the two outputs of the homodyne detector [27]. Before going to the derivation of p D (Δ f ), we notice that, while the outcome x of standard homodyne detection can assume any real value, the quantity Δ f is intrinsically discrete. However, this discretization does not represent a limitation for state reconstruction by the sampling function F nm (x, f). In particular, we notice that our method provides density matrices that are positive semidefinite within their confidence intervals. Indeed, as reported in [7], the reconstruction method also provides confidence regions evaluated from the variance of the density matrix elements by assuming that the homodyne detection is a Poissonian counting process. However, such a calculation may require large computational resources and sometimes it is better to estimate the errors by repeating the tomographic reconstruction on different experimental runs and performing the standard statistical analysis. Of course, the two methods lead to compatible results. In this work, we adopt the latter solution.
In order to obtain p D (Δ f ), let us remind that the information coming from our scheme is contained in the joint count statistics q(n, m) of the number of photons detected at each laser shot at the outputs of the BS. To find the expression of q(n, m), it is useful to consider the Glauber-Sudarshan P representation of the input state ρ is the P-function associated with ρ. After the interference of the input state with the LO (the coherent state bñ | ) at the balanced BS, the two-mode state and the corresponding joint photon-number statistics can be written as [43,44]: 2 . The effect of non unit quantum efficiency η k of the detector on the output modes k=c, d may be taken into account by replacing μ k (ζ, β) with η k μ k (ζ, β).
Starting from the joint probability (6), it is possible to calculate the distribution of the quantity D =n m, namely and the corresponding 'rescaled' version In order to investigate the performance of our reconstruction scheme, we consider some paradigmatic optical states: The coherent state, which is phase-sensitive, its phase-averaged counterpart (PHAV), which is of course phase insensitive, the single-photon Fock state, which is a nonclassical phase-insensitive state, and even and odd cat states, which are both nonclassical and phase-sensitive. For coherent and PHAV states we provide an experimental demonstration, whereas for the Fock and cat states we perform Monte Carlo simulated experiments, and analyze the effect of a non-unit quantum efficiency, η<1.

Experimental setup
For a coherent state r a a = ñá | |, the joint photon-number statistics, and the corresponding p D (Δ) distribution, can be obtained by setting P(ζ)=δ (2) (ζ−α), where δ (2) (z) is the complex Diracʼs delta function. The experimental setup is sketched in the right panel of figure 1. The second harmonic pulses (∼5 ps pulse duration) at 523 nm of a mode-locked Nd:YLF laser regeneratively amplified at 500 Hz are divided at a Mach-Zehnder interferometer into two parts in order to yield the signal and the LO. The length of one arm of the interferometer is changed in steps by means of a piezoelectric movement (Pz) in order to vary the relative phase f between the two arms in the interval [0, π]. Two variable neutral density filter wheels (ND) are used to change the balancing between the two beams, which are then recombined in a balanced BS. At the BS outputs two multi-mode 600 μm core fibers (MF) deliver the light to a pair of photon-number-resolving detectors. In detail, we employ two HPDs (mod. R10467U-40, Hamamatsu Photonics, having a quantum efficiency η ∼ 0.5 in the green spectral region), whose outputs are amplified, synchronously integrated and digitized. HPDs are commercial detectors endowed with partial photon-number-resolving capability and a good linearity up to 100 photons. As already demonstrated elsewhere (see for instance [27,25,45]), HPD response can be characterized in a self-consistent way with the same light under examination. Here we just remark that from the experimental data it is possible to calculate the gain of the detection apparatus in order to recover the distribution of detected photons at each BS output. In addition, this kind of detector allows us to obtain information on the relative phase between the two arms of the interferometer by simply monitoring the mean number of detected photons measured at each BS output as a function of the piezoelectric movement [27,46]. In the present case, 3 × 10 6 data corresponding to subsequent laser pulses are recorded and used to reconstruct the density matrix of the coherent state and of the PHAV state, which is a diagonal state obtained by randomizing the phase of a coherent state a a ñ = ñ

Coherent state
In order to process the data corresponding to the coherent state, once obtained the number of detected photons at each BS output, we calculate the shot-by-shot photon-number difference, Δ. The corresponding phase value f is determined by acquiring a suitable data sample (5×10 4 pulses) for each one of the 60 piezo positions. As the piezo moves in regular steps, the measured mean number of detected photons follows a sinusoidal trend due to the interference at the BS, from which the effective value of f for each piezo position can be evaluated. The 5×10 4 data corresponding to the same f are uniformly distributed around that value with a step of 1/(5×10 4 ).
The typical behavior of Δ f as a function of f is shown in figure 2(a) (3 × 10 5 data). Note that, due to irregularities in the piezoelectric movement, the data are not uniformly distributed in f. As mentioned in section 2, in order to estimate the uncertainties of the reconstructed states, we have chosen to perform a statistical analysis of different runs instead of following a tomographic approach, which is computationally more demanding.
In particular, the reconstruction procedure is applied to 10 sets of3 10 5 data drawn from the whole sample: the modulus of the average density matrix is shown in figure 2(b), and the typical uncertainty, obtained as standard deviation over the 10 runs, in the determination of the matrix elements is ∼10 −3 . The mean number . The fidelity [47] between the reconstructed coherent state and a coherent state with the same amplitude is F=99.4%.
As a matter of fact, the fidelity may not fully assess the quality of a reconstruction scheme [48] and therefore, in order to qualify our reconstruction scheme, we also consider the evaluation of the expectation values of some observables. To this aim, we remind that, given a generic observable, we may obtain its expectation value is a kernel, or pattern function, associated with the operatorÂ. For instance, in the following we will use the fist two moments of the quadrature operator, which may be obtained from the kernels [4]: . The last result explicitly shows the contribution to the variance due to the low-intensity LO, which becomes negligible as b a  | | | |. Notice that since we are dealing with classical states, the quantum efficiency just rescales the energy of the detected states, and we can safely set η=1: the results we obtain are thus referred to the detected states. We see that while the mean value of the quadrature is well reconstructed by the experimental data, the variance is larger than expected. Such discrepancy is likely due to phase fluctuations occurred during the measurements session.

PHAV state
In order to check whether this interpretation holds, we consider the PHAV state, which is phase insensitive and thus should not be influenced by the presence of possible phase fluctuations. As in the case of the coherent state, we save 10 sets of 3×10 5 data by calculating the shot-by-shot photon-number difference between the two BS outputs. Since the PHAV state is phase-insensitive, we randomly assign a phase value to each experimental value of Δ. A typical trace of Δ versus f is shown in figure 3(a), whereas in panel (b) the modulus of the average density matrix is presented. As expected, the off-diagonal elements are absent, and, also in this case, the typical uncertainty, obtained as standard deviation over the 10 runs, in the determination of the matrix elements is ∼10 −3 . In order to compare the reconstructed matrix to the theoretical prediction in equation (8)  , respectively. The very good agreement between theory and experiment for this phase-insensitive state confirms our considerations about the variance of the reconstructed coherent state.

Fock states
We now turn our attention to the Fock state ñ |1 , through which we explore the role played by a limited quantum efficiency for quantum-state reconstruction. Since the light source employed in the experiment is operated at 500 Hz, it does not guarantee an effective production of single-photon states. Thus, in the following we show the results obtained by using Monte Carlo simulated experiments, in which we set b = 20 .
The state ñ |1 has a highly singular P function given by: and the joint probability in equation (6) reads: where we assume that both detectors have the same quantum efficiency η. The corresponding p D (Δ) is quite clumsy and is not explicitly reported here.
In figure 4(a), we show a typical data trace, i.e.Δ f as a function of f, obtained from a Monte Carlo simulation (with 5×10 4 data) assuming ideal detection, i.e. η=1. It is clear that, as one may expect, Δ f does not depend on the relative phase f since the Fock state is phase insensitive. The modulus of the reconstructed density matrix is plotted in the panel (b) of the same figure. At variance with the experimental results, in this case the typical uncertainty in the determination of the matrix elements is ∼10 −2 , a value that in principle can be decreased either by increasing the samples or increasing the imbalance between signal and LO. Since the density matrix is diagonal, the fidelity of the state coincides with that of the photon-number distribution, for which we have F=99.9%.
Up to this point, we have shown the reliability of the reconstruction scheme and of the reconstruction strategy under ideal detection conditions. However, since photon-number-resolving detectors are real detectors, it is worth investigating whether the scheme would also work in the presence of an overall quantum efficiency η<1. When standard homodyne is employed, it is well known that a non-unit quantum efficiency rescales the mean value of the reconstructed density matrix in the case of classical states of light, such as coherent and PHAV states, whereas it deeply modifies the properties of the reconstruction in the case of nonclassical states, such as Fock states [41]. In order to check if analogous results can be achieved by means of our scheme, we consider the non-ideal detection of the Fock state ñ |1 in a simulated experiment. We assume the same value of LO adopted above, i.e. b = 20 , but now we set η = 0.4, a realistic value for many kinds of photon-numberresolving detectors [25,49,50]. In figure 4(c), we show the modulus of the reconstructed density matrix for the Fock state ñ |1 , whereas in figure 4(d) the corresponding photon-number statistics. From the reconstructed density matrix it is clear that the effect of a non-unit quantum efficiency is to add a vacuum component to the state: in this case the expected density matrix is r h and its fidelity with respect to the reconstructed one is F=99.0%. In this case the typical uncertainty in the determination of the matrix elements is ∼10 −2 .

Cat states
Now, we focus the attention on a class of states exhibiting a peculiar density operator, namely, the class of cat states [51][52][53]. These states represent a useful quantum resource for quantum teleportation, quantum computation and error correction. Here we consider the even and odd cat states defined as: In panel (a) of figures 5 and 6, we show typical data traces for the even and odd cat states with α=1, respectively, obtained from Monte Carlo simulations (with 2×10 4 data) in the case of η=1. As it is clearly evident from the two panels, the cat states are phase sensitive, at variance with Fock states. The modulus of the reconstructed density matrix for the even cat state is shown in panel (b) of figure 5, and for the odd cat state in panel (b) of figure 6. As expected, in the former case the density matrix only exhibits a non null probability only for even numbers of photons, whereas the latter one only for odd numbers. The good quality of the retrieved matrices is quantified by the high values of fidelity, which are F ∼ 99.9% and 99.8%, respectively. The situation changes if the quantum efficiency of the detector apparatus is lower than 1. For instance, in panels (c) of figures 5 and 6 we plot the reconstruction of the modulus of the density matrices for η=0.4. As it is evident, in both cases the matrices are no longer characterized by the presence of either even or odd elements. This fact is better emphasized in the panel (d) of both figures, in which the diagonal of the matrix is shown. Indeed, the non unit quantum efficiency is responsible for the appearance of odd elements in the case of the even cat state, and of even elements in the case of the odd cat state. In this respect, it is remarkable that the latter case is the worst because of the dominance of the 0-element. Nevertheless, the obtained results are in perfect agreement with the theoretical expectations, as testified by the high values of fidelity (F ∼ 99.7% for the even cat state and 99.8% for the odd one). However, we note that for the reconstruction of the cat states the use of a larger LO (b = 40 or 45 ) is needed to guarantee a good quality of reconstruction. We ascribe such a requirement to the more complex structure of the density matrix associated with these states. We will investigate the actual role of the LO amplitude in a future work.

Conclusions
In conclusion, for the first time to our knowledge we have demonstrated that photon-number-resolving detectors and a low-intensity LO can be used in a homodyne scheme instead of p-i-n photodiodes and a macroscopic LO to properly reconstruct quantum states of light. In particular, our numerical and experimental results highlight that reconstruction schemes based on the use of discrete variables may be effectively employed to retrieve the density matrix of continuous-variable systems, as well as the first two moments of the quadrature operator. At variance with standard homodyne detectors, our apparatus offers two main advantages: First of all, it is characterized by a larger experimental tolerance, since the reconstruction procedure can be applied also in the presence of a small imbalance between the two BS outputs. On the contrary, in the standard scheme an almost perfect balancing is required to avoid the damage of the circuit that amplifies the difference photocurrent. Secondly, since our apparatus is based on proportional PNR detectors to which we have direct access, each detector response can provide a self-consistent monitoring of the relative phase between signal and LO, as extensively explained in [24] and [27].
We have investigated advantages and limitations of our reconstruction scheme by performing experiments on coherent and PHAV states. The good results suggest that the homodyne-like scheme based on HPDs can find applications in discriminating among states that have the same photon-number statistics but different density matrix, such as between coherent and phase-averaged coherent states.
In addition, we have performed Monte Carlo simulated experiments with a single-photon Fock state and even and odd cat states in order to assess the role of quantum efficiency. We have found that, even in the presence of a very low value of the quantum efficiency, that further reduces the amplitude of the LO, the reconstructed density matrices corresponding to the detected states are consistent with the description model.
Our experimental scheme involves a low-intensity LO instead of a macroscopic one, thus proving that even a modest imbalance between signal and LO is sufficient to successfully perform quantum-state tomography. At variance with inversion algorithms, the reconstruction method we have applied to both experimental and simulated data properly works also in the presence of high loss. Moreover, it represents a valid alternative to maximum-likelihood methods which, being point estimators, do not naturally provide confidence intervals unless more than one data run is considered and fail in the case of a very low value of quantum detection efficiency. Overall, our results open new perspectives to quantum-state reconstruction in the mesoscopic regime, and pave the way to the use of photon-number-resolving-based detection schemes as novel hybrid schemes in Quantum Information science. For instance, they can be successfully employed also in the case in which, for some specific applications, a complete reconstruction of the states is not necessary. Moreover, the hybrid schemes could be used in all the applications in which the dual nature of light, namely wave-like and particle-like, should be exploited. For instance, the versatility of the hybrid schemes could be interesting for the codification and de-codification of information in complementary variables.