Pulsed homodyne Gaussian quantum tomography with low detection efficiency

Standard homodyne detection is an experimental method that is used for the reconstruction of quantum states of monochromatic light. In this framework, the quantum state is characterized by feeding to appropriate tomographic techniques the repeated measurement of a discrete set of field quadratures [1–3]. Quantum state reconstruction methods turn out to be of paramount importance for quantum information, for they can reveal the presence of quantum coherence and entanglement that might be used to perform tasks not possible in a classical setting [4].

For these experiments, a very high detection efficiency is typically required, along with ad-hoc designed apparatus [5]. However, new methods capable of discriminating between different quantum states of light, which demand lower detection capabilities typical of commercially available components, would facilitate the experimental realization of such apparatuses. Moreover, it would make it possible to apply quantum homodyne detection to study different physical systems where high noise conditions are unavoidable, such as out of equilibrium light matter dynamics [6]. For this purpose, non standard tomographic methods for reconstructing the state of light need to be adopted, specifically the so-called 'minimax' adaptive estimation of a quantum state [7, 8], allowing the full reconstruction of the associated Wigner function. These powerful statistical techniques are able to overcome the difficulties that arise when the more standard pattern function-based quantum tomography is used [9–16], allowing state reconstruction for arbitrary low detection efficiency, provided a sufficiently large number of data are collected.

In this article, we report on the commissioning and characterization of a time-domain homodyne detection apparatus working with coherent ultra-short light pulses, built using commercial detectors and operating in a regime of large electronic noise. By treating the shot-to-electronic-noise ratio (about 2 dB) as a detector inefficiency [17, 18], we obtain an overall detection efficiency of about 30%. By extending single-mode tomographic techniques to the pulsed regime, we show that even in these low efficiency conditions, it is still possible to discriminate between coherent and squeezed photon states with high accuracy.

The validity of the adopted minimax state reconstruction method has been tested by numerical experiments, using both coherent and squeezed states, and in real experiments for the tomographic reconstruction of coherent states with different mean photon number. These results provide concrete support for the effectiveness of the minimax adaptive state estimation even in low efficiency detection conditions, thus opening the way to wider applications of homodyne tomographic techniques to situations in quantum optics and non equilibrium spectroscopy, where the use of dedicated, highly efficient detectors is precluded.

In standard single-mode homodyne detection, the photon state under investigation, the signal, is mixed with a strong coherent reference state, the local oscillator (LO), by a $50/50$ beam splitter. The outputs are collected by two photodiodes and the difference photocurrent (homodyne photocurrent) is measured. It can be proven that when the LO is significantly more intense than the signal, the homodyne photocurrent is proportional to the signal quadrature [19]. Denoting by $\hat{a}\,$ and ${{\hat{a}\,}^{\dagger }}$ the single-mode annihilation and creation operators associated with the signal, the quadrature operator is defined as

$$\begin{eqnarray}&&{{\hat{x}\,}_{\Phi }}=\frac{\hat{a}\,{{e}^{-i\Phi }}+{{\hat{a}\,}^{\dagger }}{{e}^{i\Phi }}}{\sqrt{2}}\ ,\end{eqnarray} \tag{ 1 }$$

where $\Phi$ is the relative phase between the signal and the LO. The continuum set of quadratures with $\Phi \in \left[ 0,2\pi \right]$ provides a complete characterization of the signal state.

The standard techniques developed for the reconstruction of single-mode photon states need to be generalized when pulsed light is used [20–26]. In this case, the optical pulses are equally prepared multi mode coherent quantum states. For each single pulse a quadrature measurement is performed. In such regime, the single-mode description can still be used if the single-mode field operators are replaced by multi mode annihilation and creation operators. For this reason we make use in the following of the single-mode notation. A thorough theoretical treatment for pulsed homodyne detection is presented in the appendix.

2.1. Experimental setup

The scheme of the experimental set-up is shown in figure 1. The laser source is a mode-locked Ti:Sapphire oscillator with 80 MHz repetition rate. A beam splitter divides the incoming beam in two parts, which then interfere in a second beam splitter (NPBS in figure 1). The outputs are detected and subtracted by a commercial differential photodetector. The latter is made up of two Si/PIN photodiodes with nominal quantum efficiency ${{\eta }_{pd}}=0.85$ at 800 nm wavelength and linear response up to 0.6 mW LO power. The detector subtraction efficiency is quantified by the common mode rejection ratio (CMRR), defined as the ratio between the detector output power when both photodiodes are illuminated and the power when one of the two is screened [25]. For the present experiment CMRR $>\;36$ dB.

Zoom In Zoom Out Reset image size

The homodyne photocurrent is recorded by a digital oscilloscope with a bandwidth of 500 MHz and a sampling rate of 5 GSamples ${{\text{s}}^{-1}}$. The digitized output is numerically integrated over time intervals corresponding to the duration of the pulse. Each integral is associated with a single quadrature measurement.

In the shot noise regime with the signal beam blocked, i.e., with the signal in the vacuum state, the homodyne detector noise variance is expected to change linearly with the LO power on the top of a constant offset representing the electronic noise [27]. Figure 2 shows the detector noise variance of $8\times {{10}^{3}}$ pulses, for different values of the LO power. The noise variance grows linearly up to 0.6 mW LO power, instead, at higher powers, the photodiodes non linear effects are significant. To have the maximum shot-to-electronic noise ratio ($\approx 2$ dB) achievable in the linear regime, the experiments have been performed at 0.6 mW LO power.

Zoom In Zoom Out Reset image size

Conversely, when the signal beam is not blocked, it can be attenuated with respect to the LO by neutral density filter (${{F}_{1}}$ in figure 1). The phase difference $\Phi$ between the two arms of the interferometer can be modulated in the $\left[ 0,2\pi \right]$ range using a piezoelectric translator in the LO arm (PZT in figure 1). The homodyne photocurrent is acquired by the digital oscilloscope, and each difference pulse is integrated.

2.2. Homodyne traces

In the absence of electronic noise, the voltage V corresponding to the homodyne photocurrent, namely the experimentally accessible quantity, is proportional to the quadrature operator ${{x}_{\Phi }}$, $V=\gamma \;{{x}_{\Phi }}$, with an appropriate constant γ. The electronic noise can be generically described by a classical stochastic process δ that can be assumed to be Gaussian distributed, with zero average and variance $T/2$. The value of δ must be experimentally measured. Under these conditions the electronic noise is independent from the quadrature, and it contributes to the homodyne voltage as

$$\begin{eqnarray}&&V=\gamma \;{{x}_{\Phi }}+\delta \ .\end{eqnarray} \tag{ 2 }$$

A homodyne trace is obtained by collecting a set of homodyne voltage values ${{V}_{i}}$, corresponding to different phase values ${{\Phi }_{i}}$, associated with a large number M of piezo positions.

In the case of the vacuum state (absence of signal), these considerations are outlined in figure 2. The total vacuum homodyne variance for the chosen LO power is indeed composed of two independent contributions. An intrinsic contribution ($A={{\gamma }^{2}}/2$) is actually the shot noise, while an extrinsic contribution ($B=T/2$) is due to the electronic noise. Following [18], the electronic noise can be treated as an optical loss channel with an equivalent transmission efficiency given by

$$\begin{eqnarray}&&{{\eta }_{eq}}=\frac{A}{A+B}=\frac{{{\gamma }^{2}}}{\gamma {{\prime }^{2}}}\ ,\end{eqnarray} \tag{ 3 }$$

where $\gamma \prime =\sqrt{{{\gamma }^{2}}+T}$. Here we estimate ${{\eta }_{eq}}=0.37$. Thus the total apparatus efficiency is $\eta ={{\eta }_{eq}}\;{{\eta }_{pd}}=0.31$.

Considering that the variances of independent stochastic variables are additive, $\gamma \prime$ can be determined by using the vacuum state as reference and by assuming the variance of the quadrature operator in the vacuum to be $1/2$,

$$\begin{eqnarray}&&\gamma \prime =\sqrt{2{\left\langle{V}^{2}\right\rangle}_{0}}\ ,\end{eqnarray} \tag{ 4 }$$

where ${\left\langle{V}^{2}\right\rangle}_{0}$ is the experimental voltage variance for the vacuum state.

To consistently apply the quantum state reconstruction to the collected experimental homodyne data, it is convenient to rescale the raw data ${{V}_{i}}$ by the constant $\gamma \prime$,

$$\begin{eqnarray}&&{{Y}_{i}}={{V}_{i}}/\gamma \prime \ ,\end{eqnarray} \tag{ 5 }$$

so that the new, calibrated quantities ${{Y}_{i}}$ have variance $1/2$. Figure 3 shows the case of two optical coherent states with different mean photon numbers, each homodyne trace consisting of $M=8\times {{10}^{4}}$ experimental data. The data were collected using two different optical density for the signal attenuator (see ${{F}_{1}}$ in figure 1). Clearly even with a very low mean photon number signal (see figure 3(b)), the phase modulation of the homodyne trace is detectable.

Zoom In Zoom Out Reset image size

In the following, we prove that a minimax estimation of the Wigner function [7] can discriminate between quantum coherent and squeezed states even under the low efficiency condition of our experiment.

3.1. Minimax estimation of the Wigner function

Using the annihilation and creation operators $\hat{a}\,,{{\hat{a}\,}^{\dagger }}$ one constructs position, respectively momentum-like operators, $\hat{q}\,=\left( \hat{a}\,+{{\hat{a}\,}^{\dagger }} \right)/\sqrt{2}$, respectively $\hat{p}\,=\left( \hat{a}\,-{{\hat{a}\,}^{\dagger }} \right)/\left( i\sqrt{2} \right)$. Given a density matrix $\hat{\rho}\,$ representing a generic photon state, the associated Wigner function ${{W}_{\rho }}\left( q,p \right)$ in the variables (q, p) is defined by

$$\begin{eqnarray}&&{{W}_{\rho }}\left( q,p \right)=\frac{1}{{{\left( 2\pi \right)}^{2}}}{{\mathop{\int }^{}}_{{{\mathbb{R}}^{2}}}}\text{d}u\;\text{d}v\;{{\text{e}}^{i\left( uq+vp \right)}}\;\text{Tr}\left[ \hat{\rho}\,\;{{\text{e}}^{-i\left( u\hat{q}\,+v\hat{p}\, \right)}} \right].\end{eqnarray} \tag{ 6 }$$

Then, the mean value of any observable $\hat{O}\,$, expressible as a function of $\hat{a}\,$ and ${{\hat{a}\,}^{\dagger }}$, reads

$$\begin{eqnarray}&&{{\left\langle\hat{O}\right\rangle\,}_{\rho }}=\text{Tr}\left( \hat{\rho}\,\;\hat{O}\, \right)=2\pi \;{{\mathop{\int }^{}}_{{{\mathbb{R}}^{2}}}}dq\;dp\;{{W}_{\rho }}\left( q,p \right)O\left( q,p \right),\end{eqnarray} \tag{ 7 }$$

where the function $O\left( q,p \right)$ is connected to the operator $\hat{O}\,$ as the Wigner function is to $\hat{\rho}\,$ as in (6). The Wigner function ${{W}_{\rho }}\left( q,p \right)$ is therefore sufficient to evaluate the average value of any observable $\hat{O}\,$ with respect to any state $\hat{\rho}\,$.

The minimax reconstruction algorithm precisely gives an estimator for the Wigner function, allowing its reconstruction from the collected homodyne data ${{Y}_{i}}$, provided a proper calibration is used (equation (5)). The estimator is given by

$$\begin{eqnarray}&&W_{h,M}^{\eta }\left( q,p \right)=\frac{1}{M} \mathop\sum\limits^{M}_{i=1}\,K_{{{Y}_{i}},{{\Phi }_{i}}}^{h,\eta }\left( q,p \right),\end{eqnarray} \tag{ 8 }$$

and it depends on the overall detector efficiency η, the number of collected data M and an adaptive parameter h [7]. The kernel function $K_{{{Y}_{i}},{{\Phi }_{i}}}^{h,\eta }\left( q,p \right)$ is explicitly given by

$$\begin{eqnarray}&&K_{{{Y}_{i}},{{\Phi }_{i}}}^{h,\eta }\left( q,p \right)=\mathop{\int }^{}_{-\frac{1}{h}}^{\frac{1}{h}}dt\frac{|t|}{4\pi }{{e}^{-it\left( q\cos {{\Phi }_{i}}+p\sin {{\Phi }_{i}}-\frac{{{Y}_{i}}}{\sqrt{\eta }} \right)+{{t}^{2}}\frac{1-\eta }{4\eta }}}.\end{eqnarray} \tag{ 9 }$$

The truncation of the integration due to the parameter h is necessary because of the diverging inverse Gaussian present in the integral. In [7] it is shown that the optimal adaptive estimator, i.e., the one that minimizes uncertainties, is obtained when $h={{h}_{adap}}$, with

$$\begin{eqnarray}&&{{h}_{adap}}={{\left( \frac{2\eta \log M}{1-\eta }-\sqrt{\frac{2\eta \log M}{1-\eta }} \right)}^{-1/2}}.\end{eqnarray} \tag{ 10 }$$

Notice that the optimally adapted cutoff depends both on the detector efficiency η and the number M of collected data. When the efficiency becomes small, in order to have a good estimator of the Wigner function, a larger number M of data is needed. However, we stress that no inferior bound on η exists that prevents the convergence of the algorithm ⁵ .

Once the Wigner function is reconstructed, the expectation value of any observable $\hat{O}\,$ of the system can be estimated as follows,

$$\begin{eqnarray}&&E_{\hat{O}\,}^{h}={{\mathop{\int }}_{{{\mathbb{R}}^{2}}}}dq\;dp\ O\left( q,p \right)\;W_{h,M}^{\eta }\left( q,p \right)=\frac{1}{M}\mathop\sum\limits^{M}_{i=1}\,R_{\hat{O}\,}^{h,\eta }\left( {{Y}_{i}},{{\Phi }_{i}} \right),\end{eqnarray} \tag{ 11 }$$

where the corresponding kernel function is given by

$$\begin{eqnarray}&&R_{\hat{O}\,}^{h,\eta }\left( {{Y}_{i}},{{\Phi }_{i}} \right)={{\mathop{\int }^{}}_{{{\mathbb{R}}^{2}}}}dq\;dp\ O\left( q,p \right)\;K_{{{Y}_{i}},{{\Phi }_{i}}}^{h,\eta }\left( q,p \right).\end{eqnarray} \tag{ 12 }$$

3.2. Numerical experiments

For testing the effectiveness of this approach under our experimental conditions, we first analyze sets of numerically generated data that simulate the quadratures of a known pure quantum state $\left| \psi \right\rangle$. By knowing in advance the quantum state, it is possible to verify whether the minimax technique allows a proper reconstruction.

The generation of the fictitious data starts by using the quadrature probability distribution with $\eta =1$ associated with a quantum state, and afterwards in adding to each numerically generated state quadrature a Gaussian noise that exactly simulates the electronic noise associated with the efficiency, i.e., $\eta =0.31$, of our set up. The analysis of the simulated data proceeds, as for experimental data, through a calibration using the vacuum noise as reference. In particular, we perform two numerical experiments. The first with a coherent state $\left| \psi \right\rangle=D\left( \alpha \right)\left| 0 \right\rangle$ and another one with a displaced-squeezed state $\left| \psi \right\rangle=D\left( \alpha \right)S\left( \xi \right)\left| 0 \right\rangle$, where

$$\begin{eqnarray}&&D\left( \alpha \right)={{e}^{\alpha {{\hat{a}\,}^{\dagger }}-{{\alpha }^{*}}\hat{a}\,}},\qquad \ S\left( \xi \right)={{e}^{1/2\left( \xi {{\hat{a}\,}^{\dagger 2}}-{{\xi }^{*}}{{\hat{a}\,}^{2}} \right)}},\end{eqnarray} \tag{ 13 }$$

are the displacement and squeezing operators, respectively. Each data set consists of $M=8\times {{10}^{4}}$ quadrature measurements with the relative phase $\Phi$ ranging in the interval $\left[ 0,2\pi \right]$.

Following the algorithm sketched before, we reconstruct the Wigner function of the two quantum states. They are shown in figures 4 (a) and (b) for the coherent and the squeezed state, respectively. To remove the artifacts resulting by the numerical integration, the reported Wigner functions have been filtered by an image processing algorithm (low-pass Gaussian convolution filtering). The blue and red curves are bidimensional cuts of the Wigner functions in correspondence of the expectation value of the position and momentum operator, respectively.

Zoom In Zoom Out Reset image size

In order to assess the similarity between the reconstructed state ${{\rho }_{r}}$, corresponding to the reconstructed Wigner function ${{W}_{r}}\left( q,p \right)$, and the expected state, $\left| \psi \right\rangle$, with Wigner function ${{W}_{\psi }}\left( q,p \right)$, we use as figure of merit, the fidelity $F=\left\langle \psi \right|{{\rho }_{r}}\left| \psi \right\rangle$. This quantity essentially measures the overlap between the two quantum states; indeed, $0\leqslant F \leqslant 1$, reaching the value 1 only when the two states coincide [4]. In our case, this measure is given by

$$\begin{eqnarray}&&F=2\pi \;{{\mathop{\int }^{}}_{{{\mathbb{R}}^{2}}}}dq\;dp\ {{W}_{\psi }}\left( q,p \right)\;{{W}_{r}}\left( q,p \right).\end{eqnarray} \tag{ 14 }$$

We obtained a fidelity of 0.97 and 0.92 for the coherent and the squeezed state, respectively, despite the large noise affecting the reconstruction ⁶ . These results, together with the actual values of the quadrature variances, prove that the squeezed and coherent states can be discriminated under high noise conditions. Indeed, the features of the reconstructed Wigner functions clearly reflect the different nature of the two quantum states.

From the reconstructed Wigner function, it is possible to compute the expectation values of relevant observables with respect to the coherent and the displaced-squeezed state, in particular, the number operator $\hat{n}\,={{\hat{a}\,}^{\dagger }}\hat{a}\,$ and the position $\hat{q}\,$ and momentum $\hat{p}\,$ operators. For $\hat{q}\,$ and $\hat{p}\,$ it also possible to derive the variances, ${{\sigma }^{2}}\left[ \hat{q}\, \right]$ and ${{\sigma }^{2}}\left[ \hat{p}\, \right]$, and to estimate the squeezing parameter,

$$\begin{eqnarray}&&\xi =\frac{1}{4}\ln \left( \sigma _{\hat{p}\,}^{2}/\sigma _{\hat{q}\,}^{2} \right).\end{eqnarray} \tag{ 15 }$$

The results are summarized in table 1. For the errors evaluation we have used the standard expression for the mean average error relative to a data set $\left( {{Y}_{i}},{{\Phi }_{i}} \right)$ (recall equations (11) and (12)),

$$\begin{eqnarray}&&{{\epsilon }_{\hat{O}\,}}=\sqrt{\frac{\mathop\sum\nolimits^{M}_{i=1}\,{{\left[ R_{\hat{O}\,}^{h,\eta }\left( {{Y}_{i}};{{\Phi }_{i}} \right)-E_{\hat{O}\,}^{h} \right]}^{2}}}{M\left( M-1 \right)}}\ .\end{eqnarray} \tag{ 16 }$$

The averages computed using the Wigner function estimated from the data set fully agree with those analytically calculated from the known quantum states. In particular, even the estimated value of the squeezing parameter agrees with the expected one. This explicitly shows that the minimax adaptive estimation algorithm results are highly effective in reconstructing quantum states of light, and more specifically, in discriminating different Gaussian quantum states even under low-efficiency conditions.

Table 1.  Estimate of different expectation values for $M=8\times {{10}^{4}}$ simulated quadratures associated to the coherent state ($\alpha =3$) and the displaced-squeezed state ($\alpha =3$, $\xi =0.8$)

$\hat{O}\,$	${{\langle \hat{O}\,\rangle }_{D\left( \alpha \right)|0\rangle }}$	$E_{\hat{O}\,}^{h}$	${{\langle \hat{O}\,\rangle }_{D\left( \alpha \right)S\left( \xi \right)\beta 0}}$	$E_{\hat{O}\,}^{h}$
$\hat{n}\,$	9	$9\pm 1$	9.89	$9.9\pm 0.8$
$\hat{q}\,$	4.24	$4.2\pm 0.5$	4.24	$4.3\pm 0.3$
$\hat{p}\,$	0	$0.00\pm 0.06$	0	$0.0\pm 0.2$
$\sigma _{\hat{q}\,}^{2}$	0.5	$0.4\pm 0.2$	0.1	$0.1\pm 0.2$
$\sigma _{\hat{p}\,}^{2}$	0.5	$0.5\pm 0.2$	2.48	$2.3\pm 0.4$
ξ	0	$0.0\pm 0.1$	0.8	$0.78\pm 0.02$

3.3. Real experiments

The real data in figure 3 can now be treated in the same way. In this case, the collected homodyne traces are two coherent states with very different mean photon number.

The first step in the reconstruction procedure and determination of relevant observables is to use the estimator (8) to obtain the Wigner functions associated with the two states. The results are shown in figure 5; as in the previous numerical experiments, in order to minimize artifacts related to the numerical integration in (9), a low-pass Gaussian filter has been applied to the raw images.

Zoom In Zoom Out Reset image size

With these results, one can then estimate the mean values of relevant observables, as number, position, and momentum operators together with the squeezing parameter, obtained from the computation of the variances of position and momentum. The results are summarized in table 2. The obtained expectation values for the number operator explicitly exhibit the three order of magnitude difference between the two states. Furthermore, these results show that the analyzed experimental states are indeed minimum uncertainty states, with vanishing squeezing parameter.

Table 2.  Estimate of different expectation values for the experimental coherent states with optical density OD = 4.5 and OD = 7, respectively.

$\hat{O}\,$	$E_{\hat{O}\,}^{h}$ (OD = 4.5)	$E_{\hat{O}\,}^{h}$ (OD = 7)
$\hat{n}\,$	$51\pm 7$	$0.2\pm 0.3$
$\hat{q}\,$	$10\pm 1$	$0.7\pm 0.2$
$\hat{p}\,$	$0.1\pm 0.2$	$0.0\pm 0.2$
$\sigma _{\hat{q}\,}^{2}$	$0.7\pm 0.3$	$0.5\pm 0.2$
$\sigma _{\hat{p}\,}^{2}$	$0.6\pm 0.3$	$0.5\pm 0.3$
ξ	$0.0\pm 0.2$	$0.0\pm 0.2$

The previous results demonstrate the effectiveness of the minimax statistical techniques in the reconstruction of quantum states of light in presence of large electronic noise. These methods allow us to circumvent convergence problems that arise when using standard pattern function methods in estimating observable averages with low-detection efficiency, namely $\eta <50%$ [9–15]. Nevertheless, it is interesting to investigate the possible relations between these two approaches by replacing the cutoff truncation in the integration with respect to t in (9) with a suitable Gaussian regularization. In particular, we consider the following alternate kernel functions instead of $R_{\hat{O}\,}^{h,\eta }\left( {{Y}_{i}},{{\Phi }_{i}} \right)$,

$$\begin{eqnarray}&&R_{\hat{O}\,}^{\epsilon }\left( {{Y}_{i}},{{\Phi }_{i}} \right)={{\mathop{\int }^{}}_{{{\mathbb{R}}^{2}}}}dq\;dp\;\hat{O}\,\left( q,p \right){{e}^{-\epsilon \left( {{q}^{2}}+{{p}^{2}} \right)}}\mathop{\int }^{}_{-\infty }^{+\infty }dt\frac{|t|}{4\pi }{{e}^{-\text{i}t\left( q\cos {{\Phi }_{i}}+p\sin {{\Phi }_{i}}-\frac{{{Y}_{i}}}{\sqrt{\eta }} \right)}}\ {{e}^{{{t}^{2}}\frac{1-\eta }{4\eta }}}\ ,\end{eqnarray} \tag{ 17 }$$

where ε is a positive regularization parameter. In principle, this allows for a more direct estimate of the expectation values of any observable, provided $0<\epsilon <\eta /\left( 1-\eta \right)$. For instance, in the case of position, momentum, and number operators, one finds:

$$\begin{eqnarray}&&R_{\hat{q}\,}^{\epsilon }\left( {{Y}_{i}},{{\Phi }_{i}} \right)=\mathop{\int }^{}_{-\infty }^{+\infty }dt\;\frac{\cos {{\Phi }_{i}}}{8{{\epsilon }^{2}}}\;t\;|t|\sin \left( \frac{t\;{{Y}_{i}}}{\sqrt{\eta }} \right){{e}^{-{{t}^{2}}\,\frac{\eta -\epsilon \left( 1-\eta \right)}{4\eta \epsilon }}}\ ,\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&R_{\hat{p}\,}^{\epsilon }\left( {{Y}_{i}},{{\Phi }_{i}} \right)=\mathop{\int }^{}_{-\infty }^{+\infty }dt\;\frac{\sin {{\Phi }_{i}}}{8{{\epsilon }^{2}}}\;t\;|t|\;\sin \left( \frac{t\;{{Y}_{i}}}{\sqrt{\eta }} \right){{e}^{-{{t}^{2}}\,\frac{\eta -\epsilon \left( 1-\eta \right)}{4\eta \epsilon }}}\ ,\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&R_{\hat{n}\,}^{\epsilon }\left( {{Y}_{i}},{{\Phi }_{i}} \right)=\mathop{\int }^{}_{-\infty }^{+\infty }dt\;\frac{1}{32{{\epsilon }^{3}}}\;|t|\;\cos \left( \frac{t\;{{Y}_{i}}}{\sqrt{\eta }} \right){{e}^{-{{t}^{2}}\,\frac{\eta -\epsilon \left( 1-\eta \right)}{4\eta \epsilon }}}\left[ 4\epsilon \left( 1-\epsilon \right)-{{t}^{2}} \right]\ .\end{eqnarray} \tag{ 20 }$$

Notice that these quantities converge to functions $R_{\hat{O}\,}^{0}\left( {{Y}_{i}},{{\Phi }_{i}} \right)$ in the limit $\epsilon \to 0$:

$$\begin{eqnarray}&&R_{\hat{q}\,}^{0}\left( {{Y}_{i}},{{\Phi }_{i}} \right)=\frac{2}{\sqrt{\eta }}{{Y}_{i}}\cos {{\Phi }_{i}}\ ,\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&R_{\hat{p}\,}^{0}\left( {{Y}_{i}},{{\Phi }_{i}} \right)=\frac{2}{\sqrt{\eta }}{{Y}_{i}}\sin {{\Phi }_{i}}\ ,\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&R_{\hat{n}\,}^{0}\left( {{Y}_{i}},{{\Phi }_{i}} \right)={{\left( \frac{{{Y}_{i}}}{\sqrt{\eta }} \right)}^{2}}-\frac{1}{2\eta }\ .\end{eqnarray} \tag{ 23 }$$

These expressions coincide with the so-called kernel functions of the standard pattern function-based tomographic techniques [15]. Unfortunately, the limit $\epsilon \to 0$ cannot be taken in the case of more general observables with higher powers of $\hat{q}\,$ and $\hat{p}\,$ since the integrals diverge as powers of $1/\epsilon$. This fact may be related to the mentioned convergence problems of the standard pattern function method.

It is important to note that the estimation of the mean values of $\hat{q}\,$, $\hat{p}\,$, and $\hat{n}\,$ through the use of the limit kernel functions (21)–(23) is not sufficient to completely characterize quantum states, not even Gaussian ones, like coherent and squeezed states. Instead, the minimax techniques allow the full reconstruction of the state through its Wigner function and in principle for arbitrarily low detector efficiency η. However, as already mentioned in section 3.1, as η becomes small, a larger set of experimental data points is needed in order to have a reliable estimation of the Wigner function. This condition may de facto limit the applicability of the minimax adaptive state reconstruction techniques in special experimental instances, where a very high repetition rate state preparation is unattainable.

Here we report the characterization of a time domain balanced homodyne detection apparatus operating in the presence of large electronic noise corresponding to an overall detection efficiency $\eta =0.31$. We used the detector combined with tomographic reconstruction techniques for discriminating between different quantum states of light. A minimax adaptive reconstruction of the Wigner function has been adopted. This approach allows us to circumvent possible convergence problems arising from low detector efficiency in the standard pattern function-based quantum tomography.

The effectiveness of such a method has been verified in two ways. At first, we calculated the Wigner function for simulated data of coherent and squeezed states and then for real experimental homodyne data of coherent states with different mean photon numbers. In all cases it is proved that it is possible to efficiently reconstruct the associated Wigner function, asserting the Gaussian character of the quantum states, and evaluating their relevant parameters.

The present study demonstrates that even low efficient ($\sim 30%$) homodyne detectors can be usefully employed to study the nature of quantum states of light, provided that non-standard statistical tools as the minimax methods are used to reconstruct their Wigner functions. These results may be important whenever quantum optics techniques are employed to investigate the dynamics of out-of-equilibrium states and the presence of quantum coherence in condensed matter.

Useful discussions with Simone Cialdi and Marco Barbieri are acknowledged. This work has been supported by MIUR (FIRB 'LiCHIS'-RBFR10YQ3H), University of Trieste ('FRA 2009' and 'FRA 2011') and European Union, Seventh Framework Programme, under the project GO FAST, grant agreement no. 280555.

Homodyne detection in a pulsed regime requires a formal generalization of its theoretical description with respect to the one-mode regime; indeed, the LO and the signal at the beam splitter are not monochromatic. Classically, the electric field of a pulsed laser beam can be represented as a mode-locked superposition of amplitudes:

$$\begin{eqnarray}&&E\left( t \right)=\sum\limits^{M}_{l=-M}\,|{{\alpha }_{l}}|\;{{e}^{i\,{{\Phi }_{l}}\left( t \right)}},\qquad {{\Phi }_{l}}\left( t \right)={{\omega }_{l}}\;t+{{\varphi }_{l}},\end{eqnarray} \tag{ A.1 }$$

where the phases ${{\Phi }_{l}}\left( t \right)$ are mode-locked by the condition ${{\varphi }_{l}}=l\;{{\varphi }_{0}}$, with ${{\varphi }_{0}}$ a reference phase. Each frequency ${{\omega }_{l}}$ contributes to the field with the amplitude ${{\alpha }_{l}}=\left| {{\alpha }_{l}} \right|\;{{e}^{i{{\Phi }_{l}}}}$ and the number of contributing frequencies depends on the shape of the pulse. Quantized pulsed laser light is described by associating to each monochromatic component a coherent state $\left| {{\alpha }_{l}} \right\rangle$, that is an eigenstate of the annihilation operator ${{\hat{a}\,}_{l}}$ of photons in the mode of frequency ${{\omega }_{l}}$, ${{\hat{a}\,}_{l}}\left| {{\alpha }_{l}} \right\rangle={{\alpha }_{l}}\left| {{\alpha }_{l}} \right\rangle$, and to the entire pulse the tensor product

$$\begin{eqnarray}&&\left| {\bar{\alpha }} \right\rangle=\otimes _{l=-M}^{M}\left| {{\alpha }_{l}} \right\rangle,\end{eqnarray} \tag{ A.2 }$$

where $\bar{\alpha }$ is the vector whose components are the amplitudes ${{\alpha }_{l}}$. By means of the creation and annihilation operators ${{\hat{a}\,}_{l}}$ and $\hat{a}\,_{l}^{\dagger }$ each monochromatic coherent state reads

$$\begin{eqnarray}&&\left| {{\alpha }_{l}} \right\rangle=D\left( {{\alpha }_{l}} \right)\left| 0 \right\rangle,\qquad D\left( {{\alpha }_{l}} \right)={{e}^{{{\alpha }_{l}}\hat{a}\,_{l}^{\dagger }-\alpha _{l}^{*}{{\hat{a}\,}_{l}}}},\end{eqnarray} \tag{ A.3 }$$

where $\left| 0 \right\rangle$ is the vacuum state and $D\left( \alpha \right)$ is the so-called displacement operator.

Since the creation and annihilation operators pertaining to different modes commute, the pulsed coherent state (A.2) can be conveniently recast as

$$\begin{eqnarray}&&\left| {\bar{\alpha }} \right\rangle=D\left( {\bar{\alpha }} \right)\left| 0 \right\rangle,\qquad D\left( {\bar{\alpha }} \right)={{e}^{{{\hat{A}\,}^{\dagger }}\left( {\bar{\alpha }} \right)-\hat{A}\,\left( {\bar{\alpha }} \right)}},\end{eqnarray} \tag{ A.4 }$$

by means of a displacement operator $D\left( {\bar{\alpha }} \right)$ expressed in terms of multi mode operators

$$\begin{eqnarray}&&{{\hat{A}\,}^{\dagger }}\left( {\bar{\alpha }} \right)=\underset{l}{\mathop{\mathop{\sum }^{}}}\,{{\alpha }_{l}}\;\hat{a}\,_{l}^{\dagger },\qquad \hat{A}\,\left( {\bar{\alpha }} \right)=\underset{l}{\mathop{\mathop{\sum }^{}}}\,\alpha _{l}^{*}{{\hat{a}\,}_{l}}.\end{eqnarray} \tag{ A.5 }$$

The reason for labeling the pulsed coherent state by $\left| {\bar{\alpha }} \right\rangle$ can now be easily understood. The state of one photon of frequency ${{\omega }_{l}}$ is given by $\hat{a}\,_{l}^{\dagger }\left| 0 \right\rangle=\left| {{1}_{l}} \right\rangle$, while a generic non-monochromatic superposition of frequencies ${{\omega }_{l}}$ with amplitudes ${{\alpha }_{l}}$ corresponds to the state $\left| {{1}_{{\bar{\alpha }}}} \right\rangle={{\mathop{\sum }^{}}_{l}}{{\alpha }_{l}}\left| {{1}_{l}} \right\rangle$ that results from applying the multimode operator ${{\hat{A}\,}^{\dagger }}\left( {\bar{\alpha }} \right)$ to the vacuum state:

$$\begin{eqnarray}&&{{\hat{A}\,}^{\dagger }}\left( {\bar{\alpha }} \right)\left| 0 \right\rangle=\underset{l}{\mathop{\mathop{\sum }^{}}}\,{{\alpha }_{l}}\;\hat{a}\,_{l}^{\dagger }\left| 0 \right\rangle=\underset{l}{\mathop{\mathop{\sum }^{}}}\,{{\alpha }_{l}}\;\left| {{1}_{l}} \right\rangle=\left| {{1}_{{\bar{\alpha }}}} \right\rangle .\end{eqnarray} \tag{ A.6 }$$

Therefore, ${{\hat{A}\,}^{\dagger }}\left( {\bar{\alpha }} \right)$ is the creator operator of a single photon in the (not normalized) superposition state $\left| {{1}_{{\bar{\alpha }}}} \right\rangle$, while $\hat{A}\,\left( {\bar{\alpha }} \right)$ destroys a photon in the same state; thus, the quantum state of the pulsed laser is a coherent state associated not with a single amplitude ${{\alpha }_{l}}$, but with the vector $\bar{\alpha }$ of all the amplitudes contributing to the pulse: in other words, we have a Poissonian distribution not with respect to the number of photons in a monochromatic wave, but to the number of photons in the superposition $\left| {{1}_{{\bar{\alpha }}}} \right\rangle$.

The normalized operators

$$\begin{eqnarray}&&\hat{A}\,=\frac{\hat{A}\,\left( {\bar{\alpha }} \right)}{|\bar{\alpha }|},\qquad {{\hat{A}\,}^{\dagger }}=\frac{\hat{A}\,\left( {\bar{\alpha }} \right)}{|\bar{\alpha }|},\end{eqnarray} \tag{ A.7 }$$

where $|\bar{\alpha }{{|}^{2}}=\left\langle{{1}_{{\bar{\alpha }}}}|{{1}_{{\bar{\alpha }}}}\right\rangle$, satisfy the canonical commutation relations

$$\begin{eqnarray}&&\left[ \hat{A}\,,{{\hat{A}\,}^{\dagger }} \right]=\frac{1}{|\bar{\alpha }{{|}^{2}}}\underset{ij}{\mathop{\mathop{\sum }^{}}}\,\alpha _{i}^{*}{{\alpha }_{j}}\;\left[ {{\hat{a}\,}_{i}},\hat{a}\,_{j}^{\dagger } \right]=1.\end{eqnarray} \tag{ A.8 }$$

By means of them we can now extend the state-tomography techniques to the case in which the signal and the LO are pulsed. In the monochromatic case, when the signal mode a interferes with the LO mode b on the beam splitter, the photo-current operator is $\hat{I}\,={{\hat{a}\,}^{\dagger }}\hat{b}\,+{{\hat{b}\,}^{\dagger }}\hat{a}\,$ [19].

If more frequencies ${{\omega }_{l}}$ are present in both the LO and in the signal, each one of the corresponding mode operators will be subjected to the beam splitting transformation and the detectors will ideally register photons of all involved frequencies. Then, the photo-current operator becomes $\hat{I}\,={{\mathop{\sum }^{}}_{l}}\hat{a}\,_{l}^{\dagger }{{\hat{b}\,}_{l}}+\hat{b}\,_{l}^{\dagger }{{\hat{a}\,}_{l}}$.

If the LO is in a pulsed coherent state $\left| {\bar{z}} \right\rangle={{e}^{{{\hat{B}\,}^{\dagger }}\left( {\bar{z}} \right)-\hat{B}\,\left( {\bar{z}} \right)}}\left| 0 \right\rangle$, with generalized creation and annihilation operators ${{\hat{B}\,}^{\dagger }}\left( {\bar{z}} \right)={{\mathop{\sum }^{}}_{l}}{{z}_{l}}\hat{b}\,_{l}^{\dagger }$ and $\hat{B}\,\left( {\bar{z}} \right)={{\mathop{\sum }^{}}_{l}}z_{l}^{*}{{\hat{b}\,}_{l}}$, the phase difference $\Phi$ between the LO and the signal is changed by the action of the piezoelectric translator placed in the LO arm on all the LO modes:

$$\begin{eqnarray}&&{{\hat{b}\,}_{l}}\to {{\hat{b}\,}_{l}}\;{{e}^{i\Phi }}\ ,\qquad \hat{b}\,_{l}^{\dagger }\to \hat{b}\,_{l}^{\dagger }\;{{e}^{-i\Phi }}\ .\end{eqnarray} \tag{ A.9 }$$

The photo-current operator, which is measured by the pulsed homodyne setup, is thus given by

$$\begin{eqnarray}&&{{\hat{I}\,}_{\Phi }}=\underset{l}{\mathop{\mathop{\sum }^{}}}\,\left( \hat{a}\,_{l}^{\dagger }\;{{\hat{b}\,}_{l}}\;\;{{e}^{i\Phi }}\;+\;{{\hat{a}\,}_{l}}\;\hat{b}\,_{l}^{\dagger }\;{{e}^{-i\Phi }} \right)\ .\end{eqnarray} \tag{ A.10 }$$

Let ${{\hat{\rho}\,}_{s}}$ denote the quantum state (density matrix) of the signal field and by $\left| {\bar{z}} \right\rangle\left\langle {\bar{z}} \right|$ the projector onto the (coherent) state of the incoming pulse. Then, using that ${{\hat{b}\,}_{l}}\left| {\bar{z}} \right\rangle={{z}_{l}}\left| {\bar{z}} \right\rangle$ and the expressions in (A.5), the expectation value ${{I}_{\Phi }}=\text{Tr}\left( {{\rho }_{s}}\otimes \left| {\bar{z}} \right\rangle\left\langle {\bar{z}} \right|\;{{\hat{I}\,}_{\Phi }} \right)$ of the photo-current is calculated as follows:

$$\begin{eqnarray}\begin{matrix} {{I}_{\Phi }} & = & \underset{l}{\mathop{\mathop{\sum }^{}}}\,\left( \text{Tr}\left[ {{\hat{\rho}\,}_{s}}\;\hat{a}\,_{l}^{\dagger } \right]\cdot \left\langle\bar{z}|{{\hat{b}\,}_{l}}|\bar{z}\right\rangle\;{{e}^{i\Phi }}+h.c. \right) \\ {} & = & \underset{l}{\mathop{\mathop{\sum }^{}}}\,\left( \text{Tr}\left[ {{\hat{\rho}\,}_{s}}\;\left( \hat{a}\,_{l}^{\dagger }\;{{z}_{l}}\;{{e}^{i\Phi }}+{{\hat{a}\,}_{l}}\;z_{l}^{*}\;{{e}^{-i\Phi }} \right) \right] \right) \\ {} & = & \text{Tr}\left[ {{\hat{\rho}\,}_{s}}\left( {{\hat{A}\,}^{\dagger }}\left( {\bar{z}} \right)\;{{e}^{i\Phi }}+\hat{A}\,\left( {\bar{z}} \right)\;{{e}^{-i\Phi }} \right) \right]. \\ &&\end{matrix}\end{eqnarray} \tag{ A.11 }$$

Using (A.7) and comparing these results in the standard treatment of homodyne detection for the single-mode case, realize that in the pulsed regime one measures a quantity ${{I}_{\Phi }}$,

$$\begin{eqnarray}&&{{I}_{\Phi }}=\sqrt{2}\;|\bar{z}|\;\text{Tr}\left[ {{\hat{\rho}\,}_{s}}\;{{\hat{X}\,}_{\Phi }} \right]\ ,\end{eqnarray} \tag{ A.12 }$$

proportional to a quadrature, which generalizes that in (1):

$$\begin{eqnarray}&&{{\hat{X}\,}_{\Phi }}=\frac{{{\hat{A}\,}^{\dagger }}\;{{e}^{i\Phi }}+\hat{A}\,\;{{e}^{-i\Phi }}}{\sqrt{2}}\ .\end{eqnarray} \tag{ A.13 }$$

For this reason, in the main text of the letter we have limited the discussion to the simpler case of a single-mode homodyne detection, with the proviso that whenever a quadrature is used, it actually refers to its expression (A.13) in the pulsed regime. In particular, the measured data can be used to reconstruct the expectation values of all signal observables that can be expressed as functions of the operators $\hat{A}\,$ and ${{\hat{A}\,}^{\dagger }}$. For instance, the mean value of the photo-current operator second moment with respect to the LO pulsed coherent state $\left| {\bar{z}} \right\rangle\left\langle {\bar{z}} \right|$ is

$$\begin{eqnarray}&&\left\langle {\bar{z}} \right|\;\hat{I}\,_{\Phi }^{2}\left| {\bar{z}} \right\rangle=2\;|\bar{z}{{|}^{2}}\;\hat{X}\,_{\Phi }^{2}\;+\;\underset{l}{\mathop{\mathop{\sum }^{}}}\,\hat{a}\,_{l}^{\dagger }{{\hat{a}\,}_{l}}\ ,\end{eqnarray} \tag{ A.14 }$$

and differs from the pulsed quadrature second moment $\hat{X}\,_{\Phi }^{2}$ by the number of photons in the pulsed signal, which is to be taken much smaller than the intensity $|\bar{z}{{|}^{2}}$ of the pulsed LO. A similar suppression by $|\bar{z}{{|}^{-2}}$ occurs for the correction terms appearing in higher moments so that the distribution of the outcomes of the homodyne photocurrent is equal to that of the corresponding field quadratures. As a final remark, notice that, in the case of only one frequency mode, the above setting reduces to the standard homodyne tomographic one.