Precision Bounds on Continuous-Variable State Tomography using Classical Shadows

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


I. INTRODUCTION
The ability to estimate states accurately with as few measurements as possible -the primary goal of quantum state tomography -yields an indispensable tool in quantum information processing [1].State characterization is necessary for realizing quantum technologies in finite-dimensional tensor-product quantum systems governed by discrete variables (DV), as well as in quantum systems governed by continuous variables (CV) [2], such as electromagnetic or mechanical modes.
The recent development of classical-shadow tomography yields a succinct way to learn information about a DV quantum state through randomly chosen measurements, in such a way that the learned information can later be used to predict properties of the state [3][4][5].
A key benefit of shadow tomography is that it comes with rigorously proven guarantees on the minimum number of samples required to achieve high accuracy with high probability.However, the best guarantees require structures from finite-dimensional DV spaces, such as state and unitary designs [4] or SIC-POVMs [6], thereby obscuring any practical extension to the intrinsically infinite-dimensional CV systems.
In this paper, we apply the shadow tomography framework to a large family of well-known and well-utilized CV tomographic protocols.We distill the mathematical tools behind shadow guarantees in such a way that eliminates dependence on strictly DV ingredients and allows us to reformulate established CV protocols in the shadow framework.This reformulation yields accuracy guarantees for expectation values of local observables whose required number of samples scales polynomially with both the number of participating CV modes and the maximum occupation (a.k.a.photon) number of each mode.
CV tomography is a long-standing and well-developed field [7].Focusing on established protocols allows us to boost their credibility, as opposed to developing new protocols that may be equally efficient theoretically but whose practical utility is left as an open question (e.g.[8, Sec.VI.A]).Specifically, our work focuses on tomographic methods that can be easily implemented experimentally with existing quantum optical technology, which underpins fiber-based and free-space quantum communication and key distribution [2].To illustrate the connection to optics and our commitment to the mantra of classical shadows "measure first, ask questions later" [5], we apply our framework retroactively to experimental quantum-optical data published in 2010 [9].Naturally, our theoretical guarantees apply equally well if the corresponding protocols are performed in other CV platforms such as microwave cavities coupled to superconducting qubits [10], motional degrees of freedom of trapped ions [11], as well as optomechanical [12] and nano-acoustic [13] resonators.
The first method that we are able to recast in the shadow-tomography framework is homodyne detection [7,[14][15][16][17][18][19].We show that the number of samples needs to scale at most as the fifth power of the maxi- Design based CV tomography [8] 3-design  4  -

Multimode shadows
Local shadows  4 −  (6+1⇑3) -TABLE I.A summary of the different protocols that exist, to our knowledge, studying the sample complexities of CV state tomography.For each protocol, we indicate the measurement basis and a bound on the scaling of sample complexity (with system size  and number of modes  for multimode shadows).The sample complexity of design based shadows was estimated using known variance bounds [4,8].Additionally, for single-mode measurement methods considered in this paper, we provide the scaling observed in numerical simulations.
mum occupation number (up to a logarithmic correction) to yield reliable homodyne shadow estimates of a singlemode state.Since we require a finite occupation-number cutoff, this bound only holds for portions of states supported on finite-dimensional subspaces of the infinitedimensional Fock space.Our theoretical guarantees use several important technical bounds proved in an earlier work [20].While there have been studies of the statistical efficiency of homodyne tomography (e.g., [20][21][22][23][24]), to our knowledge, our bounds are an improvement over previous results as they allow us to analyze the sample complexity for the convergence of the operator norm of the estimated state to its ideal value.These improvements are made possible because of our use of matrix concentration inequalities reviewed in Sec.II (see Ref. [25,26]) that have not been employed in past work on CV state tomography.A recent paper has also considered the statistical efficiency of homodyne using these improved matrix concentration inequalities [27].Here, we present more explicit formulas for the bounds in the case of states with a hard photon number cutoff.Adapting our approach to the formulation in [27] also leads to explicit bounds in their case that appear to give agreement between the two papers in areas where there is overlap.On a related note, a recent paper Ref. [19] studies the sample complexity of estimating unknown Gaussian or Generalized (i.e.linear combinations of) Gaussian CV states, where a risk function/error is minimized.The distance between the outcomes of the original and estimated state, averaged over channels and measurements, was the metric of error/risk function.Our work, on the other hand, considers an absolute error via the infinity norm of the difference of the original and estimated state.Finally, we remark that our results provide upper bounds on the scaling of the number of samples required, i.e. the sample complexity, for estimating the state from random measurements.We leave open the important problem of obtaining lower bounds on the optimal sample complexity [4], which have been examined, e.g., in the case of homodyne tomography [28].
Building off the formalism in [29], we are also able to recast a large class of protocols utilizing displaced Fock states and make contact with heterodyne, photonnumber resolved (PNR), and photon-parity tomography [7].For the latter two methods, we show that the number of samples needs to scale at most as the fourth power of the maximum occupation (up to a logarithmic correction) to yield reliable PNR shadow estimates.For heterodyne measurements, our upper bound is effectively exponential in the maximum occupation.
Following our theoretical analysis of single-mode systems, we provide a generalization to multimode states.We prove that local CV shadow data can be used to reconstruct -local reduced density matrices on any subset of the modes with a sample complexity, polynomial in the total number of modes and exponential in .
We support our theoretical analysis with multiple numerical simulations and the analysis of existing experimental homodyne data.We begin by verifying that both the shadow methods indeed construct reliable estimates of a target state.We then determine how the number of samples required for both the methods scales with the maximum occupation number  .Through simulations involving the reconstruction of the vacuum state, coherent state superposition (CSS) states (i.e., cat states) and random pure states, we observe that the scaling for PNR and photon-parity measurements is similar to the calculated upper bound, whereas homodyne measurements exceed theoretical expectations and surpass the observed PNR scaling in every case.We also provide a basic demonstration of homodyne tomography on multimode separable and entangled states, finding that the latter require more samples to reach the same precision.Our results are summarized in Table I.
We begin with a general result for estimating density matrix elements of a single-mode CV state using a generic shadow-based protocol in Sec.II.Its application to homodyne measurements and photon-number resolving measurements is discussed in Sec.III.In Sec.IV, we extend our theory to multimode CV states, and support our methods by numerical results and analysis of experimental data in Sec.V. We conclude in Sec.VI.

II. SINGLE-MODE SHADOWS:
GENERAL CASE Here, we distill the necessary mathematical tools from DV shadow-based accuracy guarantees in a general lemma that we apply to specific CV protocols in later sections.We adapt the slight reformulation from [30, Appx.A] of the original shadow work, in which guarantees are shown exclusively in terms of the system state.This allows us to derive sampling guarantees that are suboptimal relative to the qubit case, but are sufficiently general to be applied to numerous CV tomographic protocols.We assume loss-free propagation and detection for our theoretical derivations.
A single-mode CV density matrix can be represented in various ways, e.g., using Wigner, ,  [31,32] or more general representations [29], moment expansions [33], homodyne measurements [7,16,17], or outer products of occupation-number (a.k.a.Fock) or other simple basis states.All of these ways effectively express the density matrix as a linear combination of some basis set of operators {()} with corresponding coefficients (), where  parameterizes some generally continuous and non-compact space  ,  = ∫   () () . ( Assuming a non-negative and normalized weight function () and sufficiently well-behaved space  allows one to re-express the above formula in a statistical fashion, i.e., as an expectation value over samples , sampled from  according to the probability distribution governed by (), We remark that classical shadow tomography in our formulation refers to tomographic protocols with the restriction that () is a probability distribution over random measurement settings and their outcomes.A general quantum state tomography protocol based solely on Eq. (1) may not have a formulation as a type of classical shadow tomography.Crucially, though, we will see below that there is a natural analysis of homodyne tomography that fits into the classical shadow framework.
The above formulation allows for a simple application of matrix concentration inequalities (loosely speaking, the matrix version of the law of large numbers; see Refs.[25,26]).If we are given a set of  independent samples  = { 1 , . . .,   }, we can construct an estimator for the state by averaging the basis operators corresponding to the outcomes in the set, We define the above estimator as a shadow of size T.Such shadows recover the state by taking their expectation value over all sample sets  , This convergence property is a necessary starting point for shadow estimation, but it says nothing about how fast the convergence happens as the number of samples increases [34].
In order to bound the rate of convergence of the above estimate, we restrict  to finite-dimensional truncated versions  ( ) , with the dimension cutoff  assumed to be sufficiently large so that the truncated operators capture the essence of the original ones.The resulting matrices have non-zero entries only for row and column indices 0 through  − 1.A "regularized" shadow is constructed with snapshots projected into the ( − 1)-photon subspace, using   = ∑ Our analysis can be readily extended to more general soft-cutoffs by replacing the projectors onto fixed photon number spaces by operators with smooth cuttoffs [35].
Note also that we do not re-normalize the above regularized density matrices or shadows.Doing so introduces additional factors that would complicate our calculations.Hence, for the sake of ease and uniformity, we chose to work with the regularized matrices as they are.
We now bound the difference between the sample average  ( )  and the actual average  ( ) =     , backing out the minimum shadow size  required for an accurate estimate.We employ the infinity norm of the difference  ( )   −  ( ) to quantify the error.
Lemma II.1.(a) Fix ,  ∈ (0, 1), and let  ( )  be an  -dimensional classical shadow of a single-mode continuous-variable quantum state .If the size of the classical shadow is at least  , such that where then the probability Proof.Applying the matrix Bernstein inequality [26][30, Appx.A] to  ( )  −  ( ) and using Eq. ( 4) yields Using the equivalence relation .
Equating the error term on the right-hand side of the inequality to  and rearranging yield the result (6).
We substituted the original variance term ∏︁ ∞ in Eq. 7a , which measures the fluctuations about the mean  ( ) , with the term ∏︁E  ( ( ) ) 2 ∏︁ ∞ for ease in calculations.We noticed that this leads to the derived bounds being higher than they would otherwise in some cases.
We see that the scaling of the required sample size in Eq. ( 6) depends logarithmically on the inverse error probability  and polynomially on the accuracy , recovering the corresponding portion of the qubit shadow result [4].To determine dependence on occupation number cutoff  , we need to know the corresponding scaling of  and  (7), which are sometimes referred to as shadow norms.This scaling depends on the detail of the protocol, and we proceed to apply this lemma to specific cases.

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

A. Homodyne Shadows
Homodyne tomography is a popular technique used in quantum optics to make quadrature measurements of CV quantum states [7].This technique allows one to measure a quadrature operator such as oscillator position, momentum, or any of their superpositions.
Let ,  † be the canonical mode lowering and raising operators satisfying (︀,  † ⌋︀ = 1.A general quadrature operator, is parameterized by  ∈ (︀0, ), where  = 0 ( = ⇑2) corresponds to the pure position (momentum) quadrature.We consider all possible values of  to be admissible and do not consider the problem of choosing an appropriate discretization.In practice, we use pseudorandom number generators to generate the values of   .Each operator admits its own set of orthogonal but nonnormalizable "eigenstates" ⋃︀  ̃︀ with   ∈ R.Even though such "eigenstates" are not quantum states but are instead distributions, they do resolve the identity for each  and thus yield valid positive operator-valued measures (POVMs) [36].For a given , a homodyne protocol yields a measurement in the corresponding set of eigenstates.
To perform a homodyne measurement, the unknown system state interferes with an ancillary mode in a coherent state ⋃︀ LO ̃︀ (coherent local oscillator or LO) via a beam splitter, as shown in Fig. 1(a).The phase of the coherent state LO is precisely the phase  which defines the quadrature that is measured, and is picked uniformly from (︀0, ).Both output amplitudes of the beam splitter,  1 and  2 respectively, are measured via homodyne detectors.Their difference  − is proportional to the quadrature amplitude of the signal: To cast homodyne tomography into shadow terminology, we first specialize the expansion in Eq. (1) to this case.Substituting parameters  → (,   ), we get the expansion (see Eq. ( 23) in Ref. [7]) where 1⇑ corresponds to picking  uniformly from (︀0, ), and ∐︀  ⋃︀⋃︀  ̃︀ to the conditional probability of the homodyne setup yielding   as the outcome.The matrix elements of the basis operators { (  , )} have standard expansions in quantum optics literature.They are defined in terms of what are known as pattern functions   () [7,37].
These pattern functions are symmetric in , , and are defined using Fock state wavefunctions   () ( ℎ energy eigen state of a harmonic oscillator) and   () ( ℎ non-normalizable solution of the Schrödinger equation of a harmonic oscillator) as Because the quadrature probability distributions are normalized for each , the above equation can automatically be interpreted in the statistical fashion of Eq. ( 2).In this interpretation, basis operators  are sampled according to the probability distribution defined by  and ∐︀  ⋃︀⋃︀  ̃︀.A regularized homodyne CV shadow of size  , constructed using the samples {  ,   }  =1 , is therefore We remark that the randomized choice of   independent of the state is what allows us to reinterpret this reconstruction method as an instance of classical shadows.We also remark that a homodyne shadow is distinct from the Wigner function projection of the state.A shadow estimates the density operator in the Fock basis from a collection of homodyne samples with randomly chosen .A Wigner function projection ∫    (  ,   ) = ∐︀  ⋃︀⋃︀  ̃︀, on the other hand, is the probability of obtaining specific measurement outcomes at a fixed values of .
We now adapt Lemma II.1 to homodyne shadows, relying on the bounds on pattern function matrix from Ref. [20].The following theorem gives an upper bound on the sample complexity of constructing shadows through homodyne samples, for an error ∏︁ ( )   −  ( ) ∏︁ ∞ ≤  with high probability (greater than 1 − ).
Theorem III.1.Using homodyne tomography, for ,  ∈ (0, 1), some positive constant  1 [20], and an N-dimensional homodyne shadow  ( )   (13) of a single-mode CV state , Pr(∏︁ ( )   −  ( ) ∏︁ ∞ ≤ ) ≥ 1 −  when the size of the shadow is at least Proof.From Theorem 3.2 of Ref. [20], we use ∏︁ 2 2 ≤  2  3 , which is valid for any value of homodyne angle  and gives us and Therefore, using Lemma II.1, the size of the shadow needs to be at least We see that the upper bound on the minimum number of samples scales as  13⇑3 log  for a photon-cutoff  .
Regarding the utility of homodyne shadows in calculating expectation values of observables, the patternfunction operators  are unequivocally more complex than the original qubit shadows [4].Nevertheless, efficient methods are known to compute them using recursive relations [38].

B. Photon Number Resolving (PNR) Shadows
Our second application is based on the  -operator formalism developed in Ref. [29].The relevant basis expansion of a single-mode CV state  can be written in terms of what are called T-operators, Here, Tr(︀ (,  * )⌋︀ plays the role of a quasiprobability distribution over the complex plane C, and the dual operators T (,  * ) are used to reconstruct the state given this distribution.
The  -operators are defined in terms of displaced Fock states, which are the result of applying a displacement operator () = exp( † −  * ) on a Fock state: ⋃︀, ̃︀ ≡ ()⋃︀̃︀.There exists a family of  (,  * ), and their corresponding duals T (,  * ), parameterized by −1 <  < 1, such that the dual operator is given by (,  * ) belonging to this family can be diagonalized as where the eigenvalues Utilizing the above expansions, Eq. ( 17) becomes which we interpret as an expectation value of samples  ()  T (,  * ) distributed according to the probabilities 1  ∐︀, ⋃︀⋃︀, ̃︀.However, in a physical interferometric measurement procedure corresponding to this expansion [see Fig. 1(b)], the displacement parameter  comes from a non-compact set.This means that we need to first truncate the sample region in order to make the above a valid (i.e., normalizable) probability distribution.
We restrict the integration of Eq. ( 21) to be over a phase-space region  of finite area .We can set this region to be, e.g., a disk of radius  max .We will then be able to account for the truncated Fock space by setting  2 max =  since that is the average occupation number of a coherent state with ⋃︀⋃︀ =  max .In that case,  = 4 scales linearly with maximum photon number cutoff.
Projecting into an  -dimensional subspace implies that we consider only states with some max photon number  − 1 and regularized  -operators T ( )  (,  * ).Combining this with the truncation of the phase-space integration yields the following approximate expression for the truncated density matrix, in which we interpret   ()  T ( )  (,  * ) as a PNR shadow sampled from the probability distribution with parameters 0 ≤  ≤  − 1 and  ∈ .

The eigenvalues 𝜆 (𝑛) 𝑟
exhibit different behaviours for different values of , and in Appx.A we analyze how many samples are required for density matrix estimation in each case.We obtain quartic scaling with  at  = 0, up to logarithmic corrections, and unfavorable exponential scaling for all other 0 < ⋃︀⋃︀ < 1.In the favorable case, the upper bound on the minimum number of samples is beating our theoretical estimates for homodyne shadows.Therefore, we see that  -operators at  = 0 have a favorable sample complexity scaling and can be used as a tomographic method to estimate an unknown state.The matrix elements of the  -operator in the Fock basis are equal to the matrix elements of the displacement operator up to a phase which can be used to efficiently compute the classical shadows.Moreover, the case of  = 0 can also be realized via photon-parity measurement tomography [39][40][41][42], providing another experimental technique to realize these shadows.
It is worth noting that in the limit  = −1, T (,  * ) =  − (,  * ) becomes a projector onto the coherent state ⋃︀̃︀, making the measurement protocol equivalent to what is known as heterodyne measurement, achieved using a similar setup as homodyne measurement [43].The basis expansion is given in terms of the Q-function () = 1  ∐︀⋃︀⋃︀̃︀.However, the sample complexity upper bound diverges in this case.It is possible to avoid this divergence by taking a limit  → −1 as the number of samples increases.Unfortunately, this method only achieves an upper bound that is exponential in  because, as noted above, for any 0 < ⋃︀⋃︀ < 1, our sample complexity upper bound diverges exponentially in  (see Appx.A).

IV. MULTIMODE SHADOWS
In this section, we show that it is possible to construct efficient representations of reduced density matrices on constant numbers of modes from local CV shadow data.
Multimode CV states are states that are made up of multiple modes of the underlying system, such as different frequencies of light in an optical system.This is analogous to having a multi-qubit state in DV.Each mode is equipped with a generally continuous and noncompact space   and M-mode, possibly entangled, quantum states reside in  = ⊗  =1   .Let the space  be parameterized by a variable  = ( 1 , . . .,   ), and spanned by the set of basis operators {()}.Extending Eq. ( 1), the density matrix of a general multimode state belonging to this space can then be written as where () are the corresponding coefficients.Boldface symbols will be used to represent multimode variables and states henceforth.All the density operators and shadows considered in this section are the regularized versions, projected into the  -photon subspace of each mode.The superscripts indicating this are omitted for the sake of brevity.
Analogous to the single-mode case in Sec.II, we require a basis expansion of  such that () ≥ 0 ∀; the nonnegativity of () allowing it be interpreted as a probability distribution over .We proceed by using an expansion of the form () = ∐︀⋃︀⋃︀̃︀, where the quantum states or distributions ⋃︀̃︀ parameterize the underlying space.For example, we can take ⋃︀̃︀ to be local position state eigenstates in the case of homodyne tomography.
We construct shadows of multimode states using local, joint single-mode measurement.For an estimator of a state  to qualify as a shadow, as introduced in Section Sec.II, it must be equal to the state  in expectation.Below we show that it can indeed be satisfied for shadows constructed through local measurements.For this calculation, we utilize another expansion of a general quantum state  in terms of some set of separable operators {  } and corresponding complex coefficients {  }, Such an expansion is possible irrespective of whether the state is separable or entangled.The estimator becomes In the following theorem, we calculate the variance and thereby the sample complexity of estimating a multimode state through local CV shadows.This result reduces to the single-mode result from II.1 when  = 1.
Theorem IV.1.Fix ,  ∈ (0, 1), and let   be an  dimensional classical shadow obtained through local measurements () = ⊗  =1 (  ) of an  -mode CV state  (27).If the size of the classical shadow is at least written in terms of the single-mode variance  2 1 (7a) and Proof.The measurement variance for an  -mode state proceeds similar to the single-mode case with the help of Eq. ( 27), Using matrix Bernstein inequality and following a calculation similar to Lemma.II.1, replacing  with   , we get the desired result in Eq. ( 29).
As a consequence, we obtain a sample complexity bound on estimating -local reduced density matrices, which are the outcome of tracing over all but k modes of a multimode state.Since we are tracing out some of the modes, we are left with a 'reduced density matrix'.
Corollary IV.1.1.If   denotes the reduced density matrix of  in  ≤  modes, and its estimator constructed through local measurements is  , , such that the size of the shadow is at least Proof.Without loss of generality, let's assume   to be the reduced density matrix over the modes 1, . . ., , with
The classical shadow of the reduced matrix  , is then the average of  such operators.The variance of estimating this reduced matrix, following earlier calculations, is One can further generalize this result to derive sample complexity bounds for determining the reduced density matrix on any choice of -local modes.The proof follows from an application of the union bound and increases the sample complexity by a factor of  [4,30].
The sum of the absolute values of the coefficients in Eq. 27 is a factor that affects the sampling complexity of multimode states.These coefficients can in general be complex.When the state is separable (i.e. can be written as a product state) this factor becomes 1, thereby resulting in a complexity that is equivalent to sampling all the modes individually, as expected of a product state.When the coefficients are all positive, they again sum up to 1 (since the trace of the state must equal 1), resulting in no increase in complexity as compared to a separable state.In the general case, we can still upper bound the sampling complexity, but the overall scaling with the number of modes has a larger prefactor in the exponent.These results are summarized in Table.I.
We would like to add that the proposed local-joint measurements are not the most general measurements for multimode CV states.Upon trying many different bases of entangled/global states, we have not found a suitable one that fulfills our criteria of a positive weight function for all states (() in Eq. 26), but have not ruled out the possibility of their existence.Studies of whether global measurements can facilitate such a statistical interpretation, and of the sample complexity when using such bases are interesting avenues for future work.
Proof.Since the 1-norm of matrix is the sum of its singular values, and using the Holder's inequality, Therefore, Using matrix Bernstein inequality, this implies  (⋃︀ (  ( −   ))⋃︀ ≥ ) , where is  and  are as defined in Theorem.IV.1.Equating the upper bound on probability to ⇑ , and using the union bound [4], a shadow size upper bounded by suffices to estimate the list of multimode observables {  }  =1 with the desired accuracy.

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

A. Single mode simulations
To analyze the applicability of our theoretical guarantees, we construct estimators of density matrices of certain target states using homodyne and PNR shadows.The measurement bases in all cases were generated using uniform pseudorandom number generators.We then analyze how close they are to the target states.Figure 2(a) shows how a reconstructed CSS state with mean-photon number 10 compares to the original density matrix elements for both homodyne and PNR tomography.A comparison in terms of the Wigner function is shown in Fig. 2(b).To achieve equally accurate shadows such that ∏︁ ( )   − ( ) ∏︁ ∞ = 0.1, we had to use 5×10 4 and 10 6 samples for homodyne and PNR, respectively.Of primary interest in the classical shadow formalism is the number of samples required to achieve a given target precision with high-probability.This minimum shadow size for a given state  was numerically obtained by repeating the estimation process for several shadow sizes and identifying the one which results in a precision , with a probability of at least 1 − .We observed that if we are looking for a fixed error in norm (i.e.fixing ), increasing the shadow size decreases the probability with which this error occurs (i.e.) as the estimation becomes more accurate.In fact, we observed that  for a fixed  is proportional to log 1⇑ and this helped us narrow down our search for minimum shadow sizes corresponding to our chosen pair (, )=(0.1,0.1).The scaling with the occupational number cutoff is then obtained by finding the minimum shadow size for different  , and the error bars indicate the statistical errors in estimating  via this process.In Fig. 2 (c)-(d), we plotted this observed minimum size for different states, for homodyne and PNR shadows respectively.
We ran our simulations for a vacuum state ⋃︀̃︀ = ⋃︀0̃︀, a CSS state ⋃︀̃︀ ∝ ⋃︀̃︀ + ⋃︀ − ̃︀ with mean photon number ⋃︀⋃︀ 2 = 10, and random pure states.The random pure states of size  had non-zero support only up to ( − 1) th Fock state.The results show that homodyne has a scaling better than  2 for the cases we tried, consistently outperforming the theoretical upper bound.For vacuum states, PNR on the other hand has a scaling much closer to the predicted upper bound of  4 .We also tested the scaling of PNR tomography of random pure states, but were computationally limited in the accessible sizes.Over the range we observed, the results are consistent with the scaling of the vacuum state.

B. Experimental data
Here, we analyze the sample complexity of experimental homodyne data from Ref. [9].The experiment consisted of creating Coherent State Superpositions (CSSs), where the created states were approximations to smallamplitude, even CSS states ⋃︀̃︀ ∝ ⋃︀̃︀ + ⋃︀ − ̃︀.Squeezed vacuum was generated through spontaneous parametric down-conversion in a KNbO 3 non-linear crystal: an upconverted laser pulse at 430 nm created the squeezed vacuum at 860 nm.After spectral filtering, the squeezed vacuum was incident on a weakly reflecting beam splitter with reflectivity .Photons detected in the reflected path herald a CSS emerging from the transmitted port of the beam splitter.The resulting CSS was then directed to a homodyne detection setup.The homodyne setup consisted of a 50/50 beam splitter, two high-efficiency photodiodes and a low-noise amplification circuit.The local oscillator was created with a strong laser field (≈ 10 9 photons per pulse) that allows access to the quadratures through balanced detection.The heralded photon detection was done using either two singlephoton avalanche diodes (SPADs) [44] or one transition edge sensor (TES) [45].The experiment analyzed here used both heralded-photon detection schemes.The CSSs were heralded upon the detection of two heralding photons, generating the approximation of a small-amplitude even cat state.
Using the homodyne data from this experiment, we constructed an estimate of the state, an example of which is shown in Fig. 3(a-b), and studied how close the reconstructed states are to the target states.Before analyzing the sample complexity, the experimental data was pre- processed to account for propagation losses and detection inefficiencies [9].
We then repeated the same numerics as in Sec.V A to find the minimum shadow size ( ) required to achieve a fixed (, ) for different occupational number ( ) cutoffs.Working with a given dataset means our shadow size is limited to the size of the dataset and that accuracy can't be improved arbitrarily.Fig. 3(c) shows how  changed with  , and we observe a scaling of  ∝  0.659±0.001for a cat CSS with  = 1.16 using the TES, and as  ∝  0.685±0.001for a cat CSS with  = 1.32 using the SPADs.This scaling of approximately  ∝  2 aligns with our numerical simulations of single-mode homodyne shadows, and is a further testament to the reliability of our simulations.

C. Multimode simulations
We now provide an example reconstruction of both separable and entangled two-mode states.
In Fig. 3(d) we can see original and reconstructed density matrices of a separable 2-mode state ⋃︀Ψ̃︀ = ⋃︀  ̃︀⊗⋃︀  ̃︀, where ⋃︀  ̃︀ is an even CSS state of amplitude  = ⌋︂ 1.5.In Fig. 3(e), we have the reconstruction of an entangled state ⋃︀Ψ̃︀ = ⋃︀̃︀ ⊗ ⋃︀̃︀ + ⋃︀ − ̃︀ ⊗ ⋃︀ − ̃︀, with the same number of measurements.The infinity-norm distance  between the original and reconstructed states is 0.03 and 0.05, respectively.The increased sample complexity needed to reach a given precision for the entangled state tomography is consistent with our general expectations from the theoretical analysis in Sec.IV.

VI. OUTLOOK AND CONCLUSIONS
We developed a formalism to derive worst-case sample complexity of a large class of state-reconstruction procedures for single-mode continuous-variable (CV) states (Lemma II.1, Theorem III.1).In an experimental context, our work provides a tool to compare the sample complexity of different measurement methods.For photon-number resolving measurements, we show through numerical simulations that our analytical derived upper bound is close to the observed complexity.In the case of homodyne tomography, we notice that the simulation results surpass our derived bound, rendering homodyne tomography more efficient than PNR and photon-parity tomography in practice.
We obtain the sample complexity of estimating multimode CV states by reconstructing shadows of the true states from local measurements on each mode, in the spirit of the qubit-based shadow tomography proposal [3][4][5].We also consider a CV manifestation of a global version of the original shadow protocol, where one uses global measurements over varying linear combinations of the constituent modes for state reconstruction.Of particular note is that, in contrast to the discrete-variable systems, our analysis does not rely on the construction of state designs (approximate or otherwise) such as in [8, Sec.VI.A].Instead, we focus on adapting existing techniques for CV tomography to fit within the framework of randomized measurements.
In summary, we recast existing homodyne and photonnumber-resolving protocols as shadow-tomography protocols and show that such protocols yield good local estimates of a multimode state using a number of samples that is polynomial in the number of modes IV.1.1.Our CV shadow framework can be extended to analyze other metrics like the total variation distance between two outcome probability distributions, which is valuable for verifying sampling experiments like Boson sampling [46].Another remaining open question is to determine robustness of CV protocols to noise, e.g., within the framework of robust shadow estimation [47].Deriving rigorous lower bounds that can match the numerically and experimentally observed scaling for homodyne tomography is also an important outstanding challenge.

FIG. 1 .
FIG. 1.(a) Balanced homodyne detection: the unknown state signal interferes with a local oscillator (a coherent state ⋃︀̃︀) and the intensities 1 and 2 at the output arms are measured.The difference in the intensities is proportional to the quadrature amplitude of the unknown state in the mode defined by the oscillator.(b)Photon-number resolving (PNR) measurement: the local oscillator is in a coherent state ⋃︀̃︀, and detecting  photon at the detector occurs with a probability ( ⋃︀) = ∐︀, ⋃︀⋃︀, ̃︀. indicates the beam splitter in both cases.