Quantitative tomography for continuous variable quantum systems

TL;DR

This paper introduces an efficient continuous variable quantum tomography method using Padua points for sampling, enabling accurate reconstruction of quantum states with fewer measurements and providing direct density matrix estimation with error bounds.

Contribution

The paper presents a novel tomography scheme utilizing Padua points and Lagrange interpolation, significantly reducing measurement requirements and allowing direct density matrix estimation with error quantification.

Findings

Reconstruction error decreases exponentially with the number of Padua points.

The method achieves quasi-linear runtime in the number of sampling points.

Density matrix elements can be estimated with linear error propagation and analytical bounds.

Abstract

We present a continuous variable tomography scheme that reconstructs the Husimi Q-function (Wigner function) by Lagrange interpolation, using measurements of the Q-function (Wigner function) at the Padua points, the optimal sampling points for two dimensional reconstruction. Our approach drastically reduces the number of measurements required compared to using equidistant points on a regular grid, although reanalysis of such experiments is possible. The reconstruction algorithm produces a reconstructed function with exponentially decreasing error and quasi-linear runtime in the number of Padua points. Moreover, using the interpolating polynomial of the Q-function, we present a technique to directly estimate the density matrix elements of the continuous variable state, with only linear propagation of input measurement error. Furthermore, we derive a state-independent analytical bound on this error, such that our estimate of the density matrix is accompanied by a measure of its uncertainty.

Figures

D of the main text. Despite Chebfun2’s non-polynomial preprocessing of the sampled data, which aims to avoid the Runge phenomenon, the reconstruction is considerably worse than when the Padua Lagrange polynomial interpolant is used, as in the main text. In particular, the relative error is no longer strictly decreasing, and appears to saturate for some density matrix elements. This specific example highlights the benefit of performing Lagrange interpolation at the Padua points in comparison to state-of-the-art polynomial approximation on an equidistant grid.

Equations

$$W(\alpha )=\frac{1}{\pi }\left[\hat{\Pi }\hat{D}(-\alpha )\rho \hat{D}(\alpha )\right],$$

$$Q(\alpha )=\frac{1}{\pi }⟨\alpha |\rho |\alpha ⟩=\frac{1}{\pi }\mathrm{Tr}\left[|0⟩⟨0|\hat{D}(-\alpha )\rho \hat{D}(\alpha )\right],$$

$$|\psi ⟩=\frac{1}{\sqrt{2}}\left(|\bar{0}⟩+i|\bar{1}⟩\right)=\frac{1}{\sqrt{2}}\left(\frac{|0⟩+|4⟩}{\sqrt{2}}+i|2⟩\right),$$

$$\forall k{\mathcal{L}}_n[f]({\alpha }_k)=v_k$$

$$||f-{\mathcal{L}}_n(f)|{|}_{\infty }\le (1+{\Lambda }_n)||f-f_n^{\ast }|{|}_{\infty }.$$

$${\rho }_{jk}=\mathrm{Tr}\left[|k⟩⟨j|\rho \right]=\pi \int P_{jk}(\alpha )Q_{\rho }(\alpha )d\alpha ,$$

$$P_{jk}(\alpha )=\frac{C_{k,j}}{2\pi r}{e}^{r^2-i(j-k)\theta }\left[{\left(-\frac{\partial }{\partial r}\right)}^{j+k}\delta (r)\right],$$

$$Q_{\rho }(r,\theta )=\sum _{\begin{array}{c}0\le m\le n\\ -n\le p\le n\end{array}}{c}_{m,p}{r}^m{e}^{ip\theta },$$

$${\rho }_{jk}=C_{k,j}\sum _{q=0}^{j+k}q!d_q^{j+k}{c}_{q,(j-k)},$$

$$\Delta {\rho }_{jk}\left[N\right]=\frac{\left|{\rho }_{jk}^{\mathrm{ideal}}-{\rho }_{jk}[N,0]\right|}{{\rho }_{jk}^{\mathrm{ideal}}},$$

$${\sigma }_{jk}\left[N,ϵ\right]\propto ϵ,$$

$${\rho }_{jk}={\rho }_{jk}[N,ϵ]\pm {\Delta }_{jk}[N]{\rho }_{jk}^{\mathrm{ideal}}\pm {\sigma }_{jk}[N,ϵ]$$

$${\rho }_{jk}=\mathrm{Tr}\left[|k⟩⟨j|\rho \right]=\pi \int P_{jk}(\alpha )Q_{\rho }(\alpha )d\alpha ,$$

$$P_{jk}(\alpha )=\frac{C_{k,j}}{2\pi r}{e}^{r^2-i(j-k)\theta }\left[{\left(-\frac{\partial }{\partial r}\right)}^{j+k}\delta (r)\right],$$

$${\rho }_{jk}=C_{k,j}{\int }_{-\pi }^{\pi }d\theta e^{-i(j-k)\theta }{\left[\frac{{\partial }^{j+k}}{\partial r^{j+k}}\left(e^{r^2}{Q}_{\rho }(r,\theta )\right)\right]}_{r=0}=C_{k,j}{\int }_{-\pi }^{\pi }d\theta e^{-i(j-k)\theta }\sum _{q=0}^{j+k}{d}_q^{j+k}{\left[\frac{{\partial }^q{Q}_{\rho }(r,\theta )}{\partial r^q}\right]}_{r=0}$$

$$\frac{{\partial }^{j+k}}{\partial r^{j+k}}\left(e^{r^2}{Q}_{\rho }(r,\theta )\right)=\sum _{q=0}^{j+k}\left(\genfrac{}{}{0ex}{}{j+k}{q}\right)\frac{{\partial }^{j+k-q}}{\partial r^{j+k-q}}\left(e^{r^2}\right)\frac{{\partial }^q{Q}_{\rho }}{\partial r^q}=e^{-r^2}\sum _{q=0}^{j+k}\left(\genfrac{}{}{0ex}{}{j+k}{q}\right)(-1{)}^{j+k-q}{H}_{j+k-q}(-r)\frac{{\partial }^q{Q}_{\rho }}{\partial r^q}$$

$$\displaystyle d_{q}^{j+k}=(-1)^{j+k-q}{{j+k}\choose{q}}H_{j+k-q}(0)=\Bigg{\{}\begin{array}[]{ll}(-1)^{\frac{3}{2}\left(j+k-q\right)}\frac{(j+k)!}{q!\left(\frac{j+k-q}{2}\right)!}&j+k-q{\rm\ even}\\ 0&j+k-q{\rm\ odd}\end{array}$$

$$Q_{\rho }(r,\theta )=\sum _{\begin{array}{c}0\le m\le n\\ -n\le p\le n\end{array}}{c}_{m,p}{r}^m{e}^{ip\theta }$$

$${\rho }_{jk}=C_{k,j}\sum _{q=0}^{j+k}q!d_q^{j+k}{c}_{q,(j-k)}.$$

$${\mathcal{N}}_{2n}(x,y)=\sum _{j=0}^{n-1}{a}_j\prod _{m=0}^{j-1}(x-x_m)(y-y_m),$$

$$f(x_t,y_t)-{\mathcal{N}}_{2n}(x_t,y_t)={\mathcal{N}}_{2(n+1)}(x_t,y_t)-{\mathcal{N}}_{2n}(x_t,y_t)=a_n\prod _{m=0}^n(x_t-x_m)(y_t-y_m),$$

$${\left|\left|\prod _{m=0}^{j-1}(x_t-x_m)(y_t-y_m)\right|\right|}_{\infty }=\frac{1}{2^{2n}},$$

$${\left|\left|f-{\mathcal{N}}_{2n}\right|\right|}_{\infty }\le \frac{{\left|\left|a_n(x_t,y_t)\right|\right|}_{\infty }}{2^{2n}}.$$

$${\left|\left|f-{\mathcal{L}}_n^{\mathrm{Pad}}\right|\right|}_{\infty }$$

$$\delta [c_{m,p}]\le 4{\left|\left|f-f_n^{\ast }\right|\right|}_{\infty }\le \frac{4{\left|\left|a_{\frac{n}{2}}(x_t,y_t)\right|\right|}_{\infty }}{2^n}.$$

$${\rho }_{jk}^{\mathrm{ideal}}{\Delta }_{jk}$$

$$\le \frac{\sqrt{k!j!}}{2}\sum _{\begin{array}{c}q=0\\ \mathrm{even}\end{array}}^{j+k}\frac{\delta [c_{m,p}]}{\left(\frac{j+k-q}{2}\right)!}\le \frac{2\sqrt{k!j!}{\left|\left|a_{\frac{n}{2}}(x_t,y_t)\right|\right|}_{\infty }}{2^n}\sum _{\begin{array}{c}q=0\\ \mathrm{even}\end{array}}^{j+k}\frac{1}{\left(\frac{j+k-q}{2}\right)!}$$

$$|\psi ⟩=\frac{1}{\sqrt{2}}\left(\frac{|0⟩+|4⟩}{\sqrt{2}}+i|2⟩\right).$$

$${\sigma }_{jk}\left[N,ϵ\right]=A_{jk}{ϵ}^p,$$

$${\sigma }_{jk}[N,ϵ]\le Kϵ\frac{\sqrt{k!j!}}{2}\sum _{\begin{array}{c}q=0\\ \mathrm{even}\end{array}}^{j+k}\frac{1}{\left(\frac{j+k-q}{2}\right)!},$$

$$A_{jk}\le K\frac{\sqrt{k!j!}}{2}\sum _{\begin{array}{c}q=0\\ \mathrm{even}\end{array}}^{j+k}\frac{1}{\left(\frac{j+k-q}{2}\right)!}.$$

Full text

Introduction—In modern implementations of quantum information protocols, quantum state tomography Hradil (1997); Altepeter et al. (2004) plays a key role. It allows for full characterization of unknown quantum states, as well as verification of prepared resource states. Continuous variable (CV) quantum systems have a wide range of applications in all areas of quantum information, ranging from quantum communication to quantum computing Braunstein and van Loock (2005). The need for fast, efficient quantum state tomography of CV systems is heightened by the continuing development of non-classical radiation sources for quantum information Ourjoumtsev et al. (2006); Kurochkin et al. (2014); Huang et al. (2015); Reimer et al. (2016); Eichler et al. (2011a); Zhong et al. (2013); Eichler et al. (2014); Toyli et al. (2016), and by recent developments in CV encodings of logical qubits Mirrahimi et al. (2014); Ofek et al. (2016); Michael et al. (2016), as well as in quantum simulation with CV systems Flurin et al. .

Quantum state tomography of CV systems typically consists of measurement of a quasi-probability distribution, such as the Wigner or Husimi Q-function Braunstein and van Loock (2005), from which the density matrix can be reconstructed Lvovsky and Raymer (2009), and this has been extensively demonstrated in experiment Kirchmair et al. (2013); Wang et al. (2009); Vlastakis et al. (2013); Bertet et al. (2002); Deleglise et al. (2008); Hofheinz et al. (2009); Wang et al. (2016). Unfortunately, full tomography of the Wigner or Q-function is inefficient, as a large number of measurements is required to sample all of phase space. To reduce the number of sampled phase space points, more advanced tomographic schemes involve displacing the state in phase space, and multiple measurements (ideally of the full photon number distribution) at each displacement Shen et al. (2016); Leibfried et al. (1996); Hofheinz et al. (2009). For these schemes there is a trade-off between the number of measurement points in phase space, and the number of operator expectation values measured at each point.

In this paper we propose an efficient method to reconstruct the full Wigner or Q-function via Lagrange interpolation, using only a small number of measured phase space points. Our method can be applied on any grid of phase space points, thus allowing for reanalysis of previous experiments. However, the ideal points to be measured are known as the Padua points Caliari et al. (2005), which are the optimal sampling points for two dimensional (2D) interpolation Bos et al. (2006); Caliari et al. (2008). Our scheme requires only a single measurement at each sample point: the Wigner or Q-function value. In addition, we show how individual density matrix elements (including off-diagonal elements) can be estimated from the interpolation reconstruction of the Q-function, such that quantum state tomography of a CV system can be performed directly, without statistical inference.

Background—It is well known that the Wigner and Q-function can be measured through a combination of coherent displacements and parity (Wigner) or ideal vacuum (Husimi-Q) measurements Lvovsky and Raymer (2009), i.e.

[TABLE]

where $W(\alpha)$ and $Q(\alpha)$ are the Wigner and Q-function, $\hat{D}(\beta)=\exp\left(\beta\hat{a}^{\dagger}-\beta^{*}\hat{a}\right)$ is the usual displacement operator, and $\hat{\Pi}=\left(-1\right)^{\hat{a}^{\dagger}\hat{a}}$ is the parity operator. Parity measurements can be implemented via interaction with a qubit Leibfried et al. (1996); Wang et al. (2009); Vlastakis et al. (2013); Bertet et al. (2002); Deleglise et al. (2008); Hofheinz et al. (2009); Wang et al. (2016), as can ideal vacuum detection Kirchmair et al. (2013); Flurin et al. , though this can also be done via photon subtraction measurements Chen et al. (2011); Govia et al. (2012); Oi et al. (2013), or heterodyne detection Kim (1997); Eichler et al. (2011b, 2012), even of an itinerant state.

Experimentally, the first step of our method consists of using one of the above techniques to measure the Wigner or Q-function at the Padua points. Using the results of these measurements, one can then reconstruct the full quasi-probability distribution using Lagrange interpolation, and, if the Q-function was measured, directly estimate the density matrix elements of the measured state.

Inspired by the binomial codes of Ref. Michael et al. (2016), we use the test state

[TABLE]

to illustrate our tomography method, where $|\bar{0}\rangle$ and $|\bar{1}\rangle$ are the logical states for the lowest order binomial code of Ref. Michael et al. (2016), with $|n\rangle$ the usual $n$-photon Fock state. This is a reasonable test state as it contain several non-zero density matrix elements of relatively large photon number, and has both real and imaginary off-diagonal elements, such that it tests our reconstruction algorithm under all possible conditions. We will now explain in detail both aspects of our tomography method, quasi-probability reconstruction and density matrix element estimation, and show the results of our method applied on the test state, including the effects of finite measurement error.

Reconstruction of quasi-probability distributions—The most naïve way to obtain a quasi-probability distribution from measured data is to simply plot the measured values, and increase the number of sampled points until an accurate reconstruction is obtained. In Fig. 1A we show the exact Q-function for our test state, and in Fig. 1B we show the result of plotting 256 sampled points (from an equidistant $16\times 16$ grid) of the Q-function. As can be seen, Fig. 1B is a very poor reconstruction of the Q-function, and increasing the number of sampled points to improve the reconstruction requires a considerable increase in experimental resources.

Instead, we propose that the Q-function (or Wigner function) should be reconstructed by Lagrange interpolation. The Q-function in particular is ideal because of its smoothness, as it is known to be an analytic function Appleby (1999). Generally, interpolation is the problem of reconstructing a function $f$ from the knowledge of its values $\{v_{k}\}_{k=1}^{N}$ at $N$ sampling points $\{\alpha_{k}\}_{k=1}^{N}$. Using a finite number $N$ of such sampling points allows one to reconstruct the Lagrange polynomial $\mathcal{L}_{n}[f]$ which obeys

[TABLE]

where $n$ is the order of the polynomial, with $N={{n+2}\choose{2}}$ for a 2D function.

However, in general, any interpolated function will differ from the polynomial $f^{*}_{n}$ (of order $n$) closest to $f$ in uniform norm $||f||_{\infty}=\sup_{x}|f(x)|$. The choice of sampling points can strongly affect the quality of the reconstruction Berrut and Trefethen (2004), and the quality of the sampling points (more generally of the interpolation scheme $\mathcal{L}_{n}$) is captured by its operator norm, known as the Lebesgue constant $\Lambda_{n}$, such that

[TABLE]

It is well known Cheney and Light (2009) that in 1D the optimal sampling points are the Chebyshev nodes, whose Lebesgue constant scales logarithmically ($\Lambda_{n}\sim\log(n)$). On the contrary, equidistant points have a Lebesgue constant exponentially large in $n$, leading to artifact oscillations known as Runge’s phenomenon.

The Padua points are the 2D equivalent of the Chebyshev points Caliari et al. (2005), and are currently the best known choice of sampling points for 2D interpolation, since their Lebesgue constant scales as $\log^{2}n$, which is expected to be optimal Bos et al. (2006). Moreover, the Padua points have a simple generating curve, such that their locations can be efficiently calculated, and allow for efficient computation of the Lagrange polynomial in time $\mathcal{O}(n^{2}\log(n))$, using a stable numerical scheme whose open-source Matlab implementation is available 111http://profs.sci.univr.it/ caliari/software.htm. Note that this running time is linear in the number of sampling points $N$ (up to logarithmic factors). This Padua interpolation technique is now integrated into standard 2D interpolation packages such as Chebfun2 Townsend and Trefethen (2013).

In Fig. 1C, we present the reconstructed Q-function for our test state using 231 Padua points (indicated by the white dots), which is an accurate reconstruction of the exact state. It even captures the features of the Q-function close to the origin, despite few sampling points in that region, which we believe is due to the analyticity of the Q-function: points away from the origin constrain the higher-order derivatives of the function, leading to a highly accurate reconstruction.

To compare to the Padua points, in Fig. 1D we consider interpolation on a square grid of 256 points using Chebfun2 Townsend and Trefethen (2013), which performs non-polynomial pre-processing on the sampled values before producing an interpolating polynomial, in order to reduce the error inherent to equidistant interpolation 222The function chebfun2 from Townsend and Trefethen (2013) for equidistant points first generates a Floater-Hormann interpolant based on rational functions Floater and Hormann (2007), and then produces a polynomial approximation of this rational interpolant. This is done to minimize the error inherent to interpolation on an equidistant grid (due to the Runge phenomenon). We compare normal Padua Lagrange interpolation to this hybrid scheme to highlight the power of the Padua points.. Fig. 1C and Fig. 1D are qualitatively similar, highlighting that any interpolation is a powerful tool. However, there are practical reasons why the Padua points are the better sampling set, as we will now discuss.

Although the theoretical Q-function we aim to reconstruct has only significant values on a rectangle whose main axis is diagonal, we chose to use interpolation points distributed on a square in order to avoid biasing the reconstruction method towards the test state. An alternative scenario is where one is reasonably certain of the state beforehand and require only verification. This can be done by thresholding, where we create a large grid of Padua points, but then only measure at points where we expect (based on our prior knowledge) the quasi-probability to be above some threshold. To illustrate this, we consider a grid of 231 Padua points, and manually set to exactly zero all sampled Q-function values whose absolute value is below $10^{-2}$. The resulting Fig. 1E is very similar to Fig. 1C, but only requires 65 Padua point measurements, which is a significant reduction of experimental effort.

Performing a similar thresholding procedure on the equidistant grid (keeping only the points with the 65 largest values) leads to Fig. 1F, which is clearly a bad reconstruction of the Q-function. The thresholded Padua points cover more of phase space than the thresholded equidistant points, allowing them to more accurately constrain the higher order derivatives of the function, which explains their advantage in verification by thresholding.

The main advantage of the Padua points is that their interpolated function has a quasi-exponentially decreasing error as the number of points is increased, which is not generally true for other point sets (including an equidistant grid). This greatly benefits the method to estimate the density matrix elements that we now introduce, as it leads to very favorable error scaling, and an analytic bound on the remaining error.

Direct estimation of density matrix elements—While quasi-probability distributions give good qualitative descriptions of CV states, for quantitative information one needs the density matrix, which can be calculated using either a linear inversion method, or by statistical inference Lvovsky and Raymer (2009). Linear inversion methods are prone to error accumulation, and thus statistical inference methods, such as maximum-likelihood reconstruction, are more commonly employed Lvovsky and Raymer (2009). However, with statistical inference it is difficult to assign a measure of confidence or error to the calculated density matrix.

Our method to directly estimate the density matrix elements uses the Lagrange interpolation reconstruction of the Q-function (calculated from measured data), without the need for statistical inference. It does not aim to return a trace one positive semidefinite matrix, but rather estimates the density matrix elements, with calculated bounds on the error in the estimate, based on both input measurement error and the reconstruction error of our algorithm. Our method is related to pattern function reconstruction using the Wigner function Lvovsky and Raymer (2009), but by using the Q-function, it avoids many of the difficulties inherent to pattern functions.

To calculate the density matrix elements of a state $\rho$ using our method, we first write them as

[TABLE]

where $Q_{\rho}(\alpha)$ is the Q-function of the measured state, and $P_{jk}(\alpha)$ is the Glauber-Sudarshan P-function representation of the operator $|k\rangle\langle j|$, which in the Fock basis is given by Sudarshan (1963)

[TABLE]

where $\alpha=re^{i\theta}$ in polar coordinates, and $C_{k,j}=\sqrt{k!\,j!}/\left(2\left(k+j\right)!\right)$, is a combinatorial factor. Eqs. (6)-(7) transform the estimation of the density matrix elements into estimating the derivatives of the sampled Q-function. Inspired from a well-known procedure to optimally estimate finite-differences on an arbitrary grid Fornberg (1981, 1988), we now show how to relate the Lagrange interpolating polynomial to these derivatives.

The interpolation polynomial can be written in polar coordinates

[TABLE]

where $c_{m,p}$ are the Q-function coefficients calculated from measured data, and as before, $n$ is the polynomial order of the reconstruction. Using the reconstructed Q-function of Eq. (8), and the P-function of Eq. (7) we obtain an expression for the calculated density matrix elements

[TABLE]

where $d_{q}^{j+k}$ are state-independent constants. Evaluating this expression requires only simple algebra on the coefficients of the reconstructed Q-function, $c_{q,(j-k)}$, which are efficiently calculated from the experimental measurements. We stress that the other factors in Eq. (9), including $d_{q}^{j+k}$ and $C_{k,j}$, are state independent constants, and can be efficiently computed once a priori.

To illustrate our method, we estimate the density matrix of our test state $|\psi\rangle$ from Q-function measurements at the Padua points. To quantify the error intrinsic to Padua reconstruction, we define the relative reconstruction error

[TABLE]

where $\rho\left[N,0\right]$ is the estimated density matrix element obtained from $N$ Q-function measurements. The relative reconstruction error is plotted in Fig. 2 for all density matrix elements of our test state $|\psi\rangle$. As can be seen, the estimated state becomes a better approximation to the actual state as $N$ increases, with the general trend a quasi-exponential decrease in error for large $N$. While our method works for any polynomial approximation of the Q-function, this quasi-exponential scaling is a result of the fact that the Padua points are the optimal reconstruction points in 2D Sup . In Sec. IV of the EPAPS Sup , we show the relative error versus $N$ for an equidistant grid and Chebfun2 interpolation; the error is appreciably worse than for Padua interpolation.

A higher degree polynomial reconstruction (more Padua points) is required to accurately estimate $\rho_{jk}$ for larger $j,k$, as is seen in Fig. 2, where the components with larger $j,k$ have larger relative reconstruction error. Thus, our method is most effective for finite superpositions of Fock states, such as binomial code states Michael et al. (2016), or states with low average photon number. Alternatively, the center of high quasi-probability can be found using efficient adaptive schemes Stenberg et al. (2015), and the state displaced to the origin. Our method can then efficiently reconstruct this new displaced state, from which the original state can be calculated.

To accurately reflect imperfect experimental measurements we introduce noise by adding to the sampled Q-function values Gaussian random noise with zero mean and standard deviation $\epsilon$. As this noise is stochastic, each reconstructed element $\rho_{jk}[N,\epsilon]$ will have a mean value $\bar{\rho}_{jk}[N,\epsilon]=\rho_{jk}[N,0]$ and a standard deviation $\sigma_{jk}\left[N,\epsilon\right]$, which are functions of both the number of Padua points $N$, and the noise level $\epsilon$. To examine how input measurement error propagates through our reconstruction algorithm we can use the standard deviation, which we calculate for a range of values of $N$, and input measurement error ranging from $\epsilon=10^{-5}$ to $\epsilon=10^{-1}$ (see Sup for further details). We find linear error scaling

[TABLE]

which could also be surmised from the fact that Eq. (9) is linear in the Q-function coefficients $c_{q,(j-k)}$ reconstructed by interpolation which we expect to propagate input measurement errors linearly. As such, our reconstruction algorithm does not amplify the input measurement error by more than a constant factor.

Input noise enters Eq. (9) of our reconstruction algorithm as error in the coefficients $c_{q,(j-k)}$. As all other terms in Eq. (9) are state-independent, propagation of this error through our algorithm will be state-independent. Therefore, the calculated scaling law of $\sigma_{jk}\left[N,\epsilon\right]$ for our test state is indicative of the scaling law for all states.

Combining the previous results, our reconstruction scheme gives an estimate $\rho_{jk}[N,\epsilon]$ for the element $\rho_{jk}$ such that

[TABLE]

where the first term in Eq. (12) accounts for the systematic bias in the estimate due to the absolute reconstruction error, and the second term describes the measurement error due to input noise (assuming one standard deviation error). As this expression shows, our scheme not only calculates the estimated value $\rho_{jk}[N,\epsilon]$, but also estimates the error, given by a linear sum of the absolute reconstruction error and the measurement error.

The input noise, $\epsilon$, can be measured experimentally by benchmarking measurements of known states, and this gives the measurement error, $\sigma_{jk}[N,\epsilon]$, up to a constant factor that can be estimated from Eq. (9) Sup . The absolute reconstruction error will be state dependent, and can most accurately be calculated by numerical simulation. Alternatively, one can iteratively increase the number of Padua points until $\rho_{jk}[N,\epsilon]$ converges to within $\sigma_{jk}[N,\epsilon]$, which indicates that the measurement error is now dominant and the reconstruction error can be ignored.

Conclusion—In this work, we have introduced a continuous variable tomography scheme for accurate and efficient reconstruction of quasi-probability distributions using Lagrange interpolation. The only experimental requirement is the ability to measure the quasi-probability distribution, and if this is done at the Padua points, the optimal sampling points in 2D, our scheme drastically reduces the number of measurements required. Reanalysis of experiments for measurements performed on other point sets is also possible, but forfeits the error bounds guaranteed by using Padua points.

Using the reconstructed Q-function, we have further shown how to estimate the system’s density matrix elements, and bounded the error in this estimation. Remarkably, the intrinsic reconstruction error (due to interpolation) decreases exponentially with the number of Padua points, while input measurement error is not amplified by more than a constant factor.

Our scheme will see immediate application in quantum communication and computing protocols with optical or microwave fields, including cavity and circuit QED setups, where it offers significant improvement over the current state of the art. Additionally, our scheme may have application in the characterization of ultra-fast electromagnetic pulses, such as in frequency-resolved optical gating, or other autocorrelation techniques.

I Direct Estimation of Density Matrix Elements

We start from the expression for the density matrix elements $\rho_{jk}$ used in the main text

[TABLE]

where $Q_{\rho}(\alpha)$ is the Q-function of the measured state, and the expression for the Glauber-Sudarshan P-function of $|k\rangle\langle j|$ in the Fock basis Sudarshan (1963)

[TABLE]

where $\alpha=re^{i\theta}$ in polar coordinates, and $C_{k,j}=\sqrt{k!\,j!}/\left(2\left(k+j\right)!\right)$, is a combinatorial factor. Substituting the expression for $P_{jk}(\alpha)$ of Eq. (14) into Eq. (13), we obtain an expression for the density matrix elements in terms of derivatives of the Q-function evaluated at the origin

[TABLE]

where the constant factors $d_{q}^{j+k}$ can be evaluated using the product rule as now follows.

We start by applying the general Leibniz rule

[TABLE]

where $H_{j+k-q}(-r)$ is the Hermite polynomial of order $j+k-q$. Evaluating Eq. (16) at $r=0$ and comparing to Eq. (15), we obtain

[TABLE]

where in the last line we have used the series expansion for the Hermite polynomials, and the fact that odd order Hermite polynomials are zero at $r=0$.

Finally, using the expression for the reconstructed Q-function in polar coordinates

[TABLE]

where $n$ is the order of the polynomial reconstruction, we see that the only components that will be survive the integral over $\theta$ in Eq. (15) will have $p=j-k$. With this in mind, substituting the reconstructed Q-function into Eq. (15) gives the expression for the density matrix elements found in the main text

[TABLE]

II Relative Reconstruction Error Scaling

In analogy to the 1D case Stewart (1996), we consider a “Newton’s” form of the 2D polynomial reconstruction (of order $2n$) of the Q-function, described by the function

[TABLE]

where $\left\{(x_{m},y_{m})\right\}^{n}_{m=0}$ is the set of interpolation points, and the coefficients $a_{j}$ can be determined iteratively as more points are added to the interpolation set. Adding one additional arbitrary point $(x_{t},y_{t})$ to this set would produce the interpolation polynomial of order $2(n+1)$, and by definition $\mathcal{N}_{2(n+1)}(x_{t},y_{t})=f(x_{t},y_{t})$. Therefore, if we wish to calculate the difference between the function and our order $2n$ reconstruction at the arbitrary point $(x_{t},y_{t})$, we see that it is given by

[TABLE]

where $a_{n}$ is a function of the whole set $\left\{(x_{m},y_{m})\right\}^{n}_{m=0}$, and $(x_{t},y_{t})$. The interpolation error $\left|\left|f-\mathcal{N}_{2n}\right|\right|_{\infty}$ is just the maximum of the above expression over all points $(x_{t},y_{t})$. Choosing the sets $\left\{x_{m}\right\}^{n}_{m=0}$ and $\left\{y_{m}\right\}^{n}_{m=0}$ to each be the Chebyshev points, we have that

[TABLE]

and as such we can bound the interpolation error from above by

[TABLE]

The best polynomial interpolation of order $2n$, $f^{*}_{2n}$, will have smaller error than the above expression, and using the Lebesgue constant for the Padua interpolation, we can say that for the Padua interpolation of order $n$ the error will be bounded as

[TABLE]

In a 1D reconstruction of Newton’s form, the analogue of the function $a_{n}(x_{t},y_{t})$ is related to the $(n+1)$ order derivative of the function $f$, which for many functions does not grow exponentially with $n$. If we assume that $a_{n}(x_{t},y_{t})$ does not grow exponentially with $n$ in the 2D case, then we obtain a quasi-exponential decrease in the Padua interpolation error.

From Ref. Caliari et al. (2008), the error in the reconstructed Q-function coefficients, $c_{m,p}$, is bounded by

[TABLE]

Using this and Eqs. (21) and (19), we can derive a bound on the reconstruction error in the calculated density matrix element, $\Delta_{jk}$, in terms of the error in the Q-function coefficients $c_{m,p}$

[TABLE]

where, as before, $n$ is the order of the Padua reconstruction. By the subscript “even” in the summations we mean that the sum is over only values of $q$ such that $j+k-q$ is even. From this expression, we see that if we assume as before that $a_{n}(x_{t},y_{t})$ does not grow exponentially with $n$, then we recover the quasi-exponential decay of the reconstruction error seen in the simulated results of Fig. 2.

III Scaling of the Standard Deviation of the Measurement Error

We calculate the standard deviation $\sigma_{jk}[N,\epsilon]$ for values of $N=(n+1)(n+2)/2$ with $n\in[11,35]$, and $\epsilon\in\{10^{-5},10^{-4},10^{-3},10^{-2},10^{-1}\}$. For all density matrix elements of the test state

[TABLE]

We find that the standard deviation of the reconstructed state scales as

[TABLE]

with $p=1\pm 1.6\times 10^{-3}$, for all $j,k$ and $N$. As an example, we plot $\sigma_{jk}\left[N,\epsilon\right]$ for fixed $N=253$ as a function of the noise level $\epsilon$ in Fig. 3, and as can be seen, $\sigma_{jk}\left[N,\epsilon\right]$ increases linearly as $\epsilon$ increases.

Given the form of Eq. (21), a sum of terms linear in $c_{q,(j-k)}$, we would expect linear scaling of $\sigma_{jk}$ with respect to the error in the reconstructed Q-function coefficients $c_{q,(j-k)}$. These coefficients are reconstructed from noisy Q-function measurements, with input measurement error $\epsilon$. Thus, the fact that the standard deviation of $\rho_{jk}$ scales linearly with $\epsilon$ implies that the error in $c_{q,(j-k)}$ must also scale linearly with $\epsilon$. As such, the Q-function reconstruction algorithm does not amplify input measurement error by more than a constant factor.

More concretely, if we assume that the error in the reconstructed Q-function coefficients $c_{q,(j-k)}$ is linear in the error in the Q-function measurements, i.e. it is $K\epsilon$, then we can bound the standard deviation in a similar way as to what was done in section II for the relative error. Such a procedure results in the expression

[TABLE]

which confirms the linear scaling of the standard deviation with $\epsilon$, and bounds the constant factor $A_{jk}$ by

[TABLE]

It is this bound that allows us to calculate error bars on our density matrix elements in a state independent way, with knowledge of only the measurement error $\epsilon$, and the constant $K$ coming from the Q-function interpolation. In general $K$ is dependent on the number of Padua points, and is difficult to calculate analytically. Heuristically we find that it is a function of both the indices $j,k$ and the number of Padua points $N$, and generally increases as these increase. For the largest $j,k,N$ tested we find $K\simeq 40$, though in general $K$ is not this large, and more often $K\simeq 1$, such as is the case for the results shown in Fig. 3.

Note that as we are introducing measurement error artificially, our results may suffer from finite sampling error, as the standard deviation of the reconstructed state is itself a stochastic variable. However, we sample from the noisy distribution sufficient times (on the order of $10^{4}$) such that the variance of the distribution for the standard deviation of the reconstructed state is several orders of magnitude smaller than the mean value that we report, i.e. $\sigma_{jk}\left[N,\epsilon\right]$.

IV Density Matrix Reconstruction with an Equidistant Grid

In Fig. 4 we plot the relative reconstruction error for the polynomial approximation generated by Chebfun2 Townsend and Trefethen (2013) on an equidistant grid of points, such as that shown in Fig. 1D of the main text. Despite Chebfun2’s non-polynomial preprocessing of the sampled data, which aims to avoid the Runge phenomenon, the reconstruction is considerably worse than when the Padua Lagrange polynomial interpolant is used, as in the main text. In particular, the relative error is no longer strictly decreasing, and appears to saturate for some density matrix elements. This specific example highlights the benefit of performing Lagrange interpolation at the Padua points in comparison to state-of-the-art polynomial approximation on an equidistant grid.