Towards optimal quantum tomography with unbalanced homodyning

TL;DR

This paper demonstrates that localized unbalanced homodyning provides more accurate quantum state tomography than heterodyning or balanced homodyning, even with realistic detectors, advancing towards optimal continuous-variable quantum measurement.

Contribution

It establishes that unbalanced homodyning is universally more powerful for quantum tomography than heterodyning across all states and schemes, including realistic detector conditions.

Findings

Unbalanced homodyning outperforms heterodyning in all tested scenarios.

Localized sampling yields more accurate estimators than delocalized methods.

The advantage persists with subunit efficiency photodetectors.

Abstract

Balanced homodyning, heterodyning and unbalanced homodyning are the three well-known sampling techniques used in quantum optics to characterize all possible photonic sources in continuous-variable quantum information theory. We show that for all quantum states and all observable-parameter tomography schemes, which includes the reconstructions of arbitrary operator moments and phase-space quasi-distributions, localized sampling with unbalanced homodyning is always tomographically more powerful (gives more accurate estimators) than delocalized sampling with heterodyning. The latter is recently known to often give more accurate parameter reconstructions than conventional marginalized sampling with balanced homodyning. This result also holds for realistic photodetectors with subunit efficiency. With examples from first- through fourth-moment tomography, we demonstrate that unbalanced homodyning can outperform balanced homodyning when heterodyning fails to do so. This new benchmark takes us one step towards optimal continuous-variable tomography with conventional photodetectors and minimal experimental components.

Equations

$$\mathbf{q}=\int \frac{(\mathrm{d}{\alpha }^{\prime })}{\pi }W\left({\alpha }^{\prime }\right){\mathbf{v}}_{\text{textsc}w}({\alpha }^{\prime })=\int \frac{(\mathrm{d}{\alpha }^{\prime })}{\pi }Q\left({\alpha }^{\prime }\right){\mathbf{v}}_{\text{textsc}p}({\alpha }^{\prime })$$

$${\mathrm{sCRB}}_{\text{textsc}het}=$$

$${\mathrm{sCRB}}_{\text{textsc}uhom}=$$

$${\mathrm{sCRB}}_{\text{textsc}het}^{\prime }=$$

$${\mathrm{sCRB}}_{\text{textsc}uhom}^{\prime }=$$

$${\mathrm{sCRB}}_{\text{textsc}bhom}=$$

$$\times [{\mathbf{v}}_{\text{textsc}w}({\alpha }^{\prime })-\mathbf{q}]\cdot [{\mathbf{v}}_{\text{textsc}w}({\alpha }^{\prime \prime })-\mathbf{q}],$$

$$w_{{\alpha }^{\prime },{\alpha }^{\prime \prime }}=$$

$$\times \left[p(x_{\vartheta }^{\prime },\vartheta )\delta (x_{\vartheta }^{\prime }-x_{\vartheta }^{\prime \prime })-p(x_{\vartheta }^{\prime },\vartheta )p(x_{\vartheta }^{\prime \prime },\vartheta )\right].$$

$$\mathbb{E}\left[\frac{n_l}{N}\right]=$$

$$\mathbb{E}\left[\frac{n_l{n}_{l^{\prime }}}{N^2}\right]=$$

$$+\frac{3p_l{p}_{l^{\prime }}}{({\sum }_l{p}_l{)}^4}\sum _{l^{\prime \prime }}{\sigma }_{l^{\prime \prime }}^2,$$

$$\dfrac{1}{A}=-\mathrm{i}\int^{\infty}_{0}\mathrm{d}t\,\mathrm{e}^{\mbox{\footnotesize$\mathrm{i}tA-\epsilon t$}}\bigg{|}_{\epsilon=0}$$

$$\mathbb{E}\left[\frac{f_l}{{\sum }_l{f}_l}\right]=$$

$$\mathbb{E}\left[\frac{f_l{f}_{l^{\prime }}}{({\sum }_l{f}_l{)}^2}\right]=$$

$$\mathbb{E}\!\left[\mathrm{e}^{\mbox{\footnotesize$\mathrm{i}\sum_{l}\lambda_{l}f_{l}$}}\right]=\prod_{l}\mathbb{E}\!\left[\mathrm{e}^{\mbox{\footnotesize$\mathrm{i}\lambda_{l}f_{l}$}}\right]\,,$$

$$\mathbb{E}\!\left[\mathrm{e}^{\mbox{\footnotesize$\mathrm{i}\lambda_{l}f_{l}$}}\right]=\mathrm{e}^{\mbox{\footnotesize$-\frac{1}{2}\lambda_{l}^{2}\sigma_{l}^{2}+\mathrm{i}\lambda_{l}p_{l}$}}$$

$$\mathbb{E}\!\left[\dfrac{f_{l}}{\sum_{l}f_{l}}\right]=-\int^{\infty}_{0}\mathrm{d}t\,\mathrm{e}^{\mbox{\footnotesize$-at^{2}+\mathrm{i}bt$}}\,(\mathrm{i}p_{l}-t\sigma_{l}^{2})\,,$$

$$\mathbb{E}\left[\frac{f_l{f}_{l^{\prime }}}{{\left({\sum }_l{f}_l\right)}^2}\right]=$$

$$\times \left[c_1(t+t^{\prime }{)}^2-\mathrm{i}{c}_2(t+t^{\prime })-c_3\right],$$

$$\sqrt{\pi}\,\mathrm{e}^{\mbox{\footnotesize$-\frac{b^{2}}{4a}$}}\,\mathrm{erfi}\!\left(\frac{b}{2\sqrt{a}}\right)\approx\frac{24a^{5/2}}{b^{5}}+\frac{4a^{3/2}}{b^{3}}+\frac{2\sqrt{a}}{b}$$

$$\mathbb{E}\left[\frac{n_l{n}_{l^{\prime }}}{N}\right]=\frac{{\sigma }_l^2{\delta }_{l,l^{\prime }}+p_l{p}_{l^{\prime }}}{{\sum }_l{p}_l}-\frac{p_l{\sigma }_{l^{\prime }}^2+p_{l^{\prime }}{\sigma }_l^2}{({\sum }_l{p}_l{)}^2}+\frac{p_l{p}_{l^{\prime }}}{({\sum }_l{p}_l{)}^3}\sum _{l^{\prime \prime }}{\sigma }_{l^{\prime \prime }}^2,$$

$$\mathbb{E}\left[\frac{f_l{f}_{l^{\prime }}}{{\sum }_l{f}_l}\right]=$$

$$\mathbb{E}\left[\frac{n_{jk}}{\mathcal{N}}\right]=$$

$$\mathbb{E}\left[\frac{n_{jk}{n}_{j^{\prime }{k}^{\prime }}}{{\mathcal{N}}^2}\right]=$$

$$-\frac{2}{b^3}(p_{jk}{w}_{j^{\prime }{k}^{\prime }}+p_{j^{\prime }{k}^{\prime }}{w}_{jk})$$

$$\mathbb{E}\left[\frac{f_{jk}}{{\sum }_{l,j,k}{\mathcal{R}}_{ljk}^{-1}{f}_{jk}}\right]=$$

$$=$$

$$\mathbb{E}\left[\frac{f_{jk}{f}_{j^{\prime }{k}^{\prime }}}{{\left({\sum }_{l,j,k}{\mathcal{R}}_{ljk}^{-1}{f}_{jk}\right)}^2}\right]$$

$$=$$

$$\mathbb{E}\!\left[\mathrm{e}^{\mbox{\footnotesize$\mathrm{i}\sum_{l,j,k}\lambda_{ljk}f_{jk}$}}\right]=\prod^{n_{\vartheta}}_{j=1}\mathbb{E}\!\left[\mathrm{e}^{\mbox{\footnotesize$\mathrm{i}\sum_{l,k}\lambda_{ljk}f_{jk}$}}\right]$$

$$\mathbb{E}\!\left[\mathrm{e}^{\mbox{\footnotesize$\mathrm{i}\sum_{l,k}\lambda_{ljk}f_{jk}$}}\right]=\mathrm{e}^{\mbox{\footnotesize$-\frac{1}{2}\sum_{l,l^{\prime},k,k^{\prime}}\lambda_{ljk}\lambda_{l^{\prime}jk^{\prime}}\Sigma_{jkk^{\prime}}+\mathrm{i}\sum_{l,k}\lambda_{ljk}p_{jk}$}}$$

$$\mathbb{E}\!\left[\dfrac{f_{jk}}{\sum_{l,j,k}\mathcal{R}^{-1}_{ljk}f_{jk}}\right]=-\int^{\infty}_{0}\mathrm{d}t\,\mathrm{e}^{\mbox{\footnotesize$-at^{2}+\mathrm{i}bt$}}\,(\mathrm{i}p_{jk}-tw_{jk})\,,$$

$$\mathbb{E}\left[\frac{f_{jk}{f}_{j^{\prime }{k}^{\prime }}}{{\left({\sum }_{l,j,k}{\mathcal{R}}_{ljk}^{-1}{f}_{jk}\right)}^2}\right]=$$

$$\times \left[c_1(t+t^{\prime }{)}^2-\mathrm{i}{c}_2(t+t^{\prime })-c_3\right],$$

$$\mathbb{E}\left[\exp (\mathrm{i}(k_l{z}_l+k_{l^{\prime }}{z}_{l^{\prime }}))\right]=\mathbb{E}\left[\exp (\mathrm{i}{k}_l{z}_l)\right]\mathbb{E}\left[\exp (\mathrm{i}{k}_{l^{\prime }}{z}_{l^{\prime }})\right]$$

$$\frac{1}{({\sum }_l{n}_l{)}^{m+m^{\prime }}}\approx \frac{(m+m^{\prime }+1)N_0{\sum }_l{p}_l-(m+m^{\prime }){\sum }_l{n}_l}{(N_0{\sum }_l{p}_l{)}^{m+m^{\prime }+1}}.$$

Full text

Introduction.—In pursuing a secure information age, the successful implementations of state-of-the-art continuous-variable (CV) quantum information and communication protocols Braunstein and van Loock (2005); Ferraro et al. (2005); Cerf et al. (2007); Andersen et al. (2010); Adesso et al. (2014) require precise reconstructions and calibrations of important properties of photonic sources. In the language of phase-space quasi-distributions that completely characterize such sources, these properties—generally the expectation values of quantum observables—can be reconstructed with either physical probabilities of a positive distribution or those derived from some aspects of a (non-singular) quasi-distribution.

There are three sampling methods considered in quantum optics that identify these two main scenarios. The first and arguably the most popular method is marginalized phase-space sampling by balanced homodyne detection (BHOM) Yuen and Chan (1983); Abbas et al. (1983); Schumaker (1984); Vogel and Risken (1989); Banaszek and Wódkiewicz (1997), which samples the marginal distributions of the Wigner function defined by quadrature directions. This requires a balanced (1:1) beam splitter, local oscillator (LO), and two photodetectors to measure photocurrent differences at the output. The second method is delocalized phase-space sampling executed with heterodyne detection (HET) that jointly measures complementary quadrature operators Arthurs and Kelly (1965); Yuen (1982); Arthurs and Goodman (1988); Martens and de Muynck (1990, 1991); Stenholm (1992); Raymer (1994); Trifonov et al. (2001); Werner (2004). This technique randomly samples the whole phase space according to the Husimi function, and usually involves a more sophisticated setup of three balanced beam splitters, LO and four photodetectors to realize such a double-BHOM scheme. The third sampling method of focus here is localized phase-space sampling with unbalanced homodyne detection (UHOM) Wallentowitz and Vogel (1996); Banaszek and Wódkiewicz (1996); Opatrný and Welsch (1997); Opatrný et al. (1997); Wallentowitz et al. (2012); Kühn and Vogel (2016); Paris (1996), which measures displaced Fock states using a highly-transmissive beam splitter, LO and two photodetectors such that the Wigner function can be directly reconstructed through the parity-operator measurement. With common photodetectors that have no photon-number resolution capabilities, this method samples the Husimi function by counting “no-click” events at the transmission arm of the signal. The displacement operation by the unbalanced beam splitter then guarantees coherent-state measurements of specified amplitudes, which data follow a binomial distribution characterized by the Husimi function at each amplitude.

The understanding of the parameter reconstruction accuracies for all sampling methods holds a fundamental link to the tomographic power of quantum measurement schemes. There exist a plethora of articles Yuen (1982); Trifonov et al. (2001); Werner (2004); Kiesel (2012) that investigated variances and measurement uncertainties, which supply information about important statistical behaviors of parameter estimators. For the purpose of analyzing tomographic power, optimality analysis for true-parameter reconstructions is in order. Recently, the relationship between the Haar-averaged Cramér–Rao bound for state estimation and the permutation group was studied in Řeháček et al. (2015); Koutný et al. (2016). In Řeháček et al. (2015); Müller et al. (2016); Teo et al. (2017), we systematically analyzed the tomographic power of both BHOM and HET using the notions of the Fisher information and Cramér–Rao bound for moment estimation and found that the latter gives higher reconstruction accuracies for typically interesting states, with Gaussian states being one important class in CV quantum information processing Lorenz et al. (2004); Ferraro et al. (2005); Lance et al. (2005); Řeháček et al. (2009); Scarani et al. (2009); Weedbrook et al. (2012). This provided irrefutable evidence of tomographic differences in parameter reconstruction for the Wigner and Husimi representations, despite their equivalence in state representation.

It will be shown here that for every used reconstruction datum, localized sampling with UHOM is always tomographically more powerful than delocalized sampling with HET for any type of observable-parameter tomography. This benefit originates from the statistical nature of the UHOM data collected at each phase-space value. We shall demonstrate that this effect can even result in a superior tomographic power over BHOM in some cases where HET is inferior. These two main results are analyzed for first- through fourth-order moment tomography with Gaussian and Fock states.

Parameters and tomographic power.—For a more concrete concept of comparing different measurement schemes, we consider the statistical mean squared error (MSE) $\mathbb{E}\!\left[(\widehat{\boldsymbol{q}}-\boldsymbol{q})^{2}\right]$ for any column $\boldsymbol{q}$ of parameters and its estimator $\widehat{\boldsymbol{q}}$. This accuracy measure is a function of both the measurement and data for $\widehat{\boldsymbol{q}}$. We consider observable parameters of the kind $\boldsymbol{q}=\left<\boldsymbol{V}\right>$ for an arbitrary column $\boldsymbol{V}$ of observables describing some list of quantum properties (which is always a function of the position $X$ and momentum $P$ operators Englert (2006)), where equivalently Cahill and Glauber (1969)

[TABLE]

is the phase-space average of either the Glauber-Sudarshan function $\boldsymbol{v}_{\textsc{p}}(\alpha)$ or the Wigner function \boldsymbol{v}_{\textsc{w}}(\alpha)=2\int(\mathrm{d}\alpha^{\prime})\,\mathrm{e}^{\mbox{\footnotesize-2|\alpha-\alpha^{\prime}|^{2}}}\boldsymbol{v}_{\textsc{p}}(\alpha^{\prime})/\pi for $\boldsymbol{V}$, respectively with the Husimi and Wigner functions for the state $\rho$ [$0\leq Q\!\left(\alpha\right)\leq 1$, $-2\leq W\!\left(\alpha\right)\leq 2$ according to the definitions for Eq. (1)]. Each point $(x,p)$ in phase space is expressed by $\alpha=(x+\mathrm{i}p)/2$, and $(\mathrm{d}\alpha)=\mathrm{d}x\,\mathrm{d}p/2$. Equation (1) covers all the interesting tomography problems. For instance, in second-moment tomography where $\boldsymbol{V}$ consists of symmetrically-ordered products of $X$ and $P$, then $\boldsymbol{v}_{\textsc{w}}(\alpha^{\prime})={(x^{\prime 2},x^{\prime}p^{\prime},p^{\prime 2})}^{\mathrm{\,\textsc{t}}}$ and $\boldsymbol{v}_{\textsc{p}}(\alpha^{\prime})={(x^{\prime 2}-1/2,x^{\prime}p^{\prime},p^{\prime 2}-1/2)}^{\mathrm{\,\textsc{t}}}$. As another example, if one is interested in Wigner function reconstruction, then $\boldsymbol{v}_{\textsc{w}}(\alpha^{\prime})=\delta(\boldsymbol{\alpha}-\alpha^{\prime})$ and $\boldsymbol{v}_{\textsc{p}}(\alpha^{\prime})$ is the kernel for the Gaussian deconvolution. For an $s$-ordered quasi-distribution, $\boldsymbol{v}_{\textsc{w}}(\alpha^{\prime})$ and $\boldsymbol{v}_{\textsc{p}}(\alpha^{\prime})$ are the relevant $s$-ordered kernels.

The detection schemes for all three sampling methods have very different kinds of data. The data sample size $N$ for BHOM is the total number of marginalized Wigner data points defined by the sampled LO phases and real voltages. For HET and UHOM, $N$ is the total number of randomly sampled phase-space values, but since only the “no-click” data for the D1 photodetector (out of a fixed number of sampling events per $\alpha$) at the transmitted arm for the signal (see Fig. 1) are used in the reconstruction, $N$ becomes a sum of binomial random integers, and is itself random. To make a fair comparison of the three schemes, the well-known (scaled) Cramér–Rao bound $\mathrm{sCRB}=\min_{\widehat{\boldsymbol{q}}}\{\mathbb{E}\!\left[N(\widehat{\boldsymbol{q}}-\boldsymbol{q})^{2}\right]\}$ is a good measure for the tomographic power of the measurement. This scaled measure consistently weights each experiment with its total sample size to average away the data aspect, and is minimized over all conceivable reconstruction strategies for $\widehat{\boldsymbol{q}}$ of some given data type. A smaller sCRB implies a greater tomographic power. For BHOM and HET, $N$ is usually a fixed constant, so that the sCRB turns into the familiar MSE per reconstruction datum. For sufficiently large coherent-state data $N$ and densely sampled phase-space points, one can show that the MSE for UHOM goes as the average shot-noise limit $\left(\propto 1/\mathbb{E}\!\left[N\right]\right)$, which again reminds us that the accuracy of $\widehat{\boldsymbol{q}}$ varies only with the used reconstruction data sample size as always. It then follows that $\min_{\widehat{\boldsymbol{q}}}\{\mathbb{E}\!\left[N(\widehat{\boldsymbol{q}}-\boldsymbol{q})^{2}\right]\}=\mathbb{E}\!\left[N\right]\min_{\widehat{\boldsymbol{q}}}\{\mathbb{E}\!\left[(\widehat{\boldsymbol{q}}-\boldsymbol{q})^{2}\right]\}$. This means that while the comparisons of the sampling methods are made by scaling away the used reconstruction data, it should not matter whether this scaling is done for every experiment or with an overall average data cost for all the experiments. Any physically meaningful definition of the tomographic power should be invariant under such a technical variation.

Main results.—Both the delocalized (HET) and localized (UHOM) phase-space sampling methods share a common trait: their data $N=\sum_{l}n_{l}$ directly reconstruct the Husimi function $Q\!\left(\alpha\right)$: $n_{l}/\sum_{l}n_{l}\approx(\updelta\alpha)Q\!\left(\alpha_{l}\right)/\pi$ at a sampled $\alpha=\alpha_{l}$ for some small pre-chosen area $(\updelta\alpha)$ of the sampled discretized phase-space. Therefore the Husimi representation of $\boldsymbol{q}$ in (1) invites an estimator of the form $\widehat{\boldsymbol{q}}=\sum_{l}n_{l}\,\boldsymbol{v}_{\textsc{p}}(\alpha_{l})/\sum_{l}n_{l}$. Furthermore, it can be shown that such a sample average estimator, for these HET and UHOM data, follows a multivariate Gaussian distribution in the limit of large $N$ with the correct mean $\boldsymbol{q}$ and data covariance, so that $\widehat{\boldsymbol{q}}$ achieves the sCRB asymptotically. After a proper statistical analysis for $N\gg 1$ and a densely sampled $\alpha_{l}$s, we have the intuitively simple expressions (see App. A)

[TABLE]

By inspection, since $Q\!\left(\alpha\right)\left[1-Q\!\left(\alpha\right)\right]\leq Q\!\left(\alpha\right)$ for any $\rho$, we immediately find that $\mathrm{sCRB}_{\textsc{uhom}}<\mathrm{sCRB}_{\textsc{het}}$. This first general result has important physical implications. It shows that localized sampling always reduces the magnitude of the phase-space distribution via the binomial deformation $Q\!\left(\alpha\right)\rightarrow Q\!\left(\alpha\right)\left[1-Q\!\left(\alpha\right)\right]$. This leads to a smaller combined reconstruction variance per datum for any list of parameters $\boldsymbol{q}$ relative to HET. In practice, the tomographic advantages of localized binomial phase-space sampling is realized with only a replacement of the balanced BS with a highly-transmissive BS, which is a minor adjustment of the BHOM setup in Fig. 1(b). For realistic photodetectors of efficiency $0<\eta\leq 1$, using the definition p(\alpha,\eta)=\left<\bm{:}\mathrm{e}^{\mbox{\footnotesize-\eta(a^{\dagger}-\alpha)(a-\alpha)}}\bm{:}\right>\leq 1 where $a$ is the usual photonic ladder operator and $\bm{:}\,\cdot\,\bm{:}$ denotes operator normal ordering, the sampled probabilities with HET are given by $\eta p(\alpha,\eta)$ while the binomial probability for “no-click” events at photodetector D1 for UHOM is $p(\alpha,\eta)$ Wallentowitz and Vogel (1996). These give the more realistic bounds (see App. C)

[TABLE]

that satisfy $\mathrm{sCRB}^{\prime}_{\textsc{uhom}}<\mathrm{sCRB}^{\prime}_{\textsc{het}}$ 111The bounds should diverge as $\eta\rightarrow 0$, which means that one cannot swap this limit with the phase-space integration.

For an arbitrary $\boldsymbol{V}$ that is a complicated function of $X$ and $P$, the general recipe for $\widehat{\boldsymbol{q}}$ with BHOM data is to adopt the Wigner representation in (1) and estimate $W\!\left(\alpha\right)$ by an application of the inverse Radon transform (invR) to the BHOM probabilities. Upon denoting the invR kernel by $\mathcal{R}^{-1}_{\alpha}(x_{\vartheta},\vartheta)=\int\mathrm{d}k|k|\exp(\mathrm{i}k(x\cos\vartheta+p\sin\vartheta-x_{\vartheta}))$ for a given LO phase $\vartheta$ and voltage $x_{\vartheta}$, the corresponding estimator for the BHOM data is given by $\widehat{\boldsymbol{q}}=\sum_{l,j,k}\mathcal{R}^{-1}_{\alpha_{l}}(x_{j},\vartheta_{k})n_{jk}\,\boldsymbol{v}_{\textsc{w}}(\alpha_{l})/\sum_{l,j,k}\mathcal{R}^{-1}_{\alpha_{l}}(x_{j},\vartheta_{k})n_{jk}$, where $n_{jk}/\sum_{j}n_{jk}$ estimates the BHOM probability $\mathrm{d}x\,p(x_{j},\vartheta_{k})$. The tomographic power of BHOM for this general recipe with invR (only one kind of estimator considered here) is measured by

[TABLE]

In practice, the estimation of the Wigner function $W\!\left(\alpha\right)$ is done in a truncated Hilbert space. As such, a direct application of the invR on the measured BHOM probabilities typically gives rise to $W\!\left(\alpha\right)$ with truncation phase-space wriggles that are otherwise absent in the infinite-dimensional limit. Together with the high sensitivity of invR to statistical fluctuations, $\mathrm{sCRB}_{\textsc{bhom}}$ is in general greater than either $\mathrm{sCRB}_{\textsc{uhom}}$ or $\mathrm{sCRB}_{\textsc{het}}$ (see Fig. 2). Thus, for general parameters where $\boldsymbol{V}$ is a complicated function of $X$ and $P$, UHOM is the best option. Although it is known that the maximum-likelihood method can reduce such reconstruction instabilities  Banaszek (1999a, b), analytical tomographic studies of such a nonlinear numerical method still form an open problem.

Certainly, a much more expedient and trusted way to estimate $\boldsymbol{q}$ (referred to as the BHOMOPT strategy) when $\boldsymbol{V}$ is a simple function of $X$ and $P$ is a direct and optimized data-processing strategy of the measured voltage values for every LO phase $\vartheta$ such that $\boldsymbol{q}$ can be efficiently reconstructed without having to go through any formal invR. For instance, in moment tomography Řeháček et al. (2009); Müller et al. (2016); Teo et al. (2017), the entries of $\boldsymbol{q}$ are linearly related to the moments of the quadrature operator, $\left<X^{m}_{\vartheta}\right>=\left<(X\cos\vartheta+P\sin\vartheta)^{m}\right>$, sampled by BHOM. Therefore, $\mathrm{sCRB}_{\textsc{bhomopt}}$ for any state using this improved reconstruction strategy can be obtained through the Fisher information of the homodyne parameter $\left<X^{m}_{\vartheta}\right>$. The theory for this was developed in Teo et al. (2017). It shall be shown that in practice, UHOM is tomographically more powerful than all other methods for moment tomography of interesting states, which forms the second main result.

Moment-tomography analysis.—We demonstrate the tomographic power of UHOM with moment tomography of orders $m=1$ through $m=4$. In particular, we study symmetrically-ordered operator moments of $X$ and $P$ that appear naturally in high-order operator covariances. As examples, we consider two classes of quantum states. The first example is the class of (centralized) Gaussian states described by a covariance matrix of eigenvalues $\mu\lambda/2$ and $\mu/(2\lambda)$, where $\mu$ is related to the thermal mean photon number or temperature and $\lambda$ describes the squeezing strength. For simplicity, we set $\mu=\lambda$, which approximately models strongly-squeezed sources with accompanying excess noise associated to the anti-squeezed quadrature due to realistic experimental imperfections (refer for example to Ref. Müller et al. (2012)). The second example is the class of Fock states of $n$ photon numbers which are arguably the most non-Glauber-Sudarshan-representable states. Even for these states, there is in general no known explicit expressions for $\mathrm{sCRB}_{\textsc{bhom}}$ in the state parameters and numerical techniques are needed to calculate its values. The expressions for the sCRBs are listed in App. D.

Figures 3 and 4 present the findings for these states. As intuitively expected, the more direct BHOMOPT reconstruction of the moments is always (exponentially) better than estimating the Wigner function with BHOM. Even then, this improved strategy still often underperforms in comparison to HET and UHOM. For the Gaussian states, when $m=1$ or 3, marginalized sampling with BHOM and BHOMOPT give the worst tomographic performance. Localized sampling with the UHOM strategy generates the most accurate estimators per reconstruction datum, and delocalized sampling with HET is second best. When $m=2$ or 4, BHOMOPT beats HET respectively for $\mu\lesssim 1.262$ and $\mu\lesssim 1.017$, after which HET catches up in tomographic power, whereas UHOM ranks the top in the respective ranges $\mu\gtrsim 1.04$ and $\mu\gtrsim 1.004$. Likewise for the Fock states, both HET and UHOM, beat BHOM and BHOMOPT for all $n$ values and $m=1,3$. When $m=2$ or 4, BHOMOPT initially outperforms HET for the vacuum state (and also the $n=1$ state for $m=2$) and subsequently becomes inferior to HET. UHOM on the other hand is superior to all methods in tomographic power for all $n>0$.

That HET surpasses BHOMOPT for the $m=1$ case for any state is a consequence of the Heisenberg-Robertson-Schrödinger uncertainty relation Teo et al. (2017). The limiting case where the two methods give identical sCRBs is when the state is of minimum uncertainty. Yet, UHOM is able to overcome this limit owing to the binomial variances. For the vacuum ($\mu=1$ or $n=0$), both UHOM and BHOMOPT are almost identical in power ($\mathrm{sCRB}_{\textsc{uhom}}/\mathrm{sCRB}_{\textsc{bhomopt}}=33/32\approx 1.031$ and $9879/9856\approx 1.002$ for $m=2$ and $4$) within experimental error margins. This justifies the use of UHOM essentially for all these states.

Conclusion.—We have first proven that, for every used reconstruction datum, localized phase-space sampling with unbalanced homodyning always beats delocalized phase-space sampling with heterodyning in tomographic power measured by the scaled Cramér–Rao for any quantum state and general multivariate observable-parameter tomography. The reason is attributed to the binomial nature of unbalanced homodyne data, which enhances the resolution of Husimi-function reconstruction with fewer experimental components. We next demonstrated that for the Gaussian states and Fock states, localized sampling almost always beats marginalized sampling with balanced homodyning in moment tomography, except for the vacuum where both methods are practically equals. These findings shed light on the general performances of sampling methods in continuous-variable tomography.

Appendix A Derivations of Eqs. (2) and (4)

To arrive at the expression for $\mathrm{sCRB}_{\textsc{het}}$, we first note that since $N$ is fixed for HET, it is sufficient to use the standard formula $\mathbb{E}\!\left[n_{l}n_{l^{\prime}}\right]=Np_{l}\delta_{l,l^{\prime}}+N(N-1)p_{l}p_{l^{\prime}}$ for the binned multinomial HET data with $p_{l}\approx(\mathrm{d}\alpha)Q\!\left(\alpha_{l}\right)/\pi$. For UHOM, we would need the averages

[TABLE]

where here $p_{l}$ is instead equal to $Q\!\left(\alpha_{l}\right)$, $\sigma_{l}^{2}=p_{l}(1-p_{l})/N_{0}$, and $N_{0}$ is the total number of detection events for each $\alpha_{l}$.

The averages of data-ratios in Eqs. (5) and (6) can be straightforwardly derived by starting with this easy but crucial integral identity

[TABLE]

for any non-zero $A$. We can then rewrite the UHOM ratio averages as

[TABLE]

after an equivalent renormalization $f_{l}=n_{l}/N_{0}$ for notational simplicity. The central object to be evaluated is thus the characteristic function \mathbb{E}\!\left[\mathrm{e}^{\mbox{\footnotesize\mathrm{i}\sum_{l}\lambda_{l}f_{l}}}\right]. As the data for distinct $\alpha_{l}$ are statistically independent, we may proceed with the decomposition

[TABLE]

where central limit theorem gives

[TABLE]

for $N_{0}\gg 1$, and the corresponding integrals

[TABLE]

and

[TABLE]

where the parameters are now $\sigma_{l}^{2}=p_{l}(1-p_{l})$, $a=\sum_{l}\sigma^{2}_{l}/2$, $b=\sum_{l}p_{l}$, $c_{1}=\sigma^{2}_{l}\sigma^{2}_{l^{\prime}}$, $c_{2}=p_{l}\sigma^{2}_{l^{\prime}}+p_{l^{\prime}}\sigma^{2}_{l}$ and $c_{3}=p_{l}p_{l^{\prime}}+\delta_{l,l^{\prime}}\sigma^{2}_{l}$. The answers to these integrals involve the imaginary error function $\mathrm{erfi}\!\left(b/2\sqrt{a}\right)$, and in the limit of large sample size where $b\gg 2\sqrt{a}$, it turns out that the asymptotic expansion

[TABLE]

gives Eqs. (5) and (6) for sufficiently dense sampling $(0<b\gg 1)$. The integral expressions for $\mathrm{sCRB}_{\textsc{het}}$ and $\mathrm{sCRB}_{\textsc{uhom}}$ are the acquired limits of such a dense sampling.

The important technical statement $\min_{\widehat{\boldsymbol{q}}}\{\mathbb{E}\!\left[N(\widehat{\boldsymbol{q}}-\boldsymbol{q})^{2}\right]\}=\mathbb{E}\!\left[N\right]\min_{\widehat{\boldsymbol{q}}}\{\mathbb{E}\!\left[(\widehat{\boldsymbol{q}}-\boldsymbol{q})^{2}\right]\}$ is then easily proven with the additional statistical identity

[TABLE]

which can also be derived with

[TABLE]

of the same parameters defined above after a similar calculation.

To get $\mathrm{sCRB}_{\textsc{bhom}}$, we need the data-ratio averages

[TABLE]

that hold when BHOM sampling is sufficiently dense ($|b|\gg 1\implies N\gg 1$), with $\mathcal{N}=\sum_{l}\sum^{n_{\vartheta}}_{j=1}\sum^{n_{x}}_{k=1}\mathcal{R}^{-1}_{ljk}n_{jk}$, $p_{jk}=\mathrm{d}x_{k}\,p(x_{k},\vartheta_{j})$, $\widetilde{N}=\sum^{n_{x}}_{k=1}n_{jk}$, $\Sigma_{jkk^{\prime}}=(p_{jk}\delta_{k,k^{\prime}}-p_{jk}p_{jk^{\prime}})/\widetilde{N}$, $w_{jk}=\sum_{k^{\prime},l^{\prime}}\mathcal{R}^{-1}_{l^{\prime}jk^{\prime}}\Sigma_{jkk^{\prime}}$, $a=\sum_{j,k,l}\mathcal{R}^{-1}_{ljk}w_{jk}/(4n_{\vartheta})$ and $b=\sum_{j,k,l}\mathcal{R}^{-1}_{ljk}p_{jk}/(2n_{\vartheta})$. The averages in (17) and (18) can be verified with (8), which yields

[TABLE]

and

[TABLE]

after a renormalization $f_{jk}=n_{jk}/\widetilde{N}$ with the constant $\widetilde{N}$.

We again realize that the data for distinct $\vartheta_{j}$ are statistically independent, which means that

[TABLE]

decomposes into the independent characteristic functions. Because the binned data $\{n_{jk}\}_{k}$ for every $j$ follows a multinomial distribution defined by the BHOM quantum probabilities $\sum_{k}p_{jk}=1$, in the limit of large $\widetilde{N}$, the column $\boldsymbol{f}_{j}=(f_{jk})$ follows a Gaussian distribution of mean $\boldsymbol{p}_{j}$ and covariance matrix $\left[\mathrm{diag}(\boldsymbol{p}_{j})-\boldsymbol{p}_{j}\boldsymbol{p}_{j}\right]/\widetilde{N}$, so that

[TABLE]

according to the central limit theorem. After the differentiations, we have

[TABLE]

and

[TABLE]

where $c_{1}=w_{jk}w_{j^{\prime}k^{\prime}}$, $c_{2}=p_{jk}w_{j^{\prime}k^{\prime}}+p_{j^{\prime}k^{\prime}}w_{jk}$ and $c_{3}=p_{jk}p_{j^{\prime}k^{\prime}}+\delta_{j,j^{\prime}}\Sigma_{jkk^{\prime}}$. From the asymptotic formula in (14) that holds for $N\gg 1$, we arrive at the results in Eqs. (17) and (17) up to $O(1/N)$ as long as phase-space sampling is sufficiently dense or $|b|\gg 1$.

It remains to show that the sMSEs are indeed sCRBs. For this, we may invoke the central limit theorem for the respective sums of random variables, the key structures of the HET and UHOM estimators, to prove that $\mathrm{sMSE}\rightarrow\mathrm{sCRB}$ in the limit of dense sampling and/or accurate sampling. The especially important verification that the summands for the UHOM estimator are asymptotically independent is given in the next section. For HET, independence is clear.

Appendix B Asymptotic independence of UHOM random variables

Let $z_{l}\equiv n_{l}/N$. Then for any two UHOM random variables $z_{l}$ and $z_{l^{\prime}}$, the one-dimensional version of Kac’s theorem Applebaum et al. (2005) states that if we can show that the two-dimensional characteristic function

[TABLE]

satisfies this decomposition rule for all real $k_{l}$ and $k_{l^{\prime}}$, then $z_{l}$ and $z_{l^{\prime}}$ are statistically independent, and vice versa. This is equivalent to showing that $\mathbb{E}\!\left[z_{l}^{m}z_{l^{\prime}}^{m^{\prime}}\right]=\mathbb{E}\!\left[z_{l}^{m}\vphantom{z_{l^{\prime}}^{m^{\prime}}}\right]\mathbb{E}\!\left[z_{l^{\prime}}^{m^{\prime}}\right]$ for any $l$, $l^{\prime}\neq l$, $m$ and $m^{\prime}$.

For sufficiently large $N_{0}$, the first-order Taylor expansion in $n_{l}$ about $N_{0}p_{l}$ well approximates

[TABLE]

Then by recalling the simple statistical fact that the $n_{l}$s of distinct $l$ are of course independent binomial random variables,

[TABLE]

where

[TABLE]

Since from App. A, we know that $n_{l}/N_{0}$ is a Gaussian random variable of mean $\mu_{l}=p_{l}$ and variance $\sigma_{l}^{2}=p_{l}(1-p_{l})/N_{0}$ for $N_{0}\gg 1$, the $m$th moment

[TABLE]

is a simple function of the $m$th-degree Hermite polynomial $\mathrm{H}_{\,m}\!\left(\,\cdot\,\right)$. Using the simple relation $\mathrm{H}_{\,m+1}\!\left(y\right)=2y\,\mathrm{H}_{\,m}\!\left(y\right)-2m\,\mathrm{H}_{\,m-1}\!\left(y\right)$ that permits changes in the polynomial degree, we obtain the useful identity

[TABLE]

that can now be applied to Eq. (30) to get

[TABLE]

For sufficiently dense sampling the sum $\sum_{l}n_{l}$ should also be independent of $n_{l}$ since the sum is contributed by very many terms that are all distinct and therefore independent from $n_{l}$. Using the asymptotic relation $\mathrm{H}_{\,m}\!\left(y\right)\approx(2y)^{m}$ for $y\gg 1$, we find that

[TABLE]

where $f_{m,m^{\prime}}(x,y)=(m+1)(1-x)+(m^{\prime}+1)(1-y)$. The next order in the asymptotic expansion of $\mathrm{H}_{\,m}\!\left(y\right)$ gives a smaller correction to $\mathbb{E}\!\left[z_{l}^{m}z_{l^{\prime}}^{m^{\prime}}\right]$.

We can repeat the exercise and obtain the asymptotic formulas

[TABLE]

and

[TABLE]

the latter is obtained from the second-order Taylor expansion of $1/(\sum_{l}n_{l})^{m+m^{\prime}}$ in $n_{l}$ about $N_{0}p_{l}$,

[TABLE]

in which the first-order term vanishes since $\mathbb{E}\!\left[n_{l}\right]=N_{0}p_{l}$. These relations inform us that all statistical bias are asymptotic in nature. Combining all elements and keeping terms up to first order in $1/\mathbb{E}\!\left[N\right]$ gives

[TABLE]

Finally, invoking Kac’s theorem confirms asymptotic independence between $z_{l}$ and $z_{l^{\prime}}$, and thereafter for the whole set $\{z_{l}\}$ of these UHOM random variables.

Appendix C Realistic detections

It is a simple matter to show that the first main result remains unchanged for realistic detections. Suppose that all photodetectors now have the efficiency $0\leq\eta\leq 1$. Then standard characteristic-function treatment (see for instance Lai and Haus (1989)) allows us to find that the more realistic measured outcomes for HET are, instead of the usual coherent states $\left|{\alpha}\right>\left<{\alpha^{*}}\right|$, given by the full-rank statistical mixtures

[TABLE]

Alternatively, Born’s rule dictates that the realistic HET setup is equivalently the perfect HET setup with the quantum state $\rho$ transformed to $\rho^{\prime}$ by a corresponding Gaussian twirling operation. Then, the expression $\mathrm{sCRB}_{\textsc{het}}^{\prime}$ in (3) of the main text can be obtained by the simple replacement $\rho\rightarrow\rho^{\prime}$.

For UHOM, the results in Wallentowitz and Vogel (1996) show that the binomial probability for “no-click” detections is transformed to p(\alpha,\eta)=\left<\bm{:}\mathrm{e}^{\mbox{\footnotesize-\eta(a^{\dagger}-\alpha^{*})(a-\alpha)}}\bm{:}\right>. Furthermore, in going from the discretized sum to the continuous integral limit (review Sec. A of this SM), we note that $\sum_{l}p^{\prime}_{l}=\sum_{l}p(\alpha_{l},\eta)\rightarrow\pi/[(\mathrm{d}\alpha)\eta]$, which contributes the multiplicative factor $\eta$ in the expression for $\mathrm{sCRB}_{\textsc{uhom}}^{\prime}$.

Appendix D Moment tomography

It is easily verified that for the $m$th operator moment ($m\geq l$) that is Weyl ordered in the position $X$ and momentum $P$, its corresponding Wigner function is given by $x^{l}p^{m-l}$ in terms of the phase-space variables $x$ and $p$. There is then a simple one-to-one relation between $\boldsymbol{v}_{\textsc{p}}(\alpha)$ and $\boldsymbol{v}_{\textsc{w}}(\alpha)$ as a consequence of the Gauss transform. These are given by

[TABLE]

Then the evaluation of $\mathrm{sCRB}_{\textsc{bhomopt}}$, $\mathrm{sCRB}_{\textsc{het}}$ and $\mathrm{sCRB}_{\textsc{uhom}}$ amounts to the evaluation of all integrals involving $\boldsymbol{v}_{\textsc{p}}(\alpha)$ and $\boldsymbol{v}_{\textsc{w}}(\alpha)$ using the identities in (39).

For $\mathrm{sCRB}_{\textsc{het}}$ and $\mathrm{sCRB}_{\textsc{uhom}}$, this can be easily accomplished with the help of characteristic functions \displaystyle\chi_{1}=\overline{\mathrm{e}^{\mbox{\footnotesizeg^{}\alpha+g\alpha^{}}}}\,\,\left[g=(u+\mathrm{i}v)/\sqrt{2}\right] and \displaystyle\chi_{2}=\overline{\overline{\mathrm{e}^{\mbox{\footnotesizeg^{}\alpha+g\alpha^{}}}}}, where the single and double overlines respectively denote the phase-space integrals with respect to $Q\!\left(\alpha\right)$ and (unnormalized) $Q\!\left(\alpha\right)\left[1-Q\!\left(\alpha\right)\right]$. Then

[TABLE]

supply the required quantities.

The centralized Gaussian states of the Husimi-function covariance matrix $\boldsymbol{G}_{\textsc{het}}$ have the characteristic functions

[TABLE]

where $\boldsymbol{g}\,\widehat{=}\,{(-g\,\,\,g^{*})}^{\mathrm{\,\textsc{t}}}$, $\boldsymbol{M}=\boldsymbol{H}^{\dagger}\boldsymbol{G}_{\textsc{het}}^{-1}\,\boldsymbol{H}$, and $\boldsymbol{H}\,\widehat{=}\dfrac{1}{\sqrt{2}}\begin{pmatrix}1&1\\ -\mathrm{i}&\mathrm{i}\end{pmatrix}$. Those of the Fock states read

[TABLE]

For the setting $\mu=\lambda$, the $\mathrm{sCRB}$ expressions are catalogued as follows:

[TABLE]

The corresponding expressions for the Fock states are given by

[TABLE]

The $\mathrm{sCRB}_{\textsc{bhomopt}}$ expressions for the improved strategy of HOM can be found in an analogous way by looking at the operator moments and calculating the Fisher information matrix Teo et al. (2017). Further analysis shall be reported elsewhere but for now, we simply supply all the final analytical results that are obtainable from the theory. These are

[TABLE]

for the Gaussian states, and

[TABLE]

for the Fock states.