Reliable quantum certification of photonic state preparations

Quantum technologies promise a variety of exciting applications. Even though impressive progress has been achieved recently, a major bottleneck currently is the lack of practical certification techniques. The challenge consists of ensuring that classically intractable quantum devices perform as expected. Here we present an experimentally friendly and reliable certification tool for photonic quantum technologies: an efficient certification test for experimental preparations of multimode pure Gaussian states, pure non-Gaussian states generated by linear-optical circuits with Fock-basis states of constant boson number as inputs, and pure states generated from the latter class by post-selecting with Fock-basis measurements on ancillary modes. Only classical computing capabilities and homodyne or hetorodyne detection are required. Minimal assumptions are made on the noise or experimental capabilities of the preparation. The method constitutes a step forward in many-body quantum certification, which is ultimately about testing quantum mechanics at large scales.

where D ∈ R 2m×2m is positive-definite and diagonal, with elements D 2j−1,2j−1 := s j ≥ 1 and D 2j,2j := s −1 j , for j ∈ [m], and O ∈ R 2m×2m and O ∈ R 2m×2m are orthogonal matrices. D describes m active single-mode squeezers in parallel, each one with squeezing parameter s j along the position quadrature. The maximum single-mode squeezing is s max := max 1≤j≤m {s j }. O and O , in turn, describe passive mode transformations that can be implemented by linear-optical networks of at most m(m − 1)/2 beam-splitters and single-mode phase shifters 8 . In the two settings considered here, i.e., for any t ∈ C G ∪ C LO , the unitarŷ U in equations (3) and (4) is such that O can be taken as the identity matrix. In the first setting, i.e., for t ∈ C G , this holds becauseÛ acts on the vacuum state vector |0 and any passive mode transformation maps the vacuum into itself. For the second setting, i.e., for t ∈ C LO , this holds simply because there we assume that the total transformation itself is passive, i.e., in that case it holds also that D = 1, so that S = O.
In both cases, coupling between different modes only takes place through the linear-optical network described by O. A general circuit can couple all m modes with each other, meaning that the quadrature operators of each output mode are linear combinations of those of all m input modes. However, often, each mode is only coupled to at most d ≤ m other modes. In these situations, O is a sparse matrix with at most 4md non-zero elements. More precisely, the columns of O are given by 2m orthonormal vectors (o (k) ) k∈[2m] each having at most 2d non-zero entries. Furthermore, since the position and momentum of each mode is coupled to at most the 2d quadratures of the same d modes, each pair o (2j−1) and o (2j) shares the same sparsity property, i.e., o (2j−1) and o (2j) have at least 2(m − d) zero entries in common, for all j ∈ [m].
Gaussian case. Using that in the Gaussian case S = O D and squaring equation (26) yieldŝ where O −1 = O T has been used and the trace is taken not over the Hilbert space but over the 2m × 2m matrix with operators as entries. Combining supplementary equation (2) and equations (27) and (29) yields Now we introduce the first moment vector γ ∈ R 2m and the symmetric second moment matrix Γ (1) ∈ R 2m×2m of p , with components γ l := r l p and Γ respectively. Since the matrix OD −2 O T is symmetric, it holds that so that we can rewrite supplementary equation (4) in terms of the observables which Arthur has access to as Supplementary Box 1 (Measurement scheme M G ). First moment vector: For each 1 ≤ l ≤ 2m, Arthur uses C 1 copies of p , with C 1 given by supplementary equation (62a), to measure the observabler l , obtaining an estimate γ l of the expectation value γ l = r l p . Second moment matrix: For each 1 ≤ l ≤ l ≤ 2m for which (OD −2 O −1 ) l,l = 2m k=1 o (k) l D −2 k,k o (k) l = 0, he uses C 2 copies of p , with C 2 given by supplementary equation (62b), to measure the observable 1 2 (r lrl +r l r l ), obtaining an estimate Γ (1) l,l of the expectation values Γ (1) l,l = Γ (1) l ,l in supplementary equation (5). Classical post-processing: He obtains the estimate F (0) of F (0) by replacing in supplementary equation (7) the actual expectation values Γ (1) and γ by the estimates Γ (1) and γ , respectively.
We will show later (see Supplementary Lemma 4 in Supplementary Note 3 and the discussion immediately after its proof) that the bound in supplementary equation (7) actually depends on at most 2mκ out of the 4m 2 entries of Γ (1) , with κ = 2 min{d 2 , m}, as defined in equation (7). Thus, only the 2mκ corresponding observables, and the 2m observables necessary for γ, as indicated in Supplementary Box 1, need to be measured. All these observables can be measured by homodyne detection 6 . Furthermore, in Supplementary Note 4 we show that only m+3 different measurement settings are required. Finally, by classical post-processing, Arthur recombines his estimates according to the third step of Supplementary Box 1 and obtains the fidelity estimate F (0) . This last step is also efficient in m.
With the definition in supplementary equation (14) and the fact that each projector P (i) is a symmetric matrix, supplementary equation (13) finally becomes Note that this is an explicit expression for F (n) in terms of the correlators in supplementary equation (14) that Arthur can measure. Due to the sparsity of O, each matrix P (i) has at most (2d) 2 non-zero entries. Then, it follows (see Supplementary Lemma 7 in Supplementary Note 3 for details) that the measurement of O m 4d 2 + 1 n observables, those listed in Supplementary Box 2, suffices for the estimation of the bound of supplementary equation (16). As in the Gaussian case, all these observables can be measured by homodyne detection 6 . Furthermore, in Supplementary Note 4 we show that at most m n 2 n+1 ≤ (2m) n /n! measurement settings are sufficient. Once again, by classical post-processing, Arthur recombines his estimates according to the third step of Supplementary Box 2 and obtains the fidelity estimate F (n) . Provided that n is constant, this last step is also efficient in m.

Supplementary Note 2 -Quantum certification of locally post-selected target states.
In this note, we extend our results to locally post-selected (m − a)-mode target states S t in C PLO . The entire treatment of the class C PLO is similar to, and follows directly from, that already seen for the classes C G or C LO . Therefore, instead of repeating all the details, we simply explain the specific differences.
The fidelity bound with post-selection. The first step is to derive the fidelity bound F S (n) given by equations (32) and (33). We proceed in a similar fashion to the Methods section in the main text. Due to equations (1) and (5), the facts that S t and t are pure, and the properties of the trace, it holds that where Tr S indicates partial trace over the Fock space of the m − a modes in S. Now, due to equation (27) with β = n, it holds that withF (n) the observable of equation (31). Using supplementary equations (18) and (19), we obtain the fidelity bound F S (n) of equations (32) and (33).
The derivation of the bound in supplementary equation (19) holds exactly the same if the m-mode target states t ∈ C LO are replaced by the more general target states t =Û |n n|Û † , withÛ any Gaussian unitary and |n any Fock-basis state, and the ancilla Fock-basis state |n A A is replaced by any generic a-mode pure product state on the modes A. What is not necessarily true for the more general post-selection scenario is that the resulting fidelity bound is tight for perfect experimental preparations. That is, there are post-selected target states for which the value of the resulting bound is too low to pass the certification test even if the preparation is perfect. Luckily, this does not happen for the experimentally more relevant case of linear-optical network target states post-selected with Fock-basis measurements, as we see in the following subsection.
Tightness of the fidelity bound for ideal experimental preparations with post-selection. Here, we show that the fidelity bound for linear-optical network target states post-selected with Fock-basis measurements is tight for ideal experimental preparations. That is, we show that, for S p = S t , just as in the cases without post-selection. We begin by expressing F where, in the second equality, the assumption that S p = S t and the definition of S t of equation (5) were used, and, in the third equality, equation (6) was used and the pure normalised target state vector was introduced, i.e., such that t = |Ψ t Ψ t |. Now, due to equation (4), we know that |Ψ t is an eigenstate of the total photon-number operatorn with eigenvalue n. In addition, since |n A A is also an eigenstate ofn and P(n A | t ) > 0 by assumption, the normalised joint system-ancilla state after the projection, Supplementary Box 3 (Certification test T PLO ). Settings adjustments: Idem as in T from Box 1. State request: Arthur provides Merlin with the classical specification n, S, a, and |n A A of the target state S t and requests a sufficient number of copies of it. Quantum measurements: He measures O m(4d 2 + 1) n multi-body correlators, each one involving between 1 and 2n + 1 modes, specified by the measurement scheme M PLO (see "The measurement scheme with post-selection" below), which can be done with a single local heterodyne setting throughout. Classical post-processing: By processing the measurement outcomes, he obtains a fidelity estimate F is the lower bound to F S given by equations (32) and (33). Accept-reject decision: If F (n) S < F T + ε, he rejects. Otherwise, he accepts.
also turns out to be an eigenstate ofn with eigenvalue n. Furthermore, sinceÛ † is passive,Û † |ψ t S ⊗ |n A A is an eigenstate ofn with eigenvalue n too and, therefore, eigenstate of n + 1 −n with eigenvalue 1. Thus, we can rewrite supplementary equation (21) as where the second equality is due to the fact that, on the eigenvalue-n eigenspace ofn, the observable n j=1n j equals the projector |1 n 1 n |. Finally, using that P(n A | t ) = Tr [ n A | A t |n A A ], with t defined by equation (4), and supplementary equation (23), one obtains that the right-hand side of supplementary equation (24) is equal to 1.
The certification test with post-selection. Next, in Supplementary Box 3, we present a test T PLO that works for post-selected target states in C PLO and which is a slightly modified version of the test T from Box 1. It is, of course, possible to unify both tests so as to account for all three classes of target states in one single test. We have, however, opted for splitting the tests into the two cases with and without post-selection to avoid an excessive notational overhead in Box 1.
The measurement scheme with post-selection. The measurement scheme to estimate F S (n) is essentially a replica of the scheme M LO to estimate F (n) , already described in detail in the Supplementary Box 2. Thus, instead of repeating all the details, we simply outline the concrete differences between M LO and M PLO . There are only three specific differences: 1. The moment tensors are now defined with respect to S p ⊗ |n A A n A | A instead of p . More precisely, we now need to estimate the tensors Γ (j) S ∈ R 2m×2m ⊗j , with elements Γ (j) S k1,l1,...,kj ,lj := r k1rl1 +r l1rk1 2 · · ·r kjrlj +r ljrkj 2 2. F S (n) is obtained dividing the expression on the right-hand side of supplementary equation (16) by P(n A | t ), and with 3. The presence of the divisor P(n A | t ) in F S (n) is the reason for the third difference. As discussed in Supplementary Lemma 11 in Supplementary Note 3, this divisor makes F S (n) more unstable than F (n) by a factor of 1/P(n A | t ). As a consequence, the number of copies of S p required to estimate each relevant moment of F S (n) is 2 instead of C ≤2(n+1) , as we discuss in Supplementary Lemma 12, also in Supplementary Note 3.
As is clear from supplementary equation (26), the estimation of Γ (j) S requires only the measurement of multi-body correlators among the (m − a) system modes in S. This is due to the facts that, after the post-selection, the system is in a product state with respect to the bipartition S versus A (see supplementary equation (23)), and that the quadrature operators in supplementary equation (26) can also be correspondingly grouped into two factors, one containing exclusively operators of modes in S and the other in A. Furthermore, since |n A A is a product state known to Arthur, he can efficiently calculate the expectation vale of any product of quadrature operators of modes in A with respect to |n A A . For instance, suppose that k 1 , l 1 , k 2 ∈ A and that l 2 , k 3 , l 3 . . . , k j , l j / ∈ A. Then, the corresponding 2j-th moment decomposes as and only the measurement of the (2j − 3)-th moment given by the second factor of supplementary equation (27) is required. As another example, consider the moments containing an odd number of quadrature operators of any A j -th mode. Since |n A A is a Fock-basis state, all these moments automatically vanish and need therefore not be measured at all.
Arthur can always efficiently obtain the Γ (j) S 's as a product of an (a priori known) expectation value with respect to |n A A of a multi-body product of quadrature operators of modes in A and a (measured) expectation value with respect to S p of a multi-body product of quadrature operators of modes in S, in a way analogous to the example of supplementary equation (27).
Formal statement of quantum certification of locally post-selected target states. Since the moments to be estimated, given in supplementary equation (26), are expectation values with respect to S p ⊗ |n A A n A | A , instead of p , a simple way to extend Theorem 3 to target states in C PLO is by redefining the variance upper bounds σ i . More precisely, taking σ i as an upper bound on the variances of any product of i phase space quadratures now in the state S p , we introduce the quantities for i ∈ [2(n + 1)]. With them, we define ς ≤i := max k≤i {ς k } as the maximal i-th generalised variance of S p . The parameters ς i quantify the maximal variances of random variables defined by products of i − j quadrature-measurement outcomes on S p renormalised by the expectation value of products of j quadrature operators with respect to |n A A , therefore automatically accounting for factorisations of the type of supplementary equation (27). They constitute non-tight upper bounds to the real variances. In particular experimental situations, tighter bounds can be found. Here, we are simply interested in taking advantage of the proof of Theorem 3 without introducing too much extra notational overhead, for which the definition of supplementary equation (28) is enough. Indeed, with these redefinitions, the following corollary follows straightforwardly from Theorem 3.
Supplementary Corollary 1 (Quantum certification of locally post-selected linear-optical network states). Under the same conditions and for the same t as in Theorem 3, test T PLO from Supplementary Box 3 is a certification test for S t ∈ C PLO and requires at most copies of a preparation S p with maximal 2(n + 1)-th generalised variance ς ≤2(n+1) , where λ > 0 is the same absolute constant as in Theorem 3.
Supplementary Corollary 1 is proven in Supplementary Note 3. Supplementary equation (29) corresponds to exactly the same expression as in equation (9) with the replacements σ → ς and ε → P(n A | t ) ε. The rescaling of ε with the factor P(n A | t ) originates directly from the new expression for the fidelity given in supplementary equation (18), which renders the fidelity bound more unstable than F (n) by the factor of 1/P(n A | t ). In most interesting cases, the post-selection success probability P(n A | t ) turns out to be exponentially small in a. Moreover, one can always come up with families of target states and post selection procedures for which P(n A | t ) decreases arbitrarily fast in m. In such cases, the scaling in supplementary equation (29) is not efficient in m, inheriting the inefficiency of the state preparation by local measurements and post selection. However, the scaling is efficient in 1/P(φ A | t ). That is, in every practical situation where state preparation via post selection is feasible, so is state certification. Interestingly, even for families of target states and post selection procedures for which P(φ A | t ) decays exponentially in a, the overall scaling (of the bound in supplementary equation (29)) with a is better than the scalings (of both the bounds in equation (9) and supplementary equation (29)) with n. Indeed, the bound in supplementary equation (29) grows, just like the bound in equation (9), faster than exponentially in n. Finally, the bound in supplementary equation (29) scales polynomially with all the other relevant parameters, including 1/ε. Thus, arbitrary m-mode target states from the class C PLO , with constant n, are certified by T PLO efficiently in m, 1/P(φ A | t ), and all the other relevant parameters.
Formal statement of robust quantum certification of locally post-selected target states. Finally, we show that our certification test for the locally post-selected target states of the class C PLO , T PLO , is also robust. In a way analogous to equation (10), where S ⊥ t is a normalised state such that Tr[ S t S ⊥ t ] = 0. Then, analogously to equation (12), we introduce the quantity is the observable given in equation (33). With these definitions, the following corollary of Theorem 5 holds true.
Supplementary Corollary 2 (Robust quantum certification of locally post-selected states). Under the same conditions as in Corollary 1, test T PLO from Supplementary Box 3 is a robust certification test for S t ∈ C PLO with fidelity gap The proof is identical to the one of Theorem 5 presented in Supplementary Note 3 but with the replacements F → F S ,

Supplementary Note 3 -Proofs of the theorems and corollaries.
Before going to the proofs, we devote two sections to establish necessary notation, review some known facts, and prove a general lemma.
Norms. Here, we introduce some helpful notation used in the proofs and review a few facts about norms on finite dimensional vector spaces. The max norm · max of a tensor is the largest of the absolute values of its entries. For a matrix A, for example, we denote the vector p-norm of a vector a by a p and the Schatten p-norm of a matrix A by A p , which is defined to be the vector p-norm of the vector of its singular values. For any matrix A, we define vec(A) to be a vector containing all the entries of A (in some order). Then one can see that and For the vector and Schatten p-norm of vectors with N elements and N × N matrices, respectively, the following inequalities hold Because the Schatten ∞-norm is induced by the vector 2-norm, i.e., it follows that for any two vectors and x Reliable estimation of expectation values from samples. We continue by proving a general large-deviation bound for estimates of expectation values from a finite number of measurements on independent copies, which we need for the proofs of Theorems 2 and 3.
Supplementary Lemma 3 (Reliable estimation of multiple expectation values from samples). Let σ > 0, ρ be a state, and let A 1 , . . . ,Â N be observables with expectation values A i := Tr ρÂ i and variances bounded as Tr be the random variable given by the measurement statistics ofÂ i on state ρ; such that, in particular, the are independent random variables and the finite sample average over c measurements ofÂ i is the random variable Then, the {A i } i are independent and, for every > 0 and α ∈ [1/2, 1), it holds that Proof. The sample averages {A i } i are independent by definition. By Chebyshev's inequality, it holds that Since the {A i } i are independent random variables, this yields Finally, To finish the proof we upper bound Using that (see Supplementary Note 6) for all x ≥ 0 it follows that To simplify the right-hand side of this inequality, we use that, since α ≥ 1 2 ≥ e −1 , it holds that ln(1/α) ≤ 1 and, therefore, 2 + 2N ln(1/α) ≥ 4. So, using again that ln(1/α) ≤ 1, we finally arrive at Proof of Theorem 2 on certification of Gaussian states. We start by proving three auxiliary lemmas specific to the fidelity bound F (0) for the Gaussian case.
The first lemma provides upper-bounds on the number of elements of Γ (1) on which the fidelity bound F (0) depends.
Supplementary Lemma 4 (Sparsity of the Gaussian fidelity bound). F (0) depends on at most 2mκ elements of Γ (1) . We call these the relevant elements of Γ (1) .
Proof. Supplementary equation (7) can be written as The last term can, in turn, be expressed as Due to the sparsity of O, as described in Supplementary Note 1, each matrix s −2 ) T has at most 4d 2 non-zero elements. Hence, summing over j, we see that F (0) depends on at most 4m min{d 2 , m} = 2κm elements of Γ (1) .
Note that the counting argument following supplementary equation (50) does not take into account the fact that Γ (1) is symmetric. Taking this fact into account, we see that, from the 4d 2 relevant elements of Γ (1) that appear in each term of the trace in supplementary equation (50), only d(2d + 1) are independent. Thus, even though 2mκ entries of Γ (1) contribute to F (0) , only m min{d(2d + 1), 4m} ≤ 2mκ of them must actually be measured.
The second auxiliary lemma bounds the deviation of F (0) from F (0) in terms of the errors made in the estimation of the individual expectation values entering F (0) .
Supplementary Lemma 5 (Stability of the Gaussian fidelity bound). Let F (0) be defined like F (0) in supplementary equation (7) but with γ and Γ (1) replaced by γ and Γ (1) and let max := γ − γ max and ε (1) Then Proof. For convenience, we define the error vector and the error matrix The fidelity estimation error can then be written as Due to Hölder's inequality, According to Supplementary Lemma 4, F (0) depends on at most 2κm entries of E (1) . Without loss of generality we can hence omit all other elements of E (1) and thus take vec(E (1) ) as a vector with at most 2κm elements. Using this fact and the second inequality in supplementary equation (35) we obtain where we have used supplementary equation (34) in the last equality. Finally, putting everything together and using that, by definition, D −2 ∞ = s 2 max , we arrive at the inequality in supplementary equation (52).
The last auxiliary lemma shows that the estimate of the fidelity lower-bound for target states t ∈ C G obtained with the measurement scheme M G in Supplementary Box 1 is reliable. This lemma is potentially interesting in its own right in scenarios other than certification.
Supplementary Lemma 6 (Reliable estimation of the Gaussian fidelity bound). Let α ∈ (0, 1/2] and ε > 0. Let F (0) be defined like F (0) in supplementary equation (7) but with γ and Γ (1) replaced by γ and Γ (1) , where γ and Γ (1) are obtained as described by M G from copies of p , with C 1 and C 2 integers such that and (62b) Then, Proof. Our proof strategy is to show that, with probability at least 1 − α, the 2m elements of γ and the 2mκ relevant elements of Γ (1) are estimated within additive errors bounded as and If the inequalities in supplementary equation (64) are fulfiled, then, due to Supplementary Lemma 5, it holds that |F (0) −F (0) | ≤ ε.
Since all 2m estimates {γ l } l are sample averages over independent copies of p , the measurement outcomes to obtain the {γ l } l are all independent random variables, for each l described by the same probability distribution. Furthermore, by assumption, the variances of these variables are all upper-bounded by σ 1 . Analogously, the measurement outcomes to obtain all 2mκ relevant estimates {Γ (1) l,l } l,l are independent random variables with variances upper-bounded by σ 2 and described, for each l and l , by the same probability distribution. Hence, according to Supplementary Lemma 3, with the choice α = √ 1 − α, taking and is sufficient for both Since the {γ l } l and the {Γ l,l } l,l are independent random variables, supplementary equations (66) imply that Finally, inserting the definitions of max and ε Proof of Theorem 2. That the total number of copies of p (see supplementary equation (61)) needed for the certification test is asymptotically upper-bounded by equation (8) can be verified by straightforward calculation using supplementary equations (62). It remains to show that (i) if p = t , then T accepts with probability at least 1 − α, i.e., and (ii) if p is such that F < F T , then T rejects with probability at least 1 − α, i.e., To show (i), we first recall that, if p = t , F (0) = 1. With this, supplementary equation (63) in Supplementary Lemma 6 implies that Since, by assumption of the theorem, the total estimation error is such that ε ≤ 1−FT 2 , it holds that 1 − ε ≥ F T + ε. Substituting the latter inequality into supplementary equation (70) yields supplementary equation (68).
To show (ii), we first note that, since F (0) ≤ F for all p , if F < F T , then On the other hand, supplementary equation (63) implies also that

Inserting supplementary equation (71) into supplementary equation (72) yields supplementary equation (69).
Proof of Theorem 3 on certification of linear optical network states. We proceed analogously to the last section and present three auxiliary lemmas specific to the fidelity bound F (n) for the linear-optical case before proving Theorem 3.
To state the first lemma in a compact form we introduce the shorthand Γ := Γ (i) i=1,...,n+1 for the collection of all the moment tensors Γ (i) . Analogously, the collection of all the estimates Γ (i) of the moment tensors, defined in Supplementary Box 2, is denoted by Γ := Γ (i) i=1,...,n+1 . Supplementary Lemma 7 (Sparsity of the linear-optical fidelity bound). We write the fidelity bound F (n) given in supplementary equation (16) as where, for each j ∈ {0, . . . , n}, f j is the bilinear functional defined by For each j, the functional f j depends on at most n j (2d) 2j elements of Γ (j) and on at most n j 2m(2d) 2j elements of Γ (j+1) . We call these the relevant elements for f j . Moreover, F (n) depends on at most N ≤2(n+1) := (1 + 2m)(4d 2 + 1) n ∈ O m 4d 2 + 1 n (75) elements of Γ. We call these the relevant elements of Γ.
The subindex "≤ 2(n + 1)" in N ≤2(n+1) makes reference to the fact that 2j-th moments with j ∈ [n + 1] are taken into account. (73) and (74) can be checked by a straightforward calculation. We use again the sparsity of O, i.e., the property that its columns o (2j−1) and o (2j) have at least 2(m − d) zero element in common. Hence, each of the symmetric matrices P (j) , defined in supplementary equation (10), has at most (2d) 2 non-zero elements. Consequently, the projectors i∈Ω (j) µ P (i) and 1 ⊗ i∈Ω (j) µ P (i) in supplementary equation (74) have at most (2d) 2j and 2m(2d) 2j non-zero elements. This implies that the first trace inside the sum in supplementary equation (74) depends on at most 2m(2d) 2j elements of Γ (j+1) and the second trace inside the sum on at most (2d) 2j elements of Γ (j) . Hence, each f j depends on at most n j (2d) 2j elements of Γ (j) and on at most n j 2m(2d) 2j elements of Γ (j+1) . This proves the statements on the sparsity of the functionals f i . From this, it follows that F (n) depends on at most It is important to mention that, as in Supplementary Lemma 4 for the Gaussian case, the symmetry in supplementary equation (15) of each Γ (j) was not taken into account. Thus, even though the lemma gives the maximal total number of relevant elements that contribute to F (n) , many of them are not independent and must therefore not be measured.

Proof. Supplementary equations
The second auxiliary lemma upper-bounds the deviation of F (n) from F (n) in terms of the errors made in the estimation of the expectation values entering F (n) .
Supplementary Lemma 8 (Stability of the linear-optical fidelity bound). Let F (n) be defined like F (n) in supplementary equation (16) but with Γ replaced by Γ and let ε max := Γ − Γ max . Then Proof. For convenience, we define, for each j ∈ [n], the error tensor Using supplementary equation (73) and the fact that f j is linear, we write the fidelity estimation error as Applying Hölder's inequality and using that the Schatten ∞-norm of a tensor product of projectors is bounded by 1 yields where the matrixẼ (j) is defined element-wise byẼ k1,l1,...,kj ,lj , where k (j) := (k 1 , . . . , k j ) and l (j) := (l 1 , . . . , l j ).

Thanks to the first bound in supplementary equation (35) and supplementary equation (33), we arrive at
According to Supplementary Lemma 7, f j depends on at most n j 2m(2d) 2j elements ofẼ (j+1) and on at most n j (2d) 2j of E (j) . Without loss of generality we can hence omit, in supplementary equation (81), all other elements inẼ (j) andẼ (j+1) and thus take vec(Ẽ (j) ) and vec(Ẽ (j+1) ) as vectors with at most n j (2d) 2j and n j 2m(2d) 2j elements, respectively. Then the second bound in supplementary equation (35) yields Next, from supplementary equation (79), it follows that Finally, using n j 1/2 ≤ n j/2 and the binomial formula, we obtain the inequality in supplementary equation (77).
The last auxiliary lemma shows that the estimate of the fidelity lower-bound for target states t ∈ C LO obtained with the measurement scheme M LO in Supplementary Box 2 is reliable. This lemma is potentially interesting in its own right in scenarios other than certification.
Supplementary Lemma 9 (Reliable estimation of the linear-optical fidelity bound). Let α ∈ (0, 1/2] and ε > 0. Let F (n) be defined like F (n) in supplementary equation (16) but with Γ replaced by Γ , where Γ is obtained as described by M LO from copies of p , with N ≤2(n+1) an integer given by supplementary equation (75) and C ≤2(n+1) an integer given by Then, Proof. Our proof strategy is similar to that of Supplementary Lemma 6. That is, we show that, with probability at least 1 − α, the N ≤2(n+1) relevant elements of Γ are estimated within additive errors bounded as If this inequality is fulfiled, then, due to Supplementary Lemma 8, it holds that |F (n) − F (n) | ≤ ε. According to Supplementary Lemma 3, with the choice α = 1 − α, taking is sufficient to get Supplementary Lemma 10 (Sparsity of the locally post-selected linear-optical fidelity bound). The fidelity bound F S (n) , defined by the same expression as F (n) in supplementary equation (16) but divided by P(n A | t ) and with Γ replaced by Γ S , can be written as where, for each j ∈ {0, . . . , n}, f j is the same linear functional as in Supplementary Lemma 7, defined by supplementary equation (74). Moreover, F S (n) depends on at most N ≤2(n+1) elements of Γ S , with N ≤2(n+1) the same as in Supplementary  Lemma 7 and given by supplementary equation (75).
Proof. The proof of the lemma is analogous to that of Supplementary Lemma 7.
Supplementary Lemma 11 (Stability of the locally post-selected linear-optical fidelity bound). Let F (n) S be defined by the same expression as F (n) in supplementary equation (16) but divided by P(n A | t ) and with Γ replaced by Γ S , and let ε max := Γ S − Γ S max . Then Proof. The proof of the lemma is similar to that of Supplementary Lemma 8, with the differences already explained in "The measurement scheme with post-selection" in Supplementary Note 2.
Supplementary Lemma 12 (Reliable estimation of the locally post-selected linear-optical fidelity bound). Let α ∈ (0, 1/2] and ε > 0. Let F (n) S be defined like F (n) in supplementary equation (16) but with Γ replaced by Γ S , where Γ S is obtained as described in "The measurement scheme with post-selection" in Supplementary Note 2 from copies of S p , with N ≤2(n+1) an integer given by supplementary equation (75) and C ≤2(n+1) an integer given by Then, Proof. The proof of the lemma is analogous to that of Supplementary Lemma 9.
Proof of Supplementary Corollary 1. The proof is analogous to the proof of Theorem 3 but with Supplementary Lemmas 10, 11, and 12 playing respectively the roles of Supplementary Lemmas 7, 8, and 9.
Proof of Supplementary Theorem 5 on robust quantum certification of post-selected target states. Crucial for the proof of this theorem is the expansion in equation (10) of p in terms of t and ⊥ t , which allows for the definition of the parameter F (n) ⊥ in equation (12).
Proof of Supplementary Theorem 5. Theorems 2 and 3 imply that p is rejected with probability at least 1−α whenever F < F T . Thus, it remains to show that if p is such that F ≥ F T + ∆, with ∆ given by equation (13), then p is accepted with probability at least 1 − α, i.e., that In turn, from supplementary equation (86), we know that supplementary equation (95) is satisfied if So, we let F ≥ F T + ∆, with ∆ given by equation (13), and prove supplementary equation (96). Using equations (6), (10), and (12), we write F (n) as where the inequality holds because F ≥ F T + ∆ by assumption and because 1 − F (n) Using equation (13) we see that Substituting supplementary equation (98) into the right-hand side of the inequality in supplementary equation (97) finishes the proof.

Supplementary Note 4 -Number of measurement settings.
In this section, we upper-bound the number of local measurement settings required for the estimation of our fidelity lower bounds. We do this explicitly only for the Gaussian and linear-optical network target states, the cases of the post-selected target states following immediately from them.
Gaussian case. Here, we show that the 2md single-quadrature and the mκ two-quadrature observables listed in Supplementary Box 1, required for the measurement scheme M G , can all be measured using m + 3 different experimental arrangements using homodyne detection. We do this by explicitly describing a homodyne measurement strategy that features such a scaling. We note that, alternatively, all the observables required for M G can also be measured with a single experimental arrangement throughout using heterodyne detection. The measurement procedure for heterodyning is nevertheless the same as in the linear-optical case below. Therefore, here we only describe the procedure with homodyning.
Any single-mode phase-space quadrature operator can be directly measured with homodyne detection 6,[11][12][13] . The two-body observablesq jqk ,q jpk , andp jpk , for j = k, can be measured by simultaneously homodyning modes j and k independently. For all possible pairs of modes, this consumes m + 2 different homodyne settings: A single setting (q 1 ,q 2 , . . . ,q m ) for all the second moments of the form q jqk p ; another single setting (p 1 ,p 2 , . . . ,p m ) for those of the form p jpk p ; and the m settings (p 1 ,q 2 , . . . ,q m ), (q 1 ,p 2 ,q 3 , . . . ,q m ), . . ., and (q 1 , . . . ,q m−1 ,p m ) for those of the form q jpk p and p jqk p with j = k. In addition, all the single-body observablesq j ,p j ,q 2 j , andp 2 j , are measured also with these same settings. With this, we have accounted, so far, for all the first moments γ l and all the second moments Γ 2j−1,2j with j ∈ [m], correspond to the single-mode observables (q jpj +p jqj )/2. To measure these, Arthur can homodyne each mode j independently in the rotated quadrature (q j +p j )/ √ 2. This requires a single setting: (q 1 +p 1 )/ √ 2, (q 2 +p 2 )/ √ 2, . . . , (q 1 +p 1 )/ √ 2 . In this setting, he can estimate all the moments of the form (q j +p j ) 2 /2 p . The latter estimates, upon subtraction of q 2 j p /2 and p 2 j p /2, whose settings have already been accounted for, finally make it possible to indirectly estimate (q jpj +p jqj )/2 p , using the equation The last setting, plus the m + 2 ones already accounted for in the previous paragraph, yields a total of m + 3 different homodyne settings, as promised. Finally, a comment on the error estimation is in order. In any measurement strategy where moments are estimated indirectly, their errors must be obtained from those of the directly measured quantities via error propagation. For instance, in the strategy just described, the error of each Γ (1) 2j−1,2j needs to be calculated from those of (q j +p j ) 2 /2 p , q 2 j p , and p 2 j p . This leads, for each indirectly estimated moment, to an increase in the number of copies of p required to attain a given error. Nevertheless, this usually has no impact on the leading terms of the total resource scaling of the protocol. For example, in the described strategy, the global scaling given in equation (8) remains unaltered.
Linear-optical case. Here, we show that the N ≤2(n+1) ∈ O m 4d 2 + 1 n observables listed in Supplementary Box 2, required for the measurement scheme M LO , can all be measured using a single experimental arrangement throughout using heterodyne detection. Before that, however, we note that the majority of the these observables can be efficiently measured with homodyne detection in a similar fashion to that described in the Gaussian case above. Namely, the scheme M LO requires the measurement of products of an even number between 2 and 2(n + 1) of quadrature operators. For arbitrary n, all products of 2n quadrature operators for which each relevant j-th mode contributes exclusively either with powers of a single quadrature (q j orp j ) or with the quadratic polynomial (q jpj +p jqj )/2 can be measured with essentially the same homodyne procedure as the one described in the Gaussian case above. One can show that this can be done with at most O( m n ) different homodyne settings. This suffices to efficiently estimate the majority of the moments Γ (n) . However, the indirect estimation of higher-order moments, such as those involving powers ofq jpj +p jqj , can be more involved. For example, using operator identities of the type of supplementary equation (99) and the canonical commutation relations, one can see that This implies that the quartic polynomials in supplementary equations (100) can be indirectly estimated through homodyne measurements in the four local settingsq j ,p j ,q j +pj √ 2 , andq j −pj √ 2 , instead of just three as in the quadratic case. In general, the number of homodyne settings required for the indirect estimation of polynomials depending on both quadratures of a same mode grows with its order 14 . In return for that, the moments involving those polynomials require homodyne measurements on fewer modes. However, it is still useful to have another alternative. Heterodyning provides such an alternative.
Independent heterodyne detection on each mode requires a single experimental setting throughout. It implementsideally-the generalised measurement 1 π m/2 |α α| α∈C m , where |α is the m-mode coherent state of amplitude α := (α 1 , α 2 , . . . , α m ) 6,11-13 . That is, for arbitrary p , heterodyne detection outputs the amplitude α ∈ C m with a probability P(α) = Tr 1 π m |α α| p |α α| =: where the Husimi Q function Q p of p has been introduced. The Q function is such that, for any p , it holds that 0 ≤ Q p (α, α * ) ≤ 1 π m for all α ∈ C m , so that it yields a well-defined classical phase-space probability distribution. It provides an experimentally convenient tool for the evaluation of expectation values of anti-normal ordered observables, i.e., observables expressed with all annihilation operators to the left of all creation operators [11][12][13] .
We emphasise that ·Γ are the same observable. The former is just an alternative expression (the anti-normal ordered one, in terms of annihilation and creation operators) of the latter (in terms of quadrature operators). Hence, in terms of ·Γ Note that for any fixed n, the anti-normal ordered expression (104) can be efficiently computed for all j ∈ [n]. Also note that normal and anti-normal ordered expressions of all powers of the number operator are explicitly given in ref. 15 . Then, by virtue of the Optical Equivalence Theorem 11-13 , the expectation value of ·Γ (j) k1,l1,...,kj ,lj · (â,â † ) with respect to any preparation p equals the expected value E Q p ·Γ (j) k1,l1,...,kj ,lj · of ·Γ (j) k1,l1,...,kj ,lj · with respect to Q p . Here, Q p is understood as a classical phase-space probability distribution and ·Γ (j) k1,l1,...,kj ,lj · as a classical random variable that takes real values ·Γ (j) k1,l1,...,kj ,lj · (α, α * ). ·Γ (j) k1,l1,...,kj ,lj · (α, α * ) is defined as ·Γ (j) k1,l1,...,kj ,lj · (â,â † ) but replacing the operator vectorâ by the complex vector α ∈ C m . That is, it holds that where the short-hand notation d 2m α := with α κ sampled from Q p . The integer C ≤2(n+1) is the number of copies of the preparation devoted to the moment Γ (j) k1,l1,...,kj ,lj , as indicated in Supplementary Box 2, and is given explicitly in supplementary equation (85). The rest of the argument follows as in the Gaussian case. Note, finally, that the definition of the maximal variances σ ≤2(n+1) remains the same too.

Supplementary Note 5 -Stability against systematic errors.
Apart from statistical errors, Arthur's measurement procedure could also have systematic errors. That is, if the characterisation of his single-mode measurement channels is erroneous, he could actually be measuring different observables from the ones he thinks he does. Theorems 2 and 3, as well as the Supplementary Corollary 1, consider only statistical errors, i.e., those that can be decreased by increasing the number of measurement repetitions (and, hence, the number of copies of p ). Since systematic errors cannot be decreased by accumulating statistics, no certification method based exclusively on the measurement statistics can rule them out. However, the stability analyses of Supplementary Lemmas 5, 8, and 11 hold regardless of the nature of errors. Thus, the experimental estimates F (0) , F (n) , and F (n) S (and, therefore, also the certification tests) turn out to be robust also against small systematic errors: The total fidelity deviations due to systematic errors scales linearly with the magnitude of the largest systematic error and polynomially in all the other relevant parameters as given in supplementary equations (52), (77), and (91).
Still, it is illustrative to consider a physically relevant example. A typical systematic error is non-unit quantum efficiency of the detectors used for homodyning. In that case, the probability density functionP of measurement outcomes r of a quadraturê r equals the ideal one P convolutioned with the normal distribution N of mean zero and squared variance (1 − η)/4η, where η is the quantum efficiency of the detectors 9 . That is,P(r) = (P * N )(r) := dr P(r )N (r − r ). Using that the first and second non-central moments of N satisfy r N := drrN (r − r ) = r (108a) and r 2 N := drr 2 N (r − r ) = r 2 + respectively, one obtains that r P = r P (109a) and r 2 P = r 2 P + 1 − η 4η . (109b) That is, the expectation value ofr is not affected by this type of systematic errors and that ofr 2 deviates from the ideal one by (1 − η)/(4η). Furthermore, the expectation values of products of quadrature operators acting on different modes are also not affected, as this type of systematic error acts independently on different modes. In the absence of statistical errors, this leads to an error vector = 0 and an error matrix E (1) that is diagonal and such that E (1) max ≤ (1 − η)/(4η), so that E (1) 1 ≤ m(1 − η)/(2η). Inserting this into supplementary equation (57), we see for instance that, for Gaussian targets, the contribution to the deviation of the fidelity estimate due to non-ideal detector efficiency in the homodyne detectors is smaller than s 2 max m 1−η 2η . This, in turn, is smaller or equal than a desired constant maximal error ε if where the approximation holds whenever s 2 max m 2ε. The scaling given by the bound in supplementary equation (110) is experimentally convenient in that, in particular, it implies that the detector inefficiency 1 − η needs to decrease only inversely proportional with the number of modes m.
Another typical systematic error is the limited power of the local oscillator field used for the homodyne detection: The homodyne (photocurrent difference) statistics, i.e., the distribution of homodyne measurement outcomes, match exactly the statistics of the corresponding quadrature only in the limit of an intense local-oscillator beam 10 . The most obvious difference is that the homodyne statistics is discrete whereas the quadrature statistics is continuous, with the former approximating the latter increasingly better as the local-oscillator power increases. However, we emphasise that our method relies on the estimation of only the expectation values of quadratures and not their full statistics. It can be seen that, provided that the local oscillator is in a coherent state, the effect of limited power is just to increase the variance of the effective quadrature without changing its expectation value with respect to the ideal case. Furthermore, in the multi-mode scenario, if the different modes are homodyned with independent local oscillators, the latter is also true for products of quadratures, as the ones considered in this work. Therefore, the effect of systematic errors due to limited homodyne local-oscillator power in our fidelity estimates is expected not to be critical either.

Supplementary Note 6 -Auxiliary mathematical relations.
Derivation of the properties of the operator valued Pochhammer-Symbol. We begin with equation (21a). The general relationship (a † j ) tn j (a j ) t = p t (n j ), for t ∈ N, can be shown by induction starting from p 0 (n j ) =n j and noting that, for all t ≥ −1, a † j p t (n)a j = a † jn j (n j − 1)(n j − 2) · · · (n j − t)a j (112) = a † jn j (n j − 1)(n j − 2) · · · (n j − (t − 1))a j (n j − (t + 1)) (113) = p t (n j )(n j − (t + 1)) (114) = p t+1 (n j ), as can be verified using the commutation relations between a j and a † j . Setting t = n j gives equation (21a). In turn, equation (21b) can be shown by noting that (a † j ) nj (a j ) nj = (a † j ) nj −1n j (a j ) nj −1 and applying supplementary equation (111), for t = n j − 1, to the right-hand side of supplementary equation (116).
Proof of supplementary equation (46). Note that for x = 0 both sides of supplementary equation (46) yield 1 and hence the bound holds in that case. We make the substitution y = 1/x and show that the bound in supplementary equation (46) holds for all x > 0 by proving the following: But this is equivalent to 2y 2 + 3y + 2 ≤ e y (2 + y).
A straightforward calculation shows that both sides and also the first derivatives of both sides coincide at y = 0, while the second derivative of the right hand side is always larger than the second derivative of the left hand side. This proves supplementary equation (117) and hence finishes the proof of the bound in supplementary equation (46).