Quantum interferences reconstruction with low homodyne detection efficiency

Standard quantum state reconstruction techniques indicate that a detection efficiency of $0.5$ is an absolute threshold below which quantum interferences cannot be measured. However, alternative statistical techniques suggest that this threshold can be overcome at the price of increasing the statistics used for the reconstruction. In the following we present numerical experiments proving that quantum interferences can be measured even with a detection efficiency smaller than $0.5$. At the same time we provide a guideline for handling the tomographic reconstruction of quantum states based on homodyne data collected by low efficiency detectors.


Introduction
Homodyne detection is an experimental method that is used to reconstruct quantum states of coherent light by repeatedly measuring a discrete set of field quadratures [1][2][3]. Usually, a very high detection efficiency and ad-hoc designed apparatuses with low electronic noise are required [4]. New methods capable of discriminating between different quantum states of light, even with low detection efficiencies, will pave the road to the application of quantum homodyne detection for studying different physical systems embedded in a high noise environment [5][6][7][8]. For this purpose, specific quantum statistical methods, based on minimax and adaptive estimation of the Wigner function, have been developed in [9][10][11]. These approaches allow for the efficient reconstruction of the Wigner function under any noise condition, at the price of acquiring larger amounts of data. Hence, they overcome the limits of more conventional pattern function quantum tomography [12][13][14][15][16][17][18][19]. The important consequence of this novel statistical approach is that the 0.5 detection efficiency threshold can be overcome and quantum tomography is still practicable when the signals are measured with appropriate statistics. The scope of this paper is to report the results of this method tested by performing numerical experiments. Indeed, we consider a linear superposition of two coherent states and numerically generate homodyne data according to the corresponding probability distribution distorted by an independent Gaussian noise simulating efficiencies lower than 0.5. By properly expanding the set of numerically generated data, we are able to reconstruct the Wigner function of the linear superposition within errors that are compatible with the theoretical bounds. Our results support the theoretical indications that homodyne reconstruction of linear superposition of quantum states is indeed possible also at efficiencies lower than 0.5.

Wigner function reconstruction
Let us consider a quantum system with one degree of freedom described by the Hilbert space L 2 (R) of square integrable functions ψ(x) over the real line. The most general states of such a system are density matricesρ, namely convex combinations of projectors |ψ j ψ j | onto normalised vector stateŝ Any density matrixρ can be completely characterised by the associated Wigner function W ρ (q, p) on the phase-space (q, p) ∈ R 2 ; namely, by the non positive-definite (pseudo) distribution defined by Hereq andp are the position and momentum operators obeying the commutation relations [q ,p] = i, = 1, and |q ± v/2 are eigenstates ofq:q|q ± v/2 = (q ± v/2)|q ± v/2 . Notice that W ρ (q, p) is a square integrable function: 2π Among the advantages of such a representation, is the possibility of expressing the mean value of any operatorÔ with respect to a stateρ as a pseudo-expectation with respect to W ρ (q, p) of an associated function O(q, p) over the phase-space, where Indeed, by direct inspection one finds In homodyne detection, a monochromatic signal photon state is mixed with a coherent reference state, a so-called local oscillator, by a 50/50 beam splitter. The output is collected by two photodiodes and the difference photocurrent is measured. It can be proved that, when the local oscillator is significantly more intense than the signal, the homodyne photocurrent is proportional to the signal quadrature [20]. Denoting byâ andâ † the single mode annihilation and creation operators associated with the signal, the quadrature operator is defined aŝ where φ is the relative phase between signal and local oscillator. The continuum set of quadratures with φ ∈ [0, π] provides a complete characterization of the signal state. Using the annihilation and creation operatorsâ,â † one constructs position and momentum-like operators,q = (â +â † )/ . With respect to the latter, the quadrature operator reads:x φ =q cos φ +p sin φ .
Quadrature operators have continuous spectrum extending over the whole real line,x φ |x = x |x ; given a generic one-mode photon state associated with a density matrixρ, its diagonal elements with respect to the (pseudo) eigenvectors represent the probability distribution over the quadrature spectrum. In homodyne detection experiments the collected data consist of n pairs of quadrature amplitudes and phases (X , Φ ): these can be considered as independent, identically distributed stochastic variables. Given the probability density p ρ (x, φ), one could reconstruct the Wigner function by substituting the integration with a sum over the pairs for a sufficiently large number of data. However, the measured values x are typically not the eigenvalues ofx φ , rather those ofx where y is a normally distributed random variable describing the possible noise that may affect the homodyne detection data and η parametrizes the detection efficiency that increases from 0 to 100% with η increasing from 0 to 1 [9]. The noise can safely be considered Gaussian and independent from the statistical properties of the quantum state, that is, y can be considered as independent fromx φ . Then, as briefly summarised in Appendix A, the Wigner function is reconstructed from a given set of n measured homodyne pairs (X , Φ ), = 1, 2, . . . , n, by means of an estimator of the form [9] W η,r h,n (q, This expression is an approximation of the Wigner function in (27) by a sum over n homodyne pairs (X , Φ ). The parameter h serves to control the divergent factor exp(γξ 2 ), while r, through the characteristic function χ r (q, p) of a circle C r (0) of radius r around the origin, restricts the reconstruction to the points (q, p) such that q 2 + p 2 ≤ r 2 . Both parameters have to be chosen in order to minimise the reconstruction error which is conveniently measured [10] by the L 2 -distance between the true Wigner function and the reconstructed one, W ρ − W η,r h,n 2 . Since such a distance is a function of the data through W η,r h,n , the L 2 -norm has to be averaged over different sets, M , of quadrature data: where E denotes the average over the M data samples, each sample consisting of n quadrature pairs (X , Φ ) corresponding to measured values of x φ with φ ∈ [0, π]. In [10], an optimal dependence of the parameters r and h upon the number of data, n, is obtained by minimizing an upper bound to ∆ η,r h,n (ρ). 1

Interfering Coherent States
Homodyne reconstruction is particularly useful to expose quantum interference effects that typically spoil positivity of the Wigner function: it is exactly these effects that are claimed not to be accessible by homodyne reconstruction in presence of efficiency lower than 50%, namely when η in (8) is smaller than 1/2 [15]. However, in [9] it is theoretically shown that η < 1/2 only requires increasingly larger data sets for achieving small reconstruction errors. However, this claim was not put to test in those studies as the values of η in the considered numerical experiments were close to 1. Instead, we here consider values η < 1/2 and reconstruct the Wigner function of the following superposition of coherent states with α any complex number α 1 + iα 2 ∈ C. The Wigner function corresponding to the pure stateρ α = |Ψ α Ψ α | is shown in Figure 1 for α ∈ R. Its general expression together with that of its Fourier transform and of the probability distributions p ρ (x, φ) and p η ρ (x, φ) are given in Appendix B.
In Appendix C, a derivation is provided of the L 2 -errors and of the optimal dependence of h and r on the number of data n and on a parameter β that takes into account the fast decay of both the Wigner function and its Fourier transform for large values of their arguments. The following upper bound to the mean square error in (11) is derived: 1 The functional relation between the parameters h and r on n also depends on an auxiliary parameter β > 0. This was introduced in [10] to characterise the localisation properties on R 2 of the Fourier transforms of the Wigner functions of the following class of density matrices addressed in that context: A β,s,L = ρ : with 0 < β < 1/4 and γ as given in (26). As explained in Appendix C, the quantities ∆ 1,2,3 do not depend on h, r and n. By taking the derivatives with respect to r and h, one finds that the upper bound to the mean square deviation is minimised, for large n, by choosing We generated M = 10 samples of n = 16 × 10 6 quadrature data (X , Φ ) distributed according to the noisy probability density p η ρ (x, φ) explicitly given in (32) of Appendix B, considering an efficiency lower than 50% (η = 0.45). Starting from each set of simulated quadrature data we reconstructed the associated Wigner function by means of (9) and (10). The averaged reconstructed Wigner functions E W η,r h,n (q, p) for η = 0.45 are shown in Figure 2 for two different values of the parameter β. The mean square error of the reconstructed Wigner functions has been computed as in (11) and compared with the mathematically predicted upper bounds ∆. The dependence of the upper bound reconstruction error on the parameter β is discussed at the end in Appendix C. In Table 1 Despite common belief, the interference features clearly appear in the reconstructed Wigner function also for efficiencies lower than 50% and the reconstruction errors are compatible with the theoretical predictions. In the next section, we make a quantitative study of the visibility of the interference effects.

A witness of interference terms
The interference effects in the state |Ψ α can be witnessed by an observableÔ α of the form With respect to an incoherent mixture of the two coherent states, its mean value is given by Therefore, from (3) it follows that the phase-space function O α (q, p) associated toÔ α is For details see (28) and (31) in Appendix B. Let us denote by W α j,rec (q, p), the estimated Wigner function W η,r h,n (q, p) in (9) for the j-th set of collected quadrature data. It yields a reconstructed mean value of which one can compute mean, Av(<Ô α > rec ), and standard deviation, Sd(<Ô α > rec ), with respect to the M sets of collected data: We computed Av( Ô α rec ) and Sd(<Ô α > rec ) with M = 10 simulated sets of noisy data with η = 0.45 for two different numbers of simulated quadrature data (see Figure 3). We repeated the procedure for different values of the parameter β. The results are presented in Figure 3, where the error bars represent the computed Sd( Ô α rec ).
In order to be compatible with the interference term present in |Ψ α , the reconstructed Wigner functions should yield an average incompatible with the incoherent mean value in (17), namely such that We thus see that the condition in (22) is verified for η = 0.45, that is the reconstructed Wigner functions are not compatible with incoherent superpositions of coherent states, if enough data are considered. We also notice that the same behavior is valid for the high efficiency η = 0.95. The dependence of the errors on β can be understood as follows: when β decreases the integration interval in (10) becomes larger and approaches the exact interval [−∞, +∞]. Nevertheless this occurs at the price of increasing the reconstruction error. This can be noted both in Figure 3 (larger error bars) and in Figure 2 (increasingly noisy effects in the reconstructed Wigner function). This problem can be overcome with a larger number of data samples M , that reduce the reconstruction noise and compensate for the effect of decreasing β.

Conclusions
We simulated quadrature data corresponding to high electronic noise and detection efficiencies lower than 0.5. Under these operating conditions the Wigner function of a linear superposition of two coherent states could be reconstructed using the tomographic techniques developed in [9]. Moreover, by taking into account the decay properties of the Wigner function along with those of its Fourier transform, we have checked that the numerical reconstruction errors are compatible with the theoretical error bounds computed there. Furthermore, the reconstruction of the quantum interference pattern of the Wigner function has been supported by a numerical study of the variance of an operatorial interference witness excluding that the reconstructed interferences might be an artefact of the reconstruction algorithm.
We have thus confirmed 1) that, as theoretically predicted in [9, 10], a 0.5 detection efficiency is not, as often stated in the quantum optics literature, an absolute threshold below which homodyne quantum state reconstruction is generically impossible and 2) that, instead, by suitably enlarging the size of the set of collected quadrature data, and using alternative techniques different from standard pattern function quantum tomography, one may indeed access quantum features even in low efficiency conditions. These results also provide the tools for quantum state reconstruction in these operating conditions and set the boundaries of the applicability of the novel statistical approach to homodyne quantum state reconstruction by checking the increase of the number of quadrature data necessary for faithful reconstruction with decreasing detector efficiency.

A Wigner function reconstruction
The quadrature probability distribution (7) can be conveniently related to the Wigner function by passing to polar coordinates u = ξ cos φ, v = ξ sin φ, such that 0 ≤ φ ≤ π and −∞ ≤ ξ ≤ +∞: |ξ| (2π) 2 e iξ(q cos φ+p sin φ) Tr ρ e −iξ(q cos φ+p sin φ) where F [p ρ (x, φ)](ξ) denotes the Fourier transform with respect to x of the probability distribution: Since y can be considered a normally distributed random variable independent ofx φ , the noise affected distribution of the eigenvalues ofx φ in (8) is given by the following convolution: Its Fourier transform is connected with that of p ρ (x, φ) according to By inserting F [p ρ ( x , φ)](ξ) into (24), one can finally write the Wigner function in terms of the noisy probability distribution p η ρ (x, φ):

C Upper bound reconstruction error estimation
Here we derive an upper bound to the mean square error of the reconstructed Wigner function. This analysis is necessary in order to find an optimal functional relation between the free parameters in the reconstruction algorithm such to minimise the reconstruction error. For this purpose we follow the techniques developed in [10] adapting them to the case on a linear superposition of coherent states. Using (9) and (10), one starts by rewriting the error in (11) as the sum of three contributions: C c r (0) denoting the region outside the circle C r (0), where q 2 + p 2 > r 2 . The first and the third term correspond to the variance and bias of the reconstructed Wigner function respectively, while the second term is the error due to restricting the reconstruction to the circle C r (0). Given a density matrixρ, the second term can be directly calculated. This is true also of the bias; indeed, because of the hypothesis that the pairs (X , Φ ) are independent identically distributed stochastic variables, it turns out that The variance contribution can be estimated as follows: firstly, by using (9) and (10), one recasts it as Then, a direct computation of the first contribution yields the upper bound with o(1) denoting a quantity which vanishes as h when h → 0. On the other hand, the second contribution can be estimated by extending the integration over the whole plane (q, p) ∈ R 2 and using (37) together with (2): Whenever β is such that the arguments of the logarithms are much smaller than the number of data n, to leading order in n the upper bound to the mean square deviation is minimised by The range of possible values of β is 0 ≤ β ≤ 1/4. However, the upper bound ∆ becomes loose when β → 1/4 and β → 0. In the first case, it is the quantity ∆ 3 (β) which diverges, in the second one, it is the variance contribution which diverges as the logarithm of the number of data. It thus follows that the range of values β ∈ [β 0 , β 1 ] where the numerical errors ∆ η,r h,n (ρ α ) are comparable with their upper bounds ∆ is roughly between β 0 = 0.04 and β 1 = 0.10 for η = 0.45 as indicated by the following Figure 4.