# Continuous-Variable Shadow Tomography

Srilekha Gandhari  Joint Center for Quantum Information and Computer Science,
NIST/University of Maryland, College Park, Maryland 20742, USA  Victor V.
Albert  Joint Center for Quantum Information and Computer Science,
NIST/University of Maryland, College Park, Maryland 20742, USA  National
Institute of Standards and Technology, 100 Bureau Dr, Gaithersburg, MD 20899,
USA  Thomas Gerrits  National Institute of Standards and Technology, 100
Bureau Dr, Gaithersburg, MD 20899, USA  Jacob M. Taylor  Joint Center for
Quantum Information and Computer Science, NIST/University of Maryland, College
Park, Maryland 20742, USA  National Institute of Standards and Technology, 100
Bureau Dr, Gaithersburg, MD 20899, USA  Joint Quantum Institute,
NIST/University of Maryland, College Park, Maryland 20742, USA  Michael J.
Gullans  Joint Center for Quantum Information and Computer Science,
NIST/University of Maryland, College Park, Maryland 20742, USA  National
Institute of Standards and Technology, 100 Bureau Dr, Gaithersburg, MD 20899,
USA

###### Abstract

Shadow tomography is a framework for constructing succinct descriptions of
quantum states, called classical shadows, with powerful methods to bound the
estimators used. We recast existing experimental protocols for continuous-
variable tomography in the classical-shadow framework, obtaining rigorous
bounds on the sample complexity for estimating density matrices from these
protocols. We analyze the efficiency of homodyne, heterodyne, photon number
resolving (PNR), and photon-parity protocols. To reach a desired precision on
the classical shadow of an $N$-photon density matrix with a high probability,
we show that homodyne detection requires an order $\mathcal{O}(N^{5})$
measurements in the worst case, whereas PNR and photon-parity detection
require $\mathcal{O}(N^{4})$ measurements in the worst case (both up to
logarithmic corrections). We benchmark these results against numerical
simulation as well as experimental data from optical homodyne experiments. We
find that numerical and experimental homodyne tomography significantly
outperforms our bounds, exhibiting a more typical scaling of the number of
measurements that is close to linear in $N$. We extend our single-mode results
to an efficient construction of multimode shadows based on local measurements.

##  I Introduction

The ability to estimate states accurately with as few measurements as possible
— the primary goal of quantum state tomography — yields an indispensable tool
in quantum information processing [1]. State characterization is necessary for
realizing quantum technologies in finite-dimensional tensor-product quantum
systems governed by discrete variables (DV), as well as in quantum systems
governed by continuous variables (CV) [2], such as electromagnetic or
mechanical modes.

The recent development of classical-shadow tomography yields a succinct way to
learn information about a DV quantum state through an optimal number of
measurements, in such a way that the learned information can later be used to
predict properties of the state [3, 4, 5]. A key benefit of shadow tomography
is that it comes with rigorously proven guarantees on the minimum number of
samples required to achieve high accuracy with high probability. However, the
best guarantees require structures from finite-dimensional DV spaces, such as
state and unitary designs [4] or SIC-POVMs [6], thereby obscuring any
practical extension to the intrinsically infinite-dimensional CV systems.

In this paper, we apply the shadow tomography framework to a large family of
well-known and well-utilized CV tomographic protocols. We distill the
mathematical tools behind shadow guarantees in such a way that eliminates
dependence on strictly DV ingredients and allows us to reformulate established
CV protocols in the shadow framework. This reformulation yields accuracy
guarantees for expectation values of local observables whose required number
of samples scales polynomially with both the number of participating CV modes
and the maximum occupation (a.k.a. photon) number of each mode.

CV tomography is a long-standing and well-developed field [7]. Focusing on
established protocols allows us to boost their credibility, as opposed to
developing new protocols that may be equally efficient theoretically but whose
practical utility is left as an open question (e.g. [8, Sec. VI.A]).
Specifically, our work focuses on tomographic methods that can be easily
implemented experimentally with existing quantum optical technology, which
underpins fiber-based and free-space quantum communication and key
distribution [2]. To illustrate the connection to optics and our commitment to
the mantra of classical shadows “measure first, ask questions later” [5], we
apply our framework retroactively to experimental quantum-optical data
published in 2010 [9]. Naturally, our theoretical guarantees apply equally
well if the corresponding protocols are performed in other CV platforms such
as microwave cavities coupled to superconducting qubits [10], motional degrees
of freedom of trapped ions [11], as well as optomechanical [12] and nano-
acoustic [13] resonators.

The first method that we are able to recast in the shadow-tomography framework
is homodyne detection [14, 15, 16, 17, 7, 18]. We show that the number of
samples needs to scale at most as the fifth power of the maximum occupation
number (up to a logarithmic correction) to yield reliable homodyne shadow
estimates of a single-mode state. Since we require a finite occupation-number
cutoff, this bound only holds for portions of states supported on finite-
dimensional subspaces of the infinite-dimensional Fock space. Our theoretical
guarantees use several important technical bounds proved in an earlier work
[19]. While there have been studies of the statistical efficiency of homodyne
tomography (e.g., [20, 21, 19, 22, 23]), to our knowledge, our bounds are new.

Building off the formalism in [24], we are also able to recast a large class
of protocols utilizing displaced Fock states and make contact with heterodyne,
photon-number resolved (PNR), and photon-parity tomography [7]. For the latter
two methods, we show that the number of samples needs to scale at most as the
fourth power of the maximum occupation (up to a logarithmic correction) to
yield reliable PNR shadow estimates. For heterodyne measurements, our upper
bound is effectively exponential in the maximum occupation.

Following our theoretical analysis of single-mode systems, we provide a
generalization to multimode states. We prove that local CV shadow data can be
used to reconstruct $k$-local reduced density matrices on any subset of the
modes with a sample complexity polynomial in the total number of modes and
exponential in $k$.

We support our theoretical analysis with multiple numerical simulations and
the analysis of existing experimental homodyne data. We begin by verifying
that both the shadow methods indeed construct reliable estimates of a target
state. We then determine how the number of samples required for both the
methods scales with the maximum occupation number $N$. Through simulations
involving the reconstruction of the vacuum state, cat states and random pure
states, we observe that the scaling for PNR and photon-parity measurements is
similar to the calculated upper bound, whereas homodyne measurements exceed
theoretical expectations and surpass the observed PNR scaling in every case.
We also provide a basic demonstration of homodyne tomography on multimode
separable and entangled states, finding that the latter require more samples
to reach the same precision.

We begin with a general result for estimating density matrix elements of a
single-mode CV state using a generic shadow-based protocol in Sec. II. Its
application to homodyne measurements and photon-number resolving measurements
is discussed in Sec. III. In Sec. IV, we extend our theory to multimode CV
states, and support our methods by numerical results and analysis of
experimental data in Sec. V. We conclude in Sec. VI.

##  II Single-mode Shadows: General case

Here, we distill the necessary mathematical tools from DV shadow-based
accuracy guarantees in a general lemma that we apply to specific CV protocols
in later sections. We adapt the slight reformulation from [25, Appx. A] of the
original shadow work, in which guarantees are shown exclusively in terms of
the system state. This allows us to derive sampling guarantees that are
suboptimal relative to the qubit case, but are sufficiently general to be
applied to numerous CV tomographic protocols.

A single-mode CV density matrix can be represented in various ways, e.g.,
using Wigner, $Q$, $P$ [26, 27] or more general representations [24], moment
expansions [28], homodyne measurements [16, 17, 7], or outer products of
occupation-number (a.k.a. Fock) or other simple basis states. All of these
ways effectively express the density matrix as a linear combination of some
basis set of operators $\\{\sigma(\mu)\\}$ with corresponding coefficients
$p(\mu)$, where $\mu$ parameterizes some generally continuous and non-compact
space $\cal X$,

$\rho=\int_{\cal X}d\mu\,p(\mu)\,\sigma(\mu)~{}.$ (1)

Assuming a non-negative and normalized weight function $p(\mu)$ and
sufficiently well-behaved space $\cal X$ allows one to re-express the above
formula in statistical fashion, i.e., as an expectation value over samples
$\sigma$, sampled from $\cal X$ according to the probability distribution
governed by $p(\mu)$,

$\rho=\operatorname*{{\mathbb{E}}}_{\mu\sim\cal
X}\sigma(\mu)=\operatorname*{{\mathbb{E}}}_{\mu}\sigma(\mu)~{}.$ (2)

The above re-interpretation allows for a simple application of matrix
concentration inequalities (loosely speaking, the matrix version of the law of
large numbers; see Refs. [29, 30]). If we are given a set of $T$ independent
samples ${\cal{T}}=\\{\mu_{1},\ldots,\mu_{T}\\}$, we can construct an
estimator for the state by averaging the basis operators corresponding to the
outcomes in the set,

$\sigma_{T}=\frac{1}{T}\sum_{i=1}^{T}\sigma(\mu_{i}).$ (3)

We define the above estimator as a shadow of size T. Such shadows recover the
state by taking their expectation value over all sample sets $\cal T$,

$\operatorname*{{\mathbb{E}}}_{{\cal
T}}~{}\sigma_{T}=\frac{1}{T}\sum_{i=1}^{T}\operatorname*{{\mathbb{E}}}_{{\mu_{i}}}\sigma(\mu_{i})=\operatorname*{{\mathbb{E}}}_{{\mu}}~{}\sigma(\mu)=\rho~{}.$
(4)

This convergence property is a necessary starting point for shadow estimation,
but it says nothing about how fast the convergence happens as the number of
samples increases.

In order to bound the rate of convergence of the above estimate, we restrict
$\rho$ to finite-dimensional truncated versions $\rho^{\scaleto{(N)}{5pt}}$,
with the dimension cutoff $N$ assumed to be sufficiently large so that the
truncated operators capture the essence of the original ones. The resulting
matrices have non-zero entries only for row and column indices $0$ through
$N-1$. A “regularized” shadow is constructed with snapshots projected into the
$(N-1)$-photon subspace, using $P_{N}=\sum_{i=0}^{N-1}|i\rangle\langle i|$,

$\sigma_{T}^{\scaleto{(N)}{5pt}}\equiv
P_{N}\,\sigma_{T}\,P_{N}=\frac{1}{T}\sum_{i=1}^{T}P_{N}\,\sigma(\mu_{i})\,P_{N}.$
(5)

Our analysis can be readily extended to more general soft-cutoffs by replacing
the projectors onto fixed photon number spaces by operators with smooth
cuttoffs [31]. Note also that we do not re-normalize the above regularized
density matrices or shadows. Doing so introduces additional factors that would
complicate our calculations. Hence, for the sake of ease and uniformity, we
chose to work with the regularized matrices as they are.

We now bound the difference between the sample average $\sigma^{(N)}_{T}$ and
the actual average $\rho^{(N)}=P_{N}\rho P_{N}$, backing out the minimum
shadow size $T$ required for an accurate estimate. We employ the infinity norm
of the difference $\sigma_{T}^{\scaleto{(N)}{5pt}}-\rho^{\scaleto{(N)}{5pt}}$
to quantify the error.

######  Lemma II.1.

(a) Fix $\epsilon,\delta\in(0,1)$, and let $\sigma_{T}^{\scaleto{(N)}{5pt}}$
be an $N$-dimensional classical shadow of a single-mode continuous-variable
quantum state $\rho$. If the size of the classical shadow is at least $T$,
such that

$T=\frac{2N^{2}\left(\nu^{2}+R\epsilon/2N\right)}{\epsilon^{2}}\left(\log\,2N+\log\,1/\delta\right),$
(6)

where

$\displaystyle\nu^{2}$
$\displaystyle=\left\|\operatorname*{{\mathbb{E}}}_{{\cal
T}}\,\left(\sigma^{\scaleto{(N)}{5pt}}_{T}\right)^{2}\right\|_{\infty}=\left\|\operatorname*{{\mathbb{E}}}_{\mu}\,\left(\sigma^{\scaleto{(N)}{5pt}}(\mu)\right)^{2}\right\|_{\infty}$
(7a) $\displaystyle R$
$\displaystyle\geq\left\|\sigma^{\scaleto{(N)}{5pt}}(\mu)-\operatorname*{{\mathbb{E}}}_{\mu}\sigma^{\scaleto{(N)}{5pt}}(\mu)\right\|_{\infty},$
(7b)

then the probability

${\rm
Pr}\left(\|\sigma_{T}^{\scaleto{(N)}{5pt}}-\rho^{\scaleto{(N)}{5pt}}\|_{\infty}\leq\epsilon\right)\geq
1-\delta.$

###### Proof.

Applying the matrix Bernstein inequality [30][25, Appx. A] to
$\sigma_{T}^{\scaleto{(N)}{5pt}}-\rho^{\scaleto{(N)}{5pt}}$ and using Eq. (4)
yields

${\rm
Pr}\left(\|\sigma_{T}^{\scaleto{(N)}{5pt}}-\rho^{\scaleto{(N)}{5pt}}\|_{\infty}\geq\tilde{\epsilon}\right)\leq
2N\,\exp\left(-\frac{T\tilde{\epsilon}^{2}/2}{\nu^{2}+R\tilde{\epsilon}/3}\right)~{}.$

Using the equivalence relation $\|X\|_{\infty}\leq\|X\|_{1}\leq
N\|X\|_{\infty}$ for an $N$ dimensional system, we have ${\rm
Pr}(\|X\|_{1}\geq\epsilon)\leq{\rm Pr}(N\|X\|_{\infty}\geq\epsilon)$.
Substituting $\tilde{\epsilon}/N\rightarrow\epsilon$ we get

${\rm
Pr}\left(\|\sigma_{T}^{\scaleto{(N)}{5pt}}-\rho^{\scaleto{(N)}{5pt}}\|_{\infty}\geq\epsilon\right)\leq
2N\,\exp\left(-\frac{T\epsilon^{2}}{2N^{2}(\nu^{2}+R\epsilon/3N)}\right).$

Equating the error term on the right-hand side of the inequality to $\delta$
and rearranging yield the result (6). ∎

We substituted the original variance term
$\|\mathbb{E}_{\mu}\left(\sigma^{\scaleto{(N)}{5pt}}-\rho^{\scaleto{(N)}{5pt}}\right)^{2}\|_{\infty}$
in the matrix Bernstein inequality, which measures the fluctuations about the
mean $\rho^{\scaleto{(N)}{5pt}}$, with the term
$\|\mathbb{E}_{\mu}\left(\sigma^{\scaleto{(N)}{5pt}}\right)^{2}\|_{\infty}$
for ease in calculations. We noticed that this leads to the derived bounds
being higher than they would otherwise.

We see that the scaling of the required sample size in Eq. (6) depends
logarithmically on the inverse error probability $\delta$ and polynomially on
the accuracy $\epsilon$, recovering the corresponding portion of the qubit
shadow result [4]. To determine dependence on occupation number cutoff $N$, we
need to know the corresponding scaling of $R$ and $\nu$ (7), which are
sometimes referred to as shadow norms. This scaling depends on the detail of
the protocol, and we proceed to apply this lemma to specific cases.

##  III Applications

We apply the general shadow bound from Sec. II to a large family of CV
tomographic protocols that include homodyne and photon-number resolving
tomography.

###  III.1 Homodyne Shadows

Homodyne tomography is a popular technique used in quantum optics to make
quadrature measurements of CV quantum states [7]. This technique allows one to
measure a quadrature operator such as oscillator position, momentum, or any of
their superpositions.

Let $a,a^{\dagger}$ be the canonical mode lowering and raising operators
satisfying $[a,a^{\dagger}]=1$. A general quadrature operator,

$X_{\theta}={\textstyle\frac{1}{\sqrt{2}}}(ae^{-i\theta}+a^{\dagger}e^{i\theta})~{},$
(8)

is parameterized by $\theta\in[0,\pi)$, where $\theta=0$ ($\theta=\pi/2$)
corresponds to the pure position (momentum) quadrature. Each operator admits
its own set of orthogonal but non-normalizable “eigenstates”
$|x_{\theta}\rangle$ with $x_{\theta}\in\mathbb{R}$. Even though such
“eigenstates” are not quantum states but are instead distributions, they do
resolve the identity for each $\theta$ and thus yield valid positive operator-
valued measures (POVMs) [32]. For a given $\theta$, a homodyne protocol yields
a measurement in the corresponding set of eigenstates.

To perform a homodyne measurement, the unknown system state interferes with an
ancillary mode in a coherent state $|\alpha_{\text{LO}}\rangle$ (coherent
local oscillator or LO) via a beam splitter, as shown in Fig. 1(a). The phase
of the coherent state LO is precisely the phase $\theta$ which defines the
quadrature that is measured, and is picked uniformly from $[0,\pi)$. Both
output amplitudes of the beam splitter, $N_{1}$ and $N_{2}$ respectively, are
measured via homodyne detectors. Their difference $N_{-}$ is proportional to
the quadrature amplitude of the signal:

$N_{-}=\sqrt{2}\,|\alpha_{LO}|\,x_{\theta}\quad\text{with
probability}\quad\langle x_{\theta}|\rho|x_{\theta}\rangle~{}.$ (9)

Figure 1: (a) Balanced homodyne detection: the unknown state signal interferes
with a local oscillator (a coherent state $|\alpha\rangle$) and the
intensities at the output arms are measured. The difference in the intensities
is proportional to the quadrature amplitude of the unknown state in the mode
defined by the oscillator. (b) Photon-number resolving (PNR) measurement: the
local oscillator is in a coherent state $|\alpha\rangle$, and detecting $N$
photon at the detector occurs with a probability $p(N|\alpha)=\langle
N,\alpha|\rho|N,\alpha\rangle$.

To cast homodyne tomography into shadow terminology, we first specialize the
expansion in Eq. (1) to this case. Substituting parameters
$\mu\rightarrow(\theta,x_{\theta})$,

$\rho=\int_{0}^{\pi}d\theta\,\int_{\mathbb{R}}dx_{\theta}\,{\textstyle\frac{1}{\pi}}\langle
x_{\theta}|\rho|x_{\theta}\rangle F(x_{\theta},\theta)\,,$ (10)

where $1/\pi$ corresponds to picking $\theta$ uniformly from $[0,\pi)$, and
$\langle x_{\theta}|\rho|x_{\theta}\rangle$ to the conditional probability of
the homodyne setup yielding $x_{\theta}$ as the outcome. The matrix elements
of the basis operators
$F_{mn}(x_{\theta},\theta)=e^{i(m-n)\theta}f_{mn}(x_{\theta})$ are defined in
terms of what are known as pattern functions $f_{mn}(x)$ [33, 7]. These
functions are symmetric in $m,n$, and are defined using Fock state
wavefunctions $\psi_{m}(x)$ ($m^{th}$ energy eigen state of a harmonic
oscillator) and $\phi_{n}(x)$ ($n^{th}$ non-normalizable solution of the
Schrödinger equation of a harmonic oscillator) as

$f_{mn}(x)=\frac{\partial}{\partial
x}\left(\psi_{m}(x)\phi_{n}(x)\right)\quad\text{for}\quad n\geq m.$ (11)

Because the quadrature probability distributions are normalized for each
$\theta$, the above equation can automatically be interpreted in the
statistical fashion of Eq. (2). In this interpretation, basis operators $F$
are sampled according to the probability distribution defined by $\theta$ and
$\langle x_{\theta}|\rho|x_{\theta}\rangle$. A regularized homodyne CV shadow
of size $T$, constructed using the samples
$\\{x_{\theta_{i}},\theta_{i}\\}_{i=1}^{T}$, is therefore

$\sigma_{T}^{\scaleto{(N)}{5pt}}=\,\frac{1}{T}\sum_{i=1}^{T}P_{N}F(x_{\theta_{i}},\theta_{i})\,P_{N}~{}.$
(12)

We now adapt Lemma II.1 to homodyne shadows, relying on the bounds on pattern
function matrix from Ref. [19]. The following theorem gives an upper bound on
the sample complexity of constructing shadows through homodyne samples, for an
error
$\|\sigma_{T}^{\scaleto{(N)}{5pt}}-\rho^{\scaleto{(N)}{5pt}}\|_{\infty}\leq\epsilon$
with high probability (greater than $1-\delta$).

######  Theorem III.1.

Using homodyne tomography, for $\epsilon,\delta\in(0,1)$, some positive
constants $C_{2}$ and $C_{3}$ [19], and an N-dimensional homodyne shadow
$\sigma_{T}^{\scaleto{(N)}{5pt}}$ (12) of a single-mode CV state $\rho$,

${\rm
Pr}(\|\sigma_{T}^{\scaleto{(N)}{5pt}}-\rho^{\scaleto{(N)}{5pt}}\|_{\infty}\leq\epsilon)\geq
1-\delta$

when the size of the shadow is at least

$T=\frac{2N^{5}C_{2}C_{3}}{\epsilon^{2}}\left(\log\,1/\delta+\log\,2N\right).$

###### Proof.

From Theorem 3.2 of Ref. [19], we use
$\|\sigma_{T}^{\scaleto{(N)}{5pt}}\|^{2}_{\infty}\leq C_{1}N^{7/3}$ and
$\|\sigma_{T}^{\scaleto{(N)}{5pt}}\|^{2}_{2}\leq C_{2}N^{3}$, which gives us

$\|\sigma_{T}^{\scaleto{(N)}{5pt}}-\rho^{\scaleto{(N)}{5pt}}\|_{\infty}\leq\|\sigma_{T}^{\scaleto{(N)}{5pt}}\|_{\infty}+1\leq\sqrt{C_{1}}N^{7/6}+1=R$
(13)

and

$\displaystyle\nu^{2}=\left\|\operatorname*{{\mathbb{E}}}_{{\cal
T}}\,\sigma_{T}^{2}\right\|_{\infty}$ $\displaystyle=\left\|\int
dx\,d\theta\,{\textstyle\frac{1}{\pi}}\langle
x_{\theta}|\rho|x_{\theta}\rangle\sigma_{T}^{2}\right\|_{\infty}$ (14a)
$\displaystyle\leq\int dx\,d\theta\,{\textstyle\frac{1}{\pi}}\langle
x_{\theta}|\rho|x_{\theta}\rangle\|\sigma_{T}^{2}\|_{\infty}$ (14b)
$\displaystyle\leq\int dx\,d\theta\,{\textstyle\frac{1}{\pi}}\langle
x_{\theta}|\rho|x_{\theta}\rangle\|\sigma_{T}\|_{2}^{2}$ (14c)
$\displaystyle\leq\int dx\,d\theta\,{\textstyle\frac{1}{\pi}}\langle
x_{\theta}|\rho|x_{\theta}\rangle C_{2}N^{3}$ (14d) $\displaystyle\leq
C_{2}N^{3}~{}.$ (14e)

Therefore, using Lemma II.1, the size of the shadow needs to be at least

$\displaystyle T$
$\displaystyle=\frac{2N^{2}\left(\sigma^{2}+R\epsilon/2N\right)}{\epsilon^{2}}\left(\log\,2N-\log\,\delta\right)$
(15)
$\displaystyle\leq\frac{2N^{2}\left(C_{2}N^{3}+1+\sqrt{C_{1}}N^{1/6}\epsilon/2\right)}{\epsilon^{2}}(\log\,2N+\log\,1/\delta).$

∎

We see that the upper bound on the minimum number of samples scales as
$N^{5}\log\,N$ for a photon-cutoff $N$.

Regarding the utility of homodyne shadows in calculating expectation values of
observables, the pattern-function operators $F$ are unequivocally more complex
than the original qubit shadows [4]. Nevertheless, efficient methods are known
to compute them using recursive relations [34].

###  III.2 Photon Number Resolving (PNR) Shadows

Our second application is based on the $T$-operator formalism developed in
Ref. [24]. The relevant basis expansion of a single-mode CV state $\rho$ can
be written in terms of what are called T-operators,

$\rho={\textstyle\frac{1}{\pi}}\int_{\mathbb{C}}d^{2}\alpha\,{\rm
Tr}[\rho\,T(\alpha,\alpha^{*})]\,\bar{T}(\alpha,\alpha^{*})~{}.$ (16)

Here, ${\rm Tr}[\rho T(\alpha,\alpha^{*})]$ plays the role of a
quasiprobability distribution over the complex plane, and the dual operators
$\bar{T}(\alpha,\alpha^{*})$ are used to reconstruct the state given this
distribution.

The $T$-operators are defined in terms of displaced Fock states (a.k.a.
generalized coherent states), which are the result of applying a displacement
operator $D(\alpha)=\exp(\alpha a^{\dagger}-\alpha^{*}a)$ on a Fock state:
$|\alpha,n\rangle\equiv D(\alpha)|n\rangle$. There exists a family of
$T(\alpha,\alpha^{*})$, and their corresponding duals
$\bar{T}(\alpha,\alpha^{*})$, parameterized by $-1<r<1$, such that the dual
operator is given by

$\bar{T}_{r}(\alpha,\alpha^{*})=T_{-r}(\alpha,\alpha^{*}).$ (17)

$T(\alpha,\alpha^{*})$ belonging to this family can be diagonalized as

$\displaystyle T_{r}(\alpha,\alpha^{*})$
$\displaystyle=\sum_{n=0}^{\infty}\lambda^{(n)}_{r}D(\alpha)|n\rangle\langle
n|D(\alpha)^{\dagger}$ (18)
$\displaystyle=\sum_{n}\lambda^{(n)}_{r}|\alpha,n\rangle\langle\alpha,n|~{},$

where the eigenvalues $\lambda^{(n)}_{r}$ are

$\lambda^{(n)}_{r}=\frac{2(-1)^{n}(1-r)^{n}}{\pi(1+r)(1+r)^{n}}~{}.$ (19)

Utilizing the above expansions, Eq. (16) becomes

$\rho=\int_{\mathbb{C}}\frac{d^{2}\alpha}{\pi}\sum_{n=0}^{\infty}\langle
n,\alpha|\rho|n,\alpha\rangle\,\lambda^{(n)}_{r}\bar{T}_{r}(\alpha,\alpha^{*})~{},$
(20)

which we would like to interpret as an expectation value of samples
$\lambda^{(n)}_{r}\bar{T}_{r}(\alpha,\alpha^{*})$ distributed according to the
probabilities $\frac{1}{\pi}\langle n,\alpha|\rho|n,\alpha\rangle$. However,
in a physical interferometric measurement procedure corresponding to this
expansion [see Fig. 1(b)], the coherent-state parameter $\alpha$ comes from a
non-compact set. This means that we need to first truncate the sample region
in order to make the above a valid (i.e., normalizable) probability
distribution.

We restrict the integration of Eq. (20) to be over a phase-space region $\cal
A$ of finite area $A$. We can set this region to be, e.g., a disk of radius
$\alpha_{\text{max}}$. We will then be able to account for the truncated Fock
space by setting $\alpha_{\text{max}}^{2}=N$ since that is the average
occupation number of a coherent state with $|\alpha|=\alpha_{\text{max}}$. In
that case, $A=4\pi N$ scales linearly with maximum photon number cutoff.

Projecting into an $N$-dimensional subspace implies that we consider only
states with some max photon number $N-1$ and regularized $T$-operators
$\bar{T}_{r}^{\scaleto{(N)}{5pt}}(\alpha,\alpha^{*})$. Combining this with the
truncation of the phase-space integration yields the following approximate
expression for the truncated density matrix,

$\rho^{\scaleto{(N)}{5pt}}\approx\int_{\cal A}\frac{d^{2}\alpha}{\pi
A}\sum_{n=0}^{N-1}\langle
n,\alpha|\rho|n,\alpha\rangle\,\,A\lambda^{(n)}_{r}\bar{T}_{r}^{\scaleto{(N)}{5pt}}(\alpha,\alpha^{*})~{},$
(21)

in which we interpret
$A\,\lambda^{(n)}_{r}\,\bar{T}_{r}^{\scaleto{(N)}{5pt}}(\alpha,\alpha^{*})$ as
a PNR shadow sampled from the probability distribution

$p(n,\alpha)=\frac{1}{A}\langle n,\alpha|\rho|n,\alpha\rangle$ (22)

with parameters $0\leq n\leq N-1$ and $\alpha\in\cal A$.

The eigenvalues $\lambda^{(n)}_{r}$ exhibit different behaviours for different
values of $r$, and in Appx. A we analyze how many samples are required for
density matrix estimation in each case. We obtain quartic scaling with $N$ at
$r=0$, up to logarithmic corrections, and unfavorable exponential scaling for
all other $0<|r|<1$. In the favorable case, the minimum number of samples is

$T\sim\frac{N^{2}\,A^{2}}{\epsilon^{2}}(\log\,2N+\log\,1/\delta)=\mathcal{O}(N^{4}\log
N)~{},$ (23)

beating our theoretical estimates for homodyne shadows.
Therefore, we see that $T$-operators at $r=0$ have a favorable sample
complexity scaling and can be used as a tomographic method to estimate an
unknown state. Moreover, the case of $r=0$ can also be realized via photon-
parity measurement tomography [35, 36, 37, 38], providing another experimental
technique to realize these shadows.

It is worth noting that in the limit $r=-1$,
$\bar{T}_{r}(\alpha,\alpha^{*})=T_{-r}(\alpha,\alpha^{*})$ becomes a projector
onto the coherent state $|\alpha\rangle$, making the measurement protocol
equivalent to what is known as heterodyne measurement, achieved using a
similar setup as homodyne measurement [39]. The basis expansion is given in
terms of the Q-function
$Q(\alpha)=\frac{1}{\pi}\langle{\alpha}|\rho|{\alpha}\rangle$. However, the
sample complexity upper bound diverges in this case. It is possible to avoid
this divergence by taking a limit $r\to-1$ as the number of samples increases.
Unfortunately, this method only achieves an upper bound that is exponential in
$N$ because, as noted above, for any $0<|r|<1$, our sample complexity upper
bound diverges exponentially in $N$ (see Appx. A).

##  IV Multimode Shadows

In this section, we show that it is possible to construct efficient
representations of reduced density matrices on constant numbers of modes from
local CV shadow data.

Multimode CV states are states that are made up of multiple modes of the
underlying system, such as different frequencies of light in an optical
system. This is analogous to having a multi-qubit state in DV. Each mode is
equipped with a generally continuous and non-compact space $\mathcal{X}_{i}$
and M-mode, possibly entangled, quantum states reside in
$\mathcal{X}=\bigotimes_{i=1}^{M}\mathcal{X}_{i}$. Let the space $\mathcal{X}$
be parameterized by a variable $\boldsymbol{\mu}=(\mu^{1},\ldots,\mu^{M})$,
and spanned by the set of basis operators
$\\{\boldsymbol{\sigma}(\boldsymbol{\mu})\\}$. Extending Eq. (1), the density
matrix of a general multimode state belonging to this space can then be
written as

$\boldsymbol{\rho}=\int_{\mathcal{X}}d\boldsymbol{\mu}~{}p(\boldsymbol{\mu})~{}\boldsymbol{\sigma}(\boldsymbol{\mu}),$
(24)

where $p(\boldsymbol{\mu})$ are the corresponding coefficients. Boldface
symbols will be used to represent multimode variables and states henceforth.
All the density operators and shadows considered in this section are the
regularized versions, projected into the $N$-photon subspace of each mode. The
superscripts indicating this are omitted for the sake of brevity.

Analogous to the single-mode case in Sec. II, we require a basis expansion of
$\boldsymbol{\rho}$ such that $p(\boldsymbol{\mu})\geq
0~{}\forall\boldsymbol{\rho}$; the non-negativity of $p(\boldsymbol{\mu})$
allowing it be interpreted as a probability distribution over
$\boldsymbol{\mu}$. We proceed by assuming that such an expansion exists and
that
$p(\boldsymbol{\mu})=\langle\boldsymbol{\mu}|\boldsymbol{\rho}|\boldsymbol{\mu}\rangle$,
implying that the underlying space is parameterized by quantum states or
distributions $|\boldsymbol{\mu}\rangle$.

We construct shadows of multimode states using local, single-mode
measurements. For an estimator of a state $\boldsymbol{\rho}$ to qualify as a
shadow, as introduced in Section Sec. II, it must be equal to the state
$\boldsymbol{\rho}$ in expectation. Below we show that it can indeed be
satisfied for shadows constructed through local measurements. For this
calculation, we utilize another expansion of a general quantum state
$\boldsymbol{\rho}$ in terms of some set of separable operators
$\\{\boldsymbol{\rho}_{\alpha}\\}$ and corresponding complex coefficients
$\\{c_{\alpha}\\}$,

$\boldsymbol{\rho}=\sum_{\alpha}c_{\alpha}\boldsymbol{\rho}_{\alpha}=\sum_{\alpha}c_{\alpha}\bigotimes_{i=1}^{M}\rho_{\alpha}^{i}~{}.$
(25)

Such an expansion is possible irrespective of whether the state is separable
or entangled. The estimator becomes

$\displaystyle\operatorname*{{\mathbb{E}}}_{\boldsymbol{\mu}}~{}\boldsymbol{\sigma}(\boldsymbol{\mu})$
$\displaystyle=\int_{\mathcal{X}}d\boldsymbol{\mu}\times\langle\boldsymbol{\mu}|\boldsymbol{\rho}|\boldsymbol{\mu}\rangle\times\bigotimes_{i=1}^{M}\sigma(\mu^{i})$
(26a)
$\displaystyle=\sum_{\alpha}c_{\alpha}\int_{\mathcal{X}}d\boldsymbol{\mu}\times\bigotimes_{i=1}^{M}\langle\mu^{i}|\rho_{\alpha}^{i}|\mu^{i}\rangle~{}\sigma(\mu^{i})$
(26b)
$\displaystyle=\sum_{\alpha}c_{\alpha}\bigotimes_{i=1}^{M}\int_{\mathcal{X}_{i}}d\mu^{i}\langle\mu^{i}|\rho_{\alpha}^{i}|\mu^{i}\rangle\sigma(\mu^{i})$
(26c)
$\displaystyle=\sum_{\alpha}c_{\alpha}\boldsymbol{\rho}_{\alpha}=\boldsymbol{\rho}~{}.$
(26d)

In the following theorem, we calculate the variance and thereby the sample
complexity of estimating a multimode state through local CV shadows. This
result reduces to the single-mode result from II.1 when $M=1$.

######  Theorem IV.1.

Fix $\epsilon,\delta\in(0,1)$, and let $\boldsymbol{\sigma}_{T}$ be an
$N^{M}$-dimensional classical shadow obtained through local measurements
$\boldsymbol{\sigma}(\boldsymbol{\mu})=\bigotimes_{i=1}^{M}\sigma(\mu^{i})$ of
an $M$-mode CV state $\boldsymbol{\rho}$ (25). If the size of the classical
shadow is at least

$T=\frac{2N^{2M}(\nu^{2M}_{1}\sum_{\alpha}|c_{\alpha}|+R\epsilon/2N)}{\epsilon^{2}}(\log
2N+\log 1/\delta),$ (27)

then

${\rm
Pr}\left(\|\boldsymbol{\sigma}_{T}-\boldsymbol{\rho}\|_{\infty}\leq\epsilon\right)\geq
1-\delta,$

written in terms of the single-mode variance $\nu^{2}_{1}$ (7a) and
$R\geq\|\boldsymbol{\sigma}(\boldsymbol{\mu})-\operatorname*{{\mathbb{E}}}_{\boldsymbol{\mu}}\boldsymbol{\sigma}(\boldsymbol{\mu})\|_{\infty}$.

###### Proof.

The measurement variance for an $M$-mode state proceeds similar to the single-
mode case with the help of Eq. (25),

$\displaystyle\nu_{M}^{2}$
$\displaystyle=\left\|\operatorname*{{\mathbb{E}}}_{\boldsymbol{\mu}}~{}\left(\boldsymbol{\sigma}(\boldsymbol{\mu})\right)^{2}\right\|_{\infty}$
(28a) $\displaystyle=\left\|\int_{{\cal
X}}\prod_{j=1}^{M}d\mu^{j}\times\langle\boldsymbol{\mu}|\boldsymbol{\rho}|\boldsymbol{\mu}\rangle\times\bigotimes_{i=1}^{M}\left(\sigma(\mu^{i})\right)^{2}\right\|_{\infty}$
(28b)
$\displaystyle=\left\|\sum_{\alpha}c_{\alpha}\bigotimes_{i}\int_{{\cal{\cal
X}}_{i}}d\mu^{i}\langle\mu^{i}|\rho_{\alpha}^{i}|\mu^{i}\rangle\left(\sigma(\mu^{i})\right)^{2}\right\|_{\infty}$
(28c) $\displaystyle\leq\sum_{\alpha}|c_{\alpha}|\,\left\|\bigotimes_{i}\int
d\mu^{i}\langle\mu^{i}|\rho_{\alpha}^{i}|\mu^{i}\rangle\left(\sigma(\mu^{i})\right)^{2}\right\|_{\infty}$
(28d) $\displaystyle=\sum_{\alpha}|c_{\alpha}|\,\prod_{i}\left\|\int
d\mu^{i}\langle\mu^{i}|\rho_{\alpha}^{i}|\mu^{i}\rangle\left(\sigma(\mu^{i})\right)^{2}\right\|_{\infty}$
(28e) $\displaystyle\leq\sum_{\alpha}|c_{\alpha}|\,\prod_{i}\nu_{1}^{2}$ (28f)
$\displaystyle\leq\nu_{1}^{2M}\sum_{\alpha}|c_{\alpha}|~{}.$ (28g)

Using matrix Bernstein inequality and following a calculation similar to
Lemma. II.1, we get the desired result in Eq. (27). ∎

Figure 2: Single-mode shadows: (a) Density matrix elements and (b) Wigner
function representation of the original state, a homodyne shadow and a PNR
shadow of an even cat state
$|\psi\rangle\propto|\alpha\rangle+|-\alpha\rangle$ with $\alpha=\sqrt{10}$.
(c-d) Numerical sample complexity of single-mode shadows; precision and
confidence parameters are $\epsilon=0.1$ and $\delta=0.1$. (c) Homodyne
Shadows: Vacuum state ($T\propto N^{1.389}$), cat state with amplitude
$\alpha=\sqrt{10}$ ($T\propto N^{1.820}$), and random pure states ($T\propto
N^{1.212}$). (d) PNR Shadows: Vacuum state ($T\propto N^{3.51}$) and pure
random states ($T\propto N^{2.49}$)

One can further generalize this result to show that it is possible to obtain
the reduced density matrix on any choice of $k$-local modes in the multimode
system. The proof follows from a simple application of the union bound and
increases the sample complexity by a factor of $k$ [4, 25].

######  Corollary IV.1.1.

If $\boldsymbol{\rho}_{k}$ denotes the reduced density matrix of
$\boldsymbol{\rho}$ in $k\leq M$ modes, and its estimator constructed through
local measurements is $\boldsymbol{\sigma}_{T,k}$, such that the size of the
shadow is at least

$T=\frac{2N^{2k}(\nu^{2k}_{1}\sum_{\alpha}|c_{\alpha}|+R\epsilon/2N)}{\epsilon^{2}}(\log
2N+\log 1/\delta)$

then

${\rm
Pr}\left(\|\boldsymbol{\sigma}_{T,k}-\boldsymbol{\rho}_{k}\|_{\infty}\leq\epsilon\right)\geq
1-\delta.$

###### Proof.

Without loss of generality, let’s assume $\boldsymbol{\rho}_{k}$ to be the
reduced density matrix over the modes $1,\ldots,k$, with $k\leq M$, i.e.
$\boldsymbol{\rho}_{k}={\rm Tr}_{k+1,..,M}(\boldsymbol{\rho})$. Define
$\boldsymbol{\mu}_{k}=(\mu^{1},\ldots,\mu^{k})$ to parameterize this $k$-mode
subspace, and the local operators

$\boldsymbol{\sigma}_{k}(\boldsymbol{\mu}_{k})={\rm
Tr}_{k+1,..,M}(\boldsymbol{\sigma}(\boldsymbol{\mu}))=\bigotimes_{i=1}^{k}\sigma(\mu^{i}).$

The classical shadow of the reduced matrix $\boldsymbol{\sigma}_{T,k}$ is then
the average of $T$ such operators. The variance of estimating this reduced
matrix, following earlier calculations, is

$\displaystyle\nu_{k}^{2}=\left\|\operatorname*{{\mathbb{E}}}_{\boldsymbol{\mu}_{k}}\left(\boldsymbol{\sigma}_{k}(\boldsymbol{\mu}_{k})\right)^{2}\right\|_{\infty}\leq\nu_{1}^{2k}\sum_{\alpha}|c_{\alpha}|~{}.$
(29)

∎

##  V Numerical and Experimental tests

In this section, we verify our shadow construction methods and present
numerical results on sample complexity scaling, with the maximum occupation
number $N$, of single-mode shadows. We also analyze experimental homodyne data
and demonstrate the reconstruction of multimode states.

###  V.1 Single mode simulations

To analyze the applicability of our theoretical guarantees, we construct
estimators of density matrices of certain target states using homodyne and PNR
shadows. We then analyze how close they are to the target states. Figure 2(a)
shows how a reconstructed cat state with mean-photon number $10$ compares to
the original density matrix elements for both homodyne and PNR tomography. A
comparison in terms of the Wigner function is shown in Fig. 2(b). To achieve
equally accurate shadows such that
$\|\sigma_{T}^{\scaleto{(N)}{5pt}}-\rho^{\scaleto{(N)}{5pt}}\|_{\infty}=0.1$,
we had to use $5\times 10^{4}$ and $10^{6}$ samples for homodyne and PNR,
respectively.

Of primary interest in the classical shadow formalism is the number of samples
required to achieve a given target precision with high-probability. Through
numerical simulations, we determined what the probability $\delta$ is for a
shadow of size $T$ to have a distance (infinity norm) greater than a fixed
error, say $\epsilon_{0}$, from the original state. On repeating this multiple
times with different shadow sizes, we observed that for a fixed error in
distance, $T\propto\log 1/\delta$. With this knowledge we can estimate the
size of a shadow required to achieve an error less than $\epsilon_{0}$ with a
probability greater than some fixed $1-\delta_{0}$.

Figure 3: Experimental test of single-mode shadows: (a) Density matrix
elements and (b) Wigner function representation of the original cat state
(left) and the reconstructed state (right) through homodyne measurements.
Original cat state: $|\psi\rangle\propto|\alpha\rangle+|-\alpha\rangle$,
$\alpha=1.32$, $T=3.10^{5}$. (c) Sample complexity scaling with maximum photon
number $N$, from the homodyne data of the experiment. The two methods ‘TES2’
and ‘APD1’ correspond to photon-detection through transition edge sensor (TES)
and avalanche photodiode (APD), respectively. The precision and confidence
parameters are $\epsilon=0.4$ and $\delta=0.1$. (d-e) Homodyne two-mode
shadows: Density matrix elements of the original states and numerically
reconstructed states. The number of samples $T=5.10^{4}$, $N=2$ and
$\alpha=\sqrt{1.5}$ in both cases. (d) 2-mode separable state:
$|\Psi\rangle=|\psi_{\alpha}\rangle\otimes|\psi_{\alpha}\rangle$, where
$|\psi_{\alpha}\rangle\propto|\alpha\rangle+|-\alpha\rangle$, and the infinity
norm distance of the reconstructed state to the ideal state is
$\epsilon=0.03$. (e) 2-mode entangled state:
$|\Psi\rangle\propto|\alpha,\alpha\rangle+|-\alpha,-\alpha\rangle$, and
$\epsilon=0.05$.

Fig. 2 (c)-(d) shows how $T$, for $\epsilon=0.1$ and $\delta=0.1$, changes
with maximum photon-number $N$ for various states using homodyne and PNR
tomography. We ran our simulations for a vacuum state
$|\psi\rangle=|0\rangle$, a cat state
$|\psi\rangle\propto|\alpha\rangle+|-\alpha\rangle$ with mean photon number
$|\alpha|^{2}=10$, and random pure states. The random pure states of size $N$
had non-zero support only up to $(N-1)^{\text{th}}$ Fock state. The results
show that homodyne has a scaling better than $N^{2}$ for the cases we tried,
consistently outperforming the theoretical upper bound. PNR on the other hand
has a scaling much closer to the predicted upper bound of $N^{4}$. However, we
are not as certain about the ultimate scaling of random pure states because
numerical precision issues prevented us from obtaining more than two points.

###  V.2 Experimental data

Here, we analyze the sample complexity of experimental homodyne data from Ref.
[9]. The experiment consisted of creating Coherent Superposition States
(CSSs), where the created states were approximations to small-amplitude, even
cat states $|\psi\rangle\propto|\alpha\rangle+|-\alpha\rangle$. Squeezed
vacuum was generated through spontaneous parametric down-conversion in a
$\text{KNbO}_{3}$ non-linear crystal: an up-converted laser pulse at $430$ nm
created the squeezed vacuum at $860$ nm. After spectral filtering, the
squeezed vacuum was incident on a weakly reflecting beam splitter with
reflectivity $R$. Photons detected in the reflected path herald a CSS emerging
from the transmitted port of the beam splitter. The resulting CSS was then
directed to a homodyne detection setup. The homodyne setup consisted of a
50/50 beam splitter, two high-efficiency photodiodes and a low-noise
amplification circuit. The heralded photon detection was done using either two
single-photon avalanche diodes (SPADs) [40] or one transition edge sensor
(TES) [41]. The experiment analyzed here used both heralded-photon detection
schemes. The CSSs were heralded upon the detection of two heralding photons,
generating the approximation of a small-amplitude even cat state.

Using the homodyne data from this experiment, we constructed shadows of the
state, an example of which is shown in Fig. 3(a-b), and studied how close the
reconstructed states are to the target states. In Fig. 3(c), we show that the
number of samples required scales as $T\propto N^{0.659\pm 0.001}$ for a cat
state of amplitude $\alpha=1.16$ using the TES, and as $T\propto N^{0.685\pm
0.001}$ for a cat state of amplitude $\alpha=1.32$ using the SPADs.

###  V.3 Multimode simulations

We now provide an example reconstruction of both separable and entangled two-
mode states.

In Fig. 3(d) we can see original and reconstructed density matrices of a
separable 2-mode state
$|\Psi\rangle=|\psi_{\alpha}\rangle\otimes|\psi_{\alpha}\rangle$, where
$|\psi_{\alpha}\rangle$ is an even cat state of amplitude $\alpha=\sqrt{1.5}$.
In Fig. 3(e), we have the reconstruction of an entangled state
$|\Psi\rangle=|\alpha\rangle\otimes|\alpha\rangle+|-\alpha\rangle\otimes|-\alpha\rangle$,
with the same number of measurements. The infinity-norm distance $\epsilon$
between the original and reconstructed states is 0.03 and 0.05, respectively.
The increased sample complexity needed to reach a given precision for the
entangled state tomography is consistent with our general expectations from
the theoretical analysis in Sec. IV.

##  VI Outlook and Conclusions

We developed a formalism to derive worst-case sample complexity of a large
class of state-reconstruction procedures for single-mode continuous-variable
(CV) states. In an experimental context, our work provides a tool to compare
the sample complexity of different measurement methods. For photon-number
resolving measurements, we show through numerical simulations that our
analytical derived upper bound is close to the observed complexity. In the
case of homodyne tomography, we notice that the simulation results surpass our
derived bound, rendering homodyne tomography more efficient than PNR and
photon-parity tomography in practice.

We obtain the sample complexity of estimating multimode CV states by
reconstructing shadows of the true states from local measurements on each
mode, in the spirit of the qubit-based shadow tomography proposal [3, 4, 5].
We also consider a CV manifestation of a global version of the original shadow
protocol, where one uses global measurements over varying linear combinations
of the constituent modes for state reconstruction. Of particular note is that,
in contrast to the discrete-variable systems, our analysis does not rely on
the construction of state designs (approximate or otherwise) such as in [8,
Sec. VI.A]. Instead, we focus on adapting existing techniques for CV
tomography to fit within the framework of randomized measurements.

In summary, we recast existing homodyne and photon-number-resolving protocols
as shadow-tomography protocols and show that such protocols yield good local
estimates of a multimode state using a number of samples that is polynomial in
the number of modes. Our CV shadow framework can be extended to analyze other
metrics like the total variation distance between two outcome probability
distributions, which is valuable for verifying sampling experiments like Boson
sampling [42]. Another remaining open question is to determine robustness of
CV protocols to noise, e.g., within the framework of robust shadow estimation
[43]. Deriving rigorous lower bounds that can match the numerically and
experimentally observed scaling for homodyne tomography is also an important
outstanding challenge.

###### Acknowledgements.

We thank Yi-Kai Liu, Dominik Hangleiter and Kunal Sharma for helpful
discussions. MJG and VVA thank Joe Iosue and Kunal Sharma for discussions and
collaborations on [8]. This work is supported in part by NIST grant
70NANB21H055_0 and NSF QLCI grant OMA-2120757. VVA thanks Olga Albert and
Ryhor Kandratsenia for providing daycare support throughout this work.

##  Appendix A Scaling with PNR Tomography

Beginning with the approximate expression for a truncated density matrix in
the displaced Fock basis (30), we want to estimate the sample complexity of
PNR shadows.

$\rho^{\scaleto{(N)}{5pt}}\approx\int_{\cal A}\frac{d^{2}\alpha}{\pi
A}\sum_{n=0}^{N-1}\langle
n,\alpha|\rho|n,\alpha\rangle\,\,A\lambda^{(n)}_{r}\bar{T}_{r}^{\scaleto{(N)}{5pt}}(\alpha,\alpha^{*})~{}$
(30)

The above expansion in terms of $T_{r}(\alpha,\alpha^{*})$ is valid for
$-1<r<1$ as mentioned in Sec. III.2. Depending on whether $r$ is positive or
negative, one of $\lambda^{(n)}_{r}$ or $\lambda^{(n)}_{-r}$ (the latter being
an eigenvalue of $T_{-r}(\alpha,\alpha^{*})$) increases with $n$ and the other
decreases with $n$. For $r=0$, the eigenvalues $\lambda^{(n)}_{0}$ are
independent of $n$. Using the triangle inequality and Eq. (19), the shadow
norm bound (7b) becomes

$\displaystyle\|A\lambda^{(n)}_{r}T_{-r}^{\scaleto{(N)}{5pt}}-\rho^{\scaleto{(N)}{5pt}}\|_{\infty}$
$\displaystyle\leq
A\|\lambda^{(n)}_{r}T_{-r}^{\scaleto{(N)}{5pt}}\|_{\infty}+A$ (31a)
$\displaystyle\leq A|\lambda^{(N)}_{-|r|}\,\lambda^{(0)}_{|r|}|+A$ (31b)
$\displaystyle\leq\frac{2A|\lambda^{(N)}_{-|r|}|}{\pi(1+|r|)}+A\equiv R,$
(31c)

with the equality holding at $r=0$. As for the variance from Eq. (7a), we have

$\displaystyle\nu^{2}$
$\displaystyle=\left\|\operatorname*{{\mathbb{E}}}_{n,\alpha}\left(A\lambda_{r}^{(n)}T_{-r}^{\scaleto{(N)}{5pt}}\right)^{2}\right\|_{\infty}$
(32a) $\displaystyle=\left\|\int_{{\cal A}}\frac{d^{2}\alpha}{\pi
A}\,\sum_{n=0}^{N-1}p(n,\alpha)\left(A\lambda_{r}^{(n)}T_{-r}^{\scaleto{(N)}{5pt}}\right)^{2}\right\|_{\infty}$
(32b) $\displaystyle\leq\int_{{\cal A}}\frac{d^{2}\alpha}{\pi
A}\,\sum_{n=0}^{N-1}p(n,\alpha)\,A^{2}\left\|(\lambda_{r}^{(n)}T_{-r}^{\scaleto{(N)}{5pt}})^{2}\right\|_{\infty}~{}$
(32c) $\displaystyle\leq\int_{{\cal A}}\frac{d^{2}\alpha}{\pi
A}\,\sum_{n=0}^{N-1}p(n,\alpha)\,A^{2}|\lambda^{(0)}_{|r|}\lambda^{(N)}_{-|r|}|^{2}$
(32d) $\displaystyle\leq A^{2}|\lambda^{(0)}_{|r|}\lambda^{(N)}_{-|r|}|^{2}$
(32e)

The remaining norm can be bounded using the formula (19) as was done in Eq.
(31), and the remaining sum and integral can be upper-bounded by a factor of
$A$, yielding

$\nu^{2}\leq\frac{16\,A^{2}}{\pi^{4}(1-|r|^{2})^{2}}\left(\frac{1+|r|}{1-|r|}\right)^{2N}~{}.$
(33)

Applying Lemma II.1, and using $A\propto N$, the number of samples is

$T\sim\frac{32\,N^{4}}{\pi^{4}\epsilon^{2}(1-|r|^{2})^{2}}\left(\frac{1+|r|}{1-|r|}\right)^{2N}\log\left(\frac{2N}{\delta}\right)~{},$
(34)

growing exponentially in $N$ for $r\neq 0$.

## References

  * Nielsen and Chuang [2011] M. A. Nielsen and I. L. Chuang, _Quantum Computation and Quantum Information: 10th Anniversary Edition_ (Cambridge University Press, 2011).
  * Braunstein and van Loock [2005] S. L. Braunstein and P. van Loock, _Quantum information with continuous variables_ , Rev. Mod. Phys. 77, 513 (2005).
  * Aaronson [2018] S. Aaronson, _Shadow tomography of quantum states_ (2018), eprint 1711.01053.
  * Huang et al. [2020] H.-Y. Huang, R. Kueng, and J. Preskill, _Predicting many properties of a quantum system from very few measurements_ , Nature Physics 16, 1050–1057 (2020).
  * Elben et al. [2022] A. Elben, S. T. Flammia, H.-Y. Huang, R. Kueng, J. Preskill, B. Vermersch, and P. Zoller, _The randomized measurement toolbox_ , arXiv:2203.11374 (2022).
  * Stricker et al. [2022] R. Stricker, M. Meth, L. Postler, C. Edmunds, C. Ferrie, R. Blatt, P. Schindler, T. Monz, R. Kueng, and M. Ringbauer, _Experimental single-setting quantum state tomography_ (2022), eprint 2206.00019.
  * Lvovsky and Raymer [2009] A. I. Lvovsky and M. G. Raymer, _Continuous-variable optical quantum-state tomography_ , Rev. Mod. Phys. 81, 299 (2009).
  * [8] J. T. Iosue, K. Sharma, M. J. Gullans, and V. V. Albert, Continuous-variable quantum state designs: theory and applications, arXiv, same day posting (2022).
  * Gerrits et al. [2010] T. Gerrits, S. Glancy, T. S. Clement, B. Calkins, A. E. Lita, A. J. Miller, A. L. Migdall, S. W. Nam, R. P. Mirin, and E. Knill, _Generation of optical coherent-state superpositions by number-resolved photon subtraction from the squeezed vacuum_ , Phys. Rev. A 82, 031802 (2010).
  * Girvin [2014] S. M. Girvin, _Circuit qed: superconducting qubits coupled to microwave photons_ , Quantum machines: measurement and control of engineered quantum systems pp. 113–256 (2014).
  * de Neeve et al. [2022] B. de Neeve, T.-L. Nguyen, T. Behrle, and J. P. Home, _Error correction of a logical grid state qubit by dissipative pumping_ , Nature Physics 18, 296 (2022).
  * Aspelmeyer et al. [2014] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, _Cavity optomechanics_ , Reviews of Modern Physics 86, 1391 (2014).
  * MacCabe et al. [2020] G. S. MacCabe, H. Ren, J. Luo, J. D. Cohen, H. Zhou, A. Sipahigil, M. Mirhosseini, and O. Painter, _Nano-acoustic resonator with ultralong phonon lifetime_ , Science 370, 840 (2020).
  * Vogel and Risken [1989] K. Vogel and H. Risken, _Determination of quasiprobability distributions in terms of probability distributions for the rotated quadrature phase_ , Phys. Rev. A 40, 2847 (1989).
  * Smithey et al. [1993] D. T. Smithey, M. Beck, M. G. Raymer, and A. Faridani, _Measurement of the wigner distribution and the density matrix of a light mode using optical homodyne tomography: Application to squeezed states and the vacuum_ , Phys. Rev. Lett. 70, 1244 (1993).
  * Leonhardt and Paul [1995] U. Leonhardt and H. Paul, _Measuring the quantum state of light_ , Progress in Quantum Electronics 19, 89 (1995).
  * Bisio et al. [2009] A. Bisio, G. Chiribella, G. M. D’Ariano, S. Facchini, and P. Perinotti, _Optimal quantum tomography_ , IEEE Journal of Selected Topics in Quantum Electronics 15, 1646 (2009).
  * D'Ariano [1995] G. M. D'Ariano, _Tomographic measurement of the density matrix of the radiation field_ , Quantum and Semiclassical Optics: Journal of the European Optical Society Part B 7, 693 (1995).
  * Artiles et al. [2005] L. M. Artiles, R. D. Gill, and M. I. Guta, _An invitation to quantum tomography_ , J. R. Stat. Soc. Ser. B (Statistical Methodol. 67, 109 (2005), eprint arXiv:math/0405595.
  * Banaszek [1998] K. Banaszek, _Maximum-likelihood estimation of photon-number distribution from homodyne statistics_ , Phys. Rev. A 57, 5013 (1998).
  * Gill and Guta [2003] R. Gill and M. Guta, _An invitation to quantum tomography_ (2003), eprint arXiv:quant-ph/0303020.
  * Albini et al. [2009] P. Albini, E. De Vito, and A. Toigo, _Quantum homodyne tomography as an informationally complete positive-operator-valued measure_ , Journal of Physics A: Mathematical and Theoretical 42, 295302 (2009).
  * Alquier et al. [2013] P. Alquier, K. Meziani, and G. Peyré, _Adaptive estimation of the density matrix in quantum homodyne tomography with noisy data_ , Inverse Problems 29, 075017 (2013).
  * Wunsche [1991] A. Wunsche, _Displaced fock states and their connection to quasiprobabilities_ , Quantum Optics: Journal of the European Optical Society Part B 3, 359 (1991).
  * Huang et al. [2021] H.-Y. Huang, R. Kueng, G. Torlai, V. V. Albert, and J. Preskill, _Provably efficient machine learning for quantum many-body problems_ (2021), eprint 2106.12627.
  * Hillery et al. [1984] M. Hillery, R. O’Connell, M. Scully, and E. Wigner, _Distribution functions in physics: Fundamentals_ , Physics Reports 106, 121 (1984).
  * Ferrie [2011] C. Ferrie, _Quasi-probability representations of quantum theory with applications to quantum information science_ , Reports on Progress in Physics 74, 116001 (2011).
  * Wünsche [1996] A. Wünsche, _Tomographic reconstruction of the density operator from its normally ordered moments_ , Phys. Rev. A 54, 5291 (1996).
  * [29] Tropp, J.A. (2016). The Expected Norm of a Sum of Independent Random Matrices: An Elementary Approach. In: Houdré, C., Mason, D., Reynaud-Bouret, P., Rosiński, J. (eds) High Dimensional Probability VII. Progress in Probability, vol 71. Birkhäuser, see pgs. 173–202.
  * Tropp [2015] J. A. Tropp, _An introduction to matrix concentration inequalities_ , Foundations and Trends in Machine Learning 8, 1 (2015).
  * [31] For example, for a regulator $R$ one can find a sufficiently large $N$ such that $||P_{N}R\rho R^{\dagger}P_{N}-R\rho R^{\dagger}||<\epsilon$ for all density matrices $\rho$ and apply our results directly based on the triangle inequality. More generally, the bounds we prove on the projected shadows could be adapted to regulators with smooth cutoffs, but we leave this approach for future work.
  * Holevo [2011] A. Holevo, _Probabilistic and Statistical Aspects of Quantum Theory_ (Edizioni della Normale, Pisa, 2011).
  * D’Ariano et al. [1994] G. M. D’Ariano, C. Macchiavello, and M. G. A. Paris, _Detection of the density matrix through optical homodyne tomography without filtered back projection_ , Phys. Rev. A 50, 4298 (1994).
  * Richter [2000] T. Richter, _Realistic pattern functions for optical homodyne tomography and determination of specific expectation values_ , Physical Review A 61 (2000).
  * Lutterbach and Davidovich [1997] L. G. Lutterbach and L. Davidovich, _Method for direct measurement of the wigner function in cavity qed and ion traps_ , Phys. Rev. Lett. 78, 2547 (1997).
  * Bertet et al. [2002] P. Bertet, A. Auffeves, P. Maioli, S. Osnaghi, T. Meunier, M. Brune, J. M. Raimond, and S. Haroche, _Direct measurement of the wigner function of a one-photon fock state in a cavity_ , Phys. Rev. Lett. 89, 200402 (2002).
  * Vlastakis et al. [2013] B. Vlastakis, G. Kirchmair, Z. Leghtas, S. E. Nigg, L. Frunzio, S. M. Girvin, M. Mirrahimi, M. H. Devoret, and R. J. Schoelkopf, _Deterministically Encoding Quantum Information Using 100-Photon Schrödinger Cat States_ , Science 342, 607 (2013).
  * Wang et al. [2016] C. Wang, Y. Y. Gao, P. Reinhold, R. W. Heeres, N. Ofek, K. Chou, C. Axline, M. Reagor, J. Blumoff, K. M. Sliwa, et al., _A Schrödinger cat living in two boxes_ , Science 352, 1087 (2016), eprint 1601.05505.
  * Cives-Esclop et al. [2000] A. Cives-Esclop, A. Luis, and L. L. Sánchez-Soto, _An eight-port detector with a local oscillator of finite intensity_ , Journal of Optics B: Quantum and Semiclassical Optics 2, 526 (2000).
  * Cova et al. [1987] S. Cova, G. Ripamonti, and A. Lacaita, _Avalanche semiconductor detector for single optical photons with a time resolution of 60 ps_ , Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 253, 482 (1987).
  * Lita et al. [2008] A. E. Lita, A. J. Miller, and S. W. Nam, _Counting near-infrared single-photons with 95% efficiency_ , Opt. Express 16, 3032 (2008).
  * Hangleiter and Eisert [2022] D. Hangleiter and J. Eisert, _Computational advantage of quantum random sampling_ , arXiv:2206.04079 (2022), eprint 2206.04079.
  * Chen et al. [2021] S. Chen, W. Yu, P. Zeng, and S. T. Flammia, _Robust shadow estimation_ , PRX Quantum 2, 030348 (2021).