Efficient verification of Boson Sampling

The demonstration of quantum speedup, also known as quantum computational supremacy, that is the ability of quantum computers to outperform dramatically their classical counterparts, is an important milestone in the field of quantum computing. While quantum speedup experiments are gradually escaping the regime of classical simulation, they still lack efficient verification protocols and rely on partial validation. Here we derive an efficient protocol for verifying with single-mode Gaussian measurements the output states of a large class of continuous-variable quantum circuits demonstrating quantum speedup, including Boson Sampling experiments, thus enabling a convincing demonstration of quantum speedup with photonic computing. Beyond the quantum speedup milestone, our results also enable the efficient and reliable certification of a large class of intractable continuous-variable multimode quantum states.


Introduction
Quantum information promises many technological applications beyond classical information [1,2,3,4].Shor's famous factoring algorithm [5] has highlighted the possibility of spectacular quantum advantages over classical computing.The ability of quantum computers to greatly outperform their classical counterparts is referred to as "quantum computational supremacy" [6,7], or quantum speedup [8].We use the latter terminology in this article.
While a universal quantum computer factoring large products of primes would provide a Ulysse Chabaud: ulysse.chabaud@gmail.comcompelling evidence of quantum speedup, the experimental requirements associated with such a demonstration preclude its realisation in the near future [9].For that reason, various sub-universal models of quantum computing, such as Boson Sampling [10], have been introduced [11,12,13,14,15].Beyond their conceptual relevance, these models, although not believed to possess the full computational power of a universal quantum computer, enable the possibility of a demonstration of quantum speedup in the near term, under certain complexity-theoretic assumptions.Each of these sub-universal models is associated with a computational problem that the model solves efficiently, but which is hard to solve for classical computers.These computational problems are sampling problems, for which the task is to draw random samples from a target probability distribution fixed by the model.In practice, any experimental implementation would suffer from noise and imperfections and would sample from an imperfect probability distribution which approximates the target one.Hence, the fact that even an approximate version of the sampling problem associated with a sub-universal model is hard for classical computers is crucial for an experimental demonstration of quantum speedup.For all models, such an approximate sampling hardness, which corresponds to the fact that a probability distribution which has a small constant total variation distance with the target probability distribution is still hard to sample classically, comes at the cost of assuming some additional reasonable-though unproven-conjectures.
The experimental demonstration of quantum speedup thus involves: (i) a quantum device solving efficiently a sampling task which is provably hard to solve for classical computers under reasonable theoretical assumptions, together with (ii) a verification that the quantum device indeed solved the hard task [7].While the former has been recently achieved with random superconducting circuits [16] and Gaussian Boson Sampling [17], the latter is still incomplete and relying on various questionable assumptions [18,19,20], and verification is a crucial missing point for a convincing demonstration of quantum speedup.
The setting of verification involves two parties, which in the language of interactive proof systems are referred to as the prover and the verifier.The prover is conceived as a powerful, untrusted party, while the verifier is a trusted party with restricted computational power.The verifier asks the prover to perform a computational task and verifies its correct behaviour, with possibly multiple rounds of interaction.While this setting corresponds to a cryptographic scenario between two parties, e.g., in the case of delegated computing, it also embeds the notion of certification of a physical experiment, where the experiment behaves as the prover and the experimenter as the verifier.Assumptions can be made on the prover, such as restricting its computational capabilities or assuming that it follows a specific behaviour.In the physical picture, this corresponds for example to choosing a noise model for an experiment.A common assumption is the so-called independently and identically distributed (i.i.d.) assumption, i.e., assuming that the output of the quantum experiment is the same at each run.Assuming i.i.d.behaviour for the prover leads to more efficient verification protocols, and this assumption can usually be removed using cryptographic techniques at the cost of an increased number of runs in the protocol [21,22].We use the language of interactive proof systems in what follows, but we emphasise that our results are not restricted to that particular context and are relevant especially in experimental scenarios.There has been much work done on verification of universal quantum computation [23,22,24], much less when the verification is limited to sub-universal models.We now briefly review the state of affairs in this direction.
A verification protocol for a quantum speedup experiment should guarantee the closeness in total variation distance between the experimental probability distribution and the target probability distribution [7].Such a verification is especially difficult because of the nature of the computational task at hand, i.e., sampling from an anti-concentrating probability distribution over an exponentially large sample space.In particular, any efficient non-interactive verification of current quantum speedup experiments with a verifier restricted to classical computations requires additional cryptographic assumptions [25].Existing verification protocols with a classical verifier based on total variation distance either rely on little-studied assumptions [12,26,27], or induce an overhead for the prover [24,28,29] which prevents a near-term use for an experimental demonstration of quantum speedup.
Weaker but more resource-efficient methods of verification with a classical verifier, which are sometimes referred to as validation [30], consist in performing a partial verification where only specific properties of the experimental probability distribution are tested, rather than closeness in total variation distance to the ideal distribution.For example, the verifier may perform statistical tests using experimental samples [14,16,31], compare the experimental probability distribution with specific mock-up distributions [32], or benchmark the individual components of the experimental setup.Ultimately, validation relies on making additional assumptions about the inner functioning of the quantum device [33].
In order to avoid relying on such undesirable assumptions, while keeping verification within nearterm experimental reach, another way for performing verification is to allow the verifier to have minimal quantum capabilities.In that context, the purpose of the verifier is to obtain a reliable bound on the trace distance between the output state produced by the prover and the ideal target state, which in turn implies a bound on the distance between the tested and target probability distributions of the samples [34].This can be done by obtaining, e.g., an estimate of the fidelity with the target state, or a fidelity witness, i.e., a tight lower bound on the fidelity [35].
In the context of quantum computing in finitedimensional Hilbert spaces [34], this minimal quantum capability corresponds to being only able to prepare single-qubit or qudit states or to perform local measurements.Protocols for verification of Instantaneous Quantum Polynomial time circuits [12] with these minimal requirements have been derived under the i.i.d.assumption with single-qubit states [36] or with local measurements [37] and more recently without the i.i.d.assumption with single-qubit states [38] or with local measurements [39].
Similarly in the context of infinite-dimensional Hilbert spaces [40], which includes Boson Sampling and related sub-universal models [10,41,42,43,44,45,46], this minimal quantum capability corresponds to being only able to prepare single-mode Gaussian states, or to perform single-mode Gaussian measurements.Gaussian state preparations and measurements are at the same time well understood theoretically [47], classically simulable [48], and routinely implemented at large scales experimentally [49,50,51].However, there is no efficient verification protocol using single-mode Gaussian state preparation nor single-mode Gaussian measurements for Boson Sampling with input single photons: current methods used for the validation of Boson sampling are either not scalable or only provide partial certificates on the tested probability distribution [30,52,53,54,55,56,57].

Results
Our contribution is three-fold: we introduce three noninteractive protocols (detailed below) for the reliable certification of continuousvariable quantum states, where the verifier only needs to perform classical computations and single-mode Gaussian measurements, namely heterodyne detection [58], in order to verify an unknown state sent by an untrusted prover.The three protocols share the same structure, depicted in Fig. 1.
For each of our three protocols, we derive one version assuming i.i.d.behaviour for the prover and another version making no assumptions whatsoever on the prover.Other than that, all protocols rely solely on the verifier having access to single-mode heterodyne detection and on the validity of quantum mechanics.Moreover, all protocols are robust and efficient, as they provide analytical confidence intervals and require a polynomial number of measurements in the number of modes and in the size of the confidence interval.
Firstly, generalising several results from [59,60], we derive a protocol for reliably estimating the fidelity of an unknown continuous-variable quantum state with any single-mode continuousvariable quantum state of bounded support over the Fock basis using heterodyne detection (Pro- f F e X c + 5 q 0 F J 5 8 5 h D 9 w P n 8 A c z e P C w = = < / l a t e x i t > ⌦N < l a t e x i t s h a 1 _ b a s e 6 4 = " If the answer is positive, the verifier may then decide to use the verified state ρ M for any quantum information processing task.In an experimental scenario, the blackbox prover corresponds to the experiment while the experimenter plays the role of the verifier; in that case, the state ρ N +M is the output quantum state over N + M runs of the experiment.

tocol 1).
Secondly, using this single-mode fidelity estimation protocol as a subroutine, we obtain a protocol for reliably estimating a fidelity witness, i.e., a tight lower bound on the fidelity, for a large class of multimode continuous-variable quantum states using heterodyne detection of an unknown state (Protocol 2).For m modes, this class of certifiable states corresponds to the m-mode states of the form: where Û is an m-mode passive linear transformation (a unitary transformation of the creation and annihilation operators of the modes), where each state |C i is a single-mode pure state with bounded support over the Fock basis, and where each operation Ĝi is a single-mode Gaussian unitary, for all i ∈ {1, . . ., m}.Thirdly, using this multimode fidelity witness estimation protocol as a subroutine, we obtain a protocol for verifying the output states of Boson Sampling experiments (Protocol 3).Our protocol provides a certificate on the total variation distance between the experimental and target probability distributions for any observable, efficiently in the number of modes and input photons, thus enabling a convincing demonstration of quantum speedup with photonic quantum computing.Moreover, the verification protocol is implemented within the same experimental setup, simply by replacing the output detectors by balanced heterodyne detection.We also optimise the efficiency of the protocol to the specific setting of Boson Sampling.
Note that other fidelity witnesses for multimode photonic state preparations have been introduced in [61].However, the number of measurements needed to estimate with constant precision the output state of a Boson Sampling interferometer with n input photons over m modes with their witnesses scales as Ω(m n+4 ), under the i.i.d.assumption, while we show that our protocol provides the same precision with O(m 2 log m) measurements.Moreover, we are able to remove the i.i.d.state preparation assumption, at the cost of an increased-though still polynomialnumber of measurements needed for the same estimate precision and confidence interval.
The rest of the paper is structured as follows: we recall heterodyne detection in section 3, we present our single-mode fidelity estimation protocol in section 4, we introduce our multimode fidelity witness estimation protocol in section 5 and we discuss the verification of Boson Sampling in section 6.Finally, we conclude in section 7.

Heterodyne detection
We write N * the set of positive integers.Hereafter, m ∈ N * denotes the number of modes and {|n } n∈N is the single-mode Fock basis.Let us define the single-mode displacement operator D(β) = exp[βâ † − β * â], for all β ∈ C, and the single-mode squeezing operator Ŝ(ξ) = exp[ 1 2 (ξ * â2 − ξâ †2 )], for all ξ ∈ C, where â † and â are the creation and annihilation operators of the mode [62], respectively.We use bold math for multimode notations.For all β, ξ ∈ C m , we write D(β) = m i=1 D(β i ) and Ŝ(ξ) = m i=1 Ŝ(ξ i ).Heterodyne detection, also called double homodyne, dual-homodyne or eight-port homodyne detection [49], is a single-mode Gaussian measure-   ment, projecting onto (unnormalised) coherent states.Mathematically, measuring a state ρ with heterodyne detection amounts to sampling from its Husimi Q phase space quasiprobability function [58]: where |α = e − |α| 2 2 n≥0 α n √ n! |n is the coherent state of amplitude α ∈ C. Heterodyne detection corresponds to a simultaneous noisy measurement of conjugate quadratures and can be implemented experimentally by mixing the state to be measured with the vacuum on a balanced beam splitter and detecting both output branches with homodyne detection.
Unbalancing the beam splitter and changing the phase of the local oscillator in the homodyne detection also allows the measurement of squeezed quadratures rotated in phase space (Fig. 2).Note that the phase of the local oscillator can be modulated with a classical postprocessing step.We refer to this Gaussian measurement as unbalanced heterodyne detection with unbalancing parameter ξ ∈ C. The positive operator-valued measure (POVM) elements of a tensor product of single-mode unbalanced heterodyne detection over m modes with unbalancing parameters ξ = (ξ 1 , . . ., ξ m ) ∈ C m are given by . Writing ξ = re iθ , the unbalalancing parameter is related to the optical setup by r = log T R , where T and R are the reflectance and transmittance of the unbalanced beam splitter, respectively, with the phase of the local oscillator being − θ 2 (see Appendix A).With these properties of heterodyne detection laid out, we detail in the following section our single-mode fidelity estimation protocol based on heterodyne detection.

Single-mode fidelity estimation
Normalised single-mode pure quantum states with bounded support over the Fock basis are referred to as core states [63,64].In what follows, we introduce a protocol for estimating the fidelity of any single-mode continuous-variable (mixed) quantum state with any target core state using heterodyne detection, by generalising heterodyne fidelity estimates from [59,60].We give two versions of our protocol, with and without the i.i.d.assumption.
We denote by E α←D [f (α)] the expected value of a function f for samples drawn from a distribution D. Let us introduce for k, l ≥ 0 the polynomials ) for z ∈ C, which are, up to a normalisation, the Laguerre 2D polynomials appearing in the expression of quasiprobability distribution of Fock states [65].Following [60], for all k, l ∈ N, we introduce with these polynomials the functions for all z ∈ C and all 0 < η < 1.Now let us define, for all p ∈ N * , all z ∈ C and all 0 < η < 1, the generalised functions: The functions z → g k,l (z, η) are polynomials multiplied by a converging Gaussian function and are thus bounded over C. Using these functions, we derive the following result: . Let ρ = +∞ i,j=0 ρ ij |i j| be a single-mode density operator.Then, where the function g k,l is defined in Eq. ( 6).
We give a proof in Appendix B. In particular, since heterodyne detection corresponds to sampling from the Husimi Q function of the measured state, this result enables estimating an element k, l of the density matrix from samples of heterodyne detection by computing the mean of the function z → g k,l (z, η) over these samples.These heterodyne estimates generalise both the ones introduced in [59], which we retrieve in the limit η → 1 and for which an assumption on the size of the support of the measured state is necessary, and the ones introduced in [60], which we retrieve by setting p = 1 and which yield less efficient estimates.
Additionnally, for any core state |C = c−1 n=0 c n |n where c ∈ N * , and for all p ∈ N * , all 0 < η < 1 and all z ∈ C, we define: Computing the mean of the function z → g (p) C (z, η) over samples from heterodyne detection of a state allows us to estimate the fidelity of the measured state with the core state |C .The estimation error may be controlled by optimising over the choice of the free parameters p and η.In particular, picking a smaller value for η decreases the error in Lemma 1 and increased the range of the function z → g (p) C (z, η), thus increasing the statistical error.
Based on this result, we now present the version of our single-mode fidelity estimation protocol under i.i.d.assumption and discuss afterwards its version without i.i.d.assumption: Protocol 1 (Single-mode fidelity estimation).Let c ∈ N * and let |C = c−1 n=0 c n |n be a core state.Let also N, M ∈ N * , and let p ∈ N * and 0 < η < 1 be free parameters.Let ρ ⊗N +M be N + M copies of an unknown single-mode (mixed) quantum state ρ.

Compute the fidelity estimate
The value F C (ρ) M obtained is an estimate of the fidelity between the M remaining copies of ρ and M copies of the target core state |C .The efficiency of the protocol is summarised by the following result: where the constant prefactor depends on the choice of the free parameters p and η in Protocol 1.
We give a detailed version of the theorem together with its proof in Appendix C, which combines Lemma 1 with Hoeffding's inequality [66].
Since the parameter p ∈ N * may vary freely, the scaling can be brought arbitrarily close to O(( M ) 2 log( 1 δ )).In practice, since the constant prefactor increases with p, an optimisation yields the optimal choice for p which minimises N when M , and δ are fixed.Moreover, the choice of the other free parameter 0 < η < 1 contributes to minimising the constant prefactor and the efficiency of Protocol 1 can also be refined by taking into account the expression of the single-mode target core state in the Fock basis (see [67] for a detailed analysis).
Hence, Protocol 1 gives a reliable and efficient way of estimating the fidelity of an unknown single-mode continuous-variable quantum state with any target core state.These core states play an important role in the characterisation of single-mode non-Gaussian states [64].Additionnally, they form a dense subset (for the trace norm) of the set of normalised single-mode pure quantum states.Hence, given any target normalised single-mode pure state |ψ , we can use our fidelity estimation protocol by targeting instead a truncation of the state |ψ in the Fock basis in order to estimate the fidelity of any single-mode continuous-variable quantum state with the state |ψ using heterodyne detection.
In an experimental scenario, one may set M = 1 in Protocol 1 and deduce information about the output state of a single run of the experiment.In a cryptographic scenario, where the state ρ is provided by an untrusted prover, the verifier may want to use some copies of the state to perform fidelity estimation of the remaining M copies, before using these remaining states for subsequent quantum information processing.In that context, the i.i.d assumption may no longer be justified: the quantum state sent by the prover ρ N +M over N + M subsystems may no longer be of the form ρ ⊗N +M , as all subsystems can possibly be entangled with each other (as well as with other quantum systems on the side of the prover).Hence, we extend Protocol 1 to a version which does not make the i.i.d.assumption, using a de Finetti reduction from [68,60].The non-i.i.d.version of Protocol 1 obtained is nearly identical, up to slight differences in the classical post-processing: a small fraction of the measured subsystems have to be discarded at random and another fraction of the samples are used for an energy test.This comes at the cost of an increased number of measurements which corresponds however to a polynomial overhead for a polynomial precision and confidence.We give a detailed analysis in Appendix D.

Multimode fidelity witness estimation
In this section, we extend our single-mode fidelity estimation protocol from the previous section to the multimode case, and we show that being able to estimate single-mode fidelities with heterodyne detection is sufficient to estimate tight fidelity witnesses, i.e., tight lower bounds on the fidelity, for the large class of multimode states in Eq. (1).
Note that the transformation in Eq. ( 1) can always be expressed as Ŝ(ξ) D(β) Û , where ξ, β ∈ C m , since single-mode Gaussian unitaries can be decomposed as a displacement, a squeezing with complex parameter, and a phase-space rotation [62].Note also that an m-mode passive linear transformation Û is associated with an m×m unitary matrix U which describes its action on the creation and annihilation operators of the modes [49].
Our extension to the multimode case is based on two observations (Lemmas 2 and 3 hereafter).
Firstly, if all the single-mode subsystems ρ i of a multimode quantum state ρ are close enough to single-mode pure states, then ρ is close to the tensor product of these single-mode pure states.Formally: Lemma 2. Let ρ be a state over m subsystems.For all i ∈ {1, . . ., m}, we denote by ρ i = Tr {1,...,m}\{i} (ρ) the reduced state of ρ over the i th subsystem.Let |ψ 1 , . . ., |ψ m be pure states.For all i ∈ {1, . . ., m}, we write where F is the fidelity, and We give a proof in Appendix E. This result shows that W is a tight lower bound on the fidelity.Importantly, this fidelity witness is only a function of the single-mode fidelities.In particular, being able to estimate single-mode fidelities with singlemode pure states using balanced heterodyne detection is enough to estimate fidelity witnesses with pure product states using balanced heterodyne detection.
At this point, we can estimate tight fidelity witnesses for pure product states using Protocol 1 for each of the single-mode subsystems in parallel.We make use of the properties of heterodyne detection in order to further extend the class of target states for which a tight fidelity witness can be efficiently obtained.
Secondly, a passive linear transformation followed by single-mode Gaussian unitary operations before unbalanced heterodyne detection is equivalent to performing balanced heterodyne detection directly, then post-processing efficiently the classical samples.Formally: We give a proof in Appendix F. A striking consequence is that a certain class of quantum operations before a balanced heterodyne measurement can be inverted efficiently by unbalancing the detection and perfoming classical operations on the samples instead.This result generalises a technique used in continuous-variable quantum key distribution protocols based on heterodyne detection, as well as in a scheme for homomorphic encryption of linear optics quantum computation [69], which allows one to apply a random unitary operation at the level of the classical samples (i.e., a multiplication of the vector of samples by a unitary matrix) rather than directly on the quantum state [70].In particular, for such a transformation V , if a multimode pure product state m i=1 |ψ i can be efficiently verified using balanced heterodyne detection, then the state V m i=1 |ψ i can be efficiently verified using unbalanced heterodyne detection by postprocessing classically the samples obtained.
This leads us to our multimode fidelity witness estimation protocol.Like for Protocol 1, we give the version of our protocol under i.i.d.assumption, and discuss afterwards its version without i.i.d.assumption: 1. Measure all m subsystems of N copies of ρ with unbalanced heterodyne detection with unbalancing parameters ξ = (ξ 1 , . . ., ξ m ), obtaining the vectors of samples γ (1) , . . ., γ (N ) ∈ C m .

Compute the fidelity witness estimate
The value W ψ obtained is an estimate of a tight lower bound on the fidelity between the M remaining copies of ρ and M copies of the m-mode target state |ψ = Ŝ(ξ) D(β) Û m i=1 |C i .The efficiency and completeness of the protocol are summarised by the following result: With the notations of Protocol 2, for all i ∈ {1, . . ., m}, let n i be the number of indices j ∈ {1, . . ., m} such that where the constant prefactor depends on the choice of the free parameters p i and η i in Protocol 2. In particular, if the tested state is perfect, i.e., We give a detailed version of the theorem together with its proof in Appendix G, which combines Theorem 1 with Lemmas 2 and 3.
Note that Protocol 2 uses Protocol 1 as a subroutine.Writing p = min i p i and c = max i c i , with In particular, compared to the single-mode case, the overhead for certifying m-mode states is only O(m which is a huge improvement compared to, e.g., the exponential scaling of tomography [71].Since the parameter p ∈ N * may vary freely, the scaling can be brought arbitrarily close to O(( M m ) 2 log(m) log( 1 δ )).The constant prefactor increases with p, and an optimisation can yield the optimal choice for p which minimises N 2 when M , m, and δ are fixed.Moreover, the choice of the other free parameters 0 < η i < 1 also contributes to minimising the constant prefactor.
Like for Protocol 1, the efficiency of Protocol 2 can also be refined by taking into account the expression of the single-mode target core states in the Fock basis.
Without the i.i.d.assumption, we obtain an equivalent protocol by using the version of our single-mode fidelity estimation protocol which does not assume i.i.d.state preparation for estimating the single-mode fidelities.Like for Protocol 1, the final protocol is nearly identical, up to slight differences in the classical post-processing.Removing the i.i.d.assumption then comes at the cost of an increased number of measurements necessary for a polynomial witness precision and confidence interval, corresponding to a polynomial overhead.We give a detailed analysis in Appendix H.
The class of states for which fidelity witnesses can be efficiently estimated using Protocol 2 includes classically intractable quantum states, such as the output states of Boson Sampling experiments, as we detail in the next section.

Verification of Boson Sampling
The setup of Boson Sampling with input singlephotons consists in n single-photon states fed into an interferometer over m modes, together with the vacuum state over m − n modes, and measured either with single-photon threshold detection [10] or with unbalanced heterodyne detection [44,45] (Fig. 3).The output probabilities are related to permanents of complex matrices, which are hard to compute [72] and even to approximate [10].A striking consequence is that classical computers cannot sample exactly from the corresponding probability distributions unless the polynomial hierarchy collapses, contradicting a widely believed conjecture from complexity theory.Aaronson and Arkhipov importantly showed that even an approximate version of their model is still hard to sample for classical computers, provided two additional conjectures on the permanent of random Gaussian matrices hold true [10].More precisely, they showed under these conjectures that sampling from a probability distribution that has small constant total variation distance with an ideal Boson Sampling distribution is classically hard in the so-called antibunching regime n = O( √ m) (in fact they used m = Ω(n 6 ) and conjectured that m = Ω(n 2 ) was sufficient, as it was later shown [73]).
The output state of a Boson Sampling experiment with n photons fed into an interferometer over m modes, described by an m-mode passive linear transformation Û reads Û |1 . . . 10 . . .0 .This state is of the form of Eq. ( 1), with Ĝi = 1 for all i ∈ {1, . . ., m}, |C i = |1 for i ∈ {1, . . ., n} and |C j = |0 for j ∈ {n+1, . . ., m}.Our fidelity witness estimation protocol from the previous section can thus be applied to certify efficiently the fidelity of the output state of a Boson Sampling experiment with the ideal target output state, using only balanced heterodyne detection.Then, our protocol yields a certificate of the total variation distance with the target probability distribution for any observable by the Fuchs-Van de Graaf inequality [35] and standard properties of the trace distance [34].Indeed, for any observable Ô and any states ρ and σ, where the left hand side is the total variation distance between the probability distributions corresponding to a measurement of the observable Ô for the states ρ and σ, D is the trace distance and F the fidelity.
We give the version of our Boson Sampling verification protocol under i.i.d.assumption and discuss afterwards its version without i.i.d.assumption. 1. Measure all m subsystems of N copies of ρ with balanced heterodyne detection, obtaining the vectors of samples γ (1) , . . ., γ (N ) ∈ C m .

Compute the fidelity witness estimate
and measure all m subsystems of the remaining M copies of ρ with unbalanced heterodyne detection (resp.single-photon threshold detection), obtaining the M vectors of samples (s (1) , . . ., The efficiency and completeness of the protocol are summarised by the following result: We use the notations of Protocol 3.With probability 1 − δ, either the protocol aborts or the samples s (1) , . . ., s (M ) are drawn from a probability distribution which has less than √ λ total variation distance with the ideal Boson Sampling probability distribution, whenever N ≥ N 3 , with where the constant prefactor depends on the choice of the free parameter η in Protocol 3.Moreover, if the tested state is perfect, i.e., ρ ⊗N +M = |ψ ψ| ⊗N +M , then Protocol 3 accepts with probability 1 − δ.
We give a proof of the theorem in Appendix I, which follows from optimising the efficiency of Protocol 2 for Boson Sampling output states, together with Eq. ( 15).The choice of the free parameters p ∈ N * even and 0 < η < 1 contributes to minimising the constant prefactor.Additionnally, the choice of between 0 and λ can be made to ensure robustness of the protocol with Eq. ( 13): if is closer to 0 then the protocol will reject a good state with less probability.On the other hand, this increases the number of samples needed with Eq. (16).     the case of photonic quantum computing, where the measurement apparatus can be easily brought apart from the rest of the experiment, allowing to treat the latter as an untrusted black-box.Moreover, we obtain a version of Protocol 3 without the i.i.d.assumption by using the version of Protocol 2 which does not assume an i.i.d.behaviour for the prover from the previous section.This slightly modifies the classical post-processing, requiring the verifier to discard at random a small fraction of the samples and perform an energy test using another fraction of the samples.This comes at the cost of a polynomial overhead for the number of measurements needed for a polynomial confidence, and we refer to Appendix H for the detailed derivation.By simply changing the unbalancing of heterodyne detection of the output modes of a Boson Sampling interferometer, one can switch between verification of Boson Sampling output states and demonstration of quantum speedup with continuous-variable measurements within the same experimental setup.Indeed, Boson Sampling interferometers with unbalanced heterodyne detection are hard to sample classically when the unbalancing of the heterodyne detec-tion is not too small [44,45], but their output can be efficiently checked simply by switching to balanced heterodyne detection with Protocol 3.This can be done within the same experimental setup using a tunable beam splitter (Fig. 3).However, showing the hardness of approximate sampling for Boson Sampling with Gaussian measurements is an important step before an experimental demonstration of quantum speedup with continuous-variable Gaussian measurements.Alternatively, by switching between balanced heterodyne detection and single-photon threshold detectors, for example with a movable mirror [75], one can switch between verification of Boson Sampling output states and demonstration of quantum speedup with discrete variable measurements, for which approximate sampling hardness is demonstrated, assuming two conjectures on the permanent of random Gaussian matrices and the conjecture that the polynomial hierarchy does not collapse [10].

Discussion
In this work, we have introduced various methods for efficiently retrieving information about single-mode and multimode continuous-variable quantum states.
We have generalised heterodyne estimates both from [59] and [60], obtaining the best of both worlds: a reliable single-mode fidelity estimation protocol with heterodyne detection which does not make an assumption of bounded support for the measured state, while remaining efficient (Protocol 1).
With this fidelity estimation protocol for single subsystems, we have derived a fidelity witness estimation protocol with heterodyne detection for a large class of quantum states over multiple subsystems (Protocol 2).This result stems from two observations: on the one hand, if all the subsystems of a quantum state are close to pure states, then this quantum state is close to the tensor product of the pure states; on the other hand, a large class of quantum operations before heterodyne detection can be reversed via classical post-processing.
Applying our fidelity witness estimation protocol to the case of Boson Sampling, we have shown how to verify efficiently Boson Sampling experiments using heterodyne detection (Proto-col 3).Our verification applies to the original Boson Sampling proposal [10], where the sampling giving rise to quantum speedup is performed using single-photon threshold detectors, as well as to related proposals [44,45], where the sampling is performed using unbalanced heterodyne detection instead.In the latter case, the verification and the computational runs only differ by the unbalancing of the verifier's heterodyne detection.Using a tunable beam splitter, this means that Boson Sampling and its verification can be implemented within the same setup.Note however that there is no rigorous approximate sampling hardness result for Boson Sampling with continuous-variable measurements, which is a crucial point before any experimental demonstration.It will likely require extending anticoncentration and average-case hardness conjectures for the permanent or the hafnian to the loop hafnian [76,77,78].On the other hand, one can already perform a demonstration of quantum speedup with Boson Sampling and discrete variable measurements, by using balanced heterodyne detection for the verification and singlephoton threshold detection for the sampling.For all our protocols, we have derived a version with i.i.d.assumption and another without this assumption.The two versions share the same structure, and given that the latter only comes at the cost of a polynomial overhead in the number of measurements and additional classical postprocessing, using the i.i.d.version of our protocol for verifying a large-scale Boson Sampling experiment would already provide a convincing evidence of quantum speedup with photonic computing.
To that end, it would be interesting to optimise further our protocols for an experimental implementation, by finding the optimal choice for the free parameters and by analysing the case where the verifier has access to imperfect rather than ideal heterodyne detection [59].As previously mentioned, it would also be interesting to show the hardness of approximate sampling for Boson Sampling with Gaussian measurements, as it would allow a verifier restricted to single-mode Gaussian measurements to perform both the verification and the hard sampling task.
Demonstration of quantum speedup is likely to be achieved not via a single experiment but rather through numerous experimental realisa-tions competing with always-improving classical simulation algorithms [79,80,19,81,20,82].In that context, but also long after the quantum speedup milestone is achieved, efficient methods for reliable certification of quantum devices will be of utmost importance [83].Beyond the quantum speedup milestone, our methods will find various applications for the reliable certification of multimode continuous-variable quantum states in existing and upcoming experiments.

Appendix
In this appendix we provide the proofs for the results stated in the main text.

A Unbalanced heterodyne detection
We relate the POVM for unbalanced heterodyne detection to the experimental detection setup of Fig. 2. Due to the tensor product structure, we only consider the single-mode setting.We use the definitions for all ξ ∈ C, and R(θ) = exp[iθâ † â] for all θ ∈ [0, 2π], for the displacement operator, the squeezing operator, and the phase operator, respectively.
When the squeezing parameter is real, the derivation can be found in Ref. [44] (note the difference in the definition of the squeezing operator with ours): the POVM elements for heterodyne detection with unbalancing parameter r ∈ R are given by for all α ∈ C, where |α, r = Ŝ(r) D(α) |0 .In this case, the squeezing parameter is related to the experimental setup by r = log T R .The general case ξ = re iθ with θ = 0 is obtained by adding a phase operator R − θ 2 in the input and changing the local oscillator amplitude from α to αe − iθ 2 (note that this may be done at the level of the classical samples, i.e., by multiplying the outcomes of the heterodyne detection by a phase e − iθ 2 ).Hence, the POVM element corresponding to the detection setup in Fig. 2 is given by where we used R θ 2 |0 = |0 .Hence Π ξ α = 1 π |α, ξ α, ξ|, which completes the proof.

B Proof of Lemma 1
We recall notations from the main text.Let us introduce for k, l ≥ 0 the polynomials for z ∈ C, which are, up to a normalisation, the Laguerre 2D polynomials.For all k, l ∈ N, we give with these polynomials the functions [60,84] for all z ∈ C and all 0 < η < 1.The function z → f k,l (z, η), being a polynomial multiplied by a converging Gaussian function, is bounded over C. Now let us define, for all p ∈ N * , all z ∈ C and all 0 < η < 1, the more general bounded functions We start by showing the following result, which generalises Lemma 3 of [84]: where the function g k,l is defined in Eq. ( 22).
Proof.We prove the lemma by induction on p ∈ N * .By Lemma 3 of [84], with E = +∞, we have, for all single-mode (mixed With Eq. ( 22), this gives hence the lemma is proven for p = 1 since q−1 0 = 1.Assuming the lemma holds for some p ∈ N * , Now by Eq. ( 24), with k = k + p and l = l + p, we have so that so the last expression in Eq. ( 28) simplifies as which completes the proof of Lemma 4.

An immediate application of Lemma 4 gives
where we used the fact that ρ is hermitian positive semidefinite in the third line, Cauchy-Schwarz inequality in the fourth line and Tr(ρ) = 1 in the last line.This last bound is tight: if q 0 ≥ p is the integer achieving the maximum in the last equation, then the inequality is an equality for ρ = |q 0 q 0 |.We have where the right hand side in the last line is an increasing function of q.Hence, for η ≤ , which corresponds to q = p in the above equation, we have for all q ≥ p, so that . With Eq. ( 31) we finally obtain Note that for η < 1 but η > , we can find the first integer q 0 such that η ≤ (q 0 −p+1)(q 0 +1) q 0 √ (k+q 0 +1)(l+q 0 +1) by incrementing q, since is an increasing function of q which goes to 1 when q goes to +∞.The bound above then gives While this can be useful to optimise further the efficiency of the protocol, we will assume for the sake of obtaining closed analytical expressions that η ≤ .

C Proof of Theorem 1
In this section, we prove the efficiency of Protocol 1, which we recall below: Protocol 1 (Single-mode fidelity estimation).Let c ∈ N * and let |C = c−1 n=0 c n |n be a core state.Let also N, M ∈ N * , and let p ∈ N * and 0 < η < 1 be free parameters.Let ρ ⊗N +M be N + M copies of an unknown single-mode (mixed) quantum state ρ.

Compute the fidelity estimate
The value F C (ρ) M obtained is an estimate of the fidelity between the M remaining copies of ρ and M copies of the target core state |C .Defining the efficiency of the protocol is summarised by the following result: With the notations of Protocol 1 and η ≤ η C,p , we have with probability greater than 1 − P iid C , where where K C,p is a constant independent of ρ given in Eq. (62).
Proof.Let us prove the theorem for M = 1.The general case is then an immediate corollary by replacing by M , since F (C ⊗M , ρ ⊗M ) = F (C, ρ) M and by Lemma 2 of [84]: for all M ∈ N * , all > 0 and all a, b ∈ [0, 1] (note that when the estimate F C (ρ) obtained is bigger than 1 we may replace it by 1 instead without loss of generality).
The proof combines Lemma 1 with Hoeffding inequality.The latter is expressed as follows: Lemma 5. (Hoeffding) Let λ > 0, let n ≥ 1, let z 1 , . . ., z n be i.i.d.complex random variables from a probability density D over R, and let g : C → R be a bounded function and let G ≥ max z g(z)−min z g(z).Then In order to apply Lemma 5, we first derive an analytical bound for the functions g k,l : Lemma 6.For all k, l ∈ N, all z ∈ C, all p ∈ N * and all η < 1, Proof.By Lemma 6 of [84], we have for all i, j ∈ N, all z ∈ C and all η < 1.Hence, with Eq. ( 22), for all k ≤ l ∈ N, for all z ∈ C and all η < 1.The same reasoning holds for l ≤ k ∈ N and we obtain for all k, l ∈ N, all z ∈ C, all p ∈ N * and all η < 1.
Let |C = c−1 n=0 c n |n be a core state.Recalling the definition from the main text, with Lemma 6 we obtain: for all z ∈ C, all p ∈ N * and all η < 1, where Let N ∈ N * , let ρ = k,l≥0 ρ kl |k l| be a density matrix and let α 1 , . . ., α N be i.i.d.samples from balanced heterodyne detection of N copies of ρ.Let us write With Lemma 5, for all λ > 0, for all z ∈ C, all p ∈ N * , where we defined the range of the function z → η c g C (z, η) and where we used Eq. ( 47) and where we used Eq. ( 46) in the first line, Lemma 1 in the third line for all 0 ≤ k, l ≤ c − 1, and where we defined For η ≤ p+1 p(p+c) , combining Eqs. ( 50) and (52) yields with probability greater than 1 − P iid C , where We now optimise over the choice of λ and η for a given precision > 0. Namely, fixing λ + η p A (p) this amounts to maximising λ 2 µ 2c p , with λ + µ = .The optimal choice is given by: This in turn gives With this choice for η, for small enough we will always have η ≤ p+1 p(p+c) .To ensure that this is always the case independently of the value of , we may fix instead Then, plugging these expressions in Eq. ( 55) we finally obtain with probability greater than 1 − P iid C , where where and where the constants C are defined in Eqs. ( 48) and ( 53), respectively.Hence, the number of samples needed for a precision > 0 and a failure probability δ > 0 scales as Note that, as indicated in Eq. ( 36), the choice of η in Eq. ( 37) is not optimal in general, and an optimisation depending on the expression of the target core state leads to a better prefactor in Eq. (63).
D Removing the i.i.d.assumption for Protocol 1 In this section, we derive a version of Protocol 1 which does not assume i.i.d.state preparation.The protocol is as follows: Protocol 4 (General single-mode fidelity estimation).Let c ∈ N * and let |C = c−1 n=0 c n |n be a core state.Let also N, M ∈ N * , and let p ∈ N * , E, S ∈ N and 0 < η < 1 be free parameters.Let ρ N +M be an unknown quantum state over N + M subsystems.We write

Compute the fidelity estimate
Note that this protocol differs from Protocol 1 by two additional classical post-processing steps, steps 2 and 3.These steps are part of a de Finetti reduction for infinite-dimensional systems [68] detailed in [84].In particular, step 3 is an energy test, and the energy parameters E and S should be chosen to guarantee completeness, i.e., if the perfect core state is sent then it passes the energy test with high probability.
The value F C (ρ) M obtained is an estimate of the fidelity between the remaining state ρ M over M subsystems and M copies of the target core state |C .The efficiency of the protocol is summarised by the following result: Theorem 5. Let > 0. With the notations of Protocol 4 and η ≤ η C,p , defined in Eq. (37), or the protocol aborts in step 3, with probability greater than 1−(P support

C
), where where G C,p is a constant independent of ρ given in Eq. (78).In particular, setting when 1 = O(poly M ), or the protocol aborts in step 3, with probability greater than 1 − P C , where for Proof.Theorem 5 is an optimised version of Theorem 5 of [84], where the heterodyne estimates f k,l (z, η) defined in Eq. ( 5) are replaced by the more efficient generalised estimates g k,l (z, η), with an additional free parameter p.The proof thus is nearly identical, and amounts to replacing occurences of the estimates f k,l (z, η) by the estimates g k,l (z, η).In particular, we use our Eq.( 52) instead of Theorem 1 of [84].For a target core state |ψ = |C , this amounts to replace by where C is defined in Eq. (53).Moreover, we use our Eq.( 47) instead of Lemma 6 of [84].For a target core state |ψ = |C , this amounts to replace where C is defined in Eq. (48).Hence, the constants used in the proof of [84] are mapped as: Note that we also adapt the notations from [84] for the global coherence of this work as: and we write c the support size of |ψ = |C instead of E + 11 .Ref. [84] sets η = mK ψ and defines for brevity With m → M , we thus set η = , or equivalently By Eq. ( 74), K ψ is mapped to , and the expression C ψ in Eq. ( 76) is thus replaced here by This expression is independent of ρ and depends on M and without loss of generality, since is upper bounded by its value for /M = 1.With these mappings, following the proof of [84] we obtain: or R > S, i.e., the protocol aborts in step 3, with probability greater than 1 ), where where G C,p is defined in Eq. (78).Setting for simplicity K = O(N ) and N = qQ, we find that all above probabilities are at least inverse polynomially small in M and 1 when and E = O(log M ) and S = o(Q).In particular, setting with Eq. ( 79) we obtain when 1 = O(poly M ), with probability 1 − P C , where for Since the parameter p ∈ N * may vary freely, the asymptotic scaling of the number of samples N needed for a precision > 0 and a polynomial confidence thus is given by N = O( M 7 4 ).Moreover, we set a choice for η in Eq. ( 77) for simplicity, but the choice of both p and η can further improve the efficiency of the protocol.

E Proof of Lemma 2
In this section, we prove a generalised version of Lemma 2 for tensor product of Hilbert spaces.Setting M = 1, we retrieve the Lemma in the main text: Proof.Since |ψ 1 , . . ., |ψ m are pure states, and for all i ∈ {1, . . ., m}, where 1 is the identity operator.The right hand side of Eq. ( 86) is obtained by writing F (ρ m×M , m i=1 ψ ⊗M i ) as a telescopic sum: (89) by linearity of the trace, where we used Tr (ρ m×M ) = 1.This gives with Eqs. ( 87) and (88).The left hand side of Eq. ( 86) is obtained using: for all i ∈ {1, . . ., m}, and summing over i.
Note that this result is not restricted to a particular protocol for fidelity estimation nor to the continuous-variable regime: if the fidelities of the single subsystems of a quantum state over m subsystems with some target pure states are higher than 1 − λ m , for some λ > 0, then by Lemma 2 the general state has fidelity at least 1 − λ with the tensor product of the m pure states.

F Proof of Lemma 3
In this section, we prove Lemma 3 from the main text which we recall below: Lemma 8. Let β, ξ ∈ C m and let V = Ŝ(ξ) D(β) Û , where Û is an m-mode passive linear transformation with m × m unitary matrix U .For all α ∈ C m , let γ = U α + β.Then, Proof.Recall that the POVM elements for a tensor product of unbalanced heterodyne detections with unbalancing parameters ξ ∈ C m are given by for all α, ξ ∈ C m , where |α, ξ = Ŝ(ξ) D(α) |0 is a tensor product of squeezed coherent states.
The POVM elements of a tensor product of single-mode balanced heterodyne detection are given by Π 0 α , for all α ∈ C m , and we have In particular, a single-mode squeezing followed by a single-mode balanced heterodyne detection can be simulated by performing directly an heterodyne detection unbalanced according to the squeezing parameter.One retrieves balanced heterodyne detection by setting the unbalancing parameter to 0 and homodyne detection by letting the modulus of the unbalancing parameter go to infinity.Passive linear transformations correspond to unitary transformations of the creation and annihilation operators of the modes.These transformations, which may be implemented by unitary optical interferometers, map coherent states to coherent states: if Û is a passive linear transformation and U is the unitary matrix describing its action on the creation and annihilation operators of the modes, an input coherent state |α is mapped to an output coherent state Û |α = |U α , where U α is obtained by multiplying the vector α by the unitary matrix U .Hence, the POVM elements of a passive linear transformation followed by a product of single-mode balanced heterodyne detections are given by for all α ∈ C m .This implies that the passive linear transformation Û † followed by a tensor product of single-mode balanced heterodyne detections can be simulated by performing the heterodyne detections first, then multiplying the vector of samples obtained by U .A similar property holds with single-mode displacements: since displacements map coherent states to coherent states, up to a global phase, by displacing their amplitude, a single-mode displacement followed by a single-mode heterodyne detection can be simulated by performing the heterodyne detection first, then translating the sample obtained according to the displacement amplitude.In particular we have for all α, β ∈ C m , where α + β = (α 1 + β 1 , . . ., α m + β m ).

Compute the fidelity witness estimate
The value W ψ obtained is an estimate of a tight lower bound on the fidelity between the M remaining copies of ρ and M copies of the m-mode target state |ψ = Ŝ(ξ) D(β) Û m i=1 |C i .The efficiency and completeness of the protocol are summarised by the following result: with probability 1 − P iid W , where where for all i ∈ {1, . . ., m}, K C i ,p i is a constant independent of ρ, defined in Eq. (62).Moreover, if the tested state is perfect, i.e., ρ ⊗N +M = |ψ ψ| ⊗N +M , then W ψ ≥ 1 − with probability 1 − P iid W .
Proof.We use the notations of the protocol above and we first prove the theorem in the case where the target state is a tensor product of single-mode pure states, i.e., Ŝ(ξ) D(β) Û = 1 and α (k) = γ (k) for all k ∈ {1, . . ., N }.We thus consider for the target state a tensor product |ψ = m i=1 |C i of core states.For all i ∈ {1, . . ., m} we write |C i = c i −1 n=0 c i,n |n , where c i ∈ N * , and ρ i = Tr {1,...,m}\{i} (ρ) the reduced state of ρ over the i th mode.Let > 0. For all i ∈ {1, . . ., m}, with η i ≤ η C i ,p i , defined in Eq. ( 37), we have by Theorem 4, with m : with probability higher than 1 − P iid C i , where where K C i ,p i is a constant independent of ρ given in Eq. ( 62).Hence, taking the union bound of the failure probabilities with probability higher than 1 − P iid W , where By Lemma 7, with ρ m×M = ρ ⊗M and With Eq. ( 103) we obtain with probability 1 − P iid W , where where for all i ∈ {1, . . ., m}, K C i ,p i is a constant independent of ρ, defined in Eq. (62).Moreover, if the tested state is perfect, i.e., ρ ⊗N +M = |ψ ψ| ⊗N +M , then F (|ψ ⊗M , ρ ⊗M ) = 1.In that case, by Eq. ( 106), with probability 1 − P iid W .We now generalise the previous proof to the case where the target state is of the form where for all i ∈ {1, . . ., m} each of the states |C i = c i −1 n=0 c i,n |n is a core state, where Û is an m-mode passive linear transformation with m × m unitary matrix U and where ξ, β ∈ C m .To that end, we make use of the properties of heterodyne detection: by Lemma 8, writing V = Ŝ(ξ) D(β) Û , the POVM { V Π 0 α V † } α∈C m can be simulated with the POVM {Π ξ γ } γ∈C m by computing α = U † (γ − β), i.e., translating the vector of samples γ by the vector of complex amplitudes −β and multiplying the vector obtained by the m × m unitary matrix U † .Formally, let ρ be an m-mode (mixed) state, then, In particular, computing the fidelity witness estimate W ψ in Protocol 2, i.e., using samples γ from unbalanced heterodyne detection of copies of ρ post-processed as α with probability 1 − P iid W by Eq. ( 106).Given that we finally obtain with probability 1 − P iid W , where P iid W is defined in Eq. ( 100).
The estimate of the fidelity witness obtained is tight when the fidelity is close to 1. Writing n i the number of indices j ∈ {1, . . ., m} such that c j p j = c i p i , for all i ∈ {1, . . ., m}, the number of samples needed for a precision > 0 of the estimate W ψ and a confidence 1 − δ scales as with Eq. (100).
H Removing the i.i.d.assumption for Protocol 2 In this section, we derive a version of Protocol 2 which does not assume i.i.d.state preparation.
For N, M ∈ N * , let H i and K k denote single-mode infinite-dimensional Hilbert spaces, for all i ∈ {1, . . ., m} and all k ∈ {1, . . ., N + M }.Let H be an m × (N + M )-mode infinite-dimensional Hilbert space.The following spaces are isomorphic: For any ρ m×(N +M ) ∈ H, we define the m-mode reduced state over K ⊗m k , for all k ∈ {1, . . ., N + M }.We also define the (N + M )-mode reduced state over H
2. Discard the last Q vectors of samples. 5.For all i ∈ {1, . . ., m}, compute the mean F C i (ρ) of the function z → g (p i ) C i (z, η i ) (defined in Eq. ( 8)) over the samples α

Compute the fidelity witness estimate W
). Akin to Protocol 4 compared to Protocol 1, this protocol differs from Protocol 2 by two additional classical post-processing steps, steps 2 and 4.These steps are part of a de Finetti reduction for infinite-dimensional systems [68] detailed in [84].In particular, step 3 is an energy test, and the energy parameters E i and S i should be chosen to guarantee completeness, i.e., if the perfect state is sent then it passes the energy test with high probability.
The value W ψ obtained is an estimate of a tight lower bound of the fidelity between the remaining state ρ m×M over m × M subsystems and M copies of the m-mode target state |ψ .The efficiency of the protocol is summarised by the following result: Theorem 7. Let > 0. With the notations of Protocol 4 and η ≤ η C,p , defined in Eq. ( 37), or the protocol aborts in step 3, with probability greater than 1 − P W , where for Proof.For the proof we only consider the case where the target state is a tensor product of core states |ψ = m i=1 |C i , i.e., Ŝ(ξ) D(β) Û = 1, since the generalisation to the multimode target states of the form Ŝ(ξ) D(β) Û m i=1 |C i using Lemma 3 is identical to the i.i.d.case, detailed in section G.Then, the proof is similar to the one of Theorem 6, using Theorem 5 instead of Theorem 4. Writing , for all i ∈ {1, . . ., m} (see Fig. 4 with N = 0), let be the reduced state over H ⊗M i of the state ρ m×M .By Theorem 5 and Eq. ( 83) in particular, for M large enough, or the protocol aborts in step 3, with probability greater than 1 − P C i , where for Replacing by m and taking the union bound of the failure probabilities we obtain or the protocol aborts in step 3, with probability greater than 1 − P W , where ).

I Optimised verification of Boson Sampling
We consider the case M = 1, the general case being retrieved by setting = M .Using Protocol 2 directly to certify the output states of a Boson Sampling experiment yields an efficiency N 2 = O(m 2+ 2 p log m log 1 δ ) for a constant error in the estimation, with a failure probability δ, where the constant prefactor scales with the free parameter p, giving an asymptotic scaling of O(m 2 log m log 1 δ ).We show in this section that a refined version of Protocol 2 using a single-mode fidelity witness as a subroutine rather than a single-mode fidelity estimate (from Protocol 1) actually yields an efficiency N 3 = O(m 2 log m log 1 δ ) already in the finite regime, thus proving Theorem 3. We use the fact that the target core states appearing in the expressions of Boson Sampling output states are Fock states, in order to control the analytical error from Lemma 4 derived in section B, which we recall below.
Let p ∈ N * , let k, l ∈ N and let 0 < η < 1.Let ρ = +∞ i,j=0 ρ ij |i j| be a density operator.Then, k,l (α, η)] = ρ kl + (−1) p+1 +∞ q=p ρ k+q,l+q η q q − 1 p − 1 where the function g k,l is defined in Eq. (22).In particular, when k = l, for a target Fock state |k k|, and for p even we obtain The sum on the right hand side is positive, so k,k (α, η)] is always a lower bound on the fidelity between the state ρ and the Fock state |k k| (with analytical error bounded independently of ρ by Eq. ( 36)).Moreover, this lower bound is tight when ρ is close to the Fock state |k k|.
With Eqs. ( 46), (47)  (137) In the case of Boson Sampling, the tensor product input is given by n Fock states |1 and m − n Fock states |0 .With the union bound of the failure probabilies above, using these single-mode fidelity witnesses in Protocol 2 instead of the fidelity estimate from Protocol 1, we obtain an estimate with precision of a tight multimode fidelity witness with failure probability when using N samples from tensor product single-mode heterodyne detection.In particular, for a constant precision and a failure probability δ, it is sufficient to use N 3 = O(m 2 log m log 1 δ ) samples from product single-mode heterodyne detection.The constant prefactor may still be optimised by playing with the choice of p even and 0 < η < 1.
Note that this optimised protocol can also be employed in the more general case of the certification of multimode quantum states of the form where each state |n i is a single-mode Fock state, where Û is an m-mode passive linear transformation and where each Ĝi is a single-mode Gaussian unitary operation, for all i ∈ {1, . . ., m}.
t e x i t s h a 1 _ b a s e 6 4 = " 9 4 w f d L R R B j C k y T C 3 Y 5 p 5 y T s 2 e g Q = " > A A A B 6 3 i c b V D L S g N B E O z 1 G e M r 6 t H L Y B A 8 h d 0 o 6 D H o x W M E 8 4 B k C b O T 2 e y Q e S w z s 0 v 3 P h a t a 1 4 x c w J / 4 H 3 + A C D 6 j k s = < / l a t e x i t > LO vacuum r < l a t e x i t s h a 1 _ b a s e 6 4 = " m c x b H M l p e x U Z O Q 8 0 I 6 5 L 7 N 8 g j U M = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m q o M e i F 4 8 t 2 F p o Q 9 l s J + 3 a z S b s b o Q S + g u 8 e F D E q z / J m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g p r 6 x u b W 8 X t 0 s 7 u 3 v 5 B + f C o r e N U M W y

Figure 2 :
Figure 2: A schematic representation of optical singlemode unblanced heterodyne detection of a state ρ, with unbalancing parameter ξ = re iθ .The dashed red lines represent beam splitters.LO stands for local oscillator, i.e., a strong coherent state.The blue cylinders are photodiode detectors.
and for all j ∈ {n + 1, . . ., m}, compute the mean F j,|0 (ρ) of the function z → g (p) |0 (z, η) over the samples α Protocol 3 allows one to verify a Boson Sampling experiment efficiently, by trusting only single-mode Gaussian measurements.Compared to other experimental platforms, the assumption of trusted measurements is especially relevant in t e x i t s h a 1 _ b a s e 6 4 = " m F b y g 4 g J W 0 Z X A v v u 8 2 L j F X T 3 N A g = " > A A A B 8 H i c b V B N S 8 N A E J 3 4 W e t X 1 a O X x S J 4 K k k V 9 F j 0 4 r G C / Z A 2 l M 1 2 0 i 7 d b M L u R i i x v 8 K L B 0 W 8 + n O 8 + W / c t j l o 6 4 O B x 3 s z z M w 1 S a x / L e j B P 0 I z q Q P O S M G i s 9 P L l d R e V A Y K 9 U d i v u D G S Z e D k p Q 4 5 6 r / T V 7 c c s j V A a J q j W H c 9 N j J 9 R Z T g T O C l 2 U 4 0 J Z S M 6 w I 6 l k k a o / W x 2 8 I S c W q V P w l j Z k o b M 1 N 8 T G Y 2 0 H k e B 7 Y y o G e p F b y r + 5 3 V S E 1 7 5 G Z d J a l C y + a I w F c T E Z P o 9 6 X O F z I i x J Z Q p b m 8 l b E g V Z c Z m V L Q h e I s v L 5 N m t e K d V 6 p 3 F + X a d R 5 H A Y 7 h B M 7 A g 0 u o w S 3 U o Q E M I n i G V 3 h z l P P i v D s f 8 9 Y V J 5 8 5 g j 9 w P n 8 A t n e Q W Q = = < / l a t e x i t > |1i < l a t e x i t s h a 1 _ b a s e 6 4 = " H 1 S a x / L e j B P 0 I z q Q P O S M G i s 9 P H l d R e V A Y K 9 U d i v u D G S Z e D k p Q 4 5 6 r / T V 7 c c s j V A a J q j W H c 9 N j J 9 R Z T g T O C l 2 U 4 0 J Z S M 6 w I 6 l k k a o / W x 2 8 I S c W q V P w l j Z k o b M 1 N 8 T G Y 2 0 H k e B 7 Y y o G e p F b y r + 5 3 V S E 1 7 5 G Z d J a l C y + a I w F c T E Z P o 9 6 X O F z I i x J Z Q p b m 8 l b E g V Z c Z m V L Q h e I s v L 5 N m t e K d V 6 p 3 F + X a d R 5 H A Y 7 h B M 7 A g 0 u o w S 3 U o Q E M I n i G V 3 h z l P P i v D s f 8 9 Y V J 5 8 5 g j 9 w P n 8 A u A K Q W g = = < / l a t e x i t > |0i < l a t e x i t s h a 1 _ b a s e 6 4 = " m F b y g 4 1 S a x / L e j B P 0 I z q Q P O S M G i s 9 P L l d R e V A Y K 9 U d i v u D G S Z e D k p Q 4 5 6 r / T V 7 c c s j V A a J q j W H c 9 N j J 9 R Z T g T O C l 2 U 4 0 J Z S M 6 w I 6 l k k a o / W x 2 8 I S c W q V P w l j Z k o b M 1 N 8 T G Y 2 0 H k e B 7 Y y o G e p F b y r + 5 3 V S E 1 7 5 G Z d J a l C y + a I w F c T E Z P o 9 6 X O F z I i x J Z Q p b m 8 l b E g V Z c Z m V L Q h e I s v L 5 N m t e K d V 6 p 3 F + X a d R 5 H A Y 7 h B M 7 A g 0 u o w S 3 U o Q E M I n i G V 3 h z l P P i v D s f 8 9 Y V J 5 8 5 g j 9 w P n 8 A t n e Q W Q = = < / l a t e x i t > |1i < l a t e x i t s h a 1 _ b a s e 6 4 = " H Y i + q G u W A D J / v O J t 8 E + C n N 5 z 0 G 4 = " > A A A B 8 H i c b V B N S 8 N A E J 3 4 W e t X 1 a O X x S J 4 K k k V 9 F j 0 4 r G C / Z A 2 l M 1 2 0 i 7 d b M L u R i i x v 8 K L B 0 W 8 + n O 8 + W / c t j l o 6 4 O B x 3 s z z M w

Figure 3 :
Figure 3: Boson Sampling with interferometer U with n photons over m modes and heterodyne detection with reconfigurable unbalancing parameter ξ.The detectors represented are homodyne detectors.Setting ξ = 0 (balanced detection) allows to certify efficiently the multimode output state with a fidelity witness.Setting ξ = 0, with |ξ| = Ω(2 − poly m )[74] allows to perform efficiently a sampling task which is hard for classical computers, unless the polynomial hierarchy collapses[44,45].
t e x i t s h a 1 _ b a s e 6 4 = " h g r b G y f U s 1 b / M r 5 + M 3 2 C k 4 D 9 z m A = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m q o M e i F 4 8 t 2 F p o Q 9 l s J + 3 a 3 S T s b o Q S + g u 8 e F D E q z / J m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g p r 6 x u b W 8 X t 0 s 7 u 3 v 5 B + f C o r e N U M W y x W M s C N 7 y y 6 u k X a t 6 F 9 V a 8 7 J S v 8 n j K M I J n M I 5 e H A F d b i D B r S A A c I z v M K b 8 + i 8 O O / O x 6 K 1 4 O Q z x / A H z u c P 1 t u M 9 Q = = < / l a t e x i t > m k < l a t e x i t s h a 1 _ b a s e 6 4 = " 8 K q b k l Z k 5 o G 1 N A g 1 i y R 6 8 e p Y K Pc = " > A A A B + 3 i c d V D L T g I x F O 3 g C / E 1 4 t J N I z F x N Z m B Q W B H d O M S E 3 k k g K R T O t D Q d i Z t x 0 g I v + L G h ca 4 9 U f c + T d 2 A B M 1 e p I m J + f c m 3 t 6 g p h R p V 3 3 w 8 q s r W 9 s b m W 3 c z u 7 e / s H 9 m G + p a J E Y t L E E Y t k J 0 C K M C p I U 1 P N S C e W B P G A k X Y w u U z 9 9 h 2 R i k b i R k 9 j 0 u d o J G h I M d J G G t j 5 H k d 6 H I Q 9 R U c c D S a 3 P D e w C 6 5 z X i y W / S p 0 n Y r r u e W a I d W K 7 9 d K 0 H P c B Q p g h c b A f u 8 N I 5 x w I j R m S K m u 5 8 a 6 P 0 N S U 8 z I P N d LF I k R n q A R 6 R o q E C e q P 1 t k n 8 N T o w x h G E n z h I Y L 9 f v G D H G l p j w w k 2 l S 9 d t L x b + 8 b q L D a n 9 G R Z x o I v D y U J g w q C O Y F g G H V B K s 2 d Q Q h C U 1 W S E e I 4 m w N n W l J X z 9 F P 5 P W k X H K z n F a 7 9 Q v 1 j V k Q X H 4 A S c A Q 9 U Q B 1 c g Q Z oA g z u w Q N 4 A s / W 3 H q 0 X q z X 5 W j G W u 0 c g R + w 3 j 4 B N 6 + U k Q = = < / l a t e x i t >

Figure 4 :
Figure 4: The notations for the subsystems of a quantum state ρ m×(N +M ) over m × (N + M ) modes, seen as N + M groups of m subsystems or as m groups of N + M subsystems.The red dots depict the single-mode subsystems of ρ m×(N +M ) .The reduced state of ρ m×(N +M ) corresponding to the i th group of N + M subsystems is denoted ρ N +M i , for all i ∈ {1, . . ., m}, and the reduced state of ρ m×(N +M ) corresponding to the k th group of m subsystems is denoted σ m k , for all k ∈ {1, . . ., N + M }.