Direct probing of the Wigner function by time-multiplexed detection of photon statistics

We investigate the capabilities of loss-tolerant quantum state characterization using a photon-number resolving, time-multiplexed detector (TMD). We employ the idea of probing the Wigner function point-by-point in phase space via photon parity measurements and displacement operations, replacing the conventional homodyne tomography. Our emphasis lies on reconstructing the Wigner function of non-Gaussian Fock states with highly negative values in a scheme that is based on a realistic experimental setup. In order to establish the concept of loss-tolerance for state characterization we show how losses can be decoupled from the impact of other experimental imperfections, i.e. the non-unity transmittance of the displacement beamsplitter and non-ideal mode overlap. We relate the experimentally accessible parameters to effective ones that are needed for an optimised state reconstruction. The feasibility of our approach is tested by Monte Carlo simulations, which provide bounds resulting from statistical errors that are due to limited data sets. Our results clearly show that high losses can be accepted for a defined parameter range, and moreover, that (in contrast to homodyne detection) mode mismatch results in a distinct signature, which can be evaluated by analysing the photon number oscillations of the displaced Fock states.


Introduction
The preparation and accurate characterization of non-Gaussian states is of paramount importance for quantum processing applications. For example the implementation of scalable continuous variable quantum information protocols require an entanglement distillation procedure, which has been proven to be impossible using only a set of Gaussian states and Gaussian operations [1,2,3]. The most fundamental non-Gaussian states are Fock states. Recent experiments have demonstrated for the first time the generation of Fock states and photon added/substracted states [4,5,6,7,8]. In all of these experiments conventional homodyne detection [9,10] followed by tomographic reconstruction was employed to prove the presence of the non-Gaussian characteristics, as well as the presence of negative values, in the Wigner functions of the quantum states [11]. These negative values are considered to be a signature of quantum behaviour. In a homodyne measurement, the detected signal consists of a beating between the quantum signal and a strong local oscillator field, which in turn has the effect that mode mismatch cannot be distinguished from losses. Another drawback in this method is the fact that homodyne tomography is always an indirect state characterization: the Wigner function can be reconstructed only after intricate computational back-projection, and the photon statistics have to be recovered by more involved means [12,13,14].
A more direct characterization of the Wigner function is possible if a photon number resolving detector is available. The Wigner function at a given position α of phase space (α ∈ C) is equal to the average value of the photon number parity operator on a quantum state displaced by −α [15]. The expectation value of the parity operator is directly inferred from a complete measurement of the photon number distribution [16,17]. As opposed to tomography, in this method the Wigner function is probed point-bypoint in phase space, resulting in a high resolution measurement. Until now, the requirement of measurements with photon number resolution has severely limited the usage of this method, although different ideas to overcome the problem have arisen for specific systems [18,19,20].
Recently, new strategies for measuring light with photon number resolution have been developed, which allow a more intensive exploration and characterization of quantum states [21]. One of them-applicable for pulsed light-is a time multiplexed detector (TMD). With the TMD, the photon number distribution of a quantum state can be revealed via an inversion procedure applied to the measured photon statistics. Experimental applications, demonstrating a reliable loss calibration, and the TMD's suitability for detecting multimode statistics have been accomplished [22,23].
In this paper, we apply for the first time to our knowledge the probing method in connection with a TMD to reconstruct the Wigner functions of single-photon and two-photon Fock states, showing the capability of loss-tolerant state characterization. This is achieved by decoupling the effects of experimental imperfections from the ideal mathematical operations. We find an equivalent model for the experimental setup, in which losses are independent from the displacement operation. In addition, mode mismatch is modeled by a convolution of a Poissonian distribution with the statistics of an ideally displaced Fock state. Via Monte Carlo simulations we evaluate the specific experimental conditions necessary to characterize reliably non-Gaussian Fock states. Using a reasonable number of data points (10 6 ) we show the extent of degradation caused by losses and mode mismatch. A simple inspection of the characteristic photon number oscillation of displaced Fock states enables us to quantify the degree of overlap. This paper is organised as follows. In Sec. 2 we briefly introduce the method for probing the Wigner function using displacement operations and photon counting. Experimental issues related to the implementation of the ideal displacement operation are discussed in Sec. 3. In Sec. 4 we recall important properties of a TMD and explain how to handle the effect of losses via an inversion algorithm. In Sec. 5 we introduce a model to study the effect of imperfect mode overlap, which is used in Sec. 6 in order to relate the degree of overlap to usual experimental observables. In the following section (Sec. 7), we exploit this theoretical framework and employ Monte Carlo simulations to analyse the reliability of the probing technique under the usual circumstances of high losses and imperfect mode overlap. We show how the Wigner function characterictics can be efficiently recovered from the response of the TMD. Finally, in Sec. 8 we conclude and highlight the main findings of our investigations.

Wigner function reconstruction by direct probing method
The probing method is based on a simple expression that relates the Wigner function at the origin of the phase space with the mean value of the photon number parity operator Π [16,17], As the Fock states |n are eigenstates of the photon number parity operator, with eigenvalues (−1) n , a natural representation for the density matrix isρ = ρ mn |m n|.
With this representation Eq. (1) simplifies to Since ρ nn is equal to the probability of counting n photons, the Wigner function can be directly measured from the photon statistics. Due to the fact that the projection onto the photon number basis always erases the information about the coherence of the state, the parity measurement returns the value of the Wigner function only at one single point of phase space. Nevertheless, the whole Wigner function can be reconstructed if appropriate displacements are applied to the state. To access its value at a given position α in phase space we need to displace the state by the quantity −α. Mathematically, this is described by the action of the displacement operatorD(−α) =D † (α). The evolution of the density operator is given byρ (−α) =D(α) †ρD (α), which results in Similarly to the expression in Eq. (2), the value of the Wigner function in Eq. (3) is determined directly by the photon statistics of the displaced state. In general terms, the reconstruction of a quasi-probability distribution function of a degraded quantum state via the direct probing technique has been analysed in the past [16,24,25]. Here, we extend the previous work by showing, how to deal with the experimental imperfections in practical terms in order to achieve a loss-tolerant state characterization. To increase the amount of information, which can be extracted from our simulations we establish a strong connection between the theory and the actual experimental parameters that constrain the reconstruction of the Wigner function. In the course of this work we evaluate the impact of imperfections in this characterization method, i.e. losses and mode mismatch between fields. In addition, an important degrading effect comes from the implementation of the displacement, which nevertheless can be handled with a loss inversion strategy similar to the inefficient detection.

Experimental implementation of the displacement operation
In an experiment the displacement operation can be well approximated by the effect of a highly asymmetric beam splitter with transmission T ≈ 1, which is used to superimpose the quantum signal of interest with a coherent reference field |β [26,27]. The beam splitter is a linear optical element that connects two input modes a and b, delivering output modes A and B, as shown in Fig. 1(I). In the Heisenberg picture, the output  annihilation operators are related to input operators by the transformationŝ with T and R = 1−T being equal to the transmission and reflection intensity coefficients. From the experimental point of view, the beam splitter will always introduce some amount of loss, since T = 1. To decouple the displacement from the loss we analyse the configuration shown in Fig. 1(II). It is straightforward to derive the conditions under which the networks (I) and (II) are equivalent, if we follow an approach introduced by Wallentowitz and Vogel [16]. Firstly, using the beam splitter transformation in Eq. (4), the output annihilation operators of each model are related to the input modes by, in which we employed the relationship D † (γ)ÂD(γ) =Â + γ. Second, we note that in a photon-counting experiment the probability of recording n photons is given by [28] P n = : e −NN n n! : with the symbol :: defining a normal ordering of the operators. The average in Eq. (8) is taken over the initial quantum state |Ψ andN =Â †Â is the number operator of the detection mode. As mode b is prepared in a coherent state in both models, and P n is evaluated by the average of normally ordered operators, the operatorb can be replaced by a complex number β, withb |β = β |β . Thus, the annihilation operators from models (I) and (II) introduced by the number operator in Eq. (8) can be substituted bŷ We are only interested in the photon number probability P n and therefore define effective photon number operatorsN eff , which are independent of losses but represent effective displacements. The number operators for networks (I) and (II) can then be substituted byN Examining the last two equations we reason that both schemes manifest equal photon number probabilities if γ = √ 1 − T β. We would like to emphasise that the loss element in Eq. (12) modelled by the beam splitter with transmission coefficient T , increases the net amount of displacement, i.e. the total amount of displacement is α = γ/ √ T . In addition to the displacement, the resulting probability is given by an attenuated signal with an the overall quantum efficiency of η = T . As a conclusion we find that even when the beam splitter transmission is far from the idealized case, the device still works as a displacing element and introduces simply some amount of extra loss that can be inverted as will be described in Sec. 4.
To simulate the possible statistics resulting from an actual measurement, we can use network (II) to calculate the photon number distributions of displaced Fock stateŝ D(γ) |m . In the absence of loss, these photon number distributions are well known and given by ρ nn (γ) = | n|D(γ)|m | 2 [29]. The effect of losses can be analysed in a very straightforward fashion. For Fock states a beam splitter with transmission coefficient T yields a statistical mixture, which is governed by binomial distribution For the Fock states |1 and |2 this results explicitly in

Loss tolerant characterization using TMD
According to the previous Sec. 3 the measurable probability P n differs from the ideal ρ nn that is required for the reconstruction of the Wigner function in Eq. (3). In this section we study how to relate these quantities and generalise the model to allow for the detection losses. In this case the quantum signal is degraded both before and after the realisation of the displacement. The network shown in Fig. 2 represents the sequence of operations undergone by the signal, where losses are modelled as beam splitters combining the signal with vacuum modes. The last beam splitter, with transmission , represents the limited quantum efficiency of the detector itself. As before, to evaluate Representation of a theoretical model in which an input quantum signal undergoes a loss, followed by a non-ideal displacement and is finally measured with inefficient detection. The displacement is realised by a high-transmission beam splitter in which the signal is superimposed with a weak reference field and the photon number resolution is achieved with a TMD.
the photon number probabilities, we introduce the effective photon number operator, which can be expressed aŝ with η = T η 1 being the overall quantum efficiency. The total amount of displacement is equal to α = γ/ √ η 1 T . Note that in contrast to losses experienced by the signal before the displacement, the detection efficiency has no impact on α.
In analogy to Eq. (8) with the detection loss η, we note that the exponential factor −ηÂ †Â can be rewritten as [1 + (1 − η)]Â †Â . Employing the series expansion of exponential function it is easy to show that where we define L n,j as an element of an upper diagonal matrix, referred as the loss matrix. Recording n events with a non-ideal detection, ensures that the signal contained at least j ≥ n photons. We note that Eq. (16) differs from Eq. (13) in the sense that the right hand term | j|D(α)|Ψ | 2 directly corresponds to the non-degraded statistics of the displaced initial stateD(α) |Ψ . If Eq. (16) is expressed in matrix form PD |Ψ = L(η) ρD |Ψ , we can easily verify that the loss inverted statistics correspond to the ideal displaced statistics. Thus, we can deduce that the effect of non-ideal detection can in principle be handled by the loss inversion method. Nevertheless, photodetectors with photon number resolution are still essential. Photon number resolution can be achieved by distributing several photons into many different temporal modes, which are then analysed with binary detectors by time multiplexing. In practice, a TMD consists of a fibre-integrated 50/50 beam splitter network with time delays implemented as fibre loops of different lengths (Fig. 2). A pulse launched into the detector is divided into a train of pulses in two separate fibre outputs [30,31]. The existence of photons in each mode (or time-bin) can be detected by avalanche photo diodes (APDs) and the so-called click-statistics are collected for ensemble measurements.
The measured click statistics p = (p 0 , p 1 , ... p n ) are related to the input photon number distribution through the relation [32] The convolution matrix C takes into account a stochastic distribution of n photons into several bins (from which follows that the number of clicks N is lower or equal to n) and the matrix L(η) describes the binomial process of loss, as introduced in Eq. (16). Therefore, by inverting equation (17) (or using a maximum likelihood technique), the photon number distribution of the quantum signal can be reconstructed from the click statistics if the overall loss η and the convolution matrix are known. Beside the numerical techniques an analytical solution can be used for inverting the losses, as presented in Appendix A. Experimentally, the calibration of the overall efficiency η = η 1 T can be reliably accomplished by taking advantage of the twin photon generation properties of a parametric down conversion process [22,23]. In addition to the overall loss, the total amount of displacement α = γ/ √ η 1 T = √ 1 − T β/ √ η 1 T needs to be properly calibrated. The transmission T of the beam splitter can be characterized by monitoring both outputs of the beam splitter and measuring the fraction of a strong input field β that is transmitted and detected by a standard photo diode. An additional measurement extracts the value of by, e.g. measuring the Poissonian statistics of the coherent reference field |β = | (1 − T ) β that is incident on the TMD. Finally, combining the results of the unbalanced homodyning of vacuum together with the knowledge of the overall losses and T , enables us to determine the amount η 1 . Thus, all losses can be experimentally characterized enabling the loss-tolerant probing of the Wigner function.

Imperfect mode overlap between quantum signal and reference beam
Mode overlap is a crucial issue for both unbalanced and balanced characterization techniques. In balanced homodyne measurements the detected photocurrentÎ is proportional to the product of the quantum signal and the local oscillator,Î ∝â β * +â † β, [14]. As only the overlapping part of the signal state is measured, imperfect mode overlap has intrinsically the same signature as a loss. Contrariwise, in unbalanced measurements such as the direct probing method, the detection is proportional to the photon number operator, represented by sum of the signal and reference field, i.e.n ∝ (â † + β * )(â + β). Each photon impinging on the APD causes a detectable photocount, and the signature of the mode mismatch becomes different from loss [33].
We model the imperfect mode overlap by the beam splitter network as shown in Fig. 3. In the following we elaborate the conditions under which this network is equivalent to the one depicted in Fig. 4. As before we relate the actual experimental parameters with the values of the equivalent network. In the network of Fig. 3, the signal and the reference beam are divided into four modes, two of them interfering. The other modes do not overlap and are just split according to the ratio T . Hereby we assume that t 1 and t 2 correspond to the fractions of the signal and reference fields, which are perfectly mode matched and the product M = t 1 t 2 defines the degree of overlap. The probability P n of detecting n photons is given by the convolution of the joint probability distribution P joint (p, m, q). The parameters p, m and q label the photon number in the three modes of the network: the signal ancilla, the displaced signal, and the reference ancilla. Using the multimode theory of photo detection [28], we can simplify the evaluation of the convolution by expressing the photon number distribution P n in the form of Eq. (8), where the photon number operatorN tot of the multimode field is equal to the sum of the number operators of each different detection modeN tot =p †p +m †m +q †q . Applying the beam splitter transformation of Eq. (4), the evolution of the different modes are given byp in whichv i , i = 1, 2, 3, are auxiliary operators that represent the modes of the vacuum, as depicted in Fig. 3. Once more, considering coherent states, the operators can be replaced by complex numbers. The number operatorN tot can be substituted bŷ where ξ = √ M 1−T T β is the amount of displacement. The parameter |ζ| = (1 − M)(1 − T )|β| can be interpreted as the amplitude of an effective coherent field.
Finally, the probability of n counts is given by P n = : : From Eq. (23), we recognise that the photon number distribution is a convolution of two different photon number distributions. The first contribution corresponds to the statistics of an ideally displaced quantum signal, detected with efficiency T , and the second corresponds to a Poissonian distribution. As a consequence, the network shown in Fig. 3 and its counterpart depicted in Fig. 4, which is conceptually simplified, exhibit equivalent photon number probabilities.  Our derivations lead to the conclusion that the detected photon number distribution can always be expressed as a convolution of the displaced state distribution and a Poissonian distribution. As first mentioned in [25], a Wigner function constructed with the distribution in Eq. (23) is composed of the product of two separate Wigner functions. The decreasing overlap between the signal and reference fields distorts the reconstructed Wigner function such that it deviates towards a Gaussian multiplied with the value of the signal state Wigner function at the origin. This results in an artificial broadening of the reconstructed Wigner functions of low-photon number Fock states.

Mode matching diagnosis
Beyond the investigation of loss-tolerant reconstruction of the Wigner function, direct probing also allows us to explore the non-trivial photon number distributions of displaced Fock states. As a result of their exact photon number, Fock states exhibit a completely undetermined phase, which can be depicted by a rotational symmetry of the field quadratures around the origin of the phase space. In other words, Fock states can be imagined as rings around the origin (uncertain phase) that have definite amplitudes (exact photon numbers). Photon number detection corresponds to evaluating the overlap of the studied signal state with Fock states. Using this concept it becomes intuitive that the photon number distributions of displaced Fock states must present oscillations as a consequence of an interference process in phase space [34,35,36]. Fig. 5 presents typical oscillations of the states |1 and |2 . Considering e.g. the displaced single photon Fock state, we clearly see that the probability of measuring one photon decreases as the displacement is increased. It is maximally suppressed when the displacement reaches the mean photon number n α = 1, and increases again for higher values of the displacement. In this section this quantum behavior is explored in order to distinguish the effects of losses from those of imperfect mode overlap between the signal and reference fields. We numerically study the reconstruction of the photon number distribution of the single photon Fock state as a function of the degree of overlap M introduced in the last section. The photon number distribution is evaluated using Eq. (23), which gives for P 0 and P 1 the following expressions We consider a beam splitter with T = 0.95 and reconstruct the photon number distributions for the cases of 100%, 50% and 0% overlap applying the loss inversion algorithm and scaling of the displacement as discussed previously. Fig. 6 illustrates the characteristic behaviour of the vacuum and one photon components of the photon number distribution with respect to the displacement applied to single photon Fock state. In the case of perfect overlap between signal and reference we expect a complete suppression of the one photon component when the displacement is α = 1 [ Fig. 6(a)], for 50% overlap, the oscillation of the vacuum and one photon components is reduced but still visible [ Fig. 6(b)]. However, for non-overlapping fields the oscillation is completely washed out and the photon number distribution is just a convolution of the single photon state statistics with a Poissonian [Fig. 6(c)]. A clear signature of the lack of mode matching becomes apparent in the vacuum component of the inverted statistics: in the case of non-overlapping fields this contribution is always zero independent of the displacement. This effect is opposed to the signal degradation caused by losses, which always increases the probability of obtaining the zeroth photon number component. These results suggest a recipe to directly deduce the degree of mode matching from the amount of oscillation of the vacuum and one photon components. In Fig. 7 we show the values of A 00 = |ρ 00 (1)−ρ 00 (0)| and A 11 = |ρ 11 (1)−ρ 11 (0)| as a function of M, where ρ ii (α) is the probability of having i photons in the photon number distribution at the displacement α. Thus, A ii represents a measure of the amplitude of the oscillation. From an experimental point of view, this method is an extremely direct way of characterizing the amount of overlap, as it solely requires the measurement of the vacuum (or one photon) component without displacement and at a single point with α = 1.

Monte Carlo Simulations
In order to characterize a quantum state experimentally it has to be prepared several times enabling an ensemble measurement, which allows us to reconstruct the photon number statistics of the state. Here, a similar scenario was achieved using a Monte Carlo simulation. The detected click statistics were simulated with the help of Eqs. (14) and (14), which describe the photon number distribution of displaced Fock states. We then replaced the parameter T with the value of the overall loss η and recalibrated the value of the displacement accordingly. By performing a Monte Carlo simulation we generated the statistics of the studied quantum states, which were then convolved with the known stochastic description of the TMD. As the states under consideration (|1 and |2 ) are symmetric in phase space, we assumed fixed values for the relative phase between signal and reference beam and changed only the amplitude of the displacement.
To recover the photon number statistics of the displaced state we employed the direct numerical inversion of Eq. (17) and evaluated the photon number parity from the statistics to yield the reconstructed value of the Wigner function. In order to estimate the statistical errors, for each sampled point of the phase space we performed ten Monte Carlo simulations, each consisting of 10 6 events. The chosen number of events was based on the expected number of experimental trials; this depends, for example, on the setup stability. This set of simulations gave at each sampled point ten simulated values for the Wigner function, from which we evaluated the mean value and the standard deviation.
There are two factors limiting the range over which the state reconstruction can be considered reliable. First, as a consequence of the TMD resolution, the mean photon number of the displaced beam is restricted. With a m-stage TMD (incorporating m 50/50 couplers) we are able to resolve at maximum 2 m clicks (equal to the number of bins into which the incident pulse is separated). The reconstruction becomes unstable as soon as the displaced state has a significant contribution of photon numbers larger than the number of the TMD bins. In our numerical investigations, we observed that the instability can be tolerated until the probability of having more than 2 m photons in the displaced state is on the order of 2 − 3%. The second limitation arises from statistical errors and can be estimated from the amount of collected events versus losses. Statistical errors are highly amplified in the loss inversion process. As the losses increase the inversion becomes unstable, yielding a singular loss matrix and negative probability components for the reconstructed photon number statistics. We define the reliable reconstruction range by a boundary value for the displacement at which the inversion becomes unstable or, similarly, we accept the values of displacement at which the photon number statistics reconstructed from all 10 7 events render components ρ nn > −0.001.
Any experimental implementation of the direct probing scheme, as envisioned with the current technology, is subjected to high losses. For example, the quantum efficiencies of APDs can reach values between 60-80%. The coupling efficiency into single-mode fibres is highly mode-sensitive and can drop down to 30-40% in case of mismatched modes. Another important source of degradation of the signal is the beam splitter used as displacing element. Typically, the transmittance of this beam splitter is around 95%. Due to the 5% reflectance, the beam splitter intrinsically attenuates the quantum signal and modifies the displacement. Effectively, in a realistic experiment for state characterization the achievable detection efficiency will lie between 20% and 30%. In addition to the effects of inefficiencies or losses, we must take into account the accuracy with which the convolution matrix of the TMD can be determined. This is caused by unbalanced 50/50 couplers, which are generally specified with 5% uncertainty.

Results
We performed the simulation for four different scenarios to clarify the effects of different experimental frames. We analysed the performance of the state characterization and estimated the limitations by considering a realisable displacement operator implemented with a beam splitter of 95% transmittance, a detection with three-stage (8-bin) TMD and perfect efficiency. Furthermore, we analysed the possible benefits associated with the use of a four-stage (16-bin) TMD taking also into account a finite detection efficiency. Next, the amount of degradation was increased, considering the experimentally achievable efficiency level. Finally, the effect of imperfect mode overlap was included.
Case 1: First we studied the idealized case of perfect detection efficiency but included the finite resolution of a TMD due to its limited number of stages. The quantum signal was displaced at a beam splitter with 95% transmittance and the TMD was assumed to be constructed from ideal 50/50 couplers. The results show that the reconstructed Wigner functions of the |1 and |2 Fock states are slightly loss-degraded around the origin due to the reflection of the signal at the beam splitter. Nevertheless, this degradation can be easily handled with the direct inversion, as shown in Fig. 8. Due to the limited number of TMD stages, a reliable reconstruction of W (α) is possible only inside a given scanning range of phase space, i.e. there is a maximum value of α that can be safely used. The ultimate boundaries for the reliable scanning ranges correspond to displacements of α = 1.5 and α = 1.2 for states |1 and |2 , respectively. Nevertheless, with lossless detection, the 8-bin TMD is adequate to scan the interesting regions of the phase space for both of these states. The simulated setup has a displacing beam splitter with 95% transmittance and ideal detection efficiency. The reconstruction is very accurate with errors smaller than used symbols. For further details see text.

Case 2:
In the second simulation we considered more realistic parameters, which might be accessible in the near future. We assumed again a beam splitter with 95% transmittance, and a detection efficiency of 60%. In this case we compared the performance of 8-and 16-bin TMDs for reconstructing the Wigner function of the single photon Fock state. As an example of an unbalanced TMD we chose couplers with splitting ratios of 45/55. However, the deconvolution was still done by assuming perfect 50/50 couplers, as the imbalance of 5% is due to uncertainty in the beam splitter specification. As the detection efficiency decreases, the effect of the loss becomes prominent restricting the reliable scanning range in phase space. Moreover, the effect of the unbalanced deconvolution becomes negligible in comparison to the effect of loss. An inspection of Fig. 9 shows that very little advantage is gained in this loss regime by choosing the larger TMD. In this case the direct inversion of the Fock state |1 breaks down when the displacement is close to 1.2 and 1.4 for 8-and 16-bin TMDs, respectively.

Case 3:
In the third simulation, we studied the reconstruction ranges for the single and two photon Fock states, assuming a currently feasible experimental situation in which the detection efficiency equals 30%. Once again, the scheme consists of a beam splitter with 95% transmittance, unbalanced TMD with 45/55 couplers but balanced deconvolution. At this high loss regime we can choose the 8-bin TMD for photon counting, since the ability to record higher photon numbers does not bring any advantage to the loss inversion. In this situation, the negative parts of the loss-degraded Wigner functions are completely washed out in the loss-degraded data. Nevertheless, the inner regions of the Wigner functions can still be reliably reconstructed up to displacements of 0.8 and 0.6 for the states |1 and |2 , respectively (Fig. 10). In the case of the Fock state |2 , the inversion causes negative probability components also close to the origin, indicating that the inversion is unstable. This is not caused by the contribution of higher photon numbers but by the accuracy at which the vacuum, one and two photon components are recorded. Nevertheless, the results appear to be accurate as the error in the parity measurement is seen to be negligible. Despite the low efficiency (30%), the simulation indicates that we are able to reconstruct the oscillatory behaviour of the photon number statistics for the Fock states |1 and |2 . The loss obscures the oscillations of the photon statistics, but similarly to the Wigner function analysis, the inversion technique enables us to reconstruct this characteristic behavior from the TMD measurement. In Fig. 11 we show this behavior, comparing the loss-degraded, inverted, and expected photon statistics at zero displacement and at the highest reliable displacement that can be applied to the quantum state. Case 4: Finally, we studied the reconstruction of the single photon Wigner function in the presence of imperfect mode overlap between the signal and reference fields. The statistics were generated according to Eq. (23). We now restricted the reconstruction range of the Wigner function to the region where probability components ρ nn > −0.003.
Other procedures, such as the error estimation, were kept the same. We utilised a beam splitter with 95% transmittance for the implementation of the displacement and a an 8-bin TMD with ideal 50/50 couplers and assumed that the single photon Fock state was affected by loss (efficiency equal to 30%) before applying a non-ideal displacement (M = 0.5). In Fig. 12 we show the photon number statistics for both the degraded and reconstructed signals at different values of the displacement. The quantum feature of photon number oscillation is washed out by the dual action of losses and imperfect overlap. Nevertheless, the loss can still be handled by applying an appropriate inversion method. In terms of displacement, the reconstruction is reliable to the boundary value of α = 1.0. As shown in Fig. 13, the Wigner function determined via the degraded signal does not display any quantum features. However, using the inversion of the statistics, the negative inner parts are recovered, although with a broadened width as a

Conclusions
We studied the reconstruction of Wigner functions of low-photon-number Fock states in a TMD-based photon counting experiment. Several parameters restrict the reliable reconstruction range of the Wigner function in phase space: the number of photons that a TMD can maximally resolve, the stability of the inversion with respect to the level of loss and the amount of statistics collected, and the imperfect mode overlap.With current technology, the dominant limitation is the efficiency achievable by the experimental realisation.
Performing Monte Carlo simulations we reconstructed the Wigner functions for single-photon and two-photon Fock states at different levels of detection efficiency. We showed that, under perfect mode-matching conditions, the Wigner function of the single-photon Fock state can be reconstructed with high confidence. The characteristic oscillatory behaviour present in the photon number distributions of displaced Fock states are extremely sensitive to losses, yet this interesting quantum feature can be reliable restored by applying a loss-inversion. A clear indication of phase space interference is given by the suppression of the one-photon component of the displaced single-photon Fock state.
Furthermore, we analysed the the situation of imperfect mode overlap. We showed how to employ the oscillatory behaviour as a tool to characterize the degree of mode overlap between the quantum signal and reference beam. Regarding the reconstruction of the Wigner function, even in a scenario of only 50% mode overlap, the simulation shows the possibility of recovering the negative values of the single-photon Fock state, although with a broader width.
Our results motivate the experimental exploration of non-Gaussian states with photon number resolving detectors. This kind of experiments offers the opportunity of loss-tolerant state characterization via the inversion of losses, along with the ability to probe the Wigner function point-by-point in phase space. Furthermore, the probing method allows advanced quantum state characterization that is intrinsically not only sensitive to a single mode but takes into account all the characteristics of the original quantum state and the coherent reference field. In our paper, we provided the tools to reliably characterize the prepared state and highlighted the different impacts of inefficient detection, convolution of statistics and mode mismatch. In contrast to standard homodyne detection, our analysis enables us to recognise the type of experimental imperfections, and gives valuable information about the extent of degradation caused by each one. Hence, our work constitutes a basis for a broad understanding and efficient manipulation of experimental data. We expect that our investigation will open a new route of detecting and studying quantum states and offer a new experimental technique for hybrid quantum information applications. statistics The result of Eq. (A.5) agrees with the independent analysis presented in Ref. [16] and corresponds to scaling the parity factor by the detection efficiency. We can apply this method for example on the loss-degraded statistics of the single photon Fock state derived in Eq. (14). The loss-inverted Wigner function can be written as