When heterodyning beats homodyning: an assessment with quadrature moments

TL;DR

This paper compares homodyne and heterodyne measurement schemes for quantum state reconstruction, demonstrating the advantages of heterodyne measurements in estimating moments, especially for non-classical states, with optimized estimators reaching maximum-likelihood efficiency.

Contribution

It introduces optimized moment estimators for both schemes and provides a comprehensive analysis showing heterodyne measurement's superiority in reconstructing quantum state moments.

Findings

Heterodyne measurement outperforms homodyne in moment reconstruction.

Optimized estimators achieve maximum-likelihood efficiency in large data limits.

Heterodyne is particularly effective for non-classical states.

Abstract

We examine the moment-reconstruction performance of both the homodyne and heterodyne (double-homodyne) measurement schemes for arbitrary quantum states and introduce moment estimators that optimize the respective schemes for any given data. In the large-data limit, these estimators are as efficient as the maximum-likelihood estimators. We then illustrate the superiority of the heterodyne measurement for the reconstruction of the first and second moments by analyzing Gaussian states and many other significant non-classical states.

Tables

Table 1. Table 1: A list of characteristic functions for all the quantum states discussed. The symbols in this table are defined as α 0 ​ e − i ​ ϑ = ( x ϑ + i ​ p ϑ ) / 2 subscript 𝛼 0 superscript e i italic-ϑ subscript 𝑥 italic-ϑ i subscript 𝑝 italic-ϑ 2 \alpha_{0}\mathrm{e}^{\mbox{\footnotesize$-\mathrm{i}\vartheta$}}=(x_{\vartheta}+\mathrm{i}p_{\vartheta})/\sqrt{2} where x 0 = x ϑ = 0 subscript 𝑥 0 subscript 𝑥 italic-ϑ 0 x_{0}=x_{\vartheta=0} and p 0 = p ϑ = 0 subscript 𝑝 0 subscript 𝑝 italic-ϑ 0 p_{0}=p_{\vartheta=0} , 𝒈 = ( − g ​ g ∗ ) t 𝒈 superscript 𝑔 superscript 𝑔 t \boldsymbol{g}={(-g\,\,\,g^{*})}^{\mathrm{\,\textsc{t}}} , 𝑴 = 𝑯 † ​ 𝑮 het − 1 ​ 𝑯 𝑴 superscript 𝑯 † superscript subscript 𝑮 het 1 𝑯 \boldsymbol{M}=\boldsymbol{H}^{\dagger}\boldsymbol{G}_{\textsc{het}}^{-1}\,\boldsymbol{H} , and 𝑯 ​ = ^ ​ 1 2 ​ ( 1 1 − i i ) 𝑯 ^ 1 2 matrix 1 1 i i \boldsymbol{H}\,\widehat{=}\dfrac{1}{\sqrt{2}}\begin{pmatrix}1&1\\ -\mathrm{i}&\mathrm{i}\end{pmatrix} .

Class of Quantum States	Quadrature Characteristic Function	Husimi Characteristic Function
Gaussian	$\exp \left(-\frac{1}{2}{\left({𝒖}_{\vartheta }^{\text{t}}𝑮{𝒖}_{\vartheta }\right)}^2{k}^2+i{𝒖}_{\vartheta }^{\text{t}}{𝒓}_{0}k\right)$ $\frac{W_{}^W{}_{}{}^W{}_{}{}^W}{W_{}^W{}_{}{}^W{}_{}{}^W}$	$e^{g^{\ast }{\alpha }_0+g{\alpha }_0^{\ast }}\exp \left(\frac{det\left\{{𝑮}_{\text{het}}\right\}}{2}{𝒈}^{†}𝑴𝒈\right)$ $\frac{W_{}^W{}_{}{}^W{}_{}{}^W}{W_{}^W{}_{}{}^W{}_{}{}^W}$
Fock	$e^{-\frac{k^2}{4}}{L}_n\left(\frac{k^2}{2}\right)$	${}_1{}^{}F_1^{}(n+1;1;{\left|g\right|}^2)$
Even/odd coherent	$e^{-\frac{k^2}{4}}\frac{\cos \left(k{x}_{\vartheta }\right)\pm e^{-2{\left|{\alpha }_0\right|}^2}\cosh \left(k{p}_{\vartheta }\right)}{1\pm e^{-2{\left|{\alpha }_0\right|}^2}}$
Displaced Fock	$e^{-\frac{k^2}{4}+ik{x}_{\vartheta }}{L}_m\left(\frac{k^2}{2}\right)$	$e^{g^{\ast }{\alpha }_0+g{\alpha }_0^{\ast }}{}_1{}^{}F_1^{}(m+1;1;{\left|g\right|}^2)$
Photon-added coherent	$e^{\frac{k^2}{4}}\frac{{}_1{}^{}F_1^{}(m+1;1;\left({\alpha }_0+\frac{ik}{\sqrt{2}}{e}^{i\vartheta }\right)\left({\alpha }_0^{\ast }+\frac{ik}{\sqrt{2}}{e}^{-i\vartheta }\right))}{{}_1{}^{}F_1^{}(m+1;1;{\left|{\alpha }_0\right|}^2)}$	$\frac{{}_1{}^{}F_1^{}(m+1;1;{\left|g+{\alpha }_0\right|}^2)}{{}_1{}^{}F_1^{}(m+1;1;{\left|{\alpha }_0\right|}^2)}$

Equations

$$\mathbf{G}=\mathrm{Re}\left\{⟨\mathbf{R}{\mathbf{R}}^{\text{textsc}t}⟩\right\}-⟨\mathbf{R}⟩{⟨\mathbf{R}⟩}^{\text{textsc}\mathrm{t}}={\mathbf{G}}_2-{\mathbf{G}}_1,$$

$$\mathcal{D}=$$

$$+\left\{\text{cross terms}\right\}.$$

$${\mathbf{F}}_{\vartheta ,m}=\frac{N}{{\sigma }^2}\frac{\partial \mu }{\partial \mathbf{a}}\frac{\partial \mu }{\partial \mathbf{a}}+\frac{1}{2{\sigma }^4}\frac{\partial {\sigma }^2}{\partial \mathbf{a}}\frac{\partial {\sigma }^2}{\partial \mathbf{a}}$$

$${\tilde{\mathbf{F}}}_{m,\text{textsc}hom}={\int }_{(\pi )}\frac{\mathrm{d}\vartheta }{\pi }\frac{{\mathbf{F}}_{\vartheta ,m}}{N}={\int }_{(\pi )}\frac{\mathrm{d}\vartheta }{\pi }\frac{1}{\sigma (\vartheta {)}^2}\frac{\partial \mu (\vartheta )}{\partial \mathbf{a}}\frac{\partial \mu (\vartheta )}{\partial \mathbf{a}}$$

$${\tilde{\mathbf{F}}}_{1,\text{textsc}hom}={\int }_{(\pi )}\frac{\mathrm{d}\theta }{\pi }\frac{{\mathbf{m}}_{\vartheta }}{⟨X_{\vartheta }^2⟩-⟨X_{\vartheta }{⟩}^2},$$

$${\mathbf{m}}_{\vartheta }=\left(\begin{array}{cc}(\cos \vartheta {)}^2& \sin \vartheta \cos \vartheta \\ \sin \vartheta \cos \vartheta & (\sin \vartheta {)}^2\end{array}\right).$$

$${\mathcal{H}}_{\text{textsc}1,hom}=\mathrm{Tr}\left\{\mathbf{G}\right\}+2\sqrt{\det \left\{\mathbf{G}\right\}}.$$

$$\widehat{⟨X_k⟩}=\frac{1}{N}\sum _{j=1}^{n_x}{n}_{jk}{x}_{jk},\sum _{k=1}^{n_{\vartheta }}\sum _{j=1}^{n_x}{n}_{jk}=\sum _{k=1}^{n_{\vartheta }}{N}_k=N.$$

$${\hat{\mathbf{r}}}_{\text{textsc}hom}^{(\text{textsc}opt)}$$

$${\mathbf{W}}_1$$

$$\widehat{⟨X_k^m⟩}=\frac{1}{N}\sum _{j=1}^{n_x}{n}_{jk}{x}_{jk}^m(m=1,2,\dots ).$$

$${\tilde{\mathbf{F}}}_{2,\text{textsc}hom}=N{\int }_{(\pi )}\frac{\mathrm{d}\theta }{\pi }\frac{{\mathbf{M}}_{\vartheta }}{⟨X_{\vartheta }^4⟩-⟨X_{\vartheta }^2{⟩}^2},$$

$${\mathbf{M}}_{\vartheta }\widehat{=}\left(\begin{array}{c}{\left(\cos \vartheta \right)}^2\\ \sqrt{2}\sin \vartheta \cos \vartheta \\ {\left(\sin \vartheta \right)}^2\end{array}\right)\left(\begin{array}{ccc}{\left(\cos \vartheta \right)}^2& \sqrt{2}\sin \vartheta \cos \vartheta & {\left(\sin \vartheta \right)}^2\end{array}\right).$$

$$\mathbf{Y}\widehat{=}\left(\begin{array}{cc}{y}_1& y_2\\ y_2& y_3\end{array}\right)\mapsto \mathrm{vec}\left(\mathbf{Y}\right)\hat{\equiv }\left(\begin{array}{c}{y}_1\\ \sqrt{2}{y}_2\\ y_3\end{array}\right)$$

$${\hat{\mathbf{G}}}_{2,\text{textsc}hom}^{(\text{textsc}opt)}$$

$${\mathbf{W}}_2$$

$${\mathbf{G}}_{\text{textsc}het}\widehat{=}\left(\begin{array}{cc}\overline{x^2}-{\overline{x}}^2& \overline{xp}-\overline{x}\overline{p}\\ \overline{xp}-\overline{x}\overline{p}& \overline{p^2}-{\overline{p}}^2\end{array}\right)$$

$${\mathbf{G}}_{\text{textsc}het}=\mathbf{G}+\frac{1}{2}.$$

$${\mathrm{Var}}_{\text{textsc}q}\left[x\right]{\mathrm{Var}}_{\text{textsc}q}\left[p\right]=\left(⟨(\Delta X{)}^2⟩+\frac{1}{2}\right)\left(⟨(\Delta P{)}^2⟩+\frac{1}{2}\right)\ge 1,$$

$${\hat{\mathbf{r}}}_{\text{textsc}het}\widehat{=}\frac{1}{N}\sum _{j=1}^N\left(\begin{array}{c}{x}_j\\ p_j\end{array}\right)$$

$${\mathcal{D}}_{1,\text{textsc}het}=\mathbb{E}\left[{\left({\hat{\mathbf{r}}}_{\text{textsc}het}-\mathbf{r}\right)}^2\right]=\frac{1}{N}\left({\mathrm{Var}}_{\text{textsc}q}\left[x\right]+{\mathrm{Var}}_{\text{textsc}q}\left[p\right]\right),$$

$${\mathcal{H}}_{1,\text{textsc}het}={\mathrm{Var}}_{\text{textsc}q}\left[x\right]+{\mathrm{Var}}_{\text{textsc}q}\left[p\right]=\mathrm{Tr}\left\{\mathbf{G}\right\}+1.$$

$${\hat{\mathbf{G}}}_{2,\text{textsc}het}\widehat{=}\frac{1}{N}\sum _{j=1}^N\left(\begin{array}{cc}{x}_j^2& x_j{p}_j\\ x_j{p}_j& p_j^2\end{array}\right),$$

$${\mathbf{G}}_{2,\text{textsc}het}={\mathbf{G}}_2+\frac{1}{2}.$$

$${\mathcal{D}}_{2,\text{textsc}het}$$

$$=\mathrm{Tr}\left\{\mathbb{E}\left[{\hat{\mathbf{G}}}_{2,\text{textsc}het}^2\right]\right\}-\mathrm{Tr}\left\{{\mathbf{G}}_{2,\text{textsc}het}^2\right\}$$

$$=\frac{1}{N}\left({\mathrm{Var}}_{\text{textsc}q}\left[x^2\right]+{\mathrm{Var}}_{\text{textsc}q}\left[p^2\right]+2{\mathrm{Var}}_{\text{textsc}q}\left[xp\right]\right).$$

$${\mathcal{H}}_{2,\text{textsc}het}={\mathrm{Var}}_{\text{textsc}q}\left[x^2\right]+{\mathrm{Var}}_{\text{textsc}q}\left[p^2\right]+2{\mathrm{Var}}_{\text{textsc}q}\left[xp\right].$$

$$\mathbf{x}=\mathrm{vec}\left({\hat{\mathbf{G}}}_{2,\text{textsc}het}\right)\hat{\equiv }\frac{1}{N}\sum _{j=1}^N\left(\begin{array}{c}{x}_j\\ \sqrt{2}{x}_j{p}_j\\ p_j\end{array}\right)$$

$$\mathbf{\Sigma }=\overline{\mathbf{x}{\mathbf{x}}^{\text{textsc}\mathrm{t}}}-\mathbf{\mu }{\mathbf{\mu }}^{\text{textsc}\mathrm{t}}\widehat{=}\frac{1}{N}\left(\begin{array}{ccc}{\mathrm{Var}}_{\text{textsc}q}\left[x^2\right]& \ast & \ast \\ \ast & 2{\mathrm{Var}}_{\text{textsc}q}\left[xp\right]& \ast \\ \ast & \ast & {\mathrm{Var}}_{\text{textsc}q}\left[p^2\right]\end{array}\right),$$

$${\mathcal{H}}_{1,\text{textsc}hom}=\mathrm{Tr}\left\{\mathbf{G}\right\}+2\sqrt{\det \left\{\mathbf{G}\right\}}\ge \mathrm{Tr}\left\{\mathbf{G}\right\}+1={\mathcal{H}}_{1,\text{textsc}het}.$$

$${\gamma }_1=\frac{{\mathcal{H}}_{1,\text{textsc}het}}{{\mathcal{H}}_{1,\text{textsc}hom}},$$

Full text

I Introduction

The next-generation quantum technologies introduce novel and innovative routes to the understanding and implementation of measurements, communication, and computation. In this respect, the manipulation of a quantum light source using continuous-variable (CV) measurements offer many advantages Braunstein and van Loock (2005); Ferraro et al. (2005); Cerf et al. (2007); Andersen et al. (2010); Adesso et al. (2014). There exist two standard $CV$ measurement schemes. The more commonly employed homodyne detection Yuen and Chan (1983); Abbas et al. (1983); Schumaker (1984), which performs an approximate measurement of rotated photonic quadratures Banaszek and Wódkiewicz (1997), probes the marginal distribution of the Wigner function of the unknown quantum state Vogel and Risken (1989). The other less widely adopted double-homodyne detection, or the heterodyne detection, involves the joint measurement of complementary observables Arthurs and Kelly (1965); Yuen (1982); Arthurs and Goodman (1988); Martens and de Muynck (1990, 1991); Raymer (1994); Trifonov et al. (2001); Werner (2004)that directly samples the phase space according to the Husimi function Stenholm (1992) and is connected to the conventional heterodyne scheme Javan et al. (1962); Read and Turner (1965); Carleton and Maloney (1968); Gerhardt et al. (1972); Yuen and Shapiro (1980); Shapiro and Wagner (1984); Shapiro (1985); Walker and Carroll (1986); Collett et al. (1987); Lai and Haus (1989).

These measurement schemes, which experimentally probe quasi-probability distributions, can also be equivalently understood as practical means to directly characterize the source in terms of the ordered moments of the quadrature operators in phase space. Gaussian states Ferraro et al. (2005) for example, which are important in analyzing CV quantum information processing Lorenz et al. (2004); Lance et al. (2005); Scarani et al. (2009); Weedbrook et al. (2012), are conveniently described by this representation since all their operator moments are functions of only the first and second moments. Therefore, estimating the first and second moments are enough to fully reconstruct the Gaussian state or verify if the reconstructed state is accurately Gaussian Řeháček et al. (2009). Higher moments come into play for general quantum states. On its own right, the topic of operator moments of quantum states draws interest in the context of generalized uncertainty relations Angulo (1993, 1994), non-classicality detection Simon et al. ; Arvind et al. (1998), entanglement detection R. (2012); Ivan et al. (2012), and cryptography Leverrier (2012); Thearle et al. (2016).

In Refs. Řeháček et al. (2015) and Müller et al. (2016), we theoretically and experimentally compared the two measurement schemes, using a polarization-squeezing setup Heersink et al. (2005); Grangier et al. (1987); Josse et al. (2004); Marquardt et al. (2007); Müller et al. (2012); Peuntinger et al. (2014) for the latter. We analyzed the physical implications of having the unavoidable Arthurs-Kelly type noise that is inherent in the joint measurement heterodyne scheme on moment reconstruction. We found that despite this additional noise, for a single-mode central-Gaussian source the heterodyne scheme still results in second-moment estimators that are more accurate than the homodyne scheme for a wide range of the squeezing strength and temperature parameter.

In this article, we extend the theory of these two CV measurement schemes to general quantum states and show that the tomographic advantage in using the heterodyne scheme carries over to other interesting and important non-Gaussian states. This message is conveyed in five main sections. Section II gives an overview of the fundamental elements in first- and second-moment tomography, as well as the concept of reconstruction accuracy. These elements are then used to discuss the general theory of moment reconstruction for the homodyne and heterodyne schemes in Sec. III. In that section, we shall also introduce optimal moment estimators that asymptotically approach the respective Cramér-Rao bounds, which are derived in Appendix A. In Sec. IV, we shall study the CV schemes in first-moment estimation where it shall be shown that heterodyne detection will always outperform homodyne detection unless the state is of minimum uncertainty, in which case the two schemes give equal reconstruction accuracy per sampling event. This result shall be discussed with some interesting classes of non-Gaussian states. Next, we study the results for second-moment estimation Sec. V with the same classes of non-Gaussian states and illustrate once again the tomographic advantages of using the heterodyne scheme in moment tomography. Finally, Sec. VI concludes the presented results in a summary.

II The covariance matrix and moment-reconstruction accuracy

In dealing with single-mode bosonic systems such as photons, for the pair of position $X$ and momentum $P$ quadrature operators obeying $[X,P]=\mathrm{i}$ (with the quantum unit $\hbar\equiv 1$) that form the column $\boldsymbol{R}={(X\,\,\,P)}^{\mathrm{\,\textsc{t}}}$, the covariance matrix can be written as

[TABLE]

where we have introduced the first- ($\boldsymbol{G}_{1}=\langle\boldsymbol{R}\rangle{\langle\boldsymbol{R}\rangle}^{\mathrm{\,\textsc{t}}}$) and second-moment ($\boldsymbol{G}_{2}=\mathrm{Re}\!\left\{\langle\boldsymbol{R}\boldsymbol{R}^{\textsc{t}}\rangle\right\}$) matrices. The two independent parameters $\left\{\langle X\rangle,\langle P\rangle\right\}$ in $\boldsymbol{G}_{1}$ and three independent parameters $\left\{\langle X^{2}\rangle,\frac{1}{2}\langle\{X,P\}\rangle,\langle P^{2}\rangle\right\}$ in $\boldsymbol{G}_{2}$ constitute the complete set of five parameters that characterize $\boldsymbol{G}$. The well-known class of Gaussian states possesses a Gaussian Wigner function or any other kind of well-behaved quasi-probability distribution. As a consequence, any Gaussian state is fully described by only $\boldsymbol{G}_{1}$ and $\boldsymbol{G}_{2}$.

The covariance matrix for any quantum state obeys the inequality $\boldsymbol{G}\geq\boldsymbol{\sigma}_{y}/2$ in terms of the Pauli matrix $\boldsymbol{\sigma}_{y}$, which is a recast of the Heisenberg-Robertson-Schrödinger (HRS) uncertainty relation for position and momentum operators. This gives the equivalent stricter inequality $\det\!\left\{\boldsymbol{G}\right\}\geq 1/4$ in addition to the standard positivity constraint for $\boldsymbol{G}$. The reconstruction of the full covariance matrix $\boldsymbol{G}$ involves the quantum tomography of all the five independent parameters that define the first and second operator moments of the state. Here, the figure of merit the reconstruction accuracy is the mean squared-error (MSE) $\mathcal{D}=\mathbb{E}\!\left[\mathrm{Tr}\!\left\{\left(\widehat{\boldsymbol{G}}-\boldsymbol{G}\right)^{2}\right\}\right]$ between $\boldsymbol{G}$ and its estimator $\widehat{\boldsymbol{G}}$. In terms of $\boldsymbol{G}_{1}$ and $\boldsymbol{G}_{2}$,

[TABLE]

To illustrate the physics behind moment reconstruction, we shall analyze both the first and second-moment reconstruction accuracy separately. In practice, these analyses are relevant to the situation where the reconstructions of $\boldsymbol{G}_{1}$ and $\boldsymbol{G}_{2}$ are carried out with independent data. For this situation, the $\left\{\text{cross terms}\right\}$ in Eq. (2) vanish so that the total MSE is the sum of the respective MSEs $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ of the reconstructed moments. From hereon, to facilitate discussions, we shall analyze the quantity $\boldsymbol{r}=\left<\boldsymbol{R}\right>$ in place of $\boldsymbol{G}_{1}$, where $\mathcal{D}_{1}=\mathbb{E}\!\left[(\widehat{\boldsymbol{r}}-\boldsymbol{r})^{2}\right]$.

In unbiased statistical estimation theory Cox (2006), the MSE $\mathcal{D}\geq\mathrm{Tr}\!\left\{\boldsymbol{F}^{-1}\right\}$ is bounded from below by the inverse of the Fisher information matrix $\boldsymbol{F}$, or the Cramér-Rao bound (CRB). Consequently, we have the respective first- and second-moment CRBs $\mathcal{D}_{1}\geq\mathrm{Tr}\!\left\{\boldsymbol{F}_{1}^{-1}\right\}$ and $\mathcal{D}_{2}\geq\mathrm{Tr}\!\left\{\boldsymbol{F}_{2}^{-1}\right\}$. Therefore, the general theory of the Fisher matrices $\boldsymbol{F}_{1}$ and $\boldsymbol{F}_{2}$ for the two CV schemes is in order.

III General theory

III.1 Homodyne detection

The homodyne detection Yuen and Chan (1983); Abbas et al. (1983); Schumaker (1984) involves a 50:50 beam splitter that introduces an interference between the optical source of an unknown state (signal) and the local oscillator (coherent-state reference source or simply LO), the latter of which is set to a much larger optical intensity than the mean intensity of the optical source of an unknown quantum state $\rho$ [see Fig. 1(a)]. A subtraction of the output photocurrents gives a distribution of voltage readouts $-\infty<x_{\vartheta}<\infty$ for the LO phase $0\leq\vartheta\leq\pi$, which essentially corresponds to the eigenvalue probability distribution of the quadrature operator $X_{\vartheta}=X\cos\vartheta+P\sin\vartheta$. It then follows that statistically, the expectation value $\left<X_{\vartheta}^{m}\right>$ for any integer value $m$ contains all measurable information about the $m$th operator moments of $X$ and $P$. Since the data acquired with this scheme are the marginals of the Wigner function, the first $(m=1)$ and second $(m=2)$ moments, or $\boldsymbol{G}$, that are reconstructed with these data may be attributed to this quasi-probability distribution function.

In a typical homodyne experiment, the value of $\vartheta$ is set to increase linearly. The data collected would then be binned for all the measured $\vartheta$ values. The data bins are mutually independent, so that the Fisher matrices $\boldsymbol{F}_{1,\textsc{hom}}$ and $\boldsymbol{F}_{2,\textsc{hom}}$ for the respective first- and second-moment CRBs can each be understood as a summation of Fisher matrices of every LO phase bin according to the additivity property of the Fisher information. In the limit of large number of sampling events $N$, the central limit theorem states that the unbiased estimator $\widehat{\left<X^{m}_{\vartheta}\right>}$ of the $m$th quadrature moment $\left<X^{m}_{\vartheta}\right>$ that is defined as an average sum of independently collected random voltage values for the phase $\vartheta$ follows a Gaussian distribution of data mean $\mu=\mu(\vartheta)=\left<X^{m}_{\vartheta}\right>$ and data variance $\sigma^{2}/N$ where $\sigma^{2}=\sigma(\vartheta)^{2}=\left<X^{2m}_{\vartheta}\right>-\left<X^{m}_{\vartheta}\right>^{2}$, so that the Fisher matrix

[TABLE]

for a given LO phase $\vartheta$ in the large-$N$ limit follows the well-known expression for Gaussian distributions, where in our case $\boldsymbol{a}$ is the column of $m$th moment parameters we are interested in reconstructing. As it is clear that only the first term of (3) would survive in this limit, we thus have the scaled homodyne Fisher matrix

[TABLE]

with respect to the number of sampling events $N$ for the complete set of homodyne quadrature-eigenstate outcomes. Scaled statistical quantities such as this one shall be the focus of this article in analyzing tomographic performances, as the scaled CRB (sCRB) represents the power-law coefficient of the MSE in this limit that determines the difficulty in obtaining an estimator of a certain pre-chosen MSE accuracy.

III.1.1 First-moment reconstruction

All information about the first moments, $\boldsymbol{a}=\boldsymbol{r}$, of the covariance matrix is completely encoded in the expectation value $\mu(\vartheta)=\langle X_{\vartheta}\rangle$. The variance for the data is then given by $\sigma(\vartheta)^{2}=\langle X_{\vartheta}^{2}\rangle-\langle X_{\vartheta}\rangle^{2}$. The scaled Fisher matrix for the first-moment estimation with homodyne data is therefore given by

[TABLE]

where

[TABLE]

The integral can be evaluated exactly, bringing us to the closed-form expression

[TABLE]

With the machinery of quantum tomography (see Appendix A.1), an observer can construct the optimal moment estimator that achieves the sCRB. Suppose that the observer collects homodyne data for $N$ sampling events and bins the voltage values into $\{x_{jk}\}$ according to a discrete number $n_{\vartheta}$ of LO phase bins $\vartheta_{k}$, where $j$ labels the $n_{x}$ real voltage values per LO phase bin $\vartheta=\vartheta_{k}$ and $k$ labels the phase bins. Then an unbiased estimator for any particular expectation value $\langle X_{k}\rangle\equiv\langle X_{\vartheta_{k}}\rangle$ would be

[TABLE]

Then upon denoting $\boldsymbol{u}_{k}\equiv\boldsymbol{u}_{\vartheta_{k}}=(\cos{\vartheta_{k}}\,\,\sin{\vartheta_{k}})^{\textsc{t}}$ where we note that $\boldsymbol{m}_{k}=\boldsymbol{u}_{k}{\boldsymbol{u}}^{\mathrm{\,\textsc{t}}}_{k}$, the optimal first-moment estimator is given by

[TABLE]

which is immediately computable given the processed data $\left\{\widehat{\left<X_{k}\right>}\right\}$ and $\left\{\widehat{\left<X^{2}_{k}\right>}\right\}$ that are defined by

[TABLE]

That this estimator achieves the sCRB asymptotically is also shown in Appendix A.1. This equivalently implies that the optimal estimator is as efficient as the maximum-likelihood (ML) estimator for the multinomially-distributed binned data $\{x_{jk}\}$.

III.1.2 Second-moment reconstruction

To estimate $\boldsymbol{G}_{2}$, it is clear that second-moment information is completely encoded in the second quadrature moment $\langle X^{2}_{\vartheta}\rangle$, which is a function of the three independent parameters $a_{1}=\langle X^{2}\rangle$, $a_{2}=\langle\frac{1}{2}\{\Delta X,\Delta P\}\rangle$ and $a_{3}=\langle P^{2}\rangle$. From Eq. (4), the corresponding $3\times 3$ Fisher matrix for these three parameters is

[TABLE]

where

[TABLE]

The analytical answer to $\widetilde{\boldsymbol{F}}_{2,\textsc{hom}}$ for an arbitrary state, and its subsequent inverse $\mathcal{H}_{2,\textsc{hom}}=\mathrm{Tr}\!\left\{\widetilde{\boldsymbol{F}}_{2,\textsc{hom}}^{-1}\right\}$ is difficult to calculate, as the denominator in the integrand generally contains trigonometric functions in a complicated manner. Nevertheless, the integral can be calculated explicitly for many interesting and important quantum sources.

The optimal second-moment estimator (see Appendix A.2) that achieves the corresponding sCRB can be cleanly expressed using the vectorization operation $\mathrm{vec}\!\left(\boldsymbol{Y}\right)$ that turns a matrix into a column according to

[TABLE]

in any pre-chosen basis, such that $\mathrm{Tr}\!\left\{\boldsymbol{Y}_{1}\boldsymbol{Y}_{2}\right\}={\mathrm{vec}\!\left(\boldsymbol{Y}_{1}\right)}^{\mathrm{\,\textsc{t}}}\mathrm{vec}\!\left(\boldsymbol{Y}_{2}\right)$ for any two $2\times 2$ symmetric matrices $\boldsymbol{Y}_{1}$ and $\boldsymbol{Y}_{2}$. Given the processed data defined in Eq. (10), the final operationally-ready expressions for this optimal estimator are given as follows:

[TABLE]

For accurate tomography, the value of $N$ is typically large enough such that $\widehat{\boldsymbol{G}}^{(\textsc{opt})}_{2,\textsc{hom}}$ is a proper covariance matrix and approaches the ML estimator that asymptotically achieves the sCRB, which is strictly speaking the correct regime where $\widehat{\boldsymbol{G}}^{(\textsc{opt})}_{2,\textsc{hom}}$ is to be used for second-moment tomography. On a separate note, optimal estimators for overcomplete quantum-state tomography of $\rho$ was developed in Zhu (2014) and later rederived in Teo (2015) with the variational principle that is also used to construct the optimal moment estimators in Appendix A.

III.2 Heterodyne detection

The heterodyne detection scheme essentially uses two homodyne setups to perform a joint measurement of two complementary observables [see Fig. 1(b)], which are in this case chosen to be the standard $X$ and $P$ quadrature pair for convenience. It is well-known (Arthurs and Kelly (1965); Yuen (1982); Arthurs and Goodman (1988); Martens and de Muynck (1990, 1991); Raymer (1994); Trifonov et al. (2001); Werner (2004)) that the product of their joint-measurement standard deviations has a larger lower bound than the usual one-half of a quantum unit given by the original Heisenberg relation owing to the additional quantum noise introduced by the joint measurement.

The outcomes for this scheme are in fact the overcomplete set of coherent states. This means that the resulting data are direct phase-space samples of the Husimi function for the statistical operator $\rho$. The technical complication of having additional measurement noise can therefore be translated completely into the phase-space language that is relevant in our subsequent analysis. Given an infinite set of the Husimi-function data, we have access to the moments $\overline{x^{k}p^{l}}$ (the overline denotes the average with respect to the Husimi function, or simply the Husimi average), with which the corresponding “$\boldsymbol{G}$” operator

[TABLE]

can be directly constructed. One can then show that for any quantum state,

[TABLE]

The corresponding Arthurs-Kelly type measurement uncertainty relation

[TABLE]

which is saturated by coherent states $[\langle(\Delta X)^{2}\rangle=\langle(\Delta P)^{2}\rangle=1/2]$, can thereafter be understood as a physical manifestation of the Gauss-Weierstrass transform [related to Eq. (16)] between the Wigner and Husimi functions if the joint-measurement data are directly used to calculate variances (here denoted by $\mathrm{Var}_{\textsc{q}}\!\left[y\right]=\overline{y^{2}}-\overline{y}^{2}$ for a complete Husimi-function data $\{y\}$). We shall show that this additional quantum noise, when combined with optimal tomography strategies, can still lead to better moment-reconstruction accuracies relative to the homodyne scheme.

III.2.1 First-moment reconstruction

From Sec. III.2, we note that the data collected from the heterodyne scheme are a scatter set of phase-space coordinates $\{(x_{j},p_{j})\}$ that are distributed according to the Husimi function. As Eq. (16) tells us that there is no difference between the state average $\boldsymbol{r}$ and Husimi average of $(x\,\,\,p)^{\textsc{t}}$, being a two-parameter estimation scheme, the first-moment sCRB with respect to the heterodyne data can again be found by taking the average of the distance between the estimator

[TABLE]

and the true column $\boldsymbol{r}^{\textsc{t}}\widehat{=}\,(\overline{x}\,\,\,\overline{p})$:

[TABLE]

so that

[TABLE]

That $N\mathcal{D}_{1,\textsc{het}}=\mathcal{H}_{1,\textsc{het}}$ follows in the limit of large $N$, where the unbiased estimator $\widehat{\boldsymbol{r}}_{\textsc{het}}$ is asymptotically optimal since in this limit, the distribution of $\widehat{\boldsymbol{r}}_{\textsc{het}}$ follows a bivariate Gaussian distribution with vanishing widths, such that $\widehat{\boldsymbol{r}}_{\textsc{het}}$ becomes the ML estimator that approaches the sCRB for this Gaussian distribution.

III.2.2 Second-moment reconstruction

Similarly, to arrive at the optimal accuracy for estimating $\boldsymbol{G}_{2}$ using heterodyne data, we define the natural second-moment estimator

[TABLE]

where $\{(x_{j},p_{j})\}$ are again the sampled Husimi-function data collected during heterodyne detection. From Eq. (16), we get

[TABLE]

The MSE $\mathcal{D}_{2,\textsc{het}}$ for heterodyne detection concerning second-moment estimation is consequently given by

[TABLE]

In the large-$N$ limit, this MSE is essentially the sCRB

[TABLE]

since $\widehat{\boldsymbol{G}}_{2,\textsc{het}}$ again becomes the ML estimator. To see this, we inspect the Fisher matrix $\boldsymbol{F}_{2,\textsc{het}}$ for the estimator $\widehat{\boldsymbol{G}}_{2,\textsc{het}}$. If we look at the random column

[TABLE]

that represents $\widehat{\boldsymbol{G}}_{2,\textsc{het}}$, we find that in the limit of large $N$, the central limit theorem again says that $\boldsymbol{x}$ follows a Gaussian distribution defined by the mean $\boldsymbol{\mu}=\overline{\boldsymbol{x}}\,\widehat{=}\,{(\overline{x^{2}}\,\,\,\sqrt{2}\,\overline{xp}\,\,\,\,\overline{p^{2}})}^{\mathrm{\,\textsc{t}}}$ and the covariance matrix

[TABLE]

so that we eventually recover the well-known result $\boldsymbol{\Sigma}=\boldsymbol{F}_{2,\textsc{het}}^{-1}$ for Gaussian scatter data that saturates the CRB as we remember that $\mathrm{Tr}\!\left\{\boldsymbol{\Sigma}\right\}=\mathcal{D}_{2,\textsc{het}}$. Equation (24) then follows tout de suite.

IV First-moment estimation

IV.1 General optimality of heterodyne tomography

As far as first-moment estimation is concerned, the general results in Eqs. (7) and (20) imply that $\mathcal{H}_{1,\textsc{het}}\leq\mathcal{H}_{1,\textsc{hom}}$ for any quantum state. This main result hinges on the physical HRS uncertainty relation, which is equivalent to the constraint $\det\!\left\{\boldsymbol{G}\right\}\geq 1/4$ for the covariance matrix $\boldsymbol{G}$. This constraint means that

[TABLE]

This implies that for all quantum states, the reconstruction accuracy of the optimal heterodyne first-moment estimator is always higher or equal to that of the optimal homodyne first-moment estimator in locating the average center of the quantum state in phase space. For minimum-uncertainty states, the accuracies of the two schemes are equal ($\mathcal{H}_{1,\textsc{hom}}=\mathcal{H}_{1,\textsc{het}}$). Subsequent well-known and interesting examples merely illustrate this fundamental fact. In terms of the first-moment performance ratio

[TABLE]

a subunit magnitude indicates that the heterodyne scheme outperforms the homodyne scheme.

IV.2 Gaussian states

For a Gaussian state where the covariance matrix $\boldsymbol{G}$ characterizes the spread of its Wigner function, the state variance of $X_{\vartheta}$ is simply

[TABLE]

From Eqs. (7) and (20), the first-moment performance ratio

[TABLE]

clearly cannot exceed one since any physical state satisfying the HRS uncertainty relation must take $\det\!\left\{\boldsymbol{G}\right\}\geq 1/4$. The maximum value of $\gamma_{1}=1$ is attained for minimum-uncertainty states.

IV.3 Fock states

A Fock state of the ket $\left|{n}\right>$ is always centered at the origin of the phase space $(\boldsymbol{r}=\boldsymbol{0})$. The circular symmetry of these states imply the fact that $(\Delta X)^{2}=(\Delta P)^{2}=n+1/2=(\Delta X_{\vartheta})^{2}$, whence

[TABLE]

since such states have zero first moments. On the other hand, for the heterodyne scheme, we get

[TABLE]

by simply using the Husimi characteristic function from Table 1 in Appendix B. Therefore, we get a

[TABLE]

that is always sub-unity unless $n=0$, a result that is again familiar from Sec. IV.2. In the limit of large photon numbers, the first-moment $\gamma_{1}$ approaches 1/2.

IV.4 Even/odd coherent states

Another popular class of non-Gaussian states with interesting phase-space quantum interference features are the even/odd coherent states characterized by the ket $\left|{\pm;\alpha_{0}}\right>=(\left|{\alpha_{0}}\right>\pm\left|{-\alpha_{0}}\right>)\mathcal{N}_{\pm}$ of appropriate normalization constants \mathcal{N}_{\pm}=1/\sqrt{2\pm 2\,\mathrm{e}^{\mbox{\footnotesize-2|\alpha_{0}|^{2}}}}, whose first moments $\boldsymbol{r}$ are all equal to zero. Using the definitions $a=\frac{1}{2}\left[\left<(\Delta X)^{2}\right>-\left<(\Delta P)^{2}\right>\right]=\alpha_{0}^{2}$ and $b_{\pm}=\frac{1}{2}\left[\left<(\Delta X)^{2}\right>+\left<(\Delta P)^{2}\right>\right]=\alpha_{0}^{2}\left[\tanh(\alpha_{0}^{2})\right]^{\pm 1}+1/2$,

[TABLE]

For the heterodyne counterpart, one finds that

[TABLE]

which contributes to the performance ratio

[TABLE]

For both types of coherent state superpositions, $\gamma_{1}\rightarrow 1$ as $\alpha_{0}\rightarrow\infty$. For even coherent states, the performance ratio $\gamma_{1}=1$ when $\alpha_{0}=0$ as it should. Otherwise, this ratio is always less than one for any positive $\alpha_{0}$. There exists a single local minimum of $\gamma_{1}\approx 0.7577$ at $\alpha_{0}\approx 1.715$. For odd coherent states, $\gamma_{1}<1$ for all $\alpha_{0}$ values, with the minimum value of $\gamma_{1}=1/3$ at $\alpha_{0}=0$. For these states, $\gamma_{1}$ increases monotonically to one as $\alpha_{0}$ tends to infinity.

IV.5 Displaced Fock states

Displacement and photon-addition are two important physical procedures that are frequently discussed in quantum physics. The different orders in which these processes are carried out on the vacuum state give output states of a different nature. Displacing an $m$-photon-added vacuum state by a complex amplitude $\alpha_{0}$ results in displaced Fock states defined by the ket $D(\alpha_{0})\left|{m}\right>$ can be effectively performed using a beam splitter with a high transmissivity and a strong coherent state Paris (1996); Banaszek and Wódkiewicz (1999).

It can be shown easily that the first-moment sCRBs are indeed given by Eqs. (31) and (32), so that the performance ratio is then completely identical to that of the usual central Fock states in Eq. (33). This reflects the physical fact that the accuracy in estimating the displacement cannot explicitly depend on where the center of the displaced Fock states is when full sets of CV measurement outcomes are considered, as the tomographic coverage of the entire phase space is then complete. This accuracy depends only on the variances, which describe the second-order symmetry and is unaffected at all by the displacement.

IV.6 Photon-added coherent states

A swap in the order of photon addition and displacement on the vacuum state gives the photon-added coherent state of $m$ added photons and reference amplitude $\alpha_{0}$ is defined by the ket $\left|{m;\alpha_{0}}\right>=\mathcal{N}_{m,|\alpha_{0}|^{2}}{A^{\dagger}}^{m}\left|{\alpha_{0}}\right>$ with the bosonic annihilation operator $A$, where the normalization constant \mathcal{N}_{m,|\alpha_{0}|^{2}}=\mathrm{e}^{\mbox{\footnotesize|\alpha_{0}|^{2}/2}}/\sqrt{m!{}_{1}\mathrm{F}_{1}\!\left({m+1};{1};|\alpha_{0}|^{2}\right)} involves the confluent hypergeometric function of the first kind ${}_{1}\mathrm{F}_{1}\!\left({a};{b};y\right)$. The integer value $m$ denotes the extent to which the mean photon number

[TABLE]

which is always larger than $|\alpha_{0}|^{2}+m$ whenever $\alpha_{0}\neq 0$, is increased nonlinearly by the operation by ${A^{\dagger}}^{m}$ on the reference coherent ket $\left|{\alpha_{0}}\right>$. This particular class of quantum states is but one of many possible kinds of photon-added states, which are of interest to the quantum community for testing some fundamental statements Parigi et al. (2007); Kim et al. (2008); Zavatta et al. (2009).

For these photon-added coherent states, the second-order symmetry is now affected by the combined action of the displacement and photon addition, so that $\left<(\Delta X)^{2}\right>$ and $\left<(\Delta P)^{2}\right>$ are functions of $m$ and $\alpha_{0}$. These expressions can be straightforwardly computed with the help of the characteristic functions given in Table 1 in Appendix B. By defining

[TABLE]

such that $b>a$, the first-moment sCRB for homodyne detection is of the same form as in Eq. (34), namely

[TABLE]

The first-moment sCRB for heterodyne detection is given by

[TABLE]

Clearly, when $\alpha_{0}=0$, the answers in Eqs. (31) and (32) for an $m$-number Fock state are reproduced exactly. With $m=0$, the respective sCRBs of a value of 2 for all $\alpha_{0}$s are furthermore consistent with Sec. IV.2. Otherwise, $\gamma_{1}$ is always sub-unity, and approaches unity as $\alpha_{0}\rightarrow\infty$.

V Second-moment estimation

V.1 Gaussian states

It seems fitting to commence the discussion of second-moment estimation with the Gaussian state, for it is natural to begin with the generalization of the results that already appeared in Refs. Řeháček et al. (2015) and Müller et al. (2016) to general noncentral Gaussian states $\left(\boldsymbol{r}\neq\boldsymbol{0}\right)$. We suppose that the Gaussian state of the covariance matrix $\boldsymbol{G}$ is centered at $\boldsymbol{r}=\boldsymbol{r}_{0}={(x_{0}\,\,\,p_{0})}^{\mathrm{\,\textsc{t}}}$. From Table 1 in Appendix B, by defining $\mu_{\vartheta}={\boldsymbol{u}}^{\mathrm{\,\textsc{t}}}_{\vartheta}\,\boldsymbol{r}_{0}$ and $\sigma_{\vartheta}^{2}={\boldsymbol{u}}^{\mathrm{\,\textsc{t}}}_{\vartheta}\,\boldsymbol{G}\,\boldsymbol{u}_{\vartheta}$, the variance for the second quadrature moment reads

[TABLE]

For central Gaussian states $(\langle X\rangle=\langle P\rangle=0)$, we have $\langle X^{4}_{\vartheta}\rangle=3\langle X^{2}_{\vartheta}\rangle^{2}$ and the scaled Fisher matrix in Eq. (11) turns into the familiar form in Řeháček et al. (2015); Müller et al. (2016). For the more general situation, one can repeat the contour-method integration in Řeháček et al. (2015) to calculate the scaled Fisher matrix in Eq. (11). The answer is given as

[TABLE]

together with the definitions

[TABLE]

When $\boldsymbol{r}_{0}=0$, we have $w_{1}=a$, $w_{2}=b$ and $w_{3}=c$ and the scaled Fisher matrix $\widetilde{\boldsymbol{F}}_{2,\textsc{hom}}$ reduces to that for the central Gaussian state in Řeháček et al. (2015). For the general setting, the full expression of $\mathcal{H}_{2,\textsc{hom}}$ is omitted here in this case due to its complexity. On the other hand, the sCRB with the heterodyne scheme for these noncentral Gaussian states can be calculated directly from Eq. (24) using the characteristic function in Table 1 and is given by

[TABLE]

where one immediately verifies the counterpart expression in Řeháček et al. (2015) for the central Gaussian states upon setting $\boldsymbol{r}_{0}=\boldsymbol{0}$.

At this stage, we reassure ourselves the physics of the problem of second-moment tomography by understanding, first, that in the case where tomography is performed on the full covariance matrix $\boldsymbol{G}$ then the sCRB, which is the minimum of the MSE, should not depend on the orientation of the two-dimensional uncertainty region (here being an ellipse for any Gaussian state) described by the eigenvectors of this matrix but only its eigenvalues owing to the form of the MSE. Additionally, the accuracy should also be independent of $\boldsymbol{r}_{0}$. When only the second-moment matrix $\boldsymbol{G}_{2}$ is reconstructed, the sCRB should also not depend on its eigenvectors but only its eigenvalues. The physics remains the same. However, there is a difference between estimating the full matrix $\boldsymbol{G}$ and estimating just $\boldsymbol{G}_{2}$. Since $\boldsymbol{G}_{2}$ is in general an increasing function of the first moments, this means that as the displacement of the center from the phase-space origin for the quantum state increases, the geometric mean of eigenvalues (GME) of $\boldsymbol{G}_{2}$ correspondingly becomes larger so that the second-order-“temperature” of the state, a terminology borrowed from Gaussian states, as described by the GME is now higher and this results in a stronger $\boldsymbol{G}_{2}$-“thermal” property much like the thermal Gaussian states. So we would expect, based on the findings in Řeháček et al. (2015), that states with large displacements give poor second-moment tomographic accuracies for both CV schemes, and yet provides a subunit

[TABLE]

performance ratio. It is also physically intuitive that the accuracies for both schemes should also be independent of the angle of displacement, but depend only on the magnitude of the displacement. For non-Gaussian states, the fourth moments arising from the structure of the MSE, which are no longer functions of the first and second moments as is the case for Gaussian states, also contribute to the sCRB, and therefore $\gamma_{2}$, as described in the general theory in Sec. III.

This physics, however, seems to be violated by the noncentral-Gaussian-state expressions in (42) and (44), namely that $\mathcal{H}_{2,\textsc{het}}$ depends on the explicit displacement vector $\boldsymbol{r}_{0}$ and covariance matrix $\boldsymbol{G}$, for instance. This mishap has nothing to do with any kind of physical violation, but has only to do with the way we specify Gaussian states. By choosing to parametrize a multivariate Gaussian distribution using the natural independent parameters $\boldsymbol{r}_{0}$ and $\boldsymbol{G}$ (the full matrix), we inadvertently change the eigenvalues of $\boldsymbol{G}_{2}$ by changing $\boldsymbol{r}_{0}$ and fixing $\boldsymbol{G}$. This becomes obvious when one finds that the two positive eigenvalues $\lambda_{\pm}$ of $\boldsymbol{G}_{2}$ is given by

[TABLE]

where $\boldsymbol{w}=\tfrac{1}{\sqrt{2}}{(1\,\,\,\,\mathrm{i})}^{\mathrm{\,\textsc{t}}}$ and $\alpha_{0}=(x_{0}+\mathrm{i}p_{0})/\sqrt{2}$. The consequence of this natural definition results in such an apparent observation. The noncentral Gaussian states so defined form the singular example in this article where this happens, and the two other noncentral non-Gaussian states which we shall soon visit do not have this technical issue.

To investigate the second-moment performance ratio $\gamma_{2}=\mathcal{H}_{\textsc{2,het}}/\mathcal{H}_{\textsc{2,hom}}$, we may reparametrize the eigenvalues of $\boldsymbol{G}$ with the squeezing strength $1\leq\lambda<\infty$ and the temperature parameter $1\leq\mu<\infty$ that is commonly adopted in describing all Gaussian states. Then $\boldsymbol{G}$ has the spectral decomposition

[TABLE]

where $\phi$ orientates the eigenvectors of $\boldsymbol{G}$. In this parametrization, we can clearly see that a large displacement magnitude contributes to a large temperature, so that a small value of $\gamma_{2}$ can be anticipated for these highly displaced or $\boldsymbol{G}_{2}$-thermal Gaussian states based on the conclusions in Řeháček et al. (2015) and Müller et al. (2016). The behavior of $\gamma_{2}$ is very similar to that for the central Gaussian states and is plotted in Fig. 2 for various values of $|\alpha_{0}|$. The lowest achievable $\gamma_{2}$ values go with the highly thermal Gaussian states ($\lambda=1$, $\mu\gg|\boldsymbol{r}_{0}|$), whose covariance matrix $\boldsymbol{G}=\mu\boldsymbol{1}/2$ is simply a multiple of the identity. Their second quadrature moment has a variance $\left<X_{\vartheta}^{4}\right>-\left<X_{\vartheta}^{2}\right>^{2}=\mu^{2}/2$, according to Eq. (41), that is of course independent of the LO phase $\vartheta$ due to the rotational symmetry. The performance ratio then takes the minimum value of $3/10$.

The maximum of $\gamma_{2}$ occurs with the coherent states ($\mu=\lambda=1$) and takes a value of $6/5$ for $\boldsymbol{r}_{0}=\boldsymbol{0}$. For larger magnitudes of $\alpha_{0}$, the value of $\gamma_{2}$ drops below unity beyond the magnitude of $|\alpha_{0}|=\sqrt{5/32}$, which can be obtained through optimization. One may verify that at this critical magnitude, $\mathcal{H}_{2,\textsc{hom}}=\mathcal{H}_{2,\textsc{het}}=63/8$. So, given a displacement magnitude larger than $\sqrt{5/32}$, the heterodyne scheme always outperforms the homodyne scheme in second-moment estimation. In the limit of large $\mu$ and $\lambda$, where we may take this limit such that $\mu=\lambda$ without loss of generality, if one considers the spectral decomposition in Eq. (47), then $\gamma_{2}$ for $\phi=0$ plotted in Fig. 3 shows the values for different $\mu$ as an indication that $\gamma_{2}\leq 1$ in this limit. Different $\phi$ values simply rotate these plots in the $x_{0}$–$p_{0}$ plane.

V.2 Fock states

Owing to the rotational symmetry of the Fock states [$\boldsymbol{G}_{2}=(n+1/2)\boldsymbol{1}$], the second and fourth quadrature moments

[TABLE]

are independent of the local-oscillator phase $\vartheta$, so that the Fisher matrix

[TABLE]

It then follows that the sCRB is given by

[TABLE]

On the other hand, the Husimi characteristic function for the Fock states in Appendix B produces the answer

[TABLE]

The performance ratio

[TABLE]

is less than one for $n\geq 2$, in which regime the Fock states are sufficiently $\boldsymbol{G}_{2}$-“thermal”. For $n=0$, we evidently obtain the familiar answer $\gamma_{2}=6/5$ for the vacuum state, whereas for $n=1$, $\gamma_{2}=16/15$. In the limit of large $n$, $\gamma_{2}\rightarrow 2/5$ (see Fig. 4).

V.3 Even/odd coherent states

Since the eigenvalues

[TABLE]

of $\boldsymbol{G}_{2}$ are simple functions of $|\alpha_{0}|^{2}$ for the even/odd $(\pm)$ coherent states, we may take $\alpha_{0}\geq 0$ without loss of generality. The quadrature moments can be easily derived with the help of Appendix B, which give the following second-moment variance

[TABLE]

By relying on the asymptotic behaviors $\coth y\approx 1/y$ and $\mathrm{cosech}\,y\approx 1/y$ of the hyperbolic trigonometric functions for small arguments, we revert to the limiting second-moment variances for $n=0$ and $n=1$, which is consistent with the fact that the even states approach the vacuum state and the odd states approach the single-photon Fock state. The Fisher matrix $\boldsymbol{F}_{2,\textsc{hom}}$ thus takes the simple form

[TABLE]

whence one obtains

[TABLE]

after carrying out the integration, matrix inversion and matrix trace. On the other hand, the Husimi-average moments of the heterodyne data contribute to the result

[TABLE]

for the heterodyne sCRB.

We once again remind the reader that the sCRBs stated in Eqs. (56) and (57) are independent of the phase of the even/odd coherent states, as this phase amounts to a rotation in phase space that is immaterial in determining the moment-estimation accuracy. For arbitrary complex values of $\alpha_{0}$, the expressions are still valid after the change $\alpha_{0}^{2}\rightarrow|\alpha_{0}|^{2}$.

The ratio $\gamma_{2}$ is greater than one for small values of $\alpha_{0}$, with the special limiting cases ($\alpha_{0}=0$) being those of the respective Fock states, and less than one for large values of $\alpha_{0}$. The crossover values for which these states become sufficiently $\boldsymbol{G}_{2}$-“thermal” such that $\gamma_{2}=1$ differ for both the even and odd states (see Fig. 5). For sufficiently large $\alpha_{0}$, $\gamma_{2}$ approaches unity from below. This can be clearly seen by taking the limit $\alpha_{0}\rightarrow\infty$. In this limit, we have $m_{\pm}\rightarrow 2\alpha_{0}^{2}=l$ so that $\mathcal{H}_{2,\textsc{hom}}\rightarrow 12\alpha_{0}^{2}\approx\mathcal{H}_{2,\textsc{het}}$. For these class of states, $\gamma_{2}$ has a stationary minimum that is again different for the two types of states, and this is elucidated in Fig. 5. At $\alpha_{0}\approx 0.631$, the $\gamma_{2}$ values for the even and odd states are equal, even though their $\boldsymbol{G}_{2}$ matrices are very different. The reason is that the combined contributions of all the second and fourth moments give an overall multiplicative factor of about 2.0694 to both $\mathcal{H}_{2,\textsc{het}}$ and $\mathcal{H}_{2,\textsc{hom}}$ for the odd state relative to the even state.

V.4 Displaced Fock states

As opposed to the previous three classes of states, the displaced Fock states (as well as the photon-added coherent states that follow) possess a nonzero quadrature first moment. As the only two parameters $\alpha_{0}=(x_{0}+\mathrm{i}p_{0})/\sqrt{2}$ and $m$ that characterize these displaced Fock states do not, in any way, restrict the covariance matrix $\boldsymbol{G}$, it is easy to show that the $\boldsymbol{G}_{2}$ geometry, and hence its reconstruction accuracy, depends only on the displacement magnitude $|\alpha_{0}|^{2}$ and not its phase. This is done by directly inspecting the eigenvalues of $\boldsymbol{G}_{2}$, namely

[TABLE]

one of which is an increasing function of $|\alpha_{0}|^{2}$. As a result, we only need to consider the case where $\alpha_{0}=x_{0}/\sqrt{2}$ is positive. As $\alpha_{0}$ increases, the GME increases, which means that the quantum state becomes more $\boldsymbol{G}_{2}$-“thermal”. We shall soon see that an increase in $|\alpha_{0}|^{2}$ results in a smaller performance ratio $\gamma_{2}$ in favor of the heterodyne scheme.

To calculate the homodyne sCRB, we first note that the relevant even-order quadrature moments (see Appendix B) supply the second-moment quadrature variance

[TABLE]

which bears striking resemblance in form with that for the even/odd coherent states, so that the sCRB also takes the same closed form as Eq. (56) inasmuch as

[TABLE]

For the heterodyne scheme, we subsequently get

[TABLE]

by again referring to Table 1.

The interplay between the discrete ($m$) and continuous $(\alpha_{0})$ parameters give rise to familiar cases that have already been analyzed previously for the Gaussian and Fock states. For $m=0$, we of course have the coherent state of amplitude $\alpha_{0}$ where the maximum $\gamma_{2}(\alpha_{0}=0)=6/5$ and the crossover point $\gamma_{2}(\alpha_{0}=\sqrt{5/32})=1$ beyond which $\gamma_{2}<1$ are reproduced by Eqs. (60) and (61). For $m=1$, we have the $m=1$ Fock state for $\alpha_{0}=0$ so that the unsurprising number $\gamma_{2}(\alpha_{0})=16/15$ comes up from the same sCRB expressions. The crossover point for $\gamma_{2}=1$ is located at $\alpha_{0}=\frac{1}{2}\sqrt{19/3-2\sqrt{87}/3}\approx 0.1696$. The performance ratio becomes subunit for all displacements $\alpha_{0}$ for $m\geq 2$, just like the Fock states. In the limit of large displacements $\alpha_{0}^{2}\gg m$, we have $m_{0}\rightarrow l$ and

[TABLE]

which approaches $1/2$ in the regime $\alpha_{0}^{2}\gg m\gg 1$.

For this two-parameter quantum state, it is interesting to look at the minimum value of $\gamma_{2}$ over all possible displacement magnitudes $\alpha_{0}$ for each $m$ [see Fig. 6(a)]. To calculate the minimum stationary points $\alpha_{0}=\widetilde{\alpha}_{0}$, we differentiate $\gamma_{2}$ with respect to $\alpha_{0}$ and set the derivative to zero. While the analytical form for the optimal $\gamma_{2}=\gamma_{2,\text{min}}$ as a complicated function of $m$ exists, the approximated forms

[TABLE]

are enough to understand the optimal-$\gamma_{2}$ curve in terms of a power law already for moderately large $m$. Interestingly, the saturation point for $\gamma_{2}$ is slightly lower than $2/5$, which is the $\gamma_{2}$ for the Fock state of an infinitely large $m$ value. This hints that the optimal center for the displaced Fock state of a large $m$ for which $\gamma_{2}=\gamma_{2,\text{min}}$ is significantly far away from the phase-space origin. This is indeed consistent with the behavior of the minimum point $\widetilde{\alpha}_{0}$, which also has a complicated closed-form expression [plotted in Fig. 6(a)], so that we only present the more useful approximated forms

[TABLE]

that highlight the main gradient features. To summarize, the minimum value of $\gamma_{2}$ essentially behaves as a power law in $m$, and the corresponding stationary minimum $\widetilde{\alpha}_{0}$ is quadratic for small $m$ and goes as a square-root curve for large $m$.

V.5 Photon-added coherent states

As in the case of the displaced Fock states, the eigenvalues of $\boldsymbol{G}_{2}$ for the photon-added coherent states,

[TABLE]

are also functions of $|\alpha_{0}|^{2}$, which correctly coincides with the physics of the second-moment estimation problem. This also means that discussing in terms of the range $\alpha_{0}\geq 0$ covers the tomography analysis sufficiently. Moreover, the eigenvalues are increasing functions of the displacement magnitude, so that the GME becomes larger with $\alpha_{0}$, thereby rendering the photon-added states more $\boldsymbol{G}_{2}$-“thermal”. This again gives a smaller performance ratio $\gamma_{2}$, or a better tomographic performance for the heterodyne scheme compared to the homodyne scheme.

Once more with the help of Table 1 in Appendix B, the quadrature moments can be written down in principle, but they are represented by bulky expressions that are hardly worth any analytical value and the Fisher-matrix integral in Eq. (11) has no known closed-form expression. However, we may still briefly discuss the important limiting cases. For $\alpha_{0}\ll\sqrt{m}$, to second order in $\alpha_{0}$, it can be shown that

[TABLE]

where the asymptotic connection with Fock states is clear. On the other hand, in the regime of large $\alpha_{0}\gg\sqrt{m}$, we find that

[TABLE]

which is the second-moment homodyne sCRB for coherent states. This is also the homodyne sCRB for large-intensity even/odd coherent states. The reason is that for large amplitudes, all these states behave like a coherent state of amplitude $\alpha_{0}$ as far as second-moment estimation is concerned since all their $\boldsymbol{G}_{2}$ eigenvalues are indistinguishable in this limit.

Upon revisiting Eq. (24), the heterodyne sCRB can be shown to have the closed form

[TABLE]

for $\alpha_{0}>0$. The behavior to leading order in $\alpha_{0}^{2}$ for $\alpha_{0}\ll\sqrt{m}$,

[TABLE]

is evidently consistent with the known result for Fock states. For $\alpha_{0}\gg\sqrt{m}$, we once again have $\mathcal{H}_{2,\textsc{het}}\approx 12\alpha_{0}^{2}$. Note also that, as expected, the first equality in Eq. (67) is equal to 5, and Eq. (67) gives a value of 6 for the vacuum state ($m=0$).

For $m\geq 2$, the ratio $\gamma_{2}<1$ for all $\alpha_{0}$. This natural extension to the result for the Fock states means that for highly nonlinear photon-“adding” operations, the performance of heterodyne detection is always better than that of homodyne detection in terms of second-moment covariance-dyadic estimation. For $m=0$, the analysis reverts to that for the coherent state, where the crossover occurs at $\alpha_{0}=\sqrt{5/32}$ after solving for $\mathcal{H}_{2,\textsc{hom}}=\mathcal{H}_{2,\textsc{het}}=2\left(3+6\alpha_{0}^{2}\right)$ so that $\gamma_{2}(\alpha_{0}>\sqrt{5/32})<1$, which is again consistent with Sec. V.1. For $m=1$, the crossover point $\alpha_{0}\approx 0.2$ may be obtained as the numerical solution. As $\alpha_{0}$ approaches infinity, previous arguments imply that $\gamma_{2}\rightarrow 1$ for any $m$.

In view of the behavior of $\gamma_{2}$, another interesting limit is the high-nonlinearity limit ($m\rightarrow\infty$). In this case, we notice that the value $\alpha_{0}=\widetilde{\alpha}_{0}$ for which $\gamma_{2}$ is minimum approaches zero. A good model to estimate this minimum point in this limit is given by

[TABLE]

which can be approximated from curve fitting. Therefore in the large-$m$ limit, the optimum performance ratio $\gamma_{2}$ is that of an intense Fock state of a large photon number and so we expect the minimum value of $\gamma_{2}$ to approach $2/5$ as discussed in Sec. V.2. In other words, for sufficiently large $m$, the minimum of $\gamma_{2}$ follows the noncentral power law

[TABLE]

Figure 6(b) succinctly highlights these observations.

VI Conclusion

We compare the moment-reconstruction performances of the homodyne and heterodyne joint-measurement measurement schemes using optimal moment estimators that minimizes the mean squared-error. We first showed that in first-moment tomography, the heterodyne scheme is always tomographically superior to, or at least as good as, the homodyne scheme for all quantum states in terms of the mean squared error of the moment estimators. The underlying physical reason is solely the Heisenberg-Robertson-Schrödinger uncertainty relation for complementary observables. For second-moment tomography, we showed that the heterodyne scheme can often outperform the homodyne scheme for Gaussian states and many other interesting and important classes of non-Gaussian states. All these states indicate a trend that a larger geometric mean of second-moment eigenvalues (second-moment “temperature”) improves the moment reconstruction accuracy with the heterodyne scheme relative to the homodyne scheme. This trend, however, is not monotonic in the second-moment “temperature”, because there is also influence from the fourth moments originating from the form of the mean squared-error, the combined contributions of both give interesting features the reconstruction accuracy, as illustrated by the examples in this article. The general theory introduced in Sec. III can be applied to higher-moment estimation that are important in general operator-moment applications and source-calibration protocols, and these shall be reported in the future.

VII Acknowledgments

We acknowledge financial support from the BK21 Plus Program (21A20131111123) funded by the Ministry of Education (MOE, Korea) and National Research Foundation of Korea (NRF), the NRF grant funded by the Korea government (MSIP) (Grant No. 2010-0018295), the Korea Institute of Science and Technology Institutional Program (Project No. 2E26680-16-P025), the European Research Council (Advanced Grant PACART), the Spanish MINECO (Grant FIS2015-67963-P), the Grant Agency of the Czech Republic (Grant No. 15-03194S), and the IGA Project of the Palacký University (Grant No. IGA PrF 2016-005).

Appendix A Optimal estimators for homodyne tomography

A.1 First-moment estimation

In this discussion, the reconstruction accuracy of the estimator $\widehat{\boldsymbol{r}}_{\textsc{hom}}$ for $\boldsymbol{r}$ shall be taken to be the usual MSE distance measure

[TABLE]

that is typically defined for columns. One straightforward way to obtain an estimator $\widehat{\boldsymbol{r}}_{\textsc{hom}}$ is to make use of $\left<X_{\vartheta}\right>=\left<X\right>\cos\vartheta+\left<P\right>\sin\vartheta$ to ascertain that

[TABLE]

for an $n_{\vartheta}\times 2$ matrix $\boldsymbol{L}_{\vartheta}$ ($n_{\vartheta}$ being the number of bins for the LO phases $\vartheta$) and a column $\boldsymbol{r}_{\vartheta}$ of $n_{\vartheta}$ true averages $\left<X_{\vartheta}\right>$. The highly overcomplete nature of the measurement thus permits us to define, for any experimentally obtained estimates of average values $\boldsymbol{r}_{\vartheta}\equiv\left(\widehat{\left<X_{1}\right>}\,\,\widehat{\left<X_{2}\right>}\,\,\ldots\,\,\widehat{\left<X_{n_{\vartheta}}\right>}\right)^{\textsc{t}}$ $\left(\overline{\widehat{\left<X_{k}\right>}}=\left<X_{\vartheta_{k}}\right>\right)$,

[TABLE]

as the linear estimator of interest using the pseudoinverse $\boldsymbol{L}_{\vartheta}^{-}$ of $\boldsymbol{L}_{\vartheta}$. This estimator, however, is suboptimal in the sense that it does not minimize the MSE $\mathcal{D}_{1,\textsc{hom}}$.

To obtain the best estimator for $\boldsymbol{r}$ [often known as the linear unbiased estimator (BLUE)] that minimizes the MSE, we resort to the linear optimization of

[TABLE]

over all possible reconstruction columns $\boldsymbol{v}_{k}$ for the estimates $\widehat{\left<X_{k}\right>}$. Data consistency according to $\widehat{\left<X_{k}\right>}=\boldsymbol{u}_{k}^{\textsc{t}}\widehat{\boldsymbol{r}}_{\textsc{hom}}$ requires these reconstruction columns, or dual columns, to satisfy the property

[TABLE]

with the measurement columns $\boldsymbol{u}_{k}=\boldsymbol{u}_{\vartheta_{k}}={(\cos\vartheta_{k}\,\,\sin\vartheta_{k})}^{\mathrm{\,\textsc{t}}}$. Logically, we must have

[TABLE]

The Lagrange function for the optimization is therefore

[TABLE]

where $\boldsymbol{\Lambda}$ is the Lagrange matrix for the dual-column constraint in (76). In terms of the dual columns,

[TABLE]

Since the unbiased estimate

[TABLE]

is an average sum of all the measured $n_{x}$ voltage readings $x_{jk}$ per LO phase that are distributed according to the multinomial distribution of random multinomial weights $\sum_{j}n_{jk}=N_{k}$, the second moment is given by

[TABLE]

The final equality is valid for sufficiently large data (bins) for all phases, as $p_{jk}\rightarrow\mathrm{d}x_{\vartheta}\,p(x_{\vartheta},\vartheta)$ and

[TABLE]

So, we finally get

[TABLE]

A simple variation of $\mathcal{L}_{\textsc{hom}}$ therefore gives

[TABLE]

or

[TABLE]

The matrix $\boldsymbol{F}\!\left(\{\left<X_{k}\right>,\left<X_{k}^{2}\right>\}\right)$ is known as the frame matrix.

The BLUE therefore depends on the true moments which are certainly unavailable in the first place, for no tomography is otherwise necessary at all. Nonetheless, one can substitute the estimated moments for them to obtain an asymptotically efficient optimal estimator that approximates the BLUE. An unbiased estimate for the second moment is given by

[TABLE]

so that the asymptotically optimal estimator is given by

[TABLE]

It is easy to see that when the estimated moments approach the true moments, this optimal estimator attains the sCRB. Directly from Eq. (83), we immediately know that its corresponding MSE is given by

[TABLE]

and all we need to realize is that for sufficiently large $N$ and uniformly-distributed quadrature outcomes, $N_{k}/N\rightarrow\mathrm{d}\vartheta/\pi$ and the frame matrix

[TABLE]

is nothing more than the Fisher matrix introduced in Eq. (5). This also means that the BLUE and the asymptotically optimal estimator are both asymptotically as efficient as the ML estimator.

This construction comes with a basic and important lesson. The simple linear estimator $\widehat{\boldsymbol{r}}^{(\textsc{lin})}_{\textsc{hom}}$ in Eq. (74), which is suboptimal, depends only on the first moments. To improve the reconstruction accuracy, more aspects of the data that are attributed to the figure of merit chosen to measure this accuracy would have to be incorporated systematically. In the case of the MSE, these are linear combinations of both the first and second moments, or at least their estimates. Put differently, we should always use the reconstruction estimator that optimize the figure of merit we choose to rank the goodness of the reconstruction.

A.2 Second-moment estimation

By the same token, we can construct the optimal estimator that approximates the BLUE for second-moment estimation by minimizing the MSE

[TABLE]

over the estimator that is of the linear form

[TABLE]

with respect to the second-moment estimates. This form is a natural extension to the column estimator $\widehat{\boldsymbol{r}}_{\textsc{hom}}$ via a generalization of the dual columns $\boldsymbol{v}_{k}$ to dual matrices $\boldsymbol{\Theta}_{k}$. Completely analogous to the discussion in Appendix A.1, consistency with $\widehat{\left<X^{2}_{k}\right>}=\boldsymbol{u}_{k}^{\textsc{t}}\widehat{\boldsymbol{G}}_{2,\textsc{hom}}\boldsymbol{u}_{k}$ implies that

[TABLE]

The above relation can be simplified by introducing the vectorization notation $\mathrm{vec}\!\left(\boldsymbol{Y}\right)$ that turns a matrix $\boldsymbol{Y}$ into a column. Since all two-dimensional matrices considered here are real and symmetric, they are essentially characterized by three real parameters. Hence in our context, given that

[TABLE]

the vectorized quantity is defined as

[TABLE]

This operation is a variant of the usual column-stacking vectorization operation to apply on $2\times 2$ real symmetric matrices for our case to make contact with the property $\mathrm{Tr}\!\left\{\boldsymbol{Y}_{1}\boldsymbol{Y}_{2}\right\}={\mathrm{vec}\!\left(\boldsymbol{Y}_{1}\right)}^{\mathrm{\,\textsc{t}}}\mathrm{vec}\!\left(\boldsymbol{Y}_{2}\right)$ between any pair of such matrices $\boldsymbol{Y}_{1}$ and $\boldsymbol{Y}_{2}$. In this notation, Eq. (92) becomes

[TABLE]

which is equivalent to the vectorized constraint

[TABLE]

As usual, to derive the expression for the optimal estimator, we first calculate $\mathcal{D}_{2,\textsc{hom}}$ in terms of the dual matrices. For this we shall need the average of the square of the estimate $\widehat{\left<X_{k}^{2}\right>}$ defined in Eq. (86):

[TABLE]

from which gives the expression

[TABLE]

for the MSE. Then, by carrying out the variation of the appropriate Lagrange function similar to the calculations in Appendix (A.1) and remembering the additional association $\boldsymbol{M}_{k}=\mathrm{vec}\!\left(\boldsymbol{m}_{k}\right)\mathrm{vec}\!\left(\boldsymbol{m}_{k}\right)^{\textsc{t}}$, we find that the optimal matrices for the BLUE are

[TABLE]

Finally, the asymptotically optimal estimator is given by

[TABLE]

That this estimator asymptotically attains the sCRB for second-moment estimation is again clear.

Appendix B List of characteristic functions

In calculating the moments for both the homodyne and heterodyne schemes, it is extremely useful to start with the relevant characteristic functions for both schemes. To facilitate the discussions in the main article, we have supplied a list of quadrature characteristic functions \left(\left<\mathrm{e}^{\mbox{\footnotesize\mathrm{i}kX_{\vartheta}}}\right>\right) for the homodyne scheme and a list of Husimi characteristic functions \displaystyle\left[\overline{\mathrm{e}^{\mbox{\footnotesizeg^{}\alpha+g\alpha^{}}}}\,,\,\,g=(u+\mathrm{i}v)/\sqrt{2}\right] for the heterodyne scheme respectively in Table 1 in this appendix section. Then the two kinds of moments can then be readily computed by the prescriptions

[TABLE]

which simply involves multiple differentiations with respect to the free variables and later setting these variables to zero. Some useful identities for the confluent hypergeometric functions and Laguerre polynomials that allow for consistency verification between two characteristic functions of different quantum states are given below:

[TABLE]