# Classical shadow tomography for continuous variables quantum systems

Simon Becker, Nilanjana Datta, Ludovico Lami, and Cambyse Rouzé Simon Becker
is with the Department of Mathematics, ETH Zürich, Rämistrasse 101, 8092
Zürich, Switzerland. Email<EMAIL_ADDRESS>Datta is with
the Department of Applied Mathematics and Theoretical Physics, Centre for
Mathematical Sciences, University of Cambridge, Cambridge CB3 0WA, United
Kingdom. Email<EMAIL_ADDRESS>Lami is with the Korteweg–de
Vries Institute for Mathematics, the Institute for Theoretical Physics, and
QuSoft, University of Amsterdam, Science Park 123, 1098 XG Amsterdam. Part of
this work was conducted while he was at the Institute for Theoretical Physics,
University of Ulm, Albert-Einstein-Allee 11, D-89069 Ulm, Germany. Email:
<EMAIL_ADDRESS>Rouzé is with the Zentrum Mathematik,
Technische Universität München, 85748 Garching, Germany. Email:
<EMAIL_ADDRESS>
###### Abstract

In this article we develop a continuous variable (CV) shadow tomography scheme
with wide ranging applications in quantum optics. Our work is motivated by the
increasing experimental and technological relevance of CV systems in quantum
information, quantum communication, quantum sensing, quantum simulations,
quantum computing and error correction. We introduce two experimentally
realisable schemes for obtaining classical shadows of CV (possibly non-
Gaussian) quantum states using only randomised Gaussian unitaries and easily
implementable Gaussian measurements such as homodyne and heterodyne detection.
For both schemes, we show that
$N=\mathcal{O}\big{(}\operatorname{poly}\big{(}\frac{1}{\epsilon},\log\big{(}\frac{1}{\delta}\big{)},M_{n}^{r+\alpha},\log(m)\big{)}\big{)}$
samples of an unknown $m$-mode state $\rho$ suffice to learn the expected
value of any $r$-local polynomial in the canonical observables of degree
$\alpha$, both with high probability $1-\delta$ and accuracy $\epsilon$, as
long as the state $\rho$ has moments of order $n>\alpha$ bounded by $M_{n}$.
By simultaneously truncating states and operators in energy and phase space,
we are able to overcome new mathematical challenges that arise due to the
infinite-dimensionality of CV systems. We also provide a scheme to learn
nonlinear functionals of the state, such as entropies over any small number of
modes, by leveraging recent energy-constrained entropic continuity bounds.
Finally, we provide numerical evidence of the efficiency of our protocols in
the case of CV states of relevance in quantum information theory, including
ground states of quadratic Hamiltonians of many-body systems and cat qubit
states. We expect our scheme to provide good recovery in learning relevant
states of 2D materials and photonic crystals.

##  I Introduction

Obtaining classical descriptions of states of quantum-mechanical systems is a
fundamental ingredient of quantum computing. It is useful for storing and
transmitting quantum information and essential for verification and
benchmarking of quantum devices. However, the underlying quantum nature of
such systems provides a huge hurdle in obtaining such a classical description:
to learn anything about a quantum state one needs to measure it, but
measurements in quantum mechanics are inherently destructive and furthermore
probabilistic, entailing that individual measurement outcomes only give
limited information about the state of the system. Consequently, in order to
obtain a classical representation of a quantum state, one requires multiple
identical copies of the state on which successive, appropriate (possibly
adaptive) single-copy measurements can be performed.

This is what is done in the traditional method of learning an unknown quantum
state, known as quantum state tomography. It is the process of inferring a
quantum state by using suitable measurements on many identical copies of the
state. There is, however, a huge practical limitation in using quantum state
tomography for many-body quantum systems. This is due to the so-called “curse
of dimensionality”: the number of parameters needed to fully specify the state
of a quantum-mechanical system grows exponentially with the system size. The
exponential number of measurements needed to infer these parameters makes
quantum state tomography infeasible for large systems. Consequently, full
quantum state tomography has only been realised in systems with few
components, in particular, in a system of ten qubits, which too required
millions of measurements. In a nutshell, obtaining a classical description of
a $d$-dimensional quantum mixed state $\rho$, given many copies of it, via
quantum state tomography, can be shown to require $\Omega(d^{2})$ copies of
the state. (More recently, it has been shown [1, 2] that $\mathcal{O}(d^{2})$
copies also suffice.) However, this number grows exponentially with the number
$n$ of qubits (since $d=2^{n}$), and the problem becomes rapidly intractable.

In 2018 Aaronson [3] pointed out that for certain concrete tasks, obtaining a
complete classical characterisation of the quantum state is unnecessary.
Instead it is often sufficient to accurately predict many useful properties of
the state. This led him to propose a novel task called shadow tomography, the
aim of which is not to learn a complete description of the unknown quantum
state but instead to simultaneously estimate the outcome probabilities
associated with a list of $M$ two-outcome measurements, $E_{1},\ldots,E_{M}$,
performed on the state, up to a desired accuracy (say, $\varepsilon$). For an
unknown $d$-dimensional quantum mixed state $\rho$, this requires the
prediction of $M$ expectation values, $\operatorname{Tr}[E_{i}\rho]$
$i=1,2,\ldots,M$, to within an additive error $\varepsilon$. Aaronson showed
that, remarkably, the number of copies of the quantum state (i.e. the sample
size) needed to make these predictions, scales polynomially in the system
size. Moreover, to predict the $M$ different expectation values, only
$\mathcal{O}({\rm{polylog}}(M))$ number of copies of the state are needed. In
spite of this advantage with respect to the sample size, implementation of
shadow tomography is impractical because it requires very expensive quantum
processing, including exponentially long quantum circuits that act
collectively across all copies of the unknown quantum state stored in a
quantum memory, as well as a lot of storage and post-processing to make the
desired predictions. Aaronson’s work was followed by a spate of papers, with
the best result on the sample complexity of shadow tomography being obtained
in [4, 5].

In 2020 Huang et al. [6] improved on the work of Aaronson, by providing an
efficient and experimentally feasible method to learn an unknown quantum state
from just a few simple, single-copy measurements. The measurement outcomes are
used to construct a minimal classical representation of the state, called its
classical shadow, which can be efficiently stored on a classical computer.
This can thereafter be used to predict many linear (and possibly polynomial)
properties of the quantum state. For example, the properties could be the
expectation values of a list of $M$ observables in the given quantum state, as
was the goal of Aaronson [3]. However, in contrast with the shadow tomography
method proposed by Aaronson, Huang et al. [6] have a strict divide between the
quantum and classical parts of their protocol: after obtaining a classical
shadow of the quantum state, all processing necessary to predict its
properties are done via classical computations. They proved that, under
certain conditions, it is possible to predict expectation values of a list of
$M$ observables for an unknown quantum state with a small constant error, with
high success probability, by using $\log(M)$ number of copies of the state.
They called this method classical shadow tomography. This novel method
provides a tractable and rigorous procedure to obtain succinct classical
descriptions of quantum many-body states, using which many useful properties
of these states can be predicted.

Classical shadow tomography was originally developed for locally finite-
dimensional systems. In contrast, recent years have seen a fast growth of the
range of applications of infinite-dimensional, continuous variable (CV)
quantum systems, e.g. collections of electromagnetic modes travelling along an
optical fibre or massive harmonic oscillators, in all areas of quantum
information [7]. Notable applications include quantum communication [8, 9, 10,
11, 12, 13, 14], quantum sensing [15, 16, 17, 18], quantum simulations [19],
quantum computing and error correction [20, 21, 22, 23, 24], and have been
enabled by a steady development of non-classical sources of radiation [25, 26,
27, 28, 29, 30, 31, 32].

Hence, CV systems are of enormous technological and experimental relevance.
This comes hand in hand with a pressing need for fast and efficient quantum
state tomography of CV systems. The aim of this paper is precisely to devise a
rigorous procedure for obtaining classical shadows of states of CV quantum
systems, thus developing an efficient and experimentally feasible procedure to
benchmark CV quantum technologies.

The traditional way of performing tomography of CV systems consists of
measuring a functional characterizing the state, typically the characteristic
function [33, 34, 35], the Wigner [36] or the Husimi Q-function [37] from
which the quantum state of the system can be reconstructed, either by inverse
linear transformations or by statistical inference techniques [38, 39, 35, 40,
41, 42, 43, 44]. This approach has also been experimentally tested in various
settings [36, 38, 45, 46, 47, 48, 49, 50, 51, 52, 53, 22, 54, 55]. However, as
for their discrete analogues, a naive full tomography of a CV quantum system
in terms of a quasi-probability distribution is highly inefficient. To remedy
this issue, more advanced tomographic schemes were proposed which involve a
displacement of the state in phase space followed by parity or multiple photon
measurements [50, 52, 53, 48, 16, 56, 47, 57]. These measurements are then
processed to reconstruct the quasi-probability distributions or their
corresponding states in the Fock basis. On the down side, these schemes
exhibit a trade-off between the number of measurement points in phase space
and the number of operator expectation values measured at each point. More
recently, a more efficient method of reconstruction of quasi-probability
distributions via Lagrange interpolation was proposed in [58]. However, to the
best of our knowledge, a rigorous analysis of the sample and computational
complexity associated to each of these methods is missing. Moreover, the
latter were mostly applied to systems of a small number of modes. Here,
instead, we propose a new scheme for building classical shadows of multi-mode
continuous variables quantum systems. Our proposal comes with rigorous
complexity bounds.

Related works: Results analogous to ours recently appeared in a concurrent and
independent work by Gandhari et al. [59]. In it too a framework generalizing
the qubit-based classical shadow tomography protocol [6] to CV systems was
developed. A key step in their method is to express the density matrix of the
reconstructed state in terms of so-called pattern functions. The latter were
originally introduced in the context of optical homodyne tomography by
D’Ariano et al. [60] and have been used extensively thereafter in quantum
tomography of CV systems (see e.g. [61] and references therein). The authors
of [59] obtained bounds on the sample complexity for estimating quantum states
for their protocol by exploiting known bounds on pattern functions [61]. The
framework of [59] is equivalent to ours in the settings of homodyne and
heterodyne detection (see Section V-C for an explanation). Even though pattern
functions do not arise explicitly in our work, they are implicit in our
results. This can be seen by a comparison of our results with known
expressions [62] for the Fourier transform of pattern functions in terms of
Laguerre polynomials.

Other recent works on learning CV quantum states (and quantum processes)
include [63, 64].

###  I-A CV classical shadow tomography

Let us start with a brief summary of our extension of the protocol of Huang et
al. [6] for obtaining a classical shadow of the quantum state of a continuous
variables quantum system. Assume that multiple (say, $N$) identical copies of
an unknown quantum state, $\rho$, are available, and one has (i) an ensemble
$\cal{U}$ of unitary operators and (ii) a quantum measurement described by the
set of measurement operators $\\{M_{x}\\}_{x}$, satisfying
$\sum_{x}M_{x}^{\dagger}M_{x}=I$, such that elements $M_{x}U$,
$U\in\mathcal{U}$, describe a tomographically complete set of measurements. In
Huang et al. [6] $\rho$ was considered to be an $n$-qubit state, and
$\\{M_{x}\\}_{x}$ was a measurement in the computational basis, in which case
$x\in\\{0,1\\}^{n}$. In a practical scheme, each ensemble ${\cal{U}}$ should
be realisable as an efficient quantum circuit, and also have a succinct
classical description.

In analogy with the discrete setting, our proposal for a CV classical shadow
tomography contains the following three main ingredients:

  1. (i)

An $m$-mode CV system in an unknown state $\rho$.

  2. (ii)

A random variable $S$ taking values in the group $\operatorname{Sp}(2m)$ of
$2m\times 2m$ symplectic matrices, and the associated unitaries $U_{S}$. For
homodyne measurements we consider e.g. random variables distributed according
to the Haar measure on $\operatorname{Sp}\cap\operatorname{SO}.$

  3. (iii)

A Gaussian positive operator valued measurement
$\big{\\{}\pazocal{D}(x)\psi\pazocal{D}(-x)\,\frac{d^{2m}x}{(2\pi)}\big{\\}}_{x\in\mathbb{R}^{2m}}$,
where $\psi\coloneqq\ket{\psi}\\!\\!\bra{\psi}$ is a fixed pure Gaussian state
and $\pazocal{D}(x)$ denotes the unitary displacement operator of phase-space
parameter $x\in\mathbb{R}^{2m}$. In practice, we will consider two types of
Gaussian measurements, namely homodyne detection along a random direction in
phase space and heterodyne detection.

A CV classical shadow of the quantum state $\rho$ is then created by using a
randomised protocol that involves repeatedly performing the following simple
steps:

  1. (a)

A symplectic matrix $S$ is selected randomly from $\operatorname{Sp}(2m)$ and
applied to a copy of $\rho$, resulting in the unitary transformation
$\rho\mapsto U_{S}\rho U_{S}^{\dagger}$.

  2. (b)

The Gaussian measurement of effect operators $M_{x}$ is performed on the
output of the previous step, yielding the post-measurement state
$\widetilde{\rho}_{x}=\psi_{x}\coloneqq\pazocal{D}(x)\psi\pazocal{D}(-x)$,
when the measurement outcome is $x$. One can attempt to partly undo the effect
of the unitary by counter-rotating $\widetilde{\rho}_{x}$, implementing the
transformation $\widetilde{\rho}_{x}\mapsto
U_{S}^{\dagger}\widetilde{\rho}_{x}U_{S}$. Naturally, due to the measurement
such transformation will not yield back the original state $\rho$ which is in
general mixed.

  3. (c)

Averaging over $S,X$, where $X=\mathbb{R}^{2m}$ is the set of measurement
outcomes equipped the probability measure derived from the measurement
outcomes, the counter-rotated state yields a quantum channel $\mathcal{M}$
that depends on the measure $\mu$:

$\displaystyle\mathcal{M}(\rho)\coloneqq{\mathbb{E}}_{S,X}\left[U_{\\!S}^{{\dagger}}\,\widetilde{\rho}_{X}U_{\\!S}^{\vphantom{{\dagger}}}\right].$
(1)

For our choice of Gaussian measurements, the effective quantum channel
$\mathcal{M}$ is a simple linear bosonic channel whose action can be
represented as a random displacement in phase space.

  4. (d)

Heuristically, one would like to define the classical shadow of $\rho$ as the
random operator

$\displaystyle\hat{\rho}_{S,x}\coloneqq$ $\displaystyle\
\mathcal{M}^{-1}\\!\left(U_{\\!S}^{{\dagger}}\widetilde{\rho}_{x}U_{\\!S}^{\vphantom{{\dagger}}}\right).$
(2)

Once the nature of the protocol, and in particular the Gaussian measurement
$\big{\\{}\pazocal{D}(x)\psi\pazocal{D}(-x)\,\frac{d^{2m}x}{(2\pi)}\big{\\}}_{x\in\mathbb{R}^{2m}}$,
has been specified, the classical shadow is simply a function of the
particular realisations $S,x$ obtained in this round of the protocol.
Importantly, it is _not_ a function of the unknown state $\rho$. The
information on $\rho$ is at this point stored _only_ in the probability
distribution associated with the measurement outcome $x$. As a matter of fact,
$\hat{\rho}_{S,x}$ is an unbiased estimator of $\rho$:

$\displaystyle\mathbb{E}_{S,X}\big{[}\hat{\rho}_{S,X}\big{]}$
$\displaystyle=\mathcal{M}^{-1}\left(\mathbb{E}_{S,X}\big{[}U_{\\!S}^{{\dagger}}\,\widetilde{\rho}_{X}U_{\\!S}^{\vphantom{{\dagger}}}\big{]}\right)$
(3)
$\displaystyle=\big{(}\mathcal{M}^{-1}\\!\\!\circ\\!\mathcal{M}\big{)}(\rho)=\rho\,.$

This is the final output of one iteration of the protocol. A purely classical
description of the operator $\hat{\rho}_{S,x}$ — which, we stress again, does
not depend on the unknown state $\rho$ — is stored in a classical memory for
future processing. Note that such a description can be produced given $S$ and
$x$ alone, by simply computing the operator defined by (2).

Let us summarise the whole procedure. The measurement yields
$\widetilde{\rho}_{x}$ as an approximation for the state
$U_{S}\rho\,U_{S}^{\dagger}$; then one attempts to find an approximation for
the original state, $\rho$, by counter-rotating $\widetilde{\rho}_{x}$ to
obtain ${U_{S}^{\dagger}}\widetilde{\rho}_{x}{U_{S}}$. Finally, by applying
$\mathcal{M}^{-1}$ one eliminates the effect of the average noise (represented
by the quantum channel $\mathcal{M}$) that the state undergoes in the
protocol. This yields the classical shadow $\hat{\rho}_{S,X}$ whose
description is stored as a string of classical data in a classical memory for
future processing. Repeating the above protocol on $N$ independent identical
copies of $\rho$ yields a string of classical shadows
$\\{\hat{\rho}^{(1)},\ldots,\hat{\rho}^{(N)}\\}$, where for simplicity we
introduced the shorthand notation
$\hat{\rho}^{(i)}\coloneqq\hat{\rho}_{S_{i},x_{i}}$. This string can be
further processed for instance by computing the empirical average

$\displaystyle\hat{\sigma}^{(N)}\coloneqq\frac{1}{N}\sum_{i=1}^{N}\hat{\rho}^{(i)}\,.$

The advantage of our CV quantum tomography protocol is most clearly visible
from the simple structure of the map $\mathcal{M}$ and its inverse at the
level of characteristic functions: for any trace class operator $Z$ with
characteristic function
$\chi_{Z}:u\in\mathbb{R}^{2m}\mapsto\operatorname{Tr}(Z\pazocal{D}(u))$ and
$|\psi\rangle\coloneqq U_{T}^{\dagger}|0\rangle$, for any
$T\in\operatorname{Sp}(2m)$,

$\displaystyle\chi_{\mathcal{M}(Z)}=\chi_{Z}\,f_{\mu,T}\,,\quad\text{ where }$
(4) $\displaystyle f_{\mu,T}(u)\coloneqq\
\int\\!\\!d\mu(S)\,e^{-\frac{1}{2}\,\|TSu\|^{2}}$ (5)

depends on the Gaussian state $\ket{\psi}$. Therefore, at least formally the
classical shadows $\hat{\rho}^{(i)}$ can be equivalently reconstructed by
their characteristic functions, which take the form

$\chi_{\hat{\rho}^{(i)}}=\chi_{U_{S_{i}}^{\dagger}\widetilde{\rho}_{x_{i}}U_{S_{i}}^{\vphantom{{\dagger}}}}\,f_{\mu,T}^{-1}\,.$
(6)

###  I-B Moment constraints

When trying to implement the above strategy, one faces however two issues that
are specific to the continuous variables setting. First, as we will see, the
quantum channel $\mathcal{M}$ is in general not surjective on the space of
trace class operators, which implies that the classical shadow
$\hat{\rho}_{S,X}$ will typically not be a trace class — and, for that matter,
not even a bounded — operator. In other words, the right-hand side of (6),
although well defined as a function on $\mathds{R}^{2m}$, is not the
characteristic function of a quantum state. Second, the parameters $S$ and $X$
that need to be stored in the classical memory at each round of the protocol
belong to continuous and unbounded sets.

To overcome the first problem, we focus on a suitable characterisation of the
classical shadow $\hat{\rho}_{S,X}$ that is well defined, namely, its
characteristic function $\chi_{\hat{\rho}_{S,X}}$, and construct the operator
$\hat{\rho}_{S,X}$ itself as defined only on a restricted domain. To obtain
effective estimates of $\rho$ from the above scheme, we will further need to
assume that $\rho$ has controlled moments of low degree. Such an assumption
allows us to show that the projection of the state $\rho$ onto the finite
subspace of Hilbert–Schmidt operators supported on the space of low energy
Fock states is sufficient for obtaining a good enough approximation of the
state via the Gaussian shadow tomography protocol.

Let us make these considerations more precise. We consider the maps

$\displaystyle\mathcal{P}_{M}(T)\coloneqq\sum_{\mathbf{n}_{1},\mathbf{n}_{2}\in\\{0,..,M\\}^{m}}\langle\mathbf{n}_{1}|T|\mathbf{n}_{2}\rangle\,\ket{\mathbf{n}_{1}}\\!\\!\bra{\mathbf{n}_{2}}\,,$
(7)

where $\ket{\mathbf{n}}$ is a multivariate Fock state, and
$\mathbf{n}\in\\{0,...,M\\}^{m}$ is a multi-index. Using a simple
approximation scheme, we also approximate the Schwartz operators
$\ket{\mathbf{n}_{1}}\\!\\!\bra{\mathbf{n}_{2}}$ by Schwartz operators
$\widetilde{Z}_{\mathbf{n}_{1}\mathbf{n}_{2}}$ with smooth compactly supported
characteristic functions and define an auxiliary map

$\widetilde{\mathcal{P}}_{M}(T)\coloneqq\sum_{\mathbf{n}_{1},\mathbf{n}_{2}\in\\{0,..,M\\}^{m}}\operatorname{Tr}(\widetilde{Z}_{\mathbf{n}_{2}\mathbf{n}_{1}}T)\,\ket{\mathbf{n}_{1}}\\!\\!\bra{\mathbf{n}_{2}}\,.$
(8)

Our goal is to approximate possibly unbounded observables of an $m$-mode CV
quantum system. For this, we introduce a norm which captures approximations of
such operators: fix an arbitrary positive integer $n$; then given $0<\alpha<n$
and two states $\rho,\rho^{\prime}$ with
$\operatorname{Tr}(\rho(I+N_{m})^{n}),\operatorname{Tr}(\rho^{\prime}(I+N_{m})^{n})<\infty$,
where $N_{m}$ stands for the $m$-mode number operator, we denote
$H_{m}\coloneqq I+N_{m}$ and

$\displaystyle\|X\|^{(\alpha)}_{1}\coloneqq\left\|H_{m}^{\alpha/2}XH_{m}^{\alpha/2}\right\|_{1}\,$

for any trace-class $X$ or which the right-hand side is finite. Then the
following bounds hold which put the above approximation scheme on a rigorous
footing:

######  Proposition 1.

Let $N_{m}$ be the number operator over $m$ modes and let $\rho$ be an
$m$-mode state such that $\operatorname{Tr}(N_{m}^{n}\rho)\eqqcolon
mE<\infty$. Then, for any $\eta^{2}\geq 2M^{2}$,

$\displaystyle\|\rho-\mathcal{P}_{M}(\rho)\|_{1}^{(\alpha)}\leq
2(1+M)^{\frac{\alpha-n}{2}}E\,,$
$\displaystyle\|\rho-\widetilde{\mathcal{P}}_{M}(\rho)\|^{(\alpha)}_{1}\leq
2(1+M)^{\frac{\alpha-n}{2}}E$ $\displaystyle\hskip
80.0831pt+(3mM)^{6mM+\alpha}e^{-\frac{\eta^{2}m}{2}}\eta^{2Mm}.$

##  II Main results

###  II-A CV classical shadows via homodyne detection

We first consider the scenario in which one performs a homodyne detection
along a random direction in phase space. More precisely, $m$ independent
random matrices $S_{1},\dots,S_{m}$ are distributed uniformly (according to
the Haar measure) on the intersection
$\operatorname{Sp}(2)\cap\operatorname{SO}(2)$ between the symplectic and the
special orthogonal group, with corresponding angles
$\theta_{1},\cdots,\theta_{m}\in[-\pi,\pi]$. The homodyne measurement is then
performed along the position axis and yields the classical outcome
$x=(x_{1},...,x_{m})\in\mathbb{R}^{2m}$. As in the qubit setting, we now
restrict ourselves to an arbitrary subset $A$ of $|A|\leq r$ modes. In that
case, the characteristic functions of the reduced shadows defined in (6) are
distributions of the form

$\displaystyle\chi_{\hat{\rho}_{A}^{(i)}}(u_{A})=\prod_{j\in
A}\sqrt{2}\pi\|u_{j}\|\,\delta\left((S_{j}u_{j})_{2}\right)e^{-iu_{j}^{\intercal}\Omega
S_{j}x_{j}}\,,$

for any $u_{A}\coloneqq\\{u_{j}\\}_{j\in A}\in\mathbb{R}^{2|A|}$, where
$\delta(x)$ is the Dirac distribution at $x\in\mathbb{R}$, where
$(S_{j}u_{j})_{2}\coloneqq\sum_{k}(S_{j})_{2k}(u_{j})_{k}$, and where $\Omega$
stands for the canonical symplectic form as defined in Equation (15). As
foreseen in the previous paragraph, this characteristic function is not square
integrable, and therefore cannot be associated to a quantum state. Indeed, we
recall that, by Plancherel’s theorem, any trace class operator, and thereby
any quantum state, gives rise to a square integrable characteristic function,
and for any two such operators $T_{1},T_{2}$,

$\displaystyle\operatorname{Tr}[T_{1}^{\dagger}T_{2}]=\int\frac{d^{2|A|}x}{(2\pi)^{|A|}}\,\chi_{T_{1}}(x)^{*}\chi_{T_{2}}(x)\,\,\equiv\braket{\chi_{T_{1}},\,\chi_{T_{2}}}\,.$
(9)

Instead, we construct a random matrix $\hat{\rho}^{(i)}_{A}(M)$ in the set
$\mathbb{M}_{(M+1)^{|A|}}(\mathbb{C})$ of matrices of size $(M+1)^{|A|}$ by
simply extending (9):

$\displaystyle\hat{\rho}^{(i)}_{A}(M)\coloneqq\sum_{\mathbf{n}_{1},\mathbf{n}_{2}\in\\{0,...,M\\}^{A}}\,\langle\chi_{\ket{\mathbf{n}_{1}}\\!\bra{\mathbf{n}_{2}}},\chi_{\hat{\rho}^{(i)}_{A}}\rangle\,\ket{\mathbf{n}_{1}}\\!\\!\bra{\mathbf{n}_{2}}\,.$

By (3), we know that the expected matrix $\mathbb{E}[\hat{\rho}^{(i)}_{A}(M)]$
coincides with the unknown state $\rho$. Moreover, by Bernstein’s matrix
concentration inequality, we have that, with high probability the empirical
average
$\sigma_{A}^{(N)}(M)\coloneqq\frac{1}{N}\sum_{i=1}^{N}\hat{\rho}^{(i)}_{A}(M)$
will well-approximate $\mathcal{P}_{M}(\rho_{A})$ on all sets $A$ of size
$|A|\leq r$ for
$N=\mathcal{O}\big{(}\operatorname{poly}\big{(}e^{r},\,\log(m)\big{)}\big{)}$.
Combining this with the approximation bounds under moment constraints derived
in Proposition (1) leads us to our first main result:

######  Theorem 2 ((CV classical shadows via homodyne detection)).

With the notation introduced above, given $0\leq\alpha<n$ and assuming that
$E^{(n)}_{r}\coloneqq\max_{|A|\leq
r}\operatorname{Tr}(\rho_{A}H_{r}^{n})<\infty$ for all $A$ of size $|A|\leq
r$, we have that for integer
$M=\Big{\lceil}\Big{(}\frac{4E_{r}^{(n)}}{\epsilon}\Big{)}^{\frac{2}{n-\alpha}}\Big{\rceil}$,
$N=\mathcal{O}\Big{(}\operatorname{poly}\Big{(}\frac{1}{\epsilon^{2}},\,M^{r+\alpha},\log\big{(}\frac{1}{\delta}\big{)},\log(m)\Big{)}\Big{)}$
and any region $A$ of size $|A|\leq r$,

$\left\|\sigma^{(N)}_{A}(M)-\rho_{A}\right\|^{(\alpha)}_{1}\leq\epsilon$ (10)

with probability at least $1-\delta$. Similarly, for
$N=\mathcal{O}\Big{(}\operatorname{poly}\Big{(}\frac{1}{\epsilon^{2}},\,M^{r+\alpha},\log\big{(}\frac{1}{\delta}\big{)},\log(L)\Big{)}\Big{)}$
we have that for any set of $L$ observables $O_{j}$ on regions $A_{j}$ of size
at most $r$ and with
$\|H_{r}^{-\frac{\alpha}{2}}O_{j}H_{r}^{-\frac{\alpha}{2}}\|_{\infty}\leq 1$,

$\displaystyle\max_{j}\big{|}\operatorname{Tr}[O_{j}\,(\sigma_{A}^{(N)}(M)-\rho_{A_{j}})]\big{|}\leq\epsilon$

with probability at least $1-\delta$.

###  II-B CV classical shadows via heterodyne detection

We also consider the case of a heterodyne detection

$\displaystyle\left\\{\frac{1}{(2\pi)^{m/2}}\ket{x}\\!\\!\bra{x}\right\\}_{x\in\mathds{R}^{2m}}\,,$
(11)

i.e. $T=I$, and all unitaries employed are passive, i.e. such that
$\big{[}U_{\\!S}^{\vphantom{{\dagger}}},\frac{1}{2}R^{\intercal}R\big{]}=0$.
In that case, we show that the function defined in (5) takes the simpler form

$f_{\mu,T}(u)=e^{-\frac{1}{2}\|u\|^{2}}.$

Therefore, the classical shadow will have improper characteristic function

$\displaystyle\chi_{\hat{\rho}^{(i)}_{A}}(u_{A})=e^{\frac{1}{4}\|u_{A}\|^{2}}\prod_{j\in
A}e^{-iu_{j}^{\intercal}\Omega x_{j}}\,.$

In this case, we consider the matrices

$\displaystyle\hat{\rho}_{A}^{(i)}(M)\coloneqq\sum_{\mathbf{n}_{1},\mathbf{n}_{2}\in\\{0,...,M\\}^{A}}\,\langle\chi_{\widetilde{Z}_{\mathbf{n}_{1}\mathbf{n}_{2}}},\chi_{\hat{\rho}^{(i)}_{A}}\rangle\,\ket{\mathbf{n}_{1}}\\!\\!\bra{\mathbf{n}_{2}}\,.$

From here, repeating essentially the same argument as in the case of homodyne
detection, defining the matrix
$\widetilde{\sigma}_{A}^{(N)}(M)\coloneqq\frac{1}{N}\sum_{i=1}^{N}\hat{\rho}^{(i)}_{A}(M)$,
we arrive at our second main result:

######  Theorem 3 ((CV classical shadows via heterodyne detection)).

With the above notation, given $0\leq\alpha<n$ and assuming that
$E^{(n)}_{r}\coloneqq\max_{|A|\leq
r}\operatorname{Tr}(\rho_{A}H_{r}^{n})<\infty$ for all $A$ of size $|A|\leq
r$, we have that for
$N=\mathcal{O}\Big{(}\operatorname{poly}\Big{(}\frac{1}{\epsilon^{2}},\,(E_{r}^{(n)})^{r+\alpha},\log\big{(}\frac{1}{\delta}\big{)},\log(m)\Big{)}\Big{)}$
and any region $A$ of size $|A|\leq r$,

$\displaystyle\|\widetilde{\sigma}^{(N)}_{A}(M)-\rho_{A}\|^{(\alpha)}_{1}\leq\epsilon$

with probability at least $1-\delta$. Similarly, for
$N=\mathcal{O}\Big{(}\operatorname{poly}\Big{(}\frac{1}{\epsilon^{2}},\,(E_{r}^{(n)})^{r+\alpha},\log\big{(}\frac{1}{\delta}\big{)},\log(L)\Big{)}\Big{)}$
we have that for any set of $L$ observables $O_{j}$ on regions $A_{j}$ of size
at most $r$ and with
$\|H_{r}^{-\frac{\alpha}{2}}O_{j}H_{r}^{-\frac{\alpha}{2}}\|_{\infty}\leq 1$,

$\displaystyle\max_{j}\big{|}\operatorname{Tr}[O_{j}\,(\widetilde{\sigma}_{A}^{(N)}(M)-\rho_{A_{j}})]\big{|}\leq\epsilon$

with probability at least $1-\delta$.

##  III Notation and basic notions

###  III-A Operators and norms

Given a separable Hilbert space ${\pazocal{H}}$, we denote by
${\pazocal{B}}({\pazocal{H}})$ the space of bounded linear operators on
${\pazocal{H}}$, and by ${\pazocal{T}}_{p}({\pazocal{H}})$ the Schatten
$p$-class, which is the Banach subspace of ${\pazocal{B}}({\pazocal{H}})$
formed by all bounded linear operators whose Schatten $p$-norm, defined as
$\|X\|_{p}=\left(\operatorname{Tr}|X|^{p}\right)^{1/p}$, is finite.
Henceforth, we refer to ${\pazocal{T}}_{1}({\pazocal{H}})$ as the set of trace
class operators. The set of quantum states (or density matrices), i.e.
positive semi-definite operators $\rho\in{\pazocal{T}}_{1}({\pazocal{H}})$ of
unit trace, is denoted by ${\pazocal{D}}({\pazocal{H}})$. The Schatten
$1$-norm, $\|\cdot\|_{1}$, is the trace norm, and the corresponding induced
distance (e.g. between quantum states) is the trace distance. Note that the
Schatten $2$-norm, $\|\cdot\|_{2}$, coincides with the Hilbert–Schmidt norm.
We also recall that a quantum channel with input system $A$ and output system
$B$ is any completely positive, trace-preserving (CPTP) linear map
$\mathcal{N}:{\pazocal{T}}_{1}({\pazocal{H}}_{A})\to{\pazocal{T}}_{1}({\pazocal{H}}_{B})$,
where ${\pazocal{H}}_{A},{\pazocal{H}}_{B}$ are the Hilbert spaces
corresponding to $A,B$, respectively.

If $\rho$ is a quantum state with spectral decomposition
$\rho=\sum_{i}p_{i}\ket{\phi_{i}}\\!\\!\bra{\phi_{i}}$, and $A$ is a positive
semi-definite operator, the expected value of $A$ on $\rho$ is defined as

$\displaystyle\operatorname{Tr}[\rho
A]\coloneqq\sum_{i:\,p_{i}>0}p_{i}\big{\|}A^{1/2}|{\phi_{i}}\rangle\big{\|}^{2}\in\mathbb{R}_{+}\cup\\{+\infty\\}\,;$
(12)

here we use the convention that $\operatorname{Tr}[\rho A]=+\infty$ if the
above series diverges or if there exists an index $i$ for which $p_{i}>0$ and
$\ket{\phi_{i}}\notin\operatorname{dom}\left(A^{1/2}\right)$. This definition
can be extended to a generic densely defined self-adjoint operator $A$ on
${\pazocal{H}}$, by considering its decomposition $A=A_{+}-A_{-}$ into
positive and negative parts, with $A_{\pm}$ being positive semi-definite
operators with mutually orthogonal supports. The operator $A$ is said to have
a finite expected value on $\rho$ if (i)
$\ket{\phi_{i}}\in\operatorname{dom}\big{(}A_{+}^{1/2}\big{)}\cap\operatorname{dom}\big{(}A_{-}^{1/2}\big{)}$
for all $i$ for which $p_{i}>0$, and (ii) the two series
$\sum_{i}p_{i}\big{\|}A_{\pm}^{1/2}\ket{\phi_{i}}\big{\|}^{2}$ both converge.
In this case, the following quantity is called the expected value of $A$ on
$\rho$:

$\displaystyle\operatorname{Tr}[\rho
A]\coloneqq\sum_{i:\,p_{i}>0}p_{i}\big{\|}A_{+}^{1/2}\ket{\phi_{i}}\big{\|}^{2}-\sum_{i:\,p_{i}>0}p_{i}\big{\|}A_{-}^{1/2}\ket{\phi_{i}}\big{\|}^{2}\,.$
(13)

Obviously, for a pair of operators $A,B$ satisfying $A\geq B$, we have that
$\operatorname{Tr}[\rho A]\geq\operatorname{Tr}[\rho B]$.

Let $A$ be an (unbounded) operator $A$ on some Banach space ${{X}}$, with
domain $\operatorname{dom}(A)$. Such an operator is called closed if its
graph, that is
$\left\\{(\ket{x},A\ket{x});\ket{x}\in\operatorname{dom}(A)\right\\}\subset
X\times X,$ is closed. The spectrum of a closed operator $A$ is defined as the
set [65, Definition 9.16]

$\operatorname{sp}(A)\coloneqq\left\\{\lambda\in\mathbb{C}:\,\lambda I-A\text{
is not continuously invertible}\right\\}.$

Henceforth, we often suppress the identity operator $I$ in the expression
$(\lambda I-A)$ for notational simplicity. We remind the reader that the
spectrum of a self-adjoint positive operator is a closed subset of the
positive real half-line [65, Proposition 9.20]. Given a possibly unbounded
operator $X$, $L_{X}$ stands for the left multiplication by $X$:
$L_{X}(Y)=XY$, whereas $R_{X}$ stands for right multiplication by $X$:
$R_{X}(Y)=YX$, whenever these products are well-defined.

###  III-B Continuous variable quantum systems

A CV system with $m$ modes is defined on the Hilbert space
${\pazocal{H}}_{m}\coloneqq L^{2}(\mathds{R}^{m})$, equipped with the multi-
mode Fock basis $\\{|\mathbf{n}\rangle\\}_{\mathbf{n}\in\mathbb{N}^{m}}$ of
eigenvectors of the number operator $N_{m}$:

$\displaystyle
N_{m}|\mathbf{n}\rangle=\Big{(}\sum\nolimits_{j=1}^{m}n_{j}\Big{)}\,|\mathbf{n}\rangle\,.$
(14)

We denote the canonical operators on each mode as $X_{j},P_{j}$
($j=1,\ldots,m$). Define the formal vector

$\displaystyle
R\coloneqq\begin{pmatrix}X_{1},\ldots,X_{n},P_{1},\ldots,P_{n}\end{pmatrix}^{\intercal}$

and the symplectic form

$\displaystyle\Omega\coloneqq\begin{pmatrix}0&\mathds{1}\\\
-\mathds{1}&0\end{pmatrix}$ (15)

(where all blocks are $m\times m$ matrices), in terms of which the canonical
commutation relations read (at least when evaluated on Schwartz functions)

$\displaystyle[R_{j},R_{k}]=i\Omega_{jk}\,.$ (16)

We also introduce also the annihilation and creation operators
$a_{j}^{\dagger},a_{j}$ ($j=1,\ldots,m$), defined by

$\displaystyle a_{j}\coloneqq\frac{X_{j}+iP_{j}}{\sqrt{2}}\,,\qquad
a_{j}^{\dagger}\coloneqq\frac{X_{j}-iP_{j}}{\sqrt{2}}\,.$ (17)

In terms of these operators, the single-mode Fock states can be constructed as

$\displaystyle\ket{n}\coloneqq\frac{(a^{\dagger})^{n}}{\sqrt{n!}}\ket{0}\,,$
(18)

with similar formulae holding for the multi-mode case. The canonical
commutation relations can also be written as
$[a_{j},a_{k}^{\dagger}]=\delta_{jk}$.

With a slight abuse of notation, we will often denote by the same symbol
$\Omega$ the symplectic form for different sets of modes. The quantum
covariance matrix $V[\rho]$ and mean vector $t[\rho]$ associated with a
generic state $\rho$ are defined by

$\displaystyle t[\rho]_{j}\coloneqq\operatorname{Tr}\left(\rho
R_{j}\right)\,,\quad
V[\rho]_{jk}\coloneqq\operatorname{Tr}\left(\rho\left\\{R_{j}-t_{j},\,R_{k}-t_{k}\right\\}\right)\,,$
(19)

provided that these expressions are well defined. For an arbitrary
$x\in\mathds{R}^{2m}$, we define the associated displacement operator by

$\displaystyle\pazocal{D}(x)\coloneqq e^{-ix^{\intercal}\Omega R}\,.$ (20)

Note that $\pazocal{D}(x)^{\dagger}=\pazocal{D}(x)^{-1}=\pazocal{D}(-x)$. By
writing $x=x^{\prime}\oplus x^{\prime\prime}$, where
$x^{\prime}\in\mathds{R}^{m}$ groups together the first $m$ components of $x$
and $x^{\prime\prime}\in\mathds{R}^{m}$ the last $m$, one can also introduce
the complex vector

$\displaystyle\alpha(x)\coloneqq\frac{1}{\sqrt{2}}\left(x^{\prime}+ix^{\prime\prime}\right),$
(21)

in terms of which we have that [66, Eq. (3.3.30)–(3.3.31)]

$\displaystyle\pazocal{D}(x)$
$\displaystyle=\exp\left[\sum\nolimits_{j=1}^{m}\left(\alpha_{j}(x)a_{j}^{\dagger}-\alpha_{j}(x)^{*}a_{j}\right)\right]$
(22)
$\displaystyle=e^{-\frac{1}{4}\|x\|^{2}}e^{\sum_{j}\alpha_{j}(x)a_{j}^{\dagger}}e^{-\sum_{j}\alpha_{j}(x)^{*}a_{j}},$
(23)

where $\|x\|^{2}\coloneqq\sum_{j}x_{j}^{2}$. In dealing with continuous
variable systems, one can stick to the real notation, employing real vectors
$x$, or move to the complex one, which uses the complex vectors $\alpha(x)$.
In this paper we will mostly follow the former convention; however, it will be
useful, occasionally, to use the latter too. In general, a prompt translation
between one set of conventions and the other can be obtained by means of (21).

Another important identity involving displacement operators is

$\displaystyle\int\frac{d^{2m}x}{(2\pi)^{m}}\,\pazocal{D}(x)Z\pazocal{D}(-x)=I\operatorname{Tr}Z\,,$
(24)

valid for all trace class operators
$Z\in{\pazocal{T}}_{1}({\pazocal{H}}_{m})$, with the integral on the left-hand
side converging in the weak sense (see [67, Proposition 3.5.1]). Using
displacement operators, we can re-write (16) in Weyl form as

$\displaystyle\pazocal{D}(x+y)=e^{\frac{i}{2}x^{\intercal}\Omega
y}\,\pazocal{D}(x)\pazocal{D}(y)\,.$ (25)

Coherent states are instead defined as

$\displaystyle\ket{x}\coloneqq\pazocal{D}(x)\ket{0}\,,$ (26)

where $\ket{0}$ is the vacuum state. We can decompose $\ket{x}$ mode-wise as
$\ket{x}=\bigotimes_{j=1}^{m}\ket{x^{(j)}}$, where
$x^{(j)}\coloneqq(x_{j},x_{m+j})^{\intercal}\in\mathds{R}^{2}$ is the sub-
vector of $x$ obtained by picking the coordinates corresponding to the
$j^{\text{th}}$ coordinate and momentum, and $\ket{x^{(j)}}$ is a single-mode
coherent state. In general, the latter can be represented in the single-mode
Fock basis according to the identity

$\displaystyle\ket{y}$
$\displaystyle=e^{-\|y\|^{2}/4}\sum_{n=0}^{\infty}\frac{(y_{1}+iy_{2})^{n}}{\sqrt{2^{n}\,n!}}\ket{n}$
(27)
$\displaystyle=e^{-\|y\|^{2}/4}\sum_{n=0}^{\infty}\frac{\alpha(y)^{n}}{\sqrt{n!}}\ket{n}\,,$

where $\alpha(y)\in\mathds{C}$ is defined in (21).

For an arbitrary trace class operator $Z$, we can construct its characteristic
function $\chi_{Z}:\mathds{R}^{2m}\to\mathds{C}$ by

$\displaystyle\chi_{Z}(x)\coloneqq\operatorname{Tr}\big{[}Z\,\pazocal{D}(x)\big{]}\,.$
(28)

Characteristic functions are always bounded and furthermore continuous,
because of the strong operator continuity of the mapping
$x\mapsto\pazocal{D}(x)$. Moreover, in the sense of weak operator convergence
it holds that

$\displaystyle
Z=\int\frac{d^{2m}x}{(2\pi)^{m}}\,\chi_{Z}(x)\,\pazocal{D}(-x)\,.$ (29)

By applying (25) and (27), one can prove that

$\displaystyle\chi_{\ket{x}\\!\bra{y}}(u)$ (30)
$\displaystyle=\operatorname{Tr}[\ket{x}\\!\bra{y}\pazocal{D}(u)]$
$\displaystyle=\braket{y}{\pazocal{D}(u)}{x}$
$\displaystyle=\braket{0}{\pazocal{D}(-y)\pazocal{D}(u)\pazocal{D}(x)}{0}$
$\displaystyle=e^{-\frac{i}{2}u^{\intercal}\Omega
x}\braket{0}{\pazocal{D}(-y)\pazocal{D}(u+x)}{0}$
$\displaystyle=e^{-\frac{i}{2}u^{\intercal}\Omega
x}e^{\frac{i}{2}y^{\intercal}\Omega(u+x)}\braket{0}{\pazocal{D}(u+x-y)}{0}$
$\displaystyle=e^{-\frac{i}{2}u^{\intercal}\Omega
x}e^{\frac{i}{2}y^{\intercal}\Omega(u+x)}e^{-\frac{1}{4}\|u+x-y\|^{2}}.$

As special cases, we conclude e.g. that

$\displaystyle\chi_{\ket{x}\\!\bra{x}}(u)$
$\displaystyle=e^{-\frac{1}{4}\|u\|^{2}-iu^{\intercal}\Omega x},$ (31)
$\displaystyle\chi_{\ket{-x}\\!\bra{x}}(u)$
$\displaystyle=e^{-\frac{1}{4}\|u-2x\|^{2}},$
$\displaystyle\chi_{\ket{x}\\!\bra{-x}}(u)$
$\displaystyle=e^{-\frac{1}{4}\|u+2x\|^{2}}.$

Given two single-mode Fock states $|k\rangle,|j\rangle$ [68] and
$x=x^{\prime}\oplus x^{\prime\prime}$:

$\displaystyle\chi_{|k\rangle\langle j|}(x)=\left\\{\begin{aligned}
&\sqrt{\frac{k!}{j!}}e^{-\frac{\pi}{2}|\omega|^{2}}(\sqrt{\pi}\omega)^{j-k}L_{k}^{(j-k)}(\pi|\omega|^{2}),~{}~{}j\geq
k\\\
&\sqrt{\frac{k!}{j!}}e^{-\frac{\pi}{2}|\omega|^{2}}(-\sqrt{\pi}\omega)^{k-j}L_{j}^{(k-j)}(\pi|\omega|^{2}),~{}~{}j\leq
k\,,\end{aligned}\right.$ (32)

where $\omega\coloneqq-x^{\prime}+i\frac{x^{\prime\prime}}{2\pi}$. Above, the
functions $L_{k}^{(j)}$ are the Laguerre polynomials, defined for any two
integers $k,j$ as

$\displaystyle
L_{k}^{(j)}(x)\coloneqq\sum_{l=0}^{k}\frac{(k+j)!}{(k-l)!(j+l)!}\,\frac{(-x)^{l}}{l!}\,.$
(33)

Interestingly, the correspondence between trace class operators and
characteristic functions is injective — even more strikingly, it can be
extended to an isometry between the space of Hilbert–Schmidt operators and
that of square integrable functions $\mathds{R}^{2m}\to\mathds{C}$ [67,
Theorem 5.3.3]. A consequence of the existence of this isometry is the quantum
Plancherel theorem, which tells us that for any two trace class operators
$Z,Z^{\prime}\in{\pazocal{T}}_{1}({\pazocal{H}}_{m})$,

$\displaystyle\operatorname{Tr}\big{[}Z^{\dagger}Z^{\prime}\big{]}=\int\frac{d^{2m}x}{(2\pi)^{m}}\,\chi_{Z}(x)^{*}\chi_{Z^{\prime}}(x)\,.$
(34)

The canonical commutation relations are invariant under so-called symplectic
unitaries, constructed as follows. A $2m\times 2m$ real matrix $S$ such that
$S\Omega S^{\intercal}=\Omega$ (or equivalently $S\Omega=\Omega
S^{-\intercal}$) is called a symplectic matrix. From the defining relation it
can be immediately seen that any symplectic matrix must satisfy $\det S=\pm
1$; however, remarkably, it turns out that in fact all symplectic matrices
have determinant $1$. To any symplectic matrix we can associate a symplectic
unitary $U_{S}$ acting on ${\pazocal{H}}_{m}$. This is defined by either of
the following relations

$\displaystyle U_{S}^{\dagger}RU_{S}^{\phantom{{\dagger}}}=SR\,,\qquad
U_{S}^{\phantom{{\dagger}}}\pazocal{D}(x)U_{S}^{\dagger}=\pazocal{D}(Sx)\,,$
(35)

where the first identity is to be understood coordinatewise:
$(SR)_{j}=\sum_{k}S_{jk}R_{k}$. Note that symplectic matrices form a group,
denoted as $\operatorname{Sp}(2m)$, and that the correspondence $S\mapsto
U_{S}$ is a group homomorphism. In particular,

$\displaystyle
U_{S}^{\dagger}=U_{S}^{-1}=U_{S^{-1}}=U_{\Omega\,S^{\text{\raisebox{0.6pt}{$\intercal$}}}\Omega^{\text{\raisebox{0.6pt}{$\intercal$}}}}.$
(36)

Also, from (35) we deduce that

$\displaystyle
t\left[U_{S}^{\phantom{{\dagger}}}\rho\,U_{S}^{\dagger}\right]=S\,t[\rho]\,,\qquad
V\left[U_{S}^{\phantom{{\dagger}}}\rho\,U_{S}^{\dagger}\right]=S\,V[\rho]S^{\intercal}\,,$
(37)

where we recall that $t[\omega]$ and $V[\omega]$ denote the mean vector and
quantum covariance matrix of the state $\omega$ as defined in (19). A generic
Gaussian unitary is obtained as the product between a symplectic unitary and a
displacement operator. States obtained by applying an arbitrary Gaussian
unitary to the vacuum state $\ket{0}=\bigotimes_{j=1}^{m}\ket{0}_{j}$ are
called pure Gaussian states. Often times, displacements can be ignored; we
will thus write an arbitrary pure Gaussian state with zero mean as

$\displaystyle\ket{\psi}=U_{S}\ket{0}\,,$ (38)

where $S\in\operatorname{Sp}(2m)$ is an arbitrary symplectic matrix.

A quantum channel that will be particularly useful to us is the Gaussian white
noise channel, defined for $\lambda>0$ by

$\displaystyle\mathcal{N}_{\lambda}(\cdot)\coloneqq\int\frac{d^{2m}x}{(2\pi\lambda)^{m}}\,e^{-\frac{\|x\|^{2}}{2\lambda}}\pazocal{D}(x)(\cdot)\pazocal{D}(-x)\,.$
(39)

Using this formula one can show that

$\displaystyle\mathcal{N}_{\lambda}:\chi_{Z}\longmapsto\chi_{\mathcal{N}_{\lambda}(Z)}(x)\coloneqq\chi_{Z}(x)\,e^{-\frac{\lambda}{2}\|x\|^{2}}\,.$
(40)

Curiously, for $\lambda\in(0,1]$ its action can be expressed alternatively as

$\displaystyle\mathcal{N}_{\lambda}(\cdot)=\int\\!\\!\frac{d^{2m}x}{(2\pi\lambda)^{m}}\,\pazocal{D}(x)\,\tau_{\\!\frac{1}{2\lambda}-\frac{1}{2}}^{\otimes
m}\pazocal{D}(x)^{\dagger}(\cdot)\pazocal{D}(x)\,\tau_{\\!\frac{1}{2\lambda}-\frac{1}{2}}^{\otimes
m}\pazocal{D}(x)^{\dagger}\\!,$ (41)

where $\tau_{\nu}$ is the single-mode thermal state with mean photon number
$\nu$, given by

$\displaystyle\tau_{\nu}=\frac{1}{\nu+1}\sum_{n=0}^{\infty}\left(\frac{\nu}{\nu+1}\right)^{n}\ket{n}\\!\\!\bra{n}\,,$
(42)

where $\ket{n}$ stands for the $n^{\text{th}}$ Fock state. Since we could not
locate a complete proof of (41) in the existing literature, we provide a self-
contained one in Appendix A. A special case of (41) is when $\lambda=1$, in
which case $\tau_{\frac{1}{2\lambda}-\frac{1}{2}}=\ket{0}\\!\\!\bra{0}$ (the
vacuum state) and

$\displaystyle\mathcal{N}_{1}(\cdot)=\int\frac{d^{2m}x}{(2\pi)^{m}}\,\ket{x}\\!\\!\bra{x}(\cdot)\ket{x}\\!\\!\bra{x}\,.$
(43)

In terms of the real mean vector and quantum covariance matrix, for all
$\lambda>0$ we have that

$\displaystyle\mathcal{N}_{\lambda}:\left\\{\begin{array}[]{lll}V&\longmapsto&V+2\lambda
I\,,\\\ t&\longmapsto&t\,.\end{array}\right.$ (46)

The above channel, $\mathcal{N}_{\lambda}$, is just an example within the
larger class of Gaussian channels. To construct the most general Gaussian
channel, take two arbitrary $2m\times 2m$ real matrices $X$ and $Y$ such that

$\displaystyle Y+i\Omega-iX\Omega X^{\intercal}\geq 0\,;$ (47)

the corresponding Gaussian channel, denoted as $\mathcal{G}_{X,Y}$, then acts
as

$\displaystyle\mathcal{G}_{X,Y}:\chi_{Z}\mapsto\chi_{\mathcal{G}_{X,Y}(Z)}(u)\coloneqq\chi_{Z}\left(\Omega^{\intercal}X^{\intercal}\Omega\,u\right)e^{-\frac{1}{4}u^{\intercal}\Omega^{\intercal}Y\Omega
u}$ (48)

and

$\displaystyle\mathcal{G}_{X,Y}:\left\\{\begin{array}[]{lll}V&\longmapsto&XVX^{\intercal}+Y\,,\\\
t&\longmapsto&Xt\,.\end{array}\right.$ (49)

It is worth observing that (47) implies that $Y\geq 0$ is positive semi-
definite.

An important class of operators on ${\pazocal{H}}_{m}$ that we consider in
this paper is the set of Schwartz operators [69]. They can be defined as those
trace class operators whose characteristic function is a Schwartz function on
$\mathbb{R}^{2m}$. We denote the set of Schwartz function as
$\mathfrak{S}(\mathbb{R}^{2m})$, and that of Schwartz operators as
$\mathfrak{S}({\pazocal{H}}_{m})$. In particular, we consider the set
$\mathfrak{S}({\pazocal{H}}_{m})_{0}$ of Schwartz operators whose
characteristic functions are compactly supported.

As spaces of Schwartz functions are usually employed as test spaces in the
rigorous theory of distributions, we can use the space
$\mathfrak{S}({\pazocal{H}}_{m})_{0}$ to formalise the definition of objects —
called symbols — that would be ill-defined as operators in the traditional
sense. For example, given _any_ smooth function $\chi$ on $\mathbb{R}^{2m}$,
we can construct a symbol $\rho_{\chi}$ with ‘characteristic function’ $\chi$.
This is defined formally as a functional
$\rho_{\chi}:\mathfrak{S}({\pazocal{H}}_{m})_{0}\to\mathds{C}$ acting as

$\displaystyle\rho_{\chi}(Z)\coloneqq\int\frac{d^{2m}x}{(2\pi)^{m}}\,\chi(x)^{*}\chi_{Z}(x)\,.$
(50)

This expression is justified by the fact that, when
$\sigma\in{\pazocal{T}}_{1}({\pazocal{H}}_{m})$, we have
$\rho_{\chi_{\sigma}}(Z)=\operatorname{Tr}(\sigma Z)$ by Plancherel’s theorem
(34). The above functional extends to the whole space
$\mathfrak{S}({\pazocal{H}}_{m})$ whenever the function
$x\mapsto\chi(x)^{*}\chi_{Z}(x)$ is integrable.

###  III-C Concentration inequalities

In this paper, we make use of Bernstein’s matrix inequality in order to prove
that the probability that, on a well-chosen finite-dimensional subspace, the
output of our shadow tomography protocol is far from the original unknown
state decays exponentially fast in the number $N$ of samples used to gather
statistics:

######  Lemma 4 ((Bernstein’s matrix inequality [70])).

Given $N$ i.i.d. random matrices
$X_{1},\dots,X_{N}\in\mathbb{M}_{n}(\mathbb{C})$ which obey
$\|X_{i}-\mathbb{E}[X_{i}]\|_{\infty}\leq R$ almost surely, for some $R>0$,
the following tail bound holds:

$\displaystyle\mathbb{P}\bigg{(}\bigg{\|}\frac{1}{N}\\!\sum_{i=1}^{N}\,(X_{i}\\!-\\!\mathbb{E}[X_{i}])\bigg{\|}_{\infty}\\!\\!\geq{\epsilon}\bigg{)}\leq
2n\,e^{-\frac{N\epsilon^{2}}{2\Sigma^{2}+2R\epsilon/3}},$ (51)

where the constant $\Sigma$ is defined as
$\Sigma^{2}\coloneqq\|\mathbb{E}[X_{1}^{2}]\|_{\infty}<\infty$.

##  IV Classical shadow tomography of a CV system

In this section, we will have a closer look at the shadow tomography scheme
sketched in Section I-A. The goal of the procedure is to construct a good
estimator of an unknown $m$-mode state $\rho$, by measuring as few i.i.d.
copies of $\rho$ as possible. To this end, we repeatedly sample symplectic
matrices from $\operatorname{Sp}(2m)$ according to some probability
distribution $\mu$, apply the corresponding symplectic unitary $U_{S}$ on one
copy of $\rho$, implementing the transformation $\rho\mapsto
U_{S}^{\phantom{{\dagger}}}\rho U_{S}^{\dagger}$, and subsequently perform a
fixed Gaussian measurement
$\big{\\{}\psi_{x}\,\frac{d^{2m}x}{(2\pi)}\big{\\}}_{x\in\mathbb{R}^{2m}}$ on
that same state. Here, $\psi_{x}=\ket{\psi_{x}}\\!\\!\bra{\psi_{x}}$ with
$\ket{\psi_{x}}\coloneqq\pazocal{D}(x)\ket{\psi}$; the Gaussian state
$\ket{\psi}$, which uniquely identifies the Gaussian measurement, is a fixed
parameter of the shadow tomography protocol. Without loss of generality, we
can take $\ket{\psi}$ to have zero mean, in which case, according to (38), we
can introduce a symplectic matrix $T\in\operatorname{Sp}(2m)$ satisfying that
$U_{T}^{\dagger}\ket{0}=\ket{\psi}$. The ${\dagger}$ here is immaterial,
thanks to (36). Now, acting with a displacement operator on the left and on
the right yields immediately

$\displaystyle\ket{\psi_{x}}=\pazocal{D}(x)\ket{\psi}=U_{\\!T}^{{\dagger}}\left(U_{\\!T}^{\vphantom{{\dagger}}}\pazocal{D}(x)U_{\\!T}^{{\dagger}}\right)\ket{0}=U_{\\!T}^{{\dagger}}\ket{Tx}\,,$
(52)

where the last step is due to the action of symplectic unitaries on
displacement operators, see (35), and to the definition (26) of coherent
states.

The measurement makes the system collapse into a random state
$\widetilde{\rho}_{x}=\psi_{x}$, where $x$ is distributed with probability
distribution

$\displaystyle
p_{\rho}(x|S)\,d^{2m}x=\braket{\psi_{x}}{U_{S}^{\phantom{{\dagger}}}\rho
U_{S}^{\dagger}}{\psi_{x}}\frac{d^{2m}x}{(2\pi)^{m}}\,.$ (53)

When combined with the probability measure $\mu$ on $\operatorname{Sp}(2m)$,
this yields a joint probability distribution on
$\operatorname{Sp}(2m)\times\mathbb{R}^{2m}$.

We then attempt to undo the effect of the symplectic unitary by applying
$U_{S}^{\dagger}$. This amounts to the mapping
$\widetilde{\rho}_{x}=\psi_{x}\mapsto
U_{\\!S}^{{\dagger}}\widetilde{\rho}_{x}U_{\\!S}^{\vphantom{{\dagger}}}$. Once
we average over the random variable $S$ and the random variable $X$ whore
realisation we denoted with $x$, the whole process yields an effective noisy
channel $\mathcal{M}$ modelled as (1). Making this more explicit, we write the
action of $\mathcal{M}$ on an arbitrary trace class operator
$Z\in\pazocal{T}(\pazocal{H}_{m})$ as

$\displaystyle\mathcal{M}(Z)=\int\\!\\!d\mu(S)\int\\!\\!\frac{d^{2m}x}{(2\pi)^{m}}\,\braket{\psi_{x}}{U_{\\!S}^{\vphantom{{\dagger}}}Z\,U_{S}^{\dagger}}{\psi_{x}}\,U_{\\!S}^{{\dagger}}\psi_{x}U_{\\!S}^{\vphantom{{\dagger}}}\,.$
(54)

Our first result allows us to express the action of $\mathcal{M}$ in a form
that is more easily amenable to investigation with phase space methods.

######  Lemma 5.

The map
$\mathcal{M}:{\pazocal{T}}_{1}({\pazocal{H}}_{m})\to{\pazocal{T}}_{1}({\pazocal{H}}_{m})$
defined in (1) can be re-expressed as

$\displaystyle\mathcal{M}(Z)=\int\\!\\!d\mu(S)\,\mathcal{G}_{I,\,2(TS)^{-1}(TS)^{-\intercal}}(Z)\,,$
(55)

where $\mathcal{G}_{X,Y}$ is the Gaussian channel defined in (48). Moreover,
the characteristic function of $\mathcal{M}(Z)$ satisfies

$\displaystyle\chi_{\mathcal{M}(Z)}(x)$
$\displaystyle=\int\\!\\!d\mu(S)\,\chi_{Z}(x)\,e^{-\frac{1}{2}\left\|TSx\right\|^{2}}=\chi_{Z}(x)\,f_{\mu,T}(x)\,,$
(56)

where

$\displaystyle
f_{\mu,T}(x)\coloneqq\int\\!\\!d\mu(S)\,e^{-\frac{1}{2}\left\|TSx\right\|^{2}}.$
(57)

###### Proof.

We start by noticing that due to (52) and thanks the fact that $S\mapsto
U_{\\!S}^{\vphantom{{\dagger}}}$ is a group homomorphism, one obtains that

$\displaystyle
U_{S}^{\dagger}\ket{\psi_{x}}=U_{S}^{\dagger}U_{T}^{\dagger}\ket{Tx}=\left(U_{T}U_{S}\right)^{\dagger}\ket{Tx}=U_{TS}^{\dagger}\ket{Tx}\,.$
(58)

Thus, for any trace class operator $Z\in{\pazocal{T}}_{1}({\pazocal{H}}_{m})$

$\displaystyle\int\\!\\!\frac{d^{2m}x}{(2\pi)^{m}}\,\braket{\psi_{x}}{U_{\\!S}^{\vphantom{{\dagger}}}ZU_{S}^{\dagger}}{\psi_{x}}\,U_{\\!S}^{{\dagger}}\psi_{x}U_{\\!S}^{\vphantom{{\dagger}}}$
$\displaystyle\quad\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(i)}}}}{{=}}\int\\!\\!\frac{d^{2m}x}{(2\pi)^{m}}\,\braket{Tx}{U_{\\!TS}^{\vphantom{{\dagger}}}ZU_{\\!TS}^{{\dagger}}}{Tx}\,U_{\\!TS}^{{\dagger}}\ket{Tx}\\!\\!\bra{Tx}U_{\\!TS}^{\vphantom{{\dagger}}}$
$\displaystyle\quad\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(ii)}}}}{{=}}\int\\!\\!\frac{d^{2m}y}{(2\pi)^{m}}\,\braket{y}{U_{\\!TS}^{\vphantom{{\dagger}}}ZU_{\\!TS}^{{\dagger}}}{y}\,U_{\\!TS}^{{\dagger}}\ket{y}\\!\\!\bra{y}U_{\\!TS}^{\vphantom{{\dagger}}}$
(59)
$\displaystyle\quad=U_{\\!TS}^{{\dagger}}\left(\int\\!\\!\frac{d^{2m}y}{(2\pi)^{m}}\,\ket{y}\\!\\!\bra{y}U_{\\!TS}^{\vphantom{{\dagger}}}ZU_{\\!TS}^{{\dagger}}\ket{y}\\!\\!\bra{y}\right)U_{\\!TS}^{\vphantom{{\dagger}}}$
$\displaystyle\quad\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(iii)}}}}{{=}}U_{\\!TS}^{{\dagger}}\,\mathcal{N}_{1}\left(U_{\\!TS}^{\vphantom{{\dagger}}}ZU_{\\!TS}^{{\dagger}}\right)U_{\\!TS}^{\vphantom{{\dagger}}}$

Here, in (i) we applied (58) twice; in (ii) we changed variable, defining
$y\coloneqq Tx$, and used the fact that $T$, being symplectic, has determinant
$1$; and in (iii) we employed the representation in (43) for the action of
$\mathcal{N}_{1}$.

Now, let us compute the action of the above transformation at the level of
covariance matrices. By applying (37) and (46), we see that

$\displaystyle
V\left[U_{\\!TS}^{{\dagger}}\,\mathcal{N}_{1}\left(U_{\\!TS}^{\vphantom{{\dagger}}}Z\,U_{\\!TS}^{{\dagger}}\right)U_{\\!TS}^{\vphantom{{\dagger}}}\right]$
$\displaystyle\qquad=(TS)^{-1}V\left[\mathcal{N}_{1}\left(U_{\\!TS}^{\vphantom{{\dagger}}}Z\,U_{\\!TS}^{{\dagger}}\right)\right](TS)^{-\intercal}$
$\displaystyle\qquad=(TS)^{-1}\left(V\left[U_{\\!TS}^{\vphantom{{\dagger}}}Z\,U_{\\!TS}^{{\dagger}}\right]+2I\right)(TS)^{-\intercal}$
$\displaystyle\qquad=(TS)^{-1}\left(TS\,V[Z](TS)^{\intercal}+2I\right)(TS)^{-\intercal}$
$\displaystyle\qquad=V[Z]+2(TS)^{-1}(TS)^{-\intercal}\,.$

Comparing the above calculation with (49), we see that

$\displaystyle
U_{\\!TS}^{{\dagger}}\,\mathcal{N}_{1}\left(U_{\\!TS}^{\vphantom{{\dagger}}}Z\,U_{\\!TS}^{{\dagger}}\right)U_{\\!TS}^{\vphantom{{\dagger}}}=\mathcal{G}_{I,\,2(TS)^{-1}(TS)^{-\intercal}}(Z)\,.$
(60)

Using this insight in (59) shows that

$\displaystyle\int\\!\\!\frac{d^{2m}x}{(2\pi)^{m}}\,\braket{\psi_{x}}{U_{\\!S}^{\vphantom{{\dagger}}}ZU_{S}^{\dagger}}{\psi_{x}}\,U_{\\!S}^{{\dagger}}\psi_{x}U_{\\!S}^{\vphantom{{\dagger}}}=\mathcal{G}_{I,\,2(TS)^{-1}(TS)^{-\intercal}}(Z)\,.$
(61)

In turn, the above identity yields (55) upon integration in $S$ with respect
to the measure $\mu$.

We conclude by computing the characteristic function (28) of both sides of
(55). For an arbitrary $x\in\mathbb{R}^{2m}$, we obtain that

$\displaystyle\chi_{\mathcal{M}(Z)}(x)$ $\displaystyle=\int
d\mu(S)\,\chi_{\mathcal{G}_{I,\,2(TS)^{-1}(TS)^{-\intercal}}(Z)}(x)$ (62)
$\displaystyle\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(iv)}}}}{{=}}\int
d\mu(S)\,\chi_{Z}(x)e^{-\frac{1}{2}x^{\intercal}\Omega^{\intercal}(TS)^{-1}(TS)^{-\intercal}\Omega
x}$
$\displaystyle\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(v)}}}}{{=}}\int
d\mu(S)\,\chi_{Z}(x)e^{-\frac{1}{2}x^{\intercal}(TS)^{\intercal}(TS)x}$
$\displaystyle=\int
d\mu(S)\,\chi_{Z}(x)e^{-\frac{1}{2}\left\|TSx\right\|^{2}}.$

Here, (iv) follows from (48), and (v) is due to the fact that since $TS$ is
symplectic, $\Omega^{\intercal}(TS)^{-1}=(TS)^{\intercal}\Omega^{\intercal}$,
or upon transposing $(TS)^{-\intercal}\Omega=\Omega(TS)$. ∎

We therefore see that the action of $\mathcal{M}$ is actually very simple,
amounting to a point-wise multiplication at the level of the characteristic
function. Channels of this form are particular examples of so-called _linear
bosonic channels_ , introduced and studied by Holevo and Werner [8]. Although
we will not use this observation in this work, it is worth noting that linear
bosonic channels are always approximable — in the strong operator sense — by
_Gaussian dilatable_ channels, i.e. channels admitting a Stinespring
representation in which the unitary is Gaussian and the ancilla is arbitrary
[71, 72].

As it turns out, $\mathcal{M}$ is a special type of linear bosonic channel
whose action is representable as a random displacement. Namely, using (25) one
sees that

$\displaystyle\mathcal{M}(Z)=$ $\displaystyle\
\int\frac{d^{2m}x}{(2\pi)^{m}}\,\widetilde{f}_{\mu,T}(x)\,\pazocal{D}(x)Z\pazocal{D}(-x)\,,$
(63) $\displaystyle\widetilde{f}_{\mu,T}(x)\coloneqq$ $\displaystyle\
\int\frac{d^{2m}u}{(2\pi)^{m}}\,f_{\mu,T}(u)\,e^{ix^{\intercal}\Omega
u}=f_{\mu,T}(x)\,,$ (64)

i.e. $f_{\mu,T}$ coincides with its own Fourier transform. The rigorous proof
of the general validity of (64) is deferred to Appendix B. In light of this
discussion, it is easy to write down, at least formally, the inverse of
$\mathcal{M}$, which acts as

$\displaystyle\mathcal{M}^{-1}:\chi_{Z}\longmapsto\chi_{\mathcal{M}^{-1}(Z)}(x)\coloneqq\frac{1}{f_{\mu,T}(x)}\,\chi_{Z}(x)\,.$
(65)

Note that $f_{\mu,T}(x)>0$ for all $x\in\mathds{R}^{2m}$, so the above
expression is always well defined. In particular, $\mathcal{M}$ is injective
as a linear map, and thus it is invertible on its range.

The problem, naturally, is that such a range is in general significantly
smaller than the space of trace class operators (we will shortly see an
example of this). This entails that the right-hand side of (65) is not always
the characteristic function of a trace class operator. This, however, does not
pose any problem since the classical shadow obtained at the output of
$\mathcal{M}^{-1}$ is, as the name suggests, just a classical object which we
will merely use as a computational tool.

Incidentally, also $\mathcal{M}^{-1}$, just like $\mathcal{M}$, acts as a mere
point-wise multiplication at the level of characteristic functions. This
implies that one can also try, as done in (63), to represent it as an affine
combination of displacement operators, i.e. by writing

$\displaystyle\mathcal{M}^{-1}(Z)=\int\frac{d^{2m}x}{(2\pi)^{m}}\,\widetilde{g}_{\mu,T}(x)\,\pazocal{D}(x)Z\pazocal{D}(-x)\,,$
(66)

where $g_{\mu,T}(x)\coloneqq\frac{1}{f_{\mu,T}(x)}$, and

$\displaystyle\widetilde{g}_{\mu,T}(x)\coloneqq\int\frac{d^{2m}u}{(2\pi)^{m}}\,g_{\mu,T}(u)\,e^{ix^{\intercal}\Omega
u}$ (67)

would be its Fourier transform. The trouble, of course, is that $g_{\mu,T}$
will not be absolutely integrable — and not even bounded — in general, so
there is little hope to define its Fourier transform unless one appeals to the
theory of distributions.

To provide a solution to this apparent issue, let us return to our original
problem. We can now formally construct the classical shadow

$\displaystyle\hat{\rho}$
$\displaystyle=\mathcal{M}^{-1}\left(U_{\\!S}^{{\dagger}}\pazocal{D}(x)\psi\pazocal{D}(-x)U_{\\!S}^{\vphantom{{\dagger}}}\right)$
$\displaystyle=\mathcal{M}^{-1}\left(U_{\\!TS}^{{\dagger}}\ket{Tx}\\!\\!\bra{Tx}U_{\\!TS}^{\vphantom{{\dagger}}}\right)\,.$

The shadow, $\hat{\rho}$, is more rigorously defined as a functional on the
set $\mathfrak{S}({\pazocal{H}}_{m})_{0}$ of Schwartz operators with compactly
supported characteristic functions via (50). Its corresponding improper
characteristic function can hence be computed as follows:

$\displaystyle\hat{\chi}(u)\equiv\chi_{\hat{\rho}}(u)$
$\displaystyle=\chi_{\mathcal{M}^{-1}\left(U_{\\!TS}^{{\dagger}}\ket{Tx}\\!\bra{Tx}U_{\\!TS}^{\vphantom{{\dagger}}}\right)}(u)$
(68)
$\displaystyle=\frac{1}{f_{\mu,T}(u)}\,\chi_{U_{\\!TS}^{{\dagger}}\ket{Tx}\\!\bra{Tx}U_{\\!TS}^{\vphantom{{\dagger}}}}(u)$
$\displaystyle=\frac{1}{f_{\mu,T}(u)}\,\chi_{\ket{Tx}\\!\bra{Tx}}(TSu)$
$\displaystyle={\frac{1}{f_{\mu,T}(u)}\,e^{-\frac{1}{4}\left\|TSu\right\|^{2}-iu^{\intercal}\Omega
S^{-1}x}}\,,$

where in the last step we used (31). In other words, for any
$Z\in\mathcal{S}({\pazocal{H}}_{m})$:

$\displaystyle\hat{\rho}(Z)\coloneqq\int\frac{d^{2m}u}{(2\pi)^{m}}\,\hat{\chi}(u)^{*}\,\chi_{Z}(u)\,,$
(69)

whenever the function $u\mapsto\hat{\chi}(u)^{*}\chi_{Z}(u)$ is integrable. In
what follows, we will also consider the reduced shadow over a subset $A$ of
$|A|=r$ modes, formally given by the partial trace
$``\hat{\rho}_{A}\coloneqq\operatorname{Tr}_{A^{c}}(\hat{\rho})"$ of the
shadow $\hat{\rho}$. Again, we will use the characteristic function to
rigorously define it: given a region $A$ of $r$ modes, it is defined for any
$u_{A}\in\mathbb{R}^{2r}$ as

$\displaystyle\hat{\chi}_{A}(u_{A})\equiv\chi_{\hat{\rho}_{A}}(u_{A})\coloneqq\chi_{\hat{\rho}}(u_{A},0)\,.$
(70)

In that case, we write for any $Z_{A}\in\mathfrak{S}({\pazocal{H}}_{r})$:

$\displaystyle\hat{\rho}_{A}(Z_{A})\coloneqq\int\frac{d^{2r}u_{A}}{(2\pi)^{r}}\,\hat{\chi}_{A}(u_{A})^{*}\,\chi_{Z_{A}}(u_{A})\,,$
(71)

whenever the function
$u_{A}\mapsto\hat{\chi}_{A}(u_{A})^{*}\,\chi_{Z_{A}}(u_{A})$ is integrable.

The following lemma further justifies the claim made in (3) that the shadow
$\hat{\rho}$ has average $\rho$ by construction.

######  Lemma 6.

For any subset $A$ of $|A|=r$ modes, and all
$Z_{A}\in\mathfrak{S}({\pazocal{H}}_{r})_{0}$ with corresponding
characteristic function $\chi_{Z_{A}}$, the random variable
$\hat{\rho}_{A}(Z_{A})$ defined via (71) is integrable and

$\displaystyle\mathbb{E}\big{[}\hat{\rho}_{A}(Z_{A})\big{]}\,=\operatorname{Tr}[\rho_{A}Z_{A}]\,,$

where the conditional expectation is taken with respect to the probability
density function
$p_{T}(S,x)\coloneqq\braket{Tx}{U_{TS}^{\phantom{{\dagger}}}\rho
U_{TS}^{\dagger}}{Tx}$ with respect to $\mu\otimes\frac{d^{2m}x}{(2\pi)^{m}}$
on $\operatorname{Sp}(2m)\times\mathbb{R}^{2m}$. The result extends to
$Z_{A}\in\mathfrak{S}({\pazocal{H}}_{r})$ under the condition of integrability
with respect to
$\mu\otimes\frac{d^{2m}x}{(2\pi)^{m}}\otimes\frac{d^{2r}u}{(2\pi)^{r}}$ of the
function

$\displaystyle(S,x,u)\mapsto\frac{\chi_{Z_{A}}(u)\,\chi_{\ket{Tx}\\!\bra{Tx}}(TS(u,0)^{\intercal})\,p_{T}(S,x)}{f_{\mu,T}(u,0)}\,.$
(72)

###### Proof.

We present the proof for $|A|=m$ since the case $|A|=r<m$ follows the exact
same strategy using that $\hat{\chi}(u_{A},0)=\hat{\chi}_{A}(u_{A})$ by
definition. Now, on the one hand, if $Z\in\mathfrak{S}({\pazocal{H}}_{m})_{0}$
then the characteristic function $\chi_{Z}$ of $Z$ is compactly supported;
therefore, owing to the boundedness of $\chi_{Z}$ we deduce that the function

$\displaystyle(S,x,u)\mapsto\frac{\chi_{Z}(u)\chi_{\ket{Tx}\\!\bra{Tx}}(TSu)\,p_{T}(S,x)}{f_{\mu,T}(u)}$
(73)

is integrable with respect to
$\mu\otimes\frac{d^{2m}x}{(2\pi)^{m}}\otimes\frac{d^{2m}u}{(2\pi)^{m}}$. On
the other hand, if only $Z\in\mathfrak{S}({\pazocal{H}}_{m})$ such
integrability is assumed by hypothesis. Therefore, in both cases thanks to
Fubini’s theorem and (68) we have that

$\displaystyle\mathbb{E}\big{[}\hat{\rho}(Z)\big{]}$ $\displaystyle\
=\\!\int\\!\\!\\!\frac{d^{2m}u}{(2\pi)^{m}}\\!\int\\!\\!\\!d\mu(S)\\!\int\\!\\!\\!\frac{d^{2m}x}{(2\pi)^{m}}\,\chi_{Z}(u)\,\frac{\chi_{\ket{Tx}\\!\bra{Tx}}(TSu)\,p_{T}(S,x)}{f_{\mu,T}(u)}$
$\displaystyle\
\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(i)}}}}{{=}}\\!\int\\!\\!\\!\frac{d^{2m}u}{(2\pi)^{m}}\\!\int\\!\\!\\!d\mu(S)\,\chi_{Z}(u)\,\frac{\chi_{\mathcal{N}_{1}\left(U_{TS}^{\phantom{{\dagger}}}\rho
U_{TS}^{\dagger}\right)}(TSu)}{f_{\mu,T}(u)}$ $\displaystyle\
\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(ii)}}}}{{=}}\\!\int\\!\\!\\!\frac{d^{2m}u}{(2\pi)^{m}}\\!\int\\!\\!\\!d\mu(S)\,\frac{\chi_{Z}(u)\,\chi_{U_{TS}^{\phantom{{\dagger}}}\rho
U_{TS}^{\dagger}}(TSu)\,e^{-\frac{1}{2}\|TSu\|^{2}}}{f_{\mu,T}(u)}$
$\displaystyle\
\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(iii)}}}}{{=}}\\!\int\\!\\!\\!\frac{d^{2m}u}{(2\pi)^{m}}\\!\int\\!\\!\\!d\mu(S)\,\frac{\chi_{Z}(u)\,\chi_{\rho}(u)\,e^{-\frac{1}{2}\|TSu\|^{2}}}{f_{\mu,T}(u)}$
$\displaystyle\
=\\!\int\\!\\!\\!\frac{d^{2m}u}{(2\pi)^{m}}\,\chi_{Z}(u)\,\chi_{\rho}(u)$
$\displaystyle\ =\operatorname{Tr}[Z\rho]\,.$

In (i) we changed variable, defining $y\coloneqq Tx$, used the fact that $T$,
being symplectic, has determinant $1$, and employed the representation in
(43); in (ii) we used (40); finally, in (iii) we leveraged (35). ∎

In this section, we have made rigorous our first intuitive notion of a CV
shadow by means of its characteristic function. We have also seen that the
latter reduces to the characteristic function of the original unknown state
$\rho$ on average when integrated against sufficiently smooth functions. The
goal of the next section is to prove that under some physically relevant
conditions such as energy boundedness of the state $\rho$, these integrals are
enough to estimate expected values of observables with respect to the state
$\rho$ to high accuracy.

###  IV-A Finite moments assumption

As we saw in the previous section, the shadow $\hat{\rho}$, which we formally
defined through its characteristic function, is in general unbounded. In this
section, we show that this issue can be fixed if we further assume that the
photon number distribution of the unknown state $\rho$ satisfies some moment
constraints. This allows us to show that the projection of $\rho$ onto a
certain finite subspace of Hilbert–Schmidt operators is enough to get a
sufficiently good approximation of it while running the shadow tomography
protocol.

More precisely, we argue that all we need is to ensure finite rank convergence
of density operators and an energy constraint on the initial state $\rho$.

Let $Z\in{\pazocal{T}}_{2}({\pazocal{H}}_{m})$ be a Hilbert Schmidt operator
and denote by $\chi_{Z}$ its characteristic function. By density, one can find
a function $\widetilde{\chi}_{Z}\in C_{c}^{\infty}(\mathbb{R}^{2m})$ such that
$\|\chi_{Z}-\widetilde{\chi}_{Z}\|_{L^{2}(\mathbb{R}^{2m})}<\varepsilon$ and
$\widetilde{\chi}_{Z}$ is the characteristic function of some operator
$\widetilde{Z}$, i.e. $\chi_{\widetilde{Z}}=\widetilde{\chi}_{Z}$. The
operator $\widetilde{Z}$ is a Schwartz operator, and is in particular trace
class, and clearly $\|Z-\widetilde{Z}\|_{2}<\varepsilon$ [69]. Thus,
$\hat{\rho}(\widetilde{Z})$ is now necessarily a well-defined quantity by the
quantum Plancherel formula (cf. (34)). In particular, we have that, for any
$Z\in{\pazocal{T}}_{2}({\pazocal{H}}_{m})$,

$\displaystyle|\operatorname{Tr}(Z\rho)-\hat{\rho}(\widetilde{Z})|$
$\displaystyle\leq|\operatorname{Tr}((Z-\widetilde{Z})\rho)+|\operatorname{Tr}(\rho\widetilde{Z})-\hat{\rho}(\widetilde{Z})|$
$\displaystyle\leq\varepsilon+|\operatorname{Tr}(\rho\widetilde{Z})-\hat{\rho}(\widetilde{Z})|\,.$

This can be used in the above context by choosing
$Z=Z_{\mathbf{n}_{1}\mathbf{n}_{2}}=\ket{\mathbf{n}_{1}}\\!\\!\bra{\mathbf{n}_{2}}$,
given two multi-mode Fock states
$\ket{\mathbf{n}_{1}}\\!\\!\bra{\mathbf{n}_{2}}\in{\pazocal{H}}_{m}$, and
defining

$\mathcal{P}_{M}(T)\coloneqq\sum_{\mathbf{n}_{1},\mathbf{n}_{2}\in\\{0,..,M\\}^{m}}\langle\mathbf{n}_{1}|T|\mathbf{n}_{2}\rangle\,\ket{\mathbf{n}_{1}}\\!\\!\bra{\mathbf{n}_{2}}\equiv
P_{M}TP_{M}\,.$ (74)

Using the above approximation scheme, we can approximate the Schwartz
operators $Z_{\mathbf{n}_{1}\mathbf{n}_{2}}$ by Schwartz operators
$\widetilde{Z}_{\mathbf{n}_{1}\mathbf{n}_{2}}$ with smooth compactly supported
characteristic functions and define an auxiliary map

$\widetilde{\mathcal{P}}_{M}(T)\coloneqq\sum_{\mathbf{n}_{1},\mathbf{n}_{2}\in\\{0,..,M\\}^{m}}\operatorname{Tr}\left[\widetilde{Z}_{\mathbf{n}_{2}\mathbf{n}_{1}}T\right]\ket{\mathbf{n}_{1}}\\!\\!\bra{\mathbf{n}_{2}}\,.$
(75)

For sake of simplicity, we will use the same notation for projected reduced
states on subsets of $r<m$ modes. We further assume that the characteristic
function of $\widetilde{Z}_{\mathbf{n}_{1}\mathbf{n}_{2}}$ is obtained from
that of $\ket{\mathbf{n}_{1}}\\!\\!\bra{\mathbf{n}_{2}}$ by point-wise
multiplication by a simple compactly supported function. Namely, we set

$\displaystyle\widetilde{\chi}_{\mathbf{n}_{1}\mathbf{n}_{2}}(u)\coloneqq\widetilde{\chi}_{Z_{\mathbf{n_{1}}\mathbf{n}_{2}}}(u)=\chi_{\ket{\mathbf{n}_{1}}\\!\bra{\mathbf{n}_{2}}}(u)\prod_{j=1}^{m}\xi_{\eta,R}\big{(}u^{(j)}\big{)}\,,$
(76)

where as usual $u^{(j)}\coloneqq(u_{j},u_{j+m})^{\intercal}\in\mathds{R}^{2}$,
and, given some parameters $0<\eta<R$, the smooth function
$\xi_{\eta,R}:\mathds{R}^{2}\to[0,1]$ satisfies that

$\displaystyle\xi_{\eta,R}(z)=\left\\{\begin{array}[]{ll}1,&\qquad\text{if
$\|z\|\leq\eta$,}\\\\[4.30554pt] 0,&\qquad\text{if $\|z\|\geq
R$.}\end{array}\right.$ (77)

In particular, this allows us to evaluate
$\rho_{\chi}\big{(}\widetilde{Z}_{\mathbf{n}_{2}\mathbf{n}_{1}}\big{)}$ for
functionals $\rho_{\chi}$ defined as in (50).

######  Lemma 7.

For all $\alpha>0$, all non-negative integers $m$ (number of modes) and $M$
(Fock truncation number), and all $\eta\geq 0$, it holds that

$\displaystyle\left\|\mathcal{L}_{(I+N_{m})^{\alpha}}\mathcal{R}_{(I+N_{m})^{\alpha}}\left(\mathcal{P}_{M}-\widetilde{\mathcal{P}}_{{M}}\right)\right\|_{2\to
1}\leq\delta_{0}(\eta,M,\alpha,m)\,,$ (78)

where the rapidly vanishing function $\delta_{0}$ is defined by

$\displaystyle\delta_{0}(\eta,M,\alpha,m)$ (79)
$\displaystyle\quad\coloneqq(mM+1)^{2\alpha}(M+1)^{m}\,3^{mM}e^{-\frac{m}{4}\eta^{2}}\left(\sum_{p=0}^{2M}\frac{\eta^{2p}}{2^{p}p!}\right)^{m/2}.$

Next, we consider an approximation of $\rho$ in the following norm: given
$0\leq\alpha<n$ and two trace class operators
$\rho,\rho^{\prime}\in{\pazocal{T}}_{1}({\pazocal{H}}_{m})$ with
$\operatorname{Tr}(\rho(I+N_{m})^{n}),\operatorname{Tr}(\rho^{\prime}(I+N_{m})^{n})<\infty$,
we denote $H_{m}\coloneqq I+N_{m}$ and

$\displaystyle\|\rho-\rho^{\prime}\|^{(\alpha)}_{1}\coloneqq\left\|H_{m}^{\alpha/2}(\rho-\rho^{\prime})H_{m}^{\alpha/2}\right\|_{1}\,.$

Similar norms were previously defined in [73] under the name of m-mode bosonic
Sobolev norms.

######  Proposition 8.

Let $N_{m}$ be the number operator on ${\pazocal{H}}_{m}$, $H_{m}=I+N_{m}$,
and $\rho$ a state such that for some $n>0$ we have
$E_{m}^{(n)}\coloneqq\operatorname{Tr}[\rho\,H_{m}^{n}]<\infty$. Define
$\mathcal{P}_{M}$ and $\widetilde{\mathcal{P}}_{M}$ as in (74) and (75),
respectively. Then for any $0\leq\alpha<n$

$\displaystyle\left\|\rho-\mathcal{P}_{M}(\rho)\right\|_{1}^{(\alpha)}\leq
2(M+2)^{-\frac{n-\alpha}{2}}\,E_{m}^{(n)},$ (80)
$\displaystyle\left\|\rho-\widetilde{\mathcal{P}}_{{M}}(\rho)\right\|^{(\alpha)}_{1}\leq
2(M+2)^{-\frac{n-\alpha}{2}}\,E_{m}^{(n)}+\delta_{0}\big{(}\eta,M,\tfrac{\alpha}{2},m\big{)}\,,$
(81)

where $\delta_{0}$ was introduced in Lemma 7.

###### Proof.

We use a duality argument for this proof. For any $O$ such that
$\big{\|}H_{m}^{-\alpha/2}OH_{m}^{-\alpha/2}\big{\|}_{\infty}\leq 1$, we find

$\displaystyle\left|\operatorname{Tr}\big{[}O(\rho-\mathcal{P}_{M}(\rho))\big{]}\right|$
$\displaystyle\quad=\left|\operatorname{Tr}\left[H_{m}^{-\frac{\alpha}{2}}OH_{m}^{-\frac{\alpha}{2}}H_{m}^{\frac{\alpha}{2}}(\rho-\mathcal{P}_{M}(\rho))H_{m}^{\frac{\alpha}{2}}\right]\right|$
$\displaystyle\quad\leq\left\|H_{m}^{\frac{\alpha}{2}}(\rho-\mathcal{P}_{M}(\rho))H_{m}^{\frac{\alpha}{2}}\right\|_{1}$
$\displaystyle\quad\leq\left\|H_{m}^{-\frac{n-\alpha}{2}}H_{m}^{\frac{n}{2}}(I-P_{M})\rho\,H_{m}^{\frac{n}{2}}H_{m}^{-\frac{n-\alpha}{2}}\right\|_{1}$
$\displaystyle\quad\quad+\left\|H_{m}^{-\frac{n-\alpha}{2}}H_{m}^{\frac{n}{2}}P_{M}\rho(I-P_{M})\,H_{m}^{\frac{n}{2}}H_{m}^{-\frac{n-\alpha}{2}}\right\|_{1}$
$\displaystyle\quad\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(i)}}}}{{\leq}}\left\|H_{m}^{-\frac{n-\alpha}{2}}(I-P_{M})\right\|_{\infty}\left\|H_{m}^{\frac{n}{2}}\rho\,H_{m}^{\frac{n}{2}}\right\|_{1}\left\|H_{m}^{-\frac{n-\alpha}{2}}\right\|_{\infty}$
$\displaystyle\quad\quad+\left\|H_{m}^{-\frac{n-\alpha}{2}}\right\|_{\infty}\|P_{M}\|_{\infty}\left\|H_{m}^{\frac{n}{2}}\rho\,H_{m}^{\frac{n}{2}}\right\|_{1}\left\|(I-P_{M})H_{m}^{-\frac{n-\alpha}{2}}\right\|_{\infty}$
$\displaystyle\quad\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(ii)}}}}{{\leq}}2(M+2)^{-\frac{n-\alpha}{2}}E_{m}^{(n)}\,.$

Here, (i) comes from a repeated application of Hölder’s inequality, and we
also observed that $P_{M}$ and $H_{m}$ commute; (ii) is because

$\displaystyle\left\|H_{m}^{-\frac{n-\alpha}{2}}(I-P_{M})\right\|_{\infty}$
$\displaystyle=\left\|\sum\nolimits_{\mathbf{n}\in\mathds{N}^{m}\setminus[M]^{m}}(1+|\mathbf{n}|)^{-\frac{n-\alpha}{2}}\ket{\mathbf{n}}\\!\\!\bra{\mathbf{n}}\right\|_{\infty}$
(82)
$\displaystyle=\max_{\mathbf{n}\in\mathds{N}^{m}\setminus[M]^{m}}(1+|\mathbf{n}|)^{-\frac{n-\alpha}{2}}$
$\displaystyle=(2+M)^{-\frac{n-\alpha}{2}}.$

This proves (80). Next, we have that

$\displaystyle\left|\operatorname{Tr}\left[O\,(\mathcal{P}_{M}-\widetilde{\mathcal{P}}_{{M}})(\rho)\right]\right|$
$\displaystyle\quad=\left|\operatorname{Tr}\left[H_{m}^{-\frac{\alpha}{2}}OH_{m}^{-\frac{\alpha}{2}}H_{m}^{\frac{\alpha}{2}}(\mathcal{P}_{M}(\rho)-\widetilde{\mathcal{P}}_{M}(\rho))H_{m}^{\frac{\alpha}{2}}\right]\right|$
$\displaystyle\quad\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(iii)}}}}{{\leq}}\left\|H_{m}^{-\frac{\alpha}{2}}OH_{m}^{-\frac{\alpha}{2}}\right\|_{\infty}\left\|H_{m}^{\frac{\alpha}{2}}\,(\mathcal{P}_{M}-\widetilde{\mathcal{P}}_{M})(\rho)H_{m}^{\frac{\alpha}{2}}\right\|_{1}$
$\displaystyle\quad\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(iv)}}}}{{\leq}}\delta_{0}\left(\eta,M,\tfrac{\alpha}{2},m\right),$

where (iii) is again Hölder’s inequality, and (iv) comes from our assumptions
together with Lemma 7. Combining this with (80) yields (81) and concludes the
proof. ∎

In order to use the approximation bounds of 8, it remains to estimate how well
the empirical average of the shadows

$\displaystyle{\hat{\sigma}}^{(N)}\coloneqq\frac{1}{N}\sum_{i=1}^{N}\hat{\rho}^{(i)}\,,$
(83)

where $\hat{\rho}^{(i)}$ are $N$ i.i.d. copies of the shadow $\hat{\rho}$
constructed in (69), approximates $\mathcal{P}_{M}(\rho)$ and
$\widetilde{\mathcal{P}}_{{M}}(\rho)$ depending on the choice of the Gaussian
shadow tomography scheme. For this, we introduce, for any set $A$ of $|A|=r$
modes, the operators

$\displaystyle{\sigma}^{(N)}_{A}(M)\coloneqq\sum_{\mathbf{n}_{1},\mathbf{n}_{2}\in\\{0,..,M\\}^{r}}\,\hat{\sigma}_{A}^{(N)}(\ket{\mathbf{n}_{2}}\\!\\!\bra{\mathbf{n}_{1}})\,\ket{\mathbf{n}_{1}}\\!\\!\bra{\mathbf{n}_{2}},$
$\displaystyle\widetilde{\sigma}_{A}^{(N)}(M)\coloneqq\sum_{\mathbf{n}_{1},\mathbf{n}_{2}\in\\{0,..,M\\}^{r}}\hat{\sigma}_{A}^{(N)}(\widetilde{Z}_{\mathbf{n}_{2}\mathbf{n}_{1}})\,\ket{\mathbf{n}_{1}}\\!\\!\bra{\mathbf{n}_{2}}\,.$
(84)

Since both operators are supported on a finite dimensional subspace, we can
resort to the matrix Bernstein inequality (LABEL:matrixBern) in order to prove
that, with high probability,
$\|{\sigma}_{A}^{(N)}(M)-{\mathcal{P}}_{{M}}(\rho_{A})\|_{1},\|\widetilde{\sigma}_{A}^{(N)}(M)-\widetilde{\mathcal{P}}_{{M}}(\rho_{A})\|_{1}\leq\epsilon$
for $N$ large enough. In the next section, we explain in more detail how we
use the matrix Bernstein inequality in the case of both the homodyne and the
heterodyne detection strategies.

##  V Homodyne and heterodyne shadow tomography

###  V-A Local homodyne detection

Let us first consider the scenario in which one performs a homodyne detection
along a random direction in phase space. Here, $m$ independent random matrices
$S_{1},\dots,S_{m}$ are distributed uniformly (according to the Haar measure)
on the intersection $\operatorname{Sp}(2)\cap\operatorname{SO}(2)$ between the
symplectic and the orthogonal group. It is useful to note that
$\operatorname{Sp}(2)\cap\operatorname{SO}(2)\simeq\operatorname{U}(1)$, where
on the right-hand side we have the unitary group of $1\times 1$ matrices, so
that the Haar measure on $\operatorname{Sp}(2)\cap\operatorname{SO}(2)$ is
essentially that on $\operatorname{U}(1)$. In this simple homodyne case the
algorithm for shadow tomography is summarised as follows:

  1. (1)

A copy of $\rho$ is loaded, and $m$ matrices
$S_{1},\dots,S_{m}\in\operatorname{U}_{1}$ are drawn at random according to
the Haar measure. In other words $U_{j}=R_{\theta_{j}}$ for some angle
$\theta_{j}\in[-\pi,\pi]$, where

$\displaystyle R_{\theta}\coloneqq\begin{pmatrix}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\end{pmatrix}\,.$

  2. (2)

The rotation $U_{S}\equiv U_{S_{1}}\otimes\dots\otimes U_{S_{m}}$ is applied
to $\rho$, where $S=S_{1}\oplus\dots\oplus S_{m}$, obtaining the state
$U_{S}\rho\,U_{S}^{\dagger}$.

  3. (3)

The output state of (2) is subjected to a homodyne measurement along the
position axis, yielding the classical outcome
$x=(x_{1},\dots,x_{m})^{\intercal}\in\mathbb{R}^{2m}$.

  4. (4)

We construct the classical shadow $\hat{\rho}$, which is the final output of
this round of the protocol.

After several rounds have been conducted, we can process the classical shadows
as we prefer. A typical method would be that of computing the empirical
average $\hat{\sigma}^{(N)}$ as defined in (83) for the above protocol. We can
then use that operator e.g. for computing expected values of observables, or
else reduced density operators, etc. To model homodyning in a rigorous way,
let us introduce a parameter $s>0$, which we will later take to infinity, and
let us set

$\displaystyle T^{\oplus m}=T^{\oplus
m}_{s}\coloneqq\begin{pmatrix}e^{-s}&0\\\ 0&e^{s}\end{pmatrix}^{\bigoplus m},$
(85)

These choices of $\mu\sim\operatorname{Haar}(\operatorname{U}_{1})$ and $T$
completely determine our shadow tomographic setting. In the limit
$s\to\infty$, the measurement will reproduce a homodyne measurement along the
position axis of each mode. To make things more concrete, given a threshold
$M\in\mathbb{N}$ and $\alpha\geq 0$ we define

$\displaystyle\hat{\rho}_{A}(M)$ $\displaystyle\coloneqq\bigotimes_{j\in
A}\,\Big{(}\sum\nolimits_{n_{1},n_{2}=0}^{M}\,\hat{\rho}_{j}(M)_{n_{1},n_{2}}\ket{n_{1}}\\!\\!\bra{n_{2}}\Big{)},$
(86) $\displaystyle\hat{\rho}_{j}(M)_{n_{1},n_{2}}$
$\displaystyle\coloneqq\int\\!\\!dy\ |y|\,e^{-i(y,0)S_{j}^{-\intercal}\Omega
S_{j}x_{j}}\,\chi_{\ket{n_{2}}\\!\bra{n_{1}}}(S_{j}(y,0)^{\intercal})\,,$

By (32), the absolute value of the characteristic function
$\chi_{|n_{2}(j)\rangle\langle n_{1}(j)|}(S_{j}(y,0)^{T})$ is upper bounded by

$\displaystyle\sqrt{\frac{n_{2}(j)!}{n_{1}(j)!}}\,e^{-\frac{\pi}{2}|\omega|^{2}}(\sqrt{\pi}|\omega|)^{M_{j}-m_{j}}\,L_{m_{j}}^{(M_{j}-m_{j})}(\pi|\omega|^{2})\,,$

where $\omega\coloneqq-x^{\prime}+i\frac{x^{\prime\prime}}{2\pi}$ with
$x=x^{\prime}\oplus x^{\prime\prime}=S_{j}(y,0)^{T}$,
$m_{j}\coloneqq\min\\{n_{1}(j),n_{2}(j)\\}$,
$M_{j}\coloneqq\max\\{n_{1}(j),n_{2}(j)\\}$, and where the Laguerre
polynomials $L_{j}^{(k)}$ were defined in (33). Since
$\frac{|y|}{2\pi}\leq|\omega|\leq|y|$, we therefore have that
$\left\|H_{r}^{\frac{\alpha}{2}}\hat{\rho}_{A}(M)H_{r}^{\frac{\alpha}{2}}\right\|_{\infty}$
is almost surely bounded by

$\displaystyle\Sigma^{(\alpha)}_{r}(M)$ (87) $\displaystyle\
\coloneqq\Bigg{\|}\,\Bigg{(}\\!\left((1+|\mathbf{n}_{1}|)(1+|\mathbf{n}_{2}|)\right)^{\frac{\alpha}{2}}\prod_{j=1}^{r}\,\sqrt{\frac{n_{2}(j)!}{n_{1}(j)!}}$
$\displaystyle\qquad\times\int\\!\\!dy\,|\sqrt{\pi}y|^{1+M_{j}-m_{j}}\,\,e^{-\frac{|y|^{2}}{8\pi}}\left|L_{M_{j}}^{M_{j}-m_{j}}\big{(}\pi|y|^{2}\big{)}\right|\Bigg{)}_{\mathbf{n}_{1},\mathbf{n}_{2}}\Bigg{\|}_{\infty},$

where $\mathbf{n}_{1},\mathbf{n}_{2}\in\\{0,\dots,M\\}^{r}$, by using the
simple fact that given two matrices $A=\\{(a_{ij})\\}_{ij}$ and
$B=\\{(b_{ij})\\}_{ij}$ with $|a_{ij}|\leq b_{ij}$ for all $i,j$, then
$\|A\|_{\infty}\leq\|B\|_{\infty}$. This can be proved as follows. For a
vector $v$ with entries $v_{i}$ and norm $\|v\|_{2}^{2}=\sum_{i}|v_{i}|^{2}$,
define the new vector $w$ with entries $w_{i}\coloneqq|v_{i}|$, so that
$\|w\|^{2}=1$. Using the triangle inequality, it is straightforward to verify
that $\|Av\|^{2}\leq\|Bw\|^{2}$; this entails that
$\|A\|_{\infty}=\sup_{\|v\|=1}\|Av\|\leq\sup_{\|w\|=1}\|Bw\|=\|B\|_{\infty}$.
Continuing, as before we then denote

$\displaystyle\sigma_{A}^{(N)}(M)\coloneqq\frac{1}{N}\sum_{i=1}^{N}\,\hat{\rho}^{(i)}_{A}(M)\,,$

where $\hat{\rho}^{(i)}_{A}(M)$ are i.i.d. random matrices of law the one of
$\hat{\rho}_{A}(M)$. With this, we are ready to state our first main result:

######  Theorem 9.

With the above notation, given $0\leq\alpha<n$ such that
$E^{(n)}_{r}\coloneqq\max_{|A|\leq
r}\operatorname{Tr}(\rho_{A}H_{r}^{n})<\infty$ and

$\displaystyle N$
$\displaystyle\geq\frac{(M+1)^{2r}}{3\epsilon^{2}}\Big{\\{}24\,\Sigma^{(\alpha)}_{r}(M)^{2}+4\Big{(}\Sigma_{r}^{(\alpha)}(M)+E_{r}^{(\alpha)}\Big{)}\epsilon\Big{\\}}$
(88)
$\displaystyle\quad\times\log\left(\frac{2\big{[}m(M+1)\big{]}^{r}}{\delta}\right),$

where

$M=\Big{\lceil}\big{(}4E_{r}^{(n)}/\epsilon\big{)}^{\frac{2}{n-\alpha}}\Big{\rceil},$
(89)

we have that for any region $A$ of size $|A|\leq r$, it holds that

$\left\|\sigma_{A}^{(N)}(M)-\rho_{A}\right\|^{(\alpha)}_{1}\leq\epsilon$ (90)

with probability at least $1-\delta$. Similarly, for

$\displaystyle N$
$\displaystyle\geq\frac{(M+1)^{2r}}{3\epsilon^{2}}\Big{\\{}24\,\Sigma^{(\alpha)}_{r}(M)^{2}+4\Big{(}\Sigma_{r}^{(\alpha)}(M)+E_{r}^{(\alpha)}\Big{)}\epsilon\Big{\\}}$
(91) $\displaystyle\quad\times\log\left(\frac{2L(M+1)^{r}}{\delta}\right)\,,$

we have that for any set of $L$ observables $O_{j}$ on regions $A_{j}$ of size
at most $r$ and with
$\left\|H_{r}^{-\frac{\alpha}{2}}O_{j}H_{r}^{-\frac{\alpha}{2}}\right\|_{\infty}\leq
1$,

$\max_{j}\big{|}\operatorname{Tr}\big{[}O_{j}\,(\sigma_{A_{j}}^{(N)}(M)-\rho_{A_{j}})\big{]}\big{|}\leq\epsilon$
(92)

with probability at least $1-\delta$.

######  Remark 10.

We essentially recover the same dependence of the number of samples in terms
of the logarithm of the number of observables/total number of modes and on the
exponential of the size of the regions $A$ as in the qubit setting of [6].

###### Proof.

The first part of the proof consists in computing the function
$f_{\mu^{\otimes m},T^{\oplus m}}$ defined as in (57) in the limit
$s\to\infty$. It is given by

$\displaystyle f_{\mu^{\otimes m},T^{\oplus m}}(u)$
$\displaystyle\coloneqq\int\\!\\!d\mu^{\otimes
m}(S_{1},\dots,S_{m})\,e^{-\frac{1}{2}\,u^{\intercal}(T^{\oplus
m}S)^{{}^{\intercal}}(T^{\oplus m}S)\,u}$
$\displaystyle=\prod_{j=1}^{m}f_{\mu,T}(u_{j})\,,$

with $u=(u_{1},\dots,u_{m})^{T}\in\mathbb{R}^{2m}$, where for each mode $j$ we
have

$\displaystyle f_{\mu,T}(u_{j})$
$\displaystyle\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(i)}}}}{{=}}\int_{-\pi}^{+\pi}\frac{d\theta}{2\pi}\,e^{-\frac{1}{2}u_{j}^{\intercal}R_{\theta}^{\intercal}D_{s}R_{\theta}u_{j}}$
(93)
$\displaystyle\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(ii)}}}}{{=}}\int_{-\pi/2}^{+\pi/2}\frac{d\theta}{\pi}\,e^{-\frac{1}{2}u_{j}^{\intercal}R_{\theta}^{\intercal}D_{s}R_{\theta}u_{j}}$
$\displaystyle\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(iii)}}}}{{=}}\int_{-\pi/2}^{+\pi/2}\frac{d\theta}{\pi}\,e^{-\frac{\|u_{j}\|^{2}}{2}\,e_{1}^{\intercal}R_{\theta}^{\intercal}D_{s}R_{\theta}e_{1}}$
$\displaystyle=\int_{-\pi/2}^{+\pi/2}\frac{d\theta}{\pi}\,e^{-\frac{\|u_{j}\|^{2}}{2}\,\left(e^{-2s}\cos^{2}\theta+e^{2s}\sin^{2}\theta\right)}$
$\displaystyle\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(iv)}}}}{{=}}e^{-\frac{\|u_{j}\|^{2}\cosh(2s)}{2}}I_{0}\Big{(}-\tfrac{\|u_{j}\|^{2}\sinh(2s)}{2}\Big{)}$
$\displaystyle=e^{-\frac{\|u_{j}\|^{2}\cosh(2s)}{2}}I_{0}\Big{(}\tfrac{\|u_{j}\|^{2}\sinh(2s)}{2}\Big{)}.$

Here, in (i) we introduced the matrix $D_{s}=T_{s}^{2}$; in (ii) we noted that
the integrand is invariant under rotations of $\pi$, in (iii) we observed that
the integral is invariant under rotations of $u_{j}$, and thus chose to
compute it for $u_{j}=\|u_{j}\|e_{1}\coloneqq\|u_{j}\|(1,0)^{\intercal}$; and
in (iv) we recognised the integral representation for the modified Bessel
function

$I_{0}(x)\coloneqq\frac{1}{\pi}\,\int_{0}^{\pi}e^{\pm x\cos\theta}d\theta\,.$

Plugging this into (68), we have found that

$\displaystyle\chi_{\hat{\rho}}(u)$
$\displaystyle=\frac{e^{-\frac{1}{4}\left\|T^{\oplus
m}Su\right\|^{2}-iu^{\intercal}\Omega Sx}}{f_{\mu^{\otimes m},T^{\oplus
m}}(u)}\,$
$\displaystyle=\prod_{j=1}^{m}\,\frac{e^{\frac{\|u_{j}\|^{2}\,\operatorname{cosh}(2s)}{2}}}{I_{0}\Big{(}\frac{\|u_{j}\|^{2}\operatorname{sinh}(2s)}{2}\Big{)}}e^{-\frac{1}{4}\left\|TS_{j}u_{j}\right\|^{2}-iu_{j}^{\intercal}\Omega
S_{j}^{-1}x_{j}}$
$\displaystyle\equiv\prod_{j=1}^{m}\,\chi_{\hat{\rho}_{j}}(u_{j})\,.$

In other words, $\chi_{\hat{\rho}}$ is formally the characteristic function of
a tensor product of forms
$\hat{\rho}=\hat{\rho}_{1}\otimes\dots\otimes\hat{\rho}_{m}$. Moreover, by the
asymptotic expansion of the modified Bessel function

$\displaystyle I_{0}\Big{(}\tfrac{\|u_{j}\|^{2}\sinh(2s)}{2}\Big{)}$
$\displaystyle=\frac{e^{\tfrac{\|u_{j}\|^{2}\sinh(2s)}{2}}}{\sqrt{\sinh(2s)\pi}\|u_{j}\|}\Big{(}1+\mathcal{O}((\|u_{j}\|^{2}\sinh(2s))^{-1})\Big{)},$

we see that

$\displaystyle\chi_{\hat{\rho}_{j}}(u_{j})=\frac{\sqrt{2\sinh(2s)\pi}\|u_{j}\|\,e^{\frac{\|u_{j}\|^{2}e^{-2s}}{2}}e^{-\frac{1}{4}\|T_{s}S_{j}u_{j}\|^{2}-iu_{j}^{\intercal}\Omega
S_{j}^{-1}x_{j}}}{1+\mathcal{O}((\|u_{j}\|^{2}\sinh(2s))^{-1})}$

which shows that $\chi_{\hat{\rho}}(u)$ can be integrated against any element
of the Fock basis as soon as $s>0$. Taking the distributional limit, we find

$\displaystyle\chi_{\hat{\rho}_{j}}(u_{j})\mathrel{\stackunder[5pt]{\xrightarrow{\makebox[22.09454pt]{$\scriptstyle\mathrm{}$}}}{\scriptstyle\\!s\to\infty\\!\,}}2\pi\|u_{j}\|\,\delta\left((S_{j}u_{j})_{2}\right)e^{-iu_{j}^{\intercal}\Omega
S_{j}^{-1}x_{j}}\,$

Next, we fix a region $A$ of size $|A|=r$ and consider for a fixed squeezing
$s$

$\displaystyle\hat{\rho}_{A}(s,M)\coloneqq\bigotimes_{j\in
A}\,\Big{(}\sum_{n_{1},n_{2}=1}^{M}\,\hat{\rho}_{j}(\ket{n_{2}}\\!\\!\bra{n_{1}}|)\,\ket{n_{1}}\\!\\!\bra{n_{2}}\Big{)}\,.$

We have from Lemma 6 that

$\displaystyle\mathbb{E}\big{[}\hat{\rho}_{A}(s,M)\big{]}=\mathcal{P}_{M}(\rho_{A})$

which holds for any squeezing parameter $s$. Therefore, we can take the limit
$s\to\infty$ above. We find

$\displaystyle\hat{\rho}_{A}(s,M)=\bigotimes_{j\in
A}\\!\bigg{(}\sum_{n_{1},n_{2}=1}^{M}\\!\\!\ket{n_{1}}\\!\\!\bra{n_{2}}\int\\!\\!\chi_{\hat{\rho}_{j}}(u_{j})\,\chi_{\ket{n_{2}}\\!\bra{n_{1}}}(u_{j})\,\frac{d^{2}u_{j}}{(2\pi)^{r}}\bigg{)}$
(94)

so that

$\displaystyle\hat{\rho}_{A}(s,M)$
$\displaystyle\mathrel{\stackunder[5pt]{\xrightarrow{\makebox[22.09454pt]{$\scriptstyle\mathrm{}$}}}{\scriptstyle\\!s\to\infty\\!\,}}\hat{\rho}_{A}(M)\coloneqq\bigotimes_{j\in
A}\hat{\rho}_{j}(M),$ (95) $\displaystyle\hat{\rho}_{j}(M)$
$\displaystyle\coloneqq\sum_{n_{1},n_{2}=1}^{M}\\!\\!\ket{n_{1}}\\!\\!\bra{n_{2}}$
$\displaystyle\qquad\times\int\\!\\!|y|\,e^{-i(y,0)S_{j}^{-\intercal}\Omega
S_{j}x_{j}}\,\chi_{\ket{n_{2}}\\!\bra{n_{1}}}(S_{j}(y,0)^{\intercal})\,dy,$

with $\mathbb{E}[\hat{\rho}_{A}(M)]=\mathcal{P}_{M}(\rho_{A})$. Then, we have
seen below Equation (86) that, almost surely,
$\|H_{r}^{\frac{\alpha}{2}}\hat{\rho}_{A}(M)H_{r}^{\frac{\alpha}{2}}\|_{\infty}\leq\Sigma^{(\alpha)}_{r}(M)$.
Moreover, we have that
$\|H_{r}^{\frac{\alpha}{2}}\mathcal{P}_{M}(\rho_{A})H_{r}^{\frac{\alpha}{2}}\|_{\infty}\leq
E_{r}^{(\alpha)}$, therefore almost surely
$\|H_{r}^{\frac{\alpha}{2}}(\hat{\rho}_{A}(M)-\mathcal{P}_{M}(\rho_{A}))H_{r}^{\frac{\alpha}{2}}\|_{\infty}\leq
E_{r}^{(\alpha)}+\Sigma^{(\alpha)}_{r}(M)$. Finally,
$\|\mathbb{E}[(H_{r}^{\frac{\alpha}{2}}\hat{\rho}_{A}(M)H_{r}^{\frac{\alpha}{2}})^{2}]\|_{\infty}\leq\Sigma_{r}^{(\alpha)}(M)^{2}$.
We can hence use the equivalence of matrix norms together with the matrix
Bernstein inequality (LABEL:matrixBern) to get

$\displaystyle\mathbb{P}\Big{(}\big{\|}\sigma_{A}^{(N)}(M)-\mathcal{P}_{M}(\rho_{A})\big{\|}^{(\alpha)}_{1}\geq{\epsilon}\Big{)}$
$\displaystyle\leq\mathbb{P}\Big{(}\big{\|}H_{r}^{\frac{\alpha}{2}}\big{(}\sigma_{A}^{(N)}(M)-\mathcal{P}_{M}(\rho_{A})\big{)}H_{r}^{\frac{\alpha}{2}}\big{\|}_{\infty}\geq{\epsilon}(M+1)^{-r}\Big{)}$
$\displaystyle\leq
2(M+1)^{r}\,\operatorname{exp}\left(-\frac{3N\epsilon^{2}(M+1)^{-2r}}{6\Sigma_{r}^{(\alpha)}(M)^{2}+2(\Sigma_{r}^{(\alpha)}(M)+E_{r}^{(\alpha)})\epsilon}\right).$
(96)

Therefore, by a union bound, we get that

$\displaystyle\mathbb{P}\Big{(}\exists A,\,|A|\leq
r:\,\big{\|}\sigma_{A}^{(N)}(M)-\mathcal{P}_{M}(\rho_{A})\big{\|}^{(\alpha)}_{1}\geq{\epsilon}\Big{)}$
$\displaystyle\leq 2\big{[}m(M+1)\big{]}^{r}$
$\displaystyle\quad\times\operatorname{exp}\left(-\frac{3N\epsilon^{2}(M+1)^{-2r}}{6\Sigma_{r}^{(\alpha)}(M)^{2}+2(\Sigma_{r}^{(\alpha)}(M)+E_{r}^{(\alpha)})\epsilon}\right)\,.$
(97)

Next, we use Proposition 8 and choose
$M=\Big{\lceil}\big{(}4E_{r}^{(n)}/\epsilon\big{)}^{\frac{2}{n-\alpha}}\Big{\rceil}$,
so that $\|\mathcal{P}_{M}(\rho_{A})-\rho_{A}\|_{1}^{(\alpha)}\leq\epsilon/2$
and

$\displaystyle\mathbb{P}\Big{(}\exists A,\,|A|\leq
r:\,\big{\|}\sigma_{A}^{(N)}(M)-\rho_{A}\big{\|}^{(\alpha)}_{1}\geq{\epsilon}\Big{)}$
$\displaystyle\leq\mathbb{P}\Big{(}\exists A,\,|A|\leq
r:\,\big{\|}\sigma_{A}^{(N)}(M)-\mathcal{P}_{M}(\rho_{A})\big{\|}^{(\alpha)}_{1}\geq\frac{\epsilon}{2}\Big{)}$
$\displaystyle\leq 2\big{[}m\big{(}M+1\big{)}\big{]}^{r}$
$\displaystyle\quad\times\operatorname{exp}\left(-\frac{3N\epsilon^{2}(M+1)^{-2r}}{24\Sigma^{(\alpha)}_{r}(M)^{2}+4(\Sigma^{(\alpha)}_{r}(M)+E_{r}^{(\alpha)})\epsilon}\right)\,.$

Therefore, choosing $N$ as in the statement of the theorem, we obtain that the
probability that on any subset $A$ of at most $r$ modes
$\|\sigma_{A}^{(N)}(M)-\rho_{A}\big{\|}^{(\alpha)}_{1}\geq{\epsilon}$ is at
least $1-\delta$.

The results for a fixed number of $L$ observables $O_{1},\dots,O_{L}$
supported on regions $A_{1},\dots,A_{L}$ of size at most $r$ follow after
replacing the above union bounds over regions $A$ by a union bound over the
observables $O_{j}$:

$\displaystyle\mathbb{P}\Big{(}\exists
j:\,|\operatorname{Tr}(O_{j}\sigma_{A_{j}}^{(N)}(M))-\operatorname{Tr}(O_{j}\,\mathcal{P}_{M}(\rho_{A_{j}}))|\geq{\epsilon}\Big{)}$
$\displaystyle\leq
2L\,\big{[}(M+1)\big{]}^{r}\,\operatorname{exp}\left(-\frac{3N\epsilon^{2}(M+1)^{-2r}}{6\Sigma_{r}^{(\alpha)}(M)^{2}+2(\Sigma_{r}^{(\alpha)}(M)+E_{r}^{(\alpha)})\epsilon}\right)$

and the result follows after replacing $m^{r}$ by $L$ in the estimate for $N$.

∎

###  V-B Local heterodyne detection

Even simpler than homodyne detection, the simplest Gaussian shadow tomographic
setting is that where the measurement employed is a heterodyne detection

$\displaystyle\left\\{\frac{1}{(2\pi)^{m/2}}\ket{x}\\!\\!\bra{x}\right\\}_{x\in\mathds{R}^{2m}}\,,$
(98)

i.e. $T=I$, and all unitaries employed are passive, i.e. such that
$\big{[}U_{\\!S}^{\vphantom{{\dagger}}},\frac{1}{2}R^{\intercal}R\big{]}=0$.
Naturally, this is the same as requiring that $S$ be (not only symplectic but
also) orthogonal. Since passive unitaries send coherent states to coherent
states, amounting to a rotation in that space, the effective measurement being
carried out on $\rho$ is the same irrespectively of $S$ — the only thing
changing is that the outcome is rotated. Mathematically, this means that the
probability distribution of the random variable
$U_{\\!S}^{{\dagger}}\pazocal{D}(x)\psi\pazocal{D}(-x)U_{\\!S}^{\vphantom{{\dagger}}}$,
where $x$ is the outcome of the heterodyne detection (98) on
$U_{\\!S}^{\vphantom{{\dagger}}}\psi U_{\\!S}^{{\dagger}}$, is the same as
that of the random variable $\pazocal{D}(y)\psi\pazocal{D}(-y)$, where $y$ is
the outcome of the heterodyne on $\psi$. This means that in the shadow
tomography protocol we can skip the unitary operation altogether without
losing any data. In what follows we will therefore set $S=I$
(deterministically) without loss of generality.

With the above simplifications, one can see that
$f_{\mu,T}(u)=e^{-\frac{1}{2}\|u\|^{2}}$, so that
$\mathcal{M}=\mathcal{N}_{1}$ (cf. (40)). Therefore, the classical shadow will
have improper characteristic function

$\displaystyle\chi_{\hat{\rho}}(u)=e^{\frac{1}{4}\|u\|^{2}-iu^{\intercal}\Omega
x}\,.$

Once again, the characteristic function tensorises: given
$u=(u_{1},\dots,u_{m})^{\intercal}$ and $x=(x_{1},\dots,x_{m})^{\intercal}$:

$\displaystyle\chi_{\hat{\rho}}(u)=\prod_{j=1}^{m}\chi_{\hat{\rho}_{j}}(u_{j})$
(99)

where

$\displaystyle\chi_{\hat{\rho}_{j}}(u_{j})\coloneqq
e^{\frac{1}{4}\|u_{j}\|^{2}-iu_{j}^{\intercal}\Omega x_{j}}\,.$ (100)

Now, if we want to use the classical shadow to compute expectation values, we
can formally use Plancherel’s relation

$\displaystyle\hat{\rho}(O)$
$\displaystyle=\int\frac{d^{2m}u}{(2\pi)^{m}}\,\chi_{\hat{\rho}}(u)\,\chi_{O}(u)$
(101)
$\displaystyle=\int\frac{d^{2m}u}{(2\pi)^{m}}\,e^{+\frac{1}{4}\|u\|^{2}-iu^{\intercal}\Omega
x}\,\chi_{O}(u)\,.$

In order for this to make sense, we should make sure that not only $O$ is
trace class (instead of bounded), but also that $\chi_{O}(u)$ decays
sufficiently rapidly, for instance like $\sim
e^{-\frac{\lambda}{4}\|u\|^{2}}$, with $\lambda>1$. This decay is too fast —
but _barely_ too fast — to be useful in practice. For instance, if $O$ has a
finite expansion in the Fock basis then $\chi_{O}(u)\sim
p(u)e^{-\frac{1}{4}\|u\|^{2}}$ as $\|u\|\to\infty$, where $p(u)$ is some
polynomial of the entries of $u$. We get rid of the diverging Gaussian in
(101), but not of the diverging polynomial. In order to take care of this
issue, we make use of the approximations
$\widetilde{Z}_{\mathbf{n}_{1}\mathbf{n}_{2}}$ of the Schwartz operators
$\ket{\mathbf{n}_{1}}\\!\\!\bra{\mathbf{n}_{2}}$ as well as of the auxiliary
map $\widetilde{\mathcal{P}}_{M}$ and the corresponding matrices
$\widetilde{\sigma}_{A}^{(N)}(M)$ introduced in (75), resp. in (84). We
consider the matrices

$\displaystyle\hat{\rho}_{A}(M)\coloneqq\sum_{\mathbf{n}_{1},\mathbf{n}_{2}\in\\{0,...,M\\}^{|A|}}\,\hat{\rho}(\widetilde{Z}_{\mathbf{n}_{2}\mathbf{n}_{1}})\,\ket{\mathbf{n}_{1}}\\!\\!\bra{\mathbf{n}_{2}},\quad\text{where}$
$\displaystyle\hat{\rho}(\widetilde{Z}_{\mathbf{n}_{2}\mathbf{n}_{1}})\coloneqq\int_{\mathbb{R}^{2r}}\widetilde{\chi}_{\mathbf{n}_{2}\mathbf{n}_{1}}(u)\,e^{\frac{1}{4}\|u\|^{2}-iu^{\intercal}\Omega
x}\,\frac{d^{2r}u}{(2\pi)^{r}}\,$

where $\widetilde{\chi}_{\mathbf{n}_{2}\mathbf{n}_{1}}$ is defined in (76).
From Lemma 6,
$\mathbb{E}\big{[}\hat{\rho}_{A}(M)\big{]}=\widetilde{\mathcal{P}}_{M}(\rho_{A})$.
For $\alpha\geq 0$, $R>0$ and $r,M\in\mathbb{N}$, we introduce the matrix norm

$\displaystyle\widetilde{\Sigma}^{(\alpha)}_{r}(M,R)\coloneqq\left\|K_{r}(M,R)\right\|_{\infty}\,,$

where the entries of the $(M+1)^{r}\times(M+1)^{r}$ matrix $K_{r}(M,R)$ are
defined by

$\displaystyle K_{r}(M,R)_{\mathbf{n_{1}},\mathbf{n_{2}}}$ $\displaystyle\
\,\coloneqq\left((1\\!+\\!|\mathbf{n}_{1}|)(1\\!+\\!|\mathbf{n}_{2}|)\right)^{\frac{\alpha}{2}}\int_{|u|\leq
R}\\!\\!|\widetilde{\chi}_{\mathbf{n}_{2}\mathbf{n}_{1}}(u)|\,e^{\frac{1}{4}\|u\|^{2}}\frac{d^{2r}u}{(2\pi)^{r}}\,,$

and $\mathbf{n_{1}},\mathbf{n_{2}}\in\\{0,\ldots,M\\}^{r}$. With these
definitions, we can write the inequality

$\displaystyle\left\|H_{r}^{\frac{\alpha}{2}}\hat{\rho}_{A}(M)\,H_{r}^{\frac{\alpha}{2}}\right\|_{\infty}\leq\widetilde{\Sigma}^{(\alpha)}_{r}(M,R)\,.$
(102)

As before, we then denote

$\displaystyle\sigma_{A}^{(N)}(M)\coloneqq\frac{1}{N}\sum_{i=1}^{N}\,\hat{\rho}^{(i)}_{A}(M)\,,$

where $\hat{\rho}^{(i)}_{A}(M)$ are i.i.d. random matrices of law the one of
$\hat{\rho}_{A}(M)$. Our second main result is stated below:

######  Theorem 11.

With the notation introduced above, given $0\leq\alpha<n$, $r\in\mathbb{N}$
such that $E^{(n)}_{r}\coloneqq\max_{|A|\leq
r}\operatorname{Tr}(\rho_{A}H_{r}^{n})<\infty$, and $0<\eta<R$, set

$\displaystyle M\\!\coloneqq\min\\!\Big{\\{}$ $\displaystyle
M^{\prime}\\!\\!\in\\!\mathbb{N}:\eta^{2}>2M^{\prime 2}$ (103) $\displaystyle
2(1\\!+\\!M^{\prime})^{\frac{\alpha-n}{2}}\\!E_{r}^{(n)}+\delta_{0}\big{(}\eta,M^{\prime}\\!\\!,\tfrac{\alpha}{2},r\big{)}\leq\tfrac{\epsilon}{2}\Big{\\}}\,,$

with $\delta_{0}$ defined as in Proposition 8. Then for any

$\displaystyle N\geq\inf_{0<\eta<R}$
$\displaystyle\frac{(M+1)^{2r}}{3\epsilon^{2}}\log\left(\frac{2\big{[}m(M+1)\big{]}^{r}}{\delta}\right)$
(104) $\displaystyle\Big{[}24\,\widetilde{\Sigma}^{(\alpha)}_{r}(M\\!,R)^{2}$
$\displaystyle\
+4\Big{(}\widetilde{\Sigma}_{r}^{(\alpha)}(M\\!,R)+E_{r}^{(\alpha)}+\delta_{0}(\eta,M,\tfrac{\alpha}{2},r)\Big{)}\epsilon\Big{]}\,,$

we have that for any region $A$ of size $|A|\leq r$, it holds that

$\displaystyle\left\|\sigma_{A}^{(N)}(M)-\rho_{A}\right\|^{(\alpha)}_{1}\leq\epsilon$
(105)

with probability at least $1-\delta$. Similarly, for

$\displaystyle N\geq\inf_{0<\eta<R}$
$\displaystyle\frac{(M+1)^{2r}}{3\epsilon^{2}}\log\left(\frac{2L(M+1)^{r}}{\delta}\right)$
(106) $\displaystyle\Big{[}24\,\widetilde{\Sigma}^{(\alpha)}_{r}(M\\!,R)^{2}$
$\displaystyle\
+4\Big{(}\widetilde{\Sigma}_{r}^{(\alpha)}(M\\!,R)+E_{r}^{(\alpha)}+\delta_{0}(\eta,M,\tfrac{\alpha}{2},r)\Big{)}\epsilon\Big{]}\,,$

we have that for any set of $L$ observables $O_{j}$ on regions $A_{j}$ of size
at most $r$ and with
$\left\|H_{r}^{-\frac{\alpha}{2}}O_{j}H_{r}^{-\frac{\alpha}{2}}\right\|_{\infty}\leq
1$, it holds that

$\displaystyle\max_{j}\big{|}\operatorname{Tr}\big{[}O_{j}\,(\sigma_{A_{j}}^{(N)}(M)-\rho_{A_{j}})\big{]}\big{|}\leq\epsilon$
(107)

with probability at least $1-\delta$.

###### Proof.

By construction, we have almost surely that
$\|H_{r}^{\frac{\alpha}{2}}\hat{\rho}_{A}(M)H_{r}^{\frac{\alpha}{2}}\|_{\infty}\leq\widetilde{\Sigma}^{(\alpha)}_{r}(M,R)$.
Moreover, we have that

$\displaystyle\|H_{r}^{\frac{\alpha}{2}}\widetilde{\mathcal{P}}_{M}(\rho_{A})H_{r}^{\frac{\alpha}{2}}\|_{\infty}$
$\displaystyle\leq\|\widetilde{\mathcal{P}}_{M}(\rho_{A})\|_{1}^{(\alpha)}$
$\displaystyle\leq\|(\widetilde{\mathcal{P}}_{M}-\mathcal{P}_{M})(\rho_{A})\|_{1}^{(\alpha)}+\|\mathcal{P}_{M}(\rho_{A})\|_{1}^{(\alpha)}$
$\displaystyle\leq\delta_{0}(\eta,M,\frac{\alpha}{2},r)+E_{r}^{(\alpha)}$

for $\eta^{2}>2M^{2}$, by Lemma 7. Therefore almost surely
$\|H_{r}^{\frac{\alpha}{2}}(\hat{\rho}_{A}(M)-\widetilde{\mathcal{P}}_{M}(\rho_{A}))H_{r}^{\frac{\alpha}{2}}\|_{\infty}\leq
E_{r}^{(\alpha)}+\delta_{0}(\eta,M,\frac{\alpha}{2},r)+\widetilde{\Sigma}^{(\alpha)}_{r}(M,R)$.
Finally,
$\|\mathbb{E}[(H_{r}^{\frac{\alpha}{2}}\hat{\rho}_{A}(M)H_{r}^{\frac{\alpha}{2}})^{2}]\|_{\infty}\leq\widetilde{\Sigma}_{r}^{(\alpha)}(M,R)^{2}$.
We hence use the matrix Bernstein inequality (LABEL:matrixBern) to get in
terms of a constant $C=C(N,\epsilon,R,\alpha,r)$ with

$C\coloneqq\frac{3N\epsilon^{2}(M+1)^{-2r}}{6\widetilde{\Sigma}_{r}^{(\alpha)}(M,R)^{2}+2(\widetilde{\Sigma}_{r}^{(\alpha)}(M,R)+\delta_{0}(\eta,M,\frac{\alpha}{2},r)+E_{r}^{(\alpha)})\epsilon}$

the estimate

$\displaystyle\mathbb{P}\Big{(}\big{\|}\widetilde{\sigma}_{A}^{(N)}(M)-\widetilde{\mathcal{P}}_{M}(\rho_{A})\big{\|}^{(\alpha)}_{1}\geq{\epsilon}\Big{)}$
(108) $\displaystyle\leq
2(M+1)^{r}\,\operatorname{exp}\left(-C(N,\epsilon,R,\alpha,r)\right)\,.$

Therefore, by a union bound, we get that

$\displaystyle\mathbb{P}\Big{(}\exists A,\,|A|\leq
r:\,\big{\|}\widetilde{\sigma}_{A}^{(N)}(M)-\widetilde{\mathcal{P}}_{M}(\rho_{A})\big{\|}^{(\alpha)}_{1}\geq{\epsilon}\Big{)}$
$\displaystyle\leq
2\big{[}m(M+1)\big{]}^{r}\,\operatorname{exp}\left(-C(N,\epsilon,R,\alpha,r)\right)\,.$

Next, we use Proposition 8 and choose
$M=\min\\{M^{\prime}\in\mathbb{N}:\,2(1+M^{\prime})^{\frac{\alpha-n}{2}}\,E_{r}^{(n)}+\delta_{0}\big{(}\eta,M^{\prime},\frac{\alpha}{2},r\big{)}\leq\frac{\epsilon}{2}\\}$,
so that
$\|\widetilde{\mathcal{P}}_{M}(\rho_{A})-\rho_{A}\|_{1}^{(\alpha)}\leq\epsilon/2$,
and

$\displaystyle\mathbb{P}\Big{(}\exists A,\,|A|\leq
r:\,\big{\|}\widetilde{\sigma}_{A}^{(N)}(M)-\rho_{A}\big{\|}^{(\alpha)}_{1}\geq{\epsilon}\Big{)}$
$\displaystyle\leq\mathbb{P}\Big{(}\exists A,\,|A|\leq
r:\,\big{\|}\widetilde{\sigma}_{A}^{(N)}(M)-\widetilde{\mathcal{P}}_{M}(\rho_{A})\big{\|}^{(\alpha)}_{1}\geq\frac{\epsilon}{2}\Big{)}\,.$

Therefore, choosing $N$ as in the statement of the theorem, we obtain that the
probability that on any subset $A$ of at most $r$ modes
$\|\widetilde{\sigma}_{A}^{(N)}(M)-\rho_{A}\big{\|}^{(\alpha)}_{1}\geq{\epsilon}$
is at least $1-\delta$.

The results for a fixed number of $L$ observables $O_{1},\dots,O_{L}$
supported on regions $A_{1},\dots,A_{L}$ of size at most $r$ follow after
replacing the above union bounds over regions $A$ by a union bound over the
observables $O_{j}$:

$\displaystyle\mathbb{P}\Big{(}\exists
j:\,|\operatorname{Tr}(O_{j}\widetilde{\sigma}_{A_{j}}^{(N)}(M))-\operatorname{Tr}(O_{j}\,\widetilde{\mathcal{P}}_{M}(\rho_{A_{j}}))|\geq{\epsilon}\Big{)}$
$\displaystyle\leq
2L\,\big{[}(M+1)\big{]}^{r}\,\operatorname{exp}\left(-C(N,\epsilon,R,\alpha,r)\right)$

and the result follows after replacing $m^{r}$ by $L$ in the estimate for $N$.

∎

###  V-C Comparison to related work

In a concurrent and independent work [23], the authors also developed a shadow
tomography protocol for bosonic systems. Their framework is equivalent to ours
in the homodyne and heterodyne settings. Let us consider the former in the
one-mode setting: the shadows in [23] are constructed in the Fock basis by
taking $\check{\rho}_{mn}\coloneqq e^{i(m-n)\theta}f_{mn}(x)$, after making a
random rotation $\theta$ and obtaining $x$ from a measurement in the position
axis. The functions $f_{mn}=f_{nm}$ are the so-called pattern functions [38,
74] and are defined using Fock state wavefunctions $\psi_{m}$ ($m^{th}$ energy
eigenstate of the harmonic oscillator) and $\phi_{n}$ ($n^{th}$ non-
normalizable solution of the Schrödinger equation of a harmonic oscillator) as

$\displaystyle f_{mn}(x)\coloneqq\frac{\partial}{\partial
x}(\psi_{m}(x)\phi_{n}(x)),\qquad n\geq m\,.$

Then, they construct the following estimator from $N$ copies of the unknown
state:

$\displaystyle\check{\sigma}^{(N)}_{mn}\coloneqq\frac{1}{N}\sum_{i=1}^{N}\,\check{\rho}^{(i)}_{mn}=\frac{1}{N}\sum_{i=1}^{N}e^{i(m-n)\theta_{i}}f_{mn}(x_{\theta_{i}})\,,$

where the angles $\theta_{i}$ are picked uniformly at random, and
$x_{\theta}=x^{\prime}_{\theta}\oplus
x^{\prime\prime}_{\theta}=\sqrt{2\pi}x^{\prime}\oplus\frac{1}{\sqrt{2\pi}}x^{\prime\prime}$
in the notations of Section V-A. It turns out that these pattern functions can
be equivalently defined in terms of their Fourier transforms (eq. (31) of
[75], see also eqs. (20)-(21) of [76]): for $n\geq m$:

$\displaystyle\widetilde{f}_{mn}(t)$ $\displaystyle\coloneqq\int
e^{-itx}\,f_{mn}(x)\,dx$
$\displaystyle=\pi(-i)^{n-m}\,\sqrt{\frac{2^{m-n}m!}{n!}}\,|t|\,t^{n-m}e^{-\frac{t^{2}}{4}}\,L_{m}^{(n-m)}\Big{(}\frac{t^{2}}{2}\Big{)}\,.$

Therefore $|\widetilde{f}_{mn}(x_{\theta})|$ corresponds to the coordinate
$\Sigma_{1}^{(0)}(M)_{mn}$ for $m,n\leq M$. As observed in [59], the pattern
functions have been already studied for tomography purposes (see e.g. Lemma
7.1 in [61]). In particular,

$\displaystyle\sum_{0\leq j\leq k\leq
M}\|f_{kj}\|^{2}_{\infty}=\mathcal{O}(M^{7/3})\,.$

These bounds can be directly used to get control over
$\|\Sigma_{1}^{(0)}(M)\|_{\infty}$. That way, we recover Theorem III.1 of
[23]. Slightly better bounds can be achieved using estimates in Lemma 4 of
[76]. It is also worth observing that, from the proof of Theorem 9, one can
extend these results to more realistic models of homodyning with finitely
squeezed resources.

##  VI Learning non-linear functionals of the states

So far, we considered properties of the quantum system which could be related
to local linear functionals of the unknown state. In finite dimensions, a
simple trick permits the estimation of non-linear functionals, e.g. the
entropy of entanglement [6, 77]. Here, we show that the technique developed in
these works combined with recent energy-constrained continuity bounds provide
us with a similar extension. For sake of conciseness, we will only consider
the entropy of a reduced CV state over $r\leq m$ modes of the $m$-mode state:

$\displaystyle S(\rho_{A})\coloneqq-\operatorname{Tr}(\rho_{A}\ln\rho_{A})\,,$

where we denote by $A$ the set of those modes, so that $|A|=r$. We also assume
that the unknown state $\rho$ has locally finite energy:
$\operatorname{Tr}((I+N^{(j)})\rho_{j})\coloneqq E<\infty$, where $N^{(j)}$
corresponds to the number operator on mode $j$, and where $\rho_{j}$ is the
reduced state on mode $j\in A$. Here, we also restrict ourselves to the
shadows constructed in our local homodyne detection scheme of (V-A), although
similar conclusions can be drawn from shadows arising from a local heterodyne
detection and states with higher moment constraints. Given the shadow
$\sigma_{A}^{(N)}(M)$ arising from the homodyne scheme of (V-A), and
$d_{p}\in\mathbb{N}$, we denote the matrix polynomial

$\displaystyle H^{(d_{p})}(\sigma_{A}^{(N)}(M))$
$\displaystyle\coloneqq-\operatorname{Tr}(\sigma_{A}^{(N)}(M)-P_{M})$
$\displaystyle\qquad\qquad+\sum_{k=2}^{d_{p}}\frac{\operatorname{Tr}\big{[}(P_{M}-\sigma_{A}^{(N)}(M))^{k}\big{]}}{k(k-1)}\,.$

######  Theorem 12.

With the notation of the previous paragraph, we have that for any $\epsilon>0$
and $M\in\mathbb{N}$,

$\displaystyle\mathbb{P}\Big{(}\exists A,\,|A|\leq
r:\,|S(\rho_{A})-H^{\big{(}\big{\lceil}\frac{3(M+1)^{r}}{\epsilon}\big{\rceil}\big{)}}(\sigma_{A}^{(N)}(M))|\geq\epsilon\Big{)}$
$\displaystyle\leq
2\big{[}m(M+1)\big{]}^{r}\,\operatorname{exp}\left(-\frac{3N\epsilon^{\prime
2}(M+1)^{-2r}}{6\Sigma_{r}^{(0)}(M)^{2}+2(\Sigma_{r}^{(0)}(M)+1)\epsilon^{\prime}}\right)\,,$

where
$\epsilon^{\prime}=\frac{\epsilon^{2}}{12(M+1)^{r}e}2^{-4(M+1)^{r}/\epsilon}$.

###### Proof.

By [78, Theorem 3], we have that for all $\gamma\in[0,rE/(1+rE)]$, and any
state $\rho_{A}^{\prime}$ on subsystem $A$ with
$\|\rho_{A}-\rho_{A}^{\prime}\|_{1}\leq 2\gamma$,

$\displaystyle|S(\rho_{A})-S(\rho_{A}^{\prime})|\leq
h(\gamma)+rE\,h\big{(}\gamma/(rE)\big{)}\,,$ (109)

where $h(.)$ denotes the binary entropy. Next, we pick
$\rho_{A}^{\prime}\coloneqq\mathcal{P}_{M}(\rho_{A})/\operatorname{Tr}(\mathcal{P}_{M}(\rho_{A}))$.
By Proposition 8, we have that, assuming $1+M>4r^{2}E^{2}$,

$\displaystyle\|\rho_{A}-\rho_{A}^{\prime}\|_{1}=\frac{\|\operatorname{Tr}(\mathcal{P}_{M}(\rho_{A}))\rho_{A}-\mathcal{P}_{M}(\rho_{A})\|_{1}}{\operatorname{Tr}(\mathcal{P}_{M}(\rho_{A}))}\,$
$\displaystyle\leq\frac{|\operatorname{Tr}(\mathcal{P}_{M}(\rho_{A}))-1|+\|\rho_{A}-\mathcal{P}_{M}(\rho_{A})\|_{1}}{1-|\operatorname{Tr}(\mathcal{P}_{M}(\rho_{A}))-1|}$
$\displaystyle\leq\frac{2\|\rho_{A}-\mathcal{P}_{M}(\rho_{A})\|_{1}}{1-\|\mathcal{P}_{M}(\rho_{A})-\rho_{A}\|_{1}}\,$
$\displaystyle\leq\frac{4{rE}}{(1+M)^{\frac{1}{2}}-2rE}\equiv 2\gamma\,.$

Next, we approximate $S(\rho_{A}^{\prime})$ in terms of
$S(\mathcal{P}_{M}(\rho_{A}))$, where the entropy of a sub-normalised positive
trace class operator $A$ is defined as $S(A)\coloneqq-\operatorname{Tr}A\ln
A$: denoting $\operatorname{Tr}(\mathcal{P}_{M}(\rho_{A}))\coloneqq\lambda\leq
1$,

$\displaystyle S(\rho_{A}^{\prime})-S(\mathcal{P}_{M}(\rho_{A}))$
$\displaystyle=-\operatorname{Tr}\frac{\mathcal{P}_{M}(\rho_{A})}{\lambda}\ln\frac{\mathcal{P}_{M}(\rho_{A})}{\lambda}+\operatorname{Tr}\mathcal{P}_{M}(\rho_{A})\ln\mathcal{P}_{M}(\rho_{A})$
$\displaystyle=\frac{1-\lambda}{1-(1-\lambda)}\,S(\mathcal{P}_{M}(\rho_{A}))+\ln\lambda$
$\displaystyle\leq\frac{\|\rho_{A}-\mathcal{P}_{M}(\rho_{A})\|_{1}}{1-\|\rho_{A}-\mathcal{P}_{M}(\rho_{A})\|_{1}}\,S(\mathcal{P}_{M}(\rho_{A}))$
$\displaystyle\leq\gamma\,S(\mathcal{P}_{M}(\rho_{A}))$ $\displaystyle\leq
r\ln(M+1)\,\gamma\,.$ (110)

Next, we use a polynomial approximation of $S(\mathcal{P}_{M}(\rho_{A}))$ as
was already done in [77]:

$\displaystyle H^{(d_{p})}(\mathcal{P}_{M}(\rho_{A}))$
$\displaystyle\coloneqq-\operatorname{Tr}(\mathcal{P}_{M}(\rho_{A})-P_{M})+\sum_{k=2}^{d_{p}}\frac{\operatorname{Tr}\big{[}(P_{M}-\mathcal{P}_{M}(\rho_{A}))^{k}\big{]}}{k(k-1)}$

where $d_{p}$ is the truncation degree which we will choose later. Since
$P_{M}$ projects onto a subspace of dimension $(M+1)^{r}$, a simple extension
of the proof leading to [77, Equation (K46)] gives us

$\displaystyle|S(\mathcal{P}_{M}(\rho_{A}))-H^{(d_{p})}(\mathcal{P}_{M}(\rho_{A}))|\leq\frac{(M+1)^{r}}{d_{p}}\,.$
(111)

Now, we re-express the function $H^{(d_{p})}(\mathcal{P}_{M}(\rho_{A}))$ as a
linear function of $(\mathcal{P}_{M}(\rho_{A}))^{\otimes d_{p}}$:

$\displaystyle H^{(d_{p})}(\mathcal{P}_{M}(\rho_{A}))$ (112)
$\displaystyle=-1+(M\\!+\\!1)^{r}+\sum_{k=2}^{d_{p}}\frac{1}{k(k-1)}\,\sum_{j=0}^{k}\binom{k}{j}\,(-1)^{j}\operatorname{Tr}((\mathcal{P}_{M}(\rho_{A}))^{j})$
$\displaystyle=-1+(M\\!+\\!1)^{r}$
$\displaystyle\qquad+\sum_{j=0}^{d_{p}}\frac{1}{j!}\,\operatorname{Tr}(S_{j}(\mathcal{P}_{M}(\rho_{A}))^{\otimes
j})\,(-1)^{j}\\!\\!\sum_{k=\max\\{2,j\\}}^{d_{p}}\frac{(k\\!-\\!2)!}{(k\\!-\\!j)!}$
$\displaystyle\equiv-1+(M\\!+\\!1)^{r}+\sum_{j=0}^{d_{p}}\frac{1}{j!}\,C_{j}\,\operatorname{Tr}(S_{j}(\mathcal{P}_{M}(\rho_{A}))^{\otimes
j})\,,$

where $S_{j}$ is a generalised swap operator over $j$ subsystems. Let us now
chose $\sigma_{A}^{(N)}(M)$ to be the classical shadow arising from the
homodyne scheme of (V-A). We have, similarly to [77, Lemma 11]

$\displaystyle|H^{(d_{p})}(\mathcal{P}_{M}(\rho_{A}))-H^{(d_{p})}(\sigma_{A}^{(N)}(M))|$
(113)
$\displaystyle\leq\sum_{j=0}^{d_{p}}\,\frac{1}{j!}|C_{j}|\,\|(\mathcal{P}_{M}(\rho_{A}))^{\otimes
j}-(\sigma_{A}^{(N)}(M))^{\otimes j}\|_{1}\,.$

Next, we estimate the trace distances in the above summand. By an arbitrary
labelling $A_{1},\dots,A_{j}$ of the subsystems so that
$(\mathcal{P}_{M}(\rho_{A}))^{\otimes
j}\coloneqq\rho_{A_{1}}^{\prime}\otimes\dots\otimes\rho_{A_{j}}^{\prime}$ and
$(\sigma_{A}^{(N)}(M))^{\otimes
j}\coloneqq\sigma_{A_{1}}\otimes\dots\otimes\sigma_{A_{j}}$, we have

$\displaystyle\|(\mathcal{P}_{M}(\rho_{A}))^{\otimes
j}-(\sigma_{A}^{(N)}(M))^{\otimes j}\|_{1}$ (114)
$\displaystyle=\|\rho_{A_{1}}^{\prime}\otimes\dots\otimes\rho_{A_{j}}^{\prime}-\sigma_{A_{1}}\otimes\dots\otimes\sigma_{A_{j}}\|_{1}$
$\displaystyle=\big{\|}(\rho_{A_{1}}^{\prime}-\sigma_{A_{1}})\otimes\rho_{A_{2}}^{\prime}\otimes\dots\otimes\rho_{A_{j}}^{\prime}$
$\displaystyle\qquad+\sigma_{A_{1}}\otimes(\rho_{A_{2}}^{\prime}\otimes\dots\otimes\rho_{A_{j}}^{\prime}-\sigma_{A_{2}}\otimes\dots\otimes\sigma_{A_{j}})\big{\|}_{1}$
$\displaystyle\leq\|\rho_{A_{1}}^{\prime}-\sigma_{A_{1}}\|_{1}\|\rho_{A_{2}}^{\prime}\|_{1}\dots\|\rho_{A_{j}}^{\prime}\|_{1}$
$\displaystyle\qquad+\|\sigma_{A_{1}}\|_{1}\,\big{\|}\rho_{A_{2}}^{\prime}\otimes\dots\otimes\rho_{A_{j}}^{\prime}-\sigma_{A_{2}}\otimes\dots\otimes\sigma_{A_{j}}\big{\|}_{1}$
$\displaystyle\leq\|\rho_{A_{1}}^{\prime}-\sigma_{A_{1}}\|_{1}$
$\displaystyle\qquad+(1+\|\rho_{A_{1}}^{\prime}\\!-\\!\sigma_{A_{1}}\|_{1})\,\big{\|}\rho_{A_{2}}^{\prime}\\!\otimes\\!\dots\\!\otimes\\!\rho_{A_{j}}^{\prime}-\sigma_{A_{2}}\\!\otimes\\!\dots\\!\otimes\\!\sigma_{A_{j}}\big{\|}_{1}.$

Now, we recall that the homodyne tomography protocol provides us with the
following concentration bound borrowed from (V-A)

$\displaystyle\mathbb{P}\Big{(}\exists A,\,|A|\leq
r:\,\big{\|}\sigma_{A}^{(N)}(M)-\mathcal{P}_{M}(\rho_{A})\big{\|}_{1}\geq{\epsilon^{\prime}}\Big{)}$
$\displaystyle\leq 2\big{[}m(M+1)\big{]}^{r}\,$
$\displaystyle\quad\times\operatorname{exp}\left(-\frac{3N\epsilon^{\prime
2}(M+1)^{-2r}}{6\Sigma_{r}^{(0)}(M)^{2}+2(\Sigma_{r}^{(0)}(M)+1)\epsilon^{\prime}}\right)\,.$
(115)

Plugging this bound into (114), we have that with high probability

$\displaystyle\|(\rho_{A}^{\prime})^{\otimes j}-(\sigma_{A}^{(N)}(M))^{\otimes
j}\|_{1}$
$\displaystyle\leq\epsilon^{\prime}+\big{(}1+\epsilon^{\prime}\big{)}\|(\mathcal{P}_{M}(\rho_{A}))^{\otimes
j-1}-(\sigma_{A}^{(N)}(M))^{\otimes j-1}\|_{1}\,.$

After iterating the procedure $j$ times, we find

$\displaystyle\|(\mathcal{P}_{M}(\rho_{A}))^{\otimes
j}-(\sigma_{A}^{(N)}(M))^{\otimes j}\|_{1}$
$\displaystyle\leq\epsilon^{\prime}\sum_{i=0}^{j-1}\big{(}1+\epsilon^{\prime}\big{)}^{i}\leq\big{(}1+\epsilon^{\prime}\big{)}^{j}-1\,.$

Into (113), we have that

$\displaystyle|H^{(d_{p})}(\mathcal{P}_{M}(\rho_{A}))-H^{(d_{p})}(\sigma_{A}^{(N)}(M))|$
$\displaystyle\leq\sum_{j=0}^{d_{p}}\frac{1}{j!}|C_{j}|\,\Big{(}\Big{(}1+\epsilon^{\prime}\Big{)}^{j}-1\Big{)}$
$\displaystyle\leq\big{(}\big{(}1+\epsilon^{\prime}\big{)}^{d_{p}}-1\big{)}\,\sum_{j=0}^{d_{p}}\,\frac{|C_{j}|}{j!}$
$\displaystyle\leq\big{(}\big{(}1+\epsilon^{\prime}\big{)}^{d_{p}}-1\big{)}\,2^{d_{p}}\,.$
(116)

Then, combining (109), (110), (111) and (116), we end up with

$\displaystyle|S(\rho_{A})-H^{(d_{p})}(\sigma_{A}^{(N)}(M))|$
$\displaystyle\leq
h(\gamma)+rEh\Big{(}\frac{\gamma}{rE}\Big{)}+2r\ln(M+1)\gamma$
$\displaystyle\quad+\frac{(M+1)^{r}}{d_{p}}+\big{(}\big{(}1+\epsilon^{\prime}\big{)}^{d_{p}}-1\big{)}\,2^{d_{p}}\,.$

Now, the first three terms on the right-hand side above are small for $M$
large enough, the fourth term is small for $d_{p}$ large enough, and the last
term is small for $N$ large enough. The result follows after choosing $M$ so
that the first three terms are smaller than $\epsilon/3$,
$d_{p}=\Big{\lceil}\frac{3(M+1)^{r}}{\epsilon}\Big{\rceil}$ so that the fourth
term is smaller than $\epsilon/3$, and finally
$\epsilon^{\prime}=\frac{\epsilon^{2}}{12(M+1)^{r}e}2^{-4(M+1)^{r}/\epsilon}$
so that the last term is smaller than $\epsilon/3$. ∎

##  VII Discretisation scheme via quasi Monte-Carlo integration

Since the objects we are manipulating in this paper are defined on continuous
metric spaces, we need to explain how to devise an efficient description
$\sigma^{\operatorname{dis}}_{N}$ of the estimator state $\sigma_{N}$ in terms
of a number of discrete parameters which scales at most polynomially with the
number $m$ of modes and such that
$\operatorname{Tr}(O(\sigma^{\operatorname{dis}}_{N}-\sigma_{N}))$ can be
controlled for the observables $O$ for which we want to learn the average
$\operatorname{Tr}(\rho O)$. One natural strategy consists e.g. in
approximating the coefficients
$\operatorname{Tr}(\widetilde{Z}_{\mathbf{n}_{1}\mathbf{n}_{2}}\sigma_{N})$ of
the matrix $\widetilde{\mathcal{P}}_{{M}}(\sigma_{N})$. Since the
characteristic function of $\widetilde{Z}_{\mathbf{n}_{1}\mathbf{n}_{2}}$ is
compactly supported in a ball $B(0,R)$ around the origin, it suffices to store
its values on a net of small enough mesh.

In order to evaluate traces
$\operatorname{Tr}(\tilde{Z}_{\mathbf{n}_{1}\mathbf{n}_{2}}\sigma_{N}),$ in an
efficient way, we shall employ quasi Monte-Carlo techniques using the isometry
(50). We start by recalling that an $L^{1}$ function $u:\Omega\to\mathbb{R}$
is of bounded variation if

$\displaystyle\operatorname{TV}(u)$
$\displaystyle\coloneqq\sup\left\\{\int_{\Omega}u(x)\operatorname{div}(\phi(x));\phi\in
C_{c}^{1}(\Omega,\mathbb{R}^{n});\|\phi\|_{\infty}\leq 1\right\\}$
$\displaystyle\quad<\infty\,$

and denote by $\operatorname{TV}(K,\mathbb{R})$ the space of functions of
bounded variation. In particular, for $u\in C^{1}(\Omega,\mathbb{R}^{n})$ one
just has

$\displaystyle\operatorname{TV}(u)=\int_{\Omega}|\nabla u(x)|\ dx.$

For a complex valued function to be of bounded variation, we accordingly
require that both its real and imaginary part are of bounded variation. Let
$K=[-l_{j},l_{j}]^{m}$ be a compact set and $\ell:[0,1]^{m}\to K$ a linear
map. We set $S=\\{t_{k}\\}_{k\in\mathbb{N}},$ where $t_{k}\in[0,1]^{m}$ is a
_Halton sequence_ , see [79] for a definition, in a pairwise prime basis
$b_{1},..,b_{m}.$ We then introduce for $f\in\operatorname{TV}(K,\mathbb{R})$,
the numerical integral

$I_{K}(f,k)=\frac{1}{k}\sum_{i=1}^{k}f(\ell(t_{i})).$

It then follows that there exists a constant $C(b_{1},...,b_{d})$ independent
of both $f$ and $N$ such that, see [79, Theorem $2.11$],

$\displaystyle\Big{|}I_{K}(f,k)-\int_{K}f(x)\ dx\Big{|}$
$\displaystyle\leq\operatorname{TV}(f)\,C(\ell(t_{1}),...,\ell(t_{m}))\,\frac{\log(k)^{m}}{k}.$

Thus, this technique allows us to approximate high-dimensional integrals with
errors that depend only very mildly on the number of modes $m.$ We can then
apply this construction to

$f(x)=\frac{\chi_{|\mathbf{n}_{1}\rangle\langle\mathbf{n}_{2}|}(x)\chi_{\sigma_{N}}(x)\prod_{j=1}^{m}\xi_{\eta,R}(x_{j})}{(2\pi)^{m}},$

with $K=[-R,R]^{m}.$

##  VIII Examples

In this section, we test our homodyne classical shadow tomography method by
means of numerical simulations.

Figure 1: The quantum characteristic function of the actual state is depicted
in the first plot on the top left, while that of the reconstructed state is
given in the plot beside it (for $N=50$) and the one below it (for $N=1000$).
For $N=50$ (above) the quantum characteristic function gets almost perfectly
approximated within the unit square but diverges outside of it. On the right,
we show the histogram of numerical measurement outcomes generating our
approximation. For $N=1000$ (below) the quantum characteristic function is
well approximated inside the square $[-2,2]^{2}$. In the center, we show the
histogram of numerical measurement outcomes generating our approximation. On
the bottom right, we show the variance of the approximated characteristic
function to the real characteristic function inside $[-2,2]^{2}.$

Figure 2: Real and imaginary parts of the quantum characteristic functions of
$\ket{0}_{\operatorname{Cat}}$(top line) and $\ket{1}_{\operatorname{Cat}}$
(center line), respectively, reconstructed with $N=200$, with reconstructed
function (left) and true one (right) on $[-2,2]$ for $\alpha=(1,1)$. At the
bottom, we illustrate the numerical measurement outcomes for the above
reconstructions ($\ket{0}_{\operatorname{Cat}}$ left,
$\ket{1}_{\operatorname{Cat}}$ right) and the variance (for
$\ket{0}_{\operatorname{Cat}}$) for fixed $\alpha=(1,1)$ and varying number of
measurements (left) and fixed number of measurements, $N=200$, and varying
$\alpha$ (right).

Figure 3: Real and imaginary part of the quantum characteristic function of
$\ket{0}_{\operatorname{Cat}}$(top) with $\alpha=10\times(1,1)$ and
$\ket{1}_{\operatorname{Cat}}$ (bottom), with
$\alpha=\frac{1}{10}\times(1,1)$, respectively reconstructed with $N=200$ with
reconstructed function (left) and true one (right)) on $[-2,2]$.

Figure 4: On the left, we see the correlation coefficient
$\log(\gamma_{1,\lceil m/2\rceil})$ (in the bottom row with $h_{XX}$ replaced
by $h_{XX}+\frac{1}{2m}Q^{2}$, where $Q$ is a GUE Matrix and $m$ is the number
of oscillators). Here $\kappa$ refers to the interaction strength between
neighboring oscillators. In the center and on the right, we see the same
correlation coefficient recovered from the shadow with 1000 and 100,000
numerical experiments sampled from (124).

Figure 5: Plot of $\chi_{1,N/2}((x_{1},0),(x_{N/2},0))$ with 1000 _numerical
measurements_ sampled from the multivariate Gaussian distribution (124). From
left to right we see simulations for
$(\kappa,N)=(0.99,1000),(0.99,10),(0.99,4)$ respectively, with recovered
characteristic functions in the top row and true ones in the bottom row. We
see that correlations between particles are almost perfectly reproduced.

###  VIII-A Vacuum state

We start with the simplest example $\rho=\ket{0}\\!\\!\bra{0}$ to fix ideas.
In this case $|\langle 0|\alpha\rangle|^{2}=e^{-|\alpha|^{2}/2}$. Measurements
of $\alpha$ are then standard normally distributed
$\alpha\sim\mathcal{N}_{\mathbb{R}^{2}}(0,1).$ Considering then $N=50$
experimental realisations, $\alpha_{1},...,\alpha_{N}\in\mathbb{R}^{2},$ the
trial quantum characteristic function is

$\chi_{\sigma_{N}}(u)=\frac{e^{\frac{|u|^{2}}{4}}}{N}\sum_{i=1}^{N}e^{-iu^{T}\Omega\alpha_{i}}.$

The function is illustrated in Figure VIII. Reduced particle characteristic
functions for up to $1000$ oscillators are illustrated in Figure VIII. This
yields an approximation of the quantum characteristic function on compact sets
around zero. Due to the exponentially increasing function
$e^{\frac{|u|^{2}}{4}}$, this approximation is only local for fixed $N$. We
also plot the variance of the true characteristic function $\chi$ vs the
reconstructed characteristic function $\tilde{\chi}$ at $N$ grid points
$(x_{n})_{i=1}^{N}$ with $x_{n}\in D$ in some domain $D$

$\mathbf{V}_{N,D}=\frac{\sum_{i=1}^{N}|\chi_{i}(x_{n})-\tilde{\chi}_{i}(x_{n})|^{2}}{\operatorname{vol}(D)N}.$

###  VIII-B Cat qubit states

For our second example, we consider a non-Gaussian, one-mode pure quantum
state called cat state. These states are used in quantum error correction [21,
24]. Given a coherent state $|\alpha\rangle$, we denote the cat states

$\displaystyle|+\rangle_{\operatorname{Cat}}=\frac{|\alpha\rangle+|-\alpha\rangle}{\mathcal{N}_{+}},\quad|-\rangle_{\operatorname{Cat}}=\frac{|\alpha\rangle-|-\alpha\rangle}{\mathcal{N}_{-}}\,,$

with normalisation constants

$\displaystyle\mathcal{N}_{\pm}\coloneqq\sqrt{(\bra{\alpha}\pm\bra{-\alpha})(\ket{\alpha}\pm\ket{-\alpha})}$
$\displaystyle=\sqrt{2\pm(\langle\alpha|-\alpha\rangle+\langle-\alpha|\alpha\rangle)}=\sqrt{2(1\pm
e^{-2|\alpha|^{2}})}.$ (117)

Here, $\langle x|\pm\rangle_{\operatorname{Cat}}=\frac{\langle
x|\alpha\rangle\pm\langle x|-\alpha\rangle}{\mathcal{N}_{\pm}},$ where
$\langle
x|\alpha\rangle=e^{-\frac{|x|^{2}+|\alpha|^{2}}{4}+\frac{\alpha\bar{x}}{2}}$
such that
$\langle-\alpha|\alpha\rangle=\langle\alpha|-\alpha\rangle=e^{-|\alpha|^{2}}.$
This way, we can define

$|0\rangle_{\operatorname{Cat}}\coloneqq\frac{|+\rangle_{\operatorname{Cat}}+|-\rangle_{\operatorname{Cat}}}{\sqrt{2}}\text{
and
}|1\rangle_{\operatorname{Cat}}\coloneqq\frac{|+\rangle_{\operatorname{Cat}}-|-\rangle_{\operatorname{Cat}}}{\sqrt{2}}.$

For $N=100$ we are able to reconstruct the characteristic function of states
$\ket{0}_{\operatorname{Cat}}$ and $\ket{1}_{\operatorname{Cat}}$ respectively
on $[-2,2]$, for $\alpha=(1,1)$, which is illustrated in Figure VIII.

The probability distribution is then

$\displaystyle\big{|}\langle
x|0\rangle_{\operatorname{Cat}}\big{|}^{2}=\frac{|\langle
x|+\rangle|^{2}+|\langle x|-\rangle|^{2}+2\Re(\overline{\langle
x|-\rangle}\langle x|+\rangle)}{2}$ (118)

and

$\displaystyle\big{|}\langle
x|1\rangle_{\operatorname{Cat}}\big{|}^{2}=\frac{|\langle
x|+\rangle|^{2}+|\langle x|-\rangle|^{2}-2\Re(\overline{\langle
x|-\rangle}\langle x|+\rangle)}{2}$ (119)

This implies that

$\displaystyle\chi_{\ket{0}_{\operatorname{Cat}}\bra{0}_{\operatorname{Cat}}}$
$\displaystyle=\frac{1}{2}\left(\chi_{\ket{+}\bra{+}}+\chi_{\ket{+}\bra{-}}+\chi_{\ket{-}\bra{+}}+\chi_{\ket{-}\bra{-}}\right)$
(120)
$\displaystyle=\frac{1}{2}(\frac{1}{\mathcal{N}_{+}^{2}}+\frac{1}{\mathcal{N}_{-}^{2}}\Bigg{)}(\chi_{\ket{\alpha}\bra{\alpha}}+\chi_{\ket{-\alpha}\bra{-\alpha}})$
$\displaystyle\quad+\frac{1}{2}\Bigg{(}\frac{1}{\mathcal{N}_{+}^{2}}-\frac{1}{\mathcal{N}_{-}^{2}}\Bigg{)}(\chi_{\ket{\alpha}\bra{-\alpha}}+\chi_{\ket{-\alpha}\bra{\alpha}})$
$\displaystyle\quad+\frac{\chi_{\ket{\alpha}\bra{\alpha}}-\chi_{\ket{-\alpha}\bra{-\alpha}}}{\mathcal{N}_{+}\mathcal{N}_{-}}$
(121)

and

$\displaystyle\chi_{\ket{1}_{\operatorname{Cat}}\bra{1}_{\operatorname{Cat}}}$
$\displaystyle=\frac{1}{2}\left(\chi_{\ket{+}\bra{+}}-\chi_{\ket{+}\bra{-}}-\chi_{\ket{-}\bra{+}}+\chi_{\ket{-}\bra{-}}\right)$
(122)
$\displaystyle=\frac{1}{2}\Bigg{(}\frac{1}{\mathcal{N}_{+}^{2}}+\frac{1}{\mathcal{N}_{-}^{2}}\Bigg{)}(\chi_{\ket{\alpha}\bra{\alpha}}+\chi_{\ket{-\alpha}\bra{-\alpha}})$
$\displaystyle\quad+\frac{1}{2}\Bigg{(}\frac{1}{\mathcal{N}_{+}^{2}}-\frac{1}{\mathcal{N}_{-}^{2}}\Bigg{)}(\chi_{\ket{\alpha}\bra{-\alpha}}+\chi_{\ket{-\alpha}\bra{\alpha}})$
$\displaystyle-\frac{\chi_{\ket{\alpha}\bra{\alpha}}-\chi_{\ket{-\alpha}\bra{-\alpha}}}{\mathcal{N}_{+}\mathcal{N}_{-}}.$
(123)

The characteristic functions on the right hand side are explicitly given by
(31).

###  VIII-C Chain of quadratic harmonic oscillators

We consider a hermitian Hamiltonian matrix $h=\begin{pmatrix}h_{XX}&h_{XP}\\\
h_{PX}&h_{PP}\end{pmatrix}$ to which we associate a quadratic chain with
Hamiltonian

$H=\langle R,hR\rangle\text{ with }R=(X,P),$

where $h_{XP}=h_{PX}=0$, $h_{PP}=\frac{\operatorname{id}}{2}$ and

$h_{XX}=\begin{pmatrix}1/2&-\frac{\kappa}{4}&\cdots&0&-\frac{\kappa}{4}\\\
-\frac{\kappa}{4}&1/2&-\frac{\kappa}{4}&&0\\\
\vdots&-\frac{\kappa}{4}&1/2&\ddots&\vdots\\\
0&&\ddots&\ddots&-\frac{\kappa}{4}\\\
-\frac{\kappa}{4}&0&\cdots&-\frac{\kappa}{4}&1/2\\\ \end{pmatrix}$

for some $\kappa\in[-1,1].$ Here, $\kappa$ refers to the interaction strength
between neighboring oscillators: small $\kappa$ means weak correlations
whereas large $\kappa$ means strong correlations. Let

$X=h_{XX}^{-1/2}\sqrt{\sqrt{h_{XX}}h_{PP}\sqrt{h_{XX}}}h_{XX}^{-1/2}.$

The ground state $|\psi\rangle$ of $H$ is a Gaussian state with covariance
matrix given by $\gamma=\operatorname{diag}(X,X^{-1})$ [80]. The decay of
spatial correlations is illustrated in Figure VIII. We will again use our
heterodyne tomography protocol to reconstruct that decay. The distribution of
the random variable arising from the heterodyne detection takes the following
expression:

$|\langle x|\psi\rangle|^{2}\propto
e^{-x^{T}\Omega^{T}\Big{(}\frac{I+\gamma}{2}\Big{)}^{-1}\Omega x}\,.$ (124)

We recall that for any two modes $i,j$, the reduced shadow characteristic
function after sampling $(x_{k})_{k=1}^{N}$ from the above distribution is
defined as

$\displaystyle\hat{\chi}(u_{i},u_{j})\coloneqq\frac{1}{N}\sum_{k=1}^{N}e^{\frac{1}{4}(\|u_{i}\|^{2}+\|u_{j}\|^{2})-i(u_{i}^{\intercal}\Omega(x_{k})_{i}+u_{j}^{\intercal}\Omega(x_{k})_{j})}\,,$

where $(x_{k})_{i},(x_{k})_{j}\in\mathbb{R}^{2}$ are the components of vector
$x_{k}\in\mathbb{R}^{2m}$ corresponding to modes $i$ and $j$. The reduced
characteristic functions of the marginal of $|\psi\rangle$ over modes $i$ and
$j$ takes the form

$\displaystyle\chi(u_{i},u_{j})=e^{-\frac{1}{4}(u_{i},u_{j})^{\intercal}\gamma(u_{i},u_{j})}$

From which the correlations coefficients corresponding to these two modes can
be read-off.

##  IX Acknowledgements

LL acknowledges financial support from the Alexander von Humboldt Foundation.
The authors are also grateful to Pembroke College for funding a workshop in
Cambridge in October 2022 where we had some useful discussions.

###  IX-A Application of shadow tomography to photonic crystals and 2D
materials

Photonic crystals are dielectric solids with a periodically modulated
refraction index. One of their many exciting features is that the photonic
bandstructure in photonic crystals, which is accessible to optical
measurements, provides insights into the electronic bandstructure in periodic
solids. Just like electrons in solids, photons are prohibited from propagating
at band-gap frequencies inside the medium. While the band structure of
photonic crystals is obtained from the classical Maxwell equations, the
filling of bands by photons relies on a quantum mechanical description of the
electromagnetic field. Such a description is also necessary to understand
effects such as spontaneous emission of photons in photonic crystals [81, 82].
We shall demonstrate now how shadow tomography can be applied to exhibit the
dispersion surface of the Bloch-Floquet bands. Our description here is
semiclassical, as we treat the photons quantum-mechanically and the crystallic
band structure classically.

In our discussion of applications of shadow tomography, we will then focus on
optical analogues of graphene called _photonic graphene_

We start by using classical electrodynamics to describe the propagation of
electromagnetic waves in photonic crystals. We will then quantise the
electromagnetic fields to convert the classical picture to a quantum picture
involving photons and illustrate how the methods developed in this article can
be used in photonic crystals. This semiclassical approach is used to simplify
the analysis and get explicit formulas for photonic states and electromagnetic
fields.

####  IX-A1 Derivation of the Helmholtz equation