Efficient estimation of multipartite quantum coherence

Quantification of coherence lies at the heart of quantum information processing and fundamental physics. Exact evaluation of coherence measures generally needs a full reconstruction of the density matrix, which becomes intractable for large-scale multipartite systems. Here, we propose a systematic theoretical approach to efficiently estimating lower and upper bounds of coherence in multipartite states. Under the stabilizer formalism, the lower bound is determined by the spectrum estimation method with a small number of measurements and the upper bound is determined by a single measurement. We verify our theory with a four-qubit optical quantum system.We experimentally implement various multi-qubit entangled states, including the Greenberger-Horne-Zeilinger state, the cluster state, and the W state, and show how their coherence are efficiently inferred from measuring few observables.

Several experiments have been reported regarding the efficient detection of robustness of coherence [34,35]. However, the experimental detection of general coherence measure, such as the relative entropy of coherence [14], is still missing. Theoretical proposals to estimate general multipartite coherence without costly state tomography have also been proposed [36][37][38]. While the initial proposals either need copies of the prepared multipartite state [36] or complicated post-processing [37], the spectrum estimation method was recently proposed [38], which only requires local measurements and easy-to-compute post-processing. Nevertheless, the performance of the spectrum estimation method highly depends on the choice of the measurements, and how it works for a general multipartite state still needs further study. Moreover, * These authors contributed equally to this work † luhe@sdu.edu.cn existing works generally focus on the lower bound of coherence. For a given quantum state, the maximal entanglement or discord it can generates are upper bounded by its coherence [15,32], which makes detecting the upper bound of coherence important as it can indicates whether the given quantum state can generate sufficient resource for a certain QIP task.
In this work, we theoretically address these issues by proposing two methods that can respectively detect the lower and upper bound of coherence for all multi-qubit stabilizer states. The lower bound detection is based on the spectrum estimation method [38] and the stabilizer theory [39,40], which only requires few local observable measurements for stabilizer states. The upper bound detection is based on the monogamy of coherence with a single local measurement. Experimentally, we prepare five stabilizer states of up to four qubits and demonstrate how few number of measurements could enable us to infer multipartite coherence.

A. Lower bound estimation
Under the computational basis {|i : i ∈ {0, 1} ⊗n } of an nqubit state, we consider the relative entropy of coherence [14] with S VN = −tr[ρ log 2 ρ] being the von Neumann entropy and ρ d = i i|ρ|i |i i| being the diagonal part of ρ. The relative entropy of coherence characterizes the asymptotic distillable coherence under different types of incoherent operations [41,42], quantifies the genuine randomness that can be extracted from measuring the quantum state in the computational basis [43][44][45], captures the deviation from thermodynamic equilibrium [46], etc. We thus focus on the estimation, in particular, the lower and upper bounds, of the relative entropy of coherence for general multipartite states. The lower bound l c (ρ) of the coherence C RE (ρ) can be obtained by spectrum estimation and the majorization theory [47] as where d = (d 1 , ..., d 2 n ) are the diagonal elements of ρ, p = (p 1 , ..., p 2 n ) is the estimated probability distribution of the measurement on a certain entangled basis {|ψ k } 2 n k=1 , ∨ is majorization joint, and ∧ p∈X p is the majorization meet of all probability distributions in X [38]. Here the majorization join and meet are defined based on majorization. Specifically, given two probability distributions a = (a 1 , a 2 , . . . , a n ) and b = (b 1 , b 2 , . . . , b n ) with a 1 ≥ a 2 ≥ . . . ≥ a n and A probability distribution c is called the majorization join (meet) of a and b if it satisfies: (i) c a, b (c ≺ a, b), and (ii) c ≺c (c c) for anyc that satisfies a, b ≺c (a, b c) [47]. Here, we consider p is selected from the set X, which satisfies X = {p|Ap ≥ α, Bp = β}. A and B are matrices and α and β are vectors. "≥" represents component-wise comparison.
To calculate l c (ρ) via Eq. 2, it is crucial to set constraints Ap ≥ α and Bp = β with experimentally collected data, then find the "largest" distribution p majorized by all probability distributions in X, i.e., ∧ p∈X p. According to the hermiticity of density matrix ρ, we can set A as {2 n -dimensional identity matrix and α = 0 for Ap ≥ α, by which p k ≥ 0 is guaranteed. {In the following section, we introduce the procedure to construct constraint Bp = β via stabilizer formalism.

B. Constructing constraint via stabilizer formalism
An observable S i stabilizes an n-qubit state |ψ if |ψ is an eigenstate with eigenvalue +1 of S i , i.e., S i |ψ = |ψ . The set S of operators S i is the stabilizer of |ψ , and |ψ is a so-called a stabilizer state [40,48]. For an n-qubit state, there are n stabilizing operators {S 1 , ..., S n } that can uniquely determine |ψ . Here S 1 , ..., S n are the generators of the set S, and we denote S = S 1 , ..., S n . Note that |ψ is not only stabilized by {S 1 , ..., S n }, but also their products. Thus, there could be in total 2 n stabilizer operators in S.
Given an n-qubit stabilizer state |ψ 1 associated with stabilizer S, there exists an orthonormal basis {|ψ k } 2 n k=1 including |ψ 1 = |ψ , where |ψ k is uniquely specified by S but with different eigenvalues, i.e., S i |ψ k = a ik |ψ k with eigenvalues a ik = ±1. For example, the Bell state |Φ + = (|00 + |11 )/ √ 2 can be specified by S = X (1) X (2) , Z (1) Z (2) . Hereafter, X (j) , Y (j) and Z (j) denote the Pauli matrices σ x , σ y and σ z acting on the j-th qubit. In the following, qubit index j may be omitted if there is no confusion. The four Bell states |Φ ± and |Ψ ± are specified by S = XX, ZZ as well, with eigenvalues of (+1, +1), (−1, +1), (+1, −1) and (−1, −1) respectively. The basis {|ψ k } 2 n k=1 associated with the same stabilizer S is also called graphdiagonal basis [49,50] as stabilizer state is equivalent to a graph state under local Clifford operations [51]. In the graph-diagonal basis, the stabilizing operator can be written as S i = k a ik |ψ k ψ k |, and its expected value on a given quantum state ρ is S i = tr(S i ρ) = k p k a ik , where the parameters p k = ψ k |ρ|ψ k form a probability distribution p = (p 1 , ..., p 2 n ) with p k ≥ 0 and k p k = 1.
The expected values of stabilizer S on ρ leads to 2 n equations, and can be represented in matrix form    a 11 . . . a 12 n . . . . . . . . .
from which we can construct the constraint Bp = β. However, in practice, there does not always exist solutions of ∧ p∈X p in Eq. 2 with constraint Eq. 3 (as reflected by our experimental results). The experimentally generated state always has a distance to target state due to the inevitable imperfections, which might lead to no solution of ∧ p∈X p with inputs of { S i }. On the other hand, experimentally obtained S i is always associated with statistical errors. We address these issues by introducing the experimental standard deviation σ i of S i and relax the constraint Eq. 3 to an inequality form of where wσ i with w ≥ 0 is the deviation to the mean value S i represented in σ i . To this end, an experimentally accessible constraint is formulated as β − ≤ Bp ≤ β + . In practice, instead of measuring all the stabilizers, which is impractical for a large quantum state, we can select a small subset of stabilizers so that the number of measurement does not scale exponentially to the number of qubits. Note that I n = 1 must be set in Eq. 4 to ensure k p k = 1. If we apply the scheme on the graph-diagonal states ρ = k λ k |ψ k ψ k | with λ = (λ 1 , ..., λ 2 n ) the spectrum of ρ, we have p = λ. Thus, d ∨ (∧ p∈X p) = λ implies l c (ρ) = C RE (ρ), which indicates that the estimated lower bound of coherence is tight for graph-diagonal states.
We emphasize that relaxing the constraint to β − ≤ Bp ≤ β + does not increasing the risk of overestimation of l c (ρ). Suppose that X 1 and X 2 are two feasible sets of probability distributions, and satisfy X 1 ⊆ X 2 . X 1 and X 2 are restricted to d ≺ ∧ p∈X1 p and d ≺ ∧ p∈X2 p, otherwise the result of Eq. 2 is 0. According to the definition of majorization meet, ∧ p∈X p is the "largest" distribution majorized by all probability distributions in X. Therefore, ∧ p∈X p becomes "smaller" when we enlarger the range of X, i.e, ∧ p∈X1 p ∧ p∈X2 p for X 1 ⊆ X 2 . According to the strict Schur concavity of the Shannon entropy S, we obtain S(∧ p∈X1 p) < S(∧ p∈X2 p), which implies S(d) − S(∧ p∈X1 p) > S(d) − S(∧ p∈X2 p). Thus, we conclude that l c (ρ) decreases with enlarging the range of X.
Similar constraints can also be formulated for multi-qubit states that do not obviously fit the stabilizer formalism. The stabilizing operators of such kind of n-qubit states |ψ could be determined by finding its unitary dynamics U ψ acting on |0 ⊗n , i.e., |ψ = U ψ |0 ⊗n [48]. As |0 ⊗n is stabilized by S |0 ⊗n i = Z (i) , ∀i ∈ {1, 2, ..., n}, the stabilizing operator of C. Upper bound estimation Theorem 1. Let M(d) be a set of states with the same diag- We first consider the case that there are no non-zero ele- We define a completely positive and trace-preserving (CPTP) According to the definition of strictly incoherent operation [52,53], K α is a strictly incoherent operator and Λ 1 (·) is a strictly incoherent operation. For the case that there are m (m < 2 n ) non-zero elements in d, we first transform |ψ d ψ d | into a block diagonal matrix via via a permutation matrix M , i.e., acting on m-dimensional Hilbert space. Thus, we can derive that Note that any permutation matrix is strictly incoherent unitary [54]. For any two strictly incoherent operations Λ α and Λ β , the operation Λ = Λ α • Λ β is also a strictly incoherent operation [55]. Thus, To upper bound the coherence of an n-qubit state ρ, we measure it in the computational basis {|i }, which yields a distribution d = (d 1 , ..., d 2 n ). Based on the monogamy of the relative entropy of coherence [14], the coherence of ρ is upper bounded by the coherence u c (ρ) of state |ψ d = Note that this lower bound is tight for pure states, and is capable for various coherence measures.
For each experimentally generated state ρ ψ expt , we measure the expected values of its stabilizers S i associated with the corresponding statistical errors σ i . We refer to Appendix A for details of stabilizing operators of |ψ and the corresponding graph-diagonal basis. The measured expected values of stabilizing operators are presented in Appendix B. With the measured S i and σ i , we construct the constraint Eq. 4. Thus, we can calculate l c ω,m (ρ ψ expt ) via solving Eq. 2, where ω represents the setting of deviations and m is the number of S i we construct the constraint. We set ω as non-negative integers from 0 to 3, and m from 1 to 2 n − 1 as I n = 1 must be set in the constraint. In our calculation, we treat the case of no solution as l stabilizers for other three states. With the increasing of m, the maximum l c ω,m (ρ ψ expt ) increases accordingly. When m gets close to 2 n − 1, there exists situations that we can not obtain l c ω,m (ρ ψ expt ) for smaller w (the values are 0 at the right lower corner of each figure in Fig. 2). As aforementioned, this is caused by the experimental imperfections, such as slight misalignment of optical elements during data collection, which introduces small variation of prepared ρ ψ expt . The issue is improved by extending the range of S ψ i , i.e., increasing ω. As shown in Fig. 2b-Fig. 2e, most l c ω,m (ρ ψ expt ) has valid solutions for large m by setting ω = 3. Moreover, the accuracy of estimated l c ω,m (ρ ψ expt ) increases along with m as well. We investigate this by calculating the normalized distance between l c ω,m (ρ ψ expt ) and C RE (ρ ψ expt ), i.e., 1 − l c ω,m (ρ ψ expt )/C RE (ρ ψ expt ). C RE (ρ ψ expt ) is calculated with Eq. 1 by reconstructing ρ ψ expt via quantum state tomographic technology (See Appendix B for the reconstructed ρ ψ expt ). The distances between maximal l c 3,m (ρ ψ expt ) and C RE (ρ ψ expt ) are shown in Fig. 3a, from which we observe that the distance drops down quickly with increasing of m and tends to converge at m = 5.
As aforementioned, the choice of selecting m stabilizers from {S ψ i } is not unique except m = 2 n − 1. An important property is the successful probability of obtaining valid l c ω,m by randomly selecting m stabilizers. We show the percentage of valid l c ω,5 (ρ ψ expt ) (l c ω,5 (ρ ψ expt ) > 0) for m = 5 in Fig. 3b. By increasing ω, the probability of getting valid l c ω,5 (ρ ψ expt ) is enhanced, especially for ρ W3 expt . However, one cannot increasing ω arbitrarily. A larger ω represents a smaller probability we could obtain S i in range , which is less than 0.3% for ω = 3. This is also reflected by normalized distance of maximal l c ω,5 (ρ ψ expt ) shown in Fig. 3c, which indicates the inaccuracy of l c ω,5 (ρ ψ expt ) increases when we extend the range of S i . Also, the results in Fig. 3c agree with our claim that relaxing the constraint decreases the estimated value of l c (ρ). We conclude that set ω = 3 is reasonable in experiment, under which we observe the probability of getting valid l c ω,5 (ρ ψ expt ) is 100% for ρ GHZ3 expt , ρ W3 expt , ρ C4 expt and ρ W4 expt , and that of 86% for ρ GHZ4 expt . Note that the estimated l c (ρ) of ρ GHZ3 expt , ρ W3 expt , ρ C4 expt and ρ GHZ4 expt is slightly more accurate than that of ρ W4 expt as shown in Fig. 3a and Fig. 3c. The main reason is that the fidelity of prepared ρ W4 expt is slightly lower than that of other four states as shown in Appendix B. The lower fidelity implies a larger distance between prepared state and target state, which indicates λ d ∨ (∧ p∈X p) so that C RE (ρ) = S(d) − S(λ) > l c (ρ) = S(d) − S(d ∨ (∧ p∈X p)). However, it does not indicate the lower fidelity always leads to the bigger gap between C RE (ρ) and l c (ρ). If the prepared state is still a graph-diagonal state after considering the experimental imperfections, the estimated lower bound on such state is tight as well, i.e., l c (ρ) = C RE (ρ). In Appendix C, we analyze our experimental imperfections and show how it affect the prepared states.
Finally, we estimate the upper bound u c of ρ ψ expt by measuring the probability distribution d ψ expt on basis of Z ⊗n . The probability distribution d ψ expt are shown in Appendix B, by which we can calculate u c (ρ ψ expt ) according to Eq. 8. The results of u c (ρ ψ expt ) are shown with blue dash line in Fig. 3d. For all the five states, we observe that C RE (ρ ψ expt ) lies within the range bounded by l c 3,5 (ρ ψ expt ) and u c (ρ ψ expt ).

IV. CONCLUSION
To conclude, we introduce an efficient and experimentally friendly estimation method for detecting coherence of multipartite states. We demonstrate that the coherence with high accuracy as well as high successful probability can be efficiently estimated with a few measurements for various multiqubit states. The procedure to obtain the lower bound is based on few measurements of the stabilizing operators, similar to multipartite entanglement detection [60] and multipartite Bell inequalities [61]. It thus indicates that coherence and other resources can be inferred by the same set of measurements, which can further benefit our understanding of the connection between coherence and other resources [15,30,31,33].
There are several open follow-up problems. The scheme to detect the lower bound of multipartite coherence is efficient and tight for graph-diagonal states as well as special types of quantum states that admit efficient classical representation (such as the W state). Whether our scheme will work for more general multi-qubit quantum state is an interesting future work. Besides, one may concern whether our scheme can be generalized to high-dimensional cases. For the multiqudit stabilizer state, its stabilizing operators constructed by generalized Pauli group are generally non-Hermitian [62][63][64], which can not be observed directly in experiment. Recent studies indicate that expected values of non-Hermitian operators could be measured via weak measurements [65,66]. However, the definite answer may require rather sophisticated analysis. Finally, it is worth noting that a fidelity-based method was recently proposed to detect the lower bound of multipartite coherence via the convex roof construction [67].
Thus, another open follow-up question is whether our scheme can be further generalized to other coherence measures while maintaining their appealing feature of efficiency and simplicity.

ACKNOWLEDGMENTS
We are grateful to anonymous referee for providing very useful comments on an earlier version of this manuscript. This work is supported by the National Natural Science Foundation The GHZ 3 -diagonal basis can be obtained from the computational basis |k 1 k 2 k 3 by acting a 3-qubit unitary operation 000 ).
operators are The GHZ 4 -diagonal basis can be obtained from the computational basis |k 1 k 2 k 3 k 4 by acting a 4-qubit unitary operation is the unitary matrix of controlled-Z (CZ) operation. Sixteen GHZ 4 -diagonal bases are determined by |ψ GHZ4 i = U GHZ4 |k 1 k 2 k 3 k 4 and shown in Table II.
IZXX and IIZZ. Other stabilizing operators are The Cluster-diagonal basis can be obtained from the computational basis |k 1 k 2 k 3 k 4 by acting a 4-qubit unitary operation . Then, the Cluster-diagonal bases are determined by |ψ C4 i = U C4 |k 1 k 2 k 3 k 4 and shown in Table III. 3. W state |W 3 can be transformed from |000 by unitary U W3 = (XZI + IXZ + ZIX)/ √ 3, i.e., |W 3 = U W3 |000 [60]. Thus the generators of |W 3 is derived by Other stabilizing operators are The W 3 -diagonal bases are shown in Table IV. Similarly, We can find a possible unitary operator U W4 = (ZZZX + ZZXI + ZXII + XIII)/2 to generate |W 4 from |0000 . Thus, the generators of |W 4 can be obtained, i.e., and other stabilizing operators are The W 4 -diagonal bases are shown in Table V.

Appendix B: Details of experimental realizations and results
In this Appendix, we provide further details about our experimental setup. It would be useful to bear in mind the following:   As shown in Fig. 4, the power of the pump light can be adjusted by a HWP and a PBS. After the PBS, horizontal polarization |H p is rotated to |+ p = 1 √ 2 (|H p + |V p ) by a HWP set at 22.5 • . Pump beam is focused into PPKTP crystal with beam waist of 74µm by two lenses L 1 , whose focal length is 75mm and 125mm respectively. PPKTP crystal, with dimensions of 10 mm (length) × 2 mm (width) × 1 mm (thickness) and poling period of Λ = 10.025µm, is held in a home-built copper oven and the temperature is controlled by a homemade temperature controller, which is set at 29 • C to realize the optimum type-II phase matching at 810 nm. Then, the pump beam is split by a dual-wavelength PBS, and coherently pumps PPKTP in the clockwise and counterclockwise direction respectively. The clockwise and counterclockwise photons are recombined at the dual-wavelength PBS to generate polarization-entangled photons with ideal form of . Two photons are filtered by narrow band filter (NBF) with full width at half maximum (FWHM) of 3nm, and coupled into single-mode fiber by lenses of focal length 200 mm (L 2 and L 3 ) and objective lenses (not shown in Fig. 4). HWP@90°v on path h&v of photon b 1 2 (|H 1 |h 2 |V 3 |h 4 + |H 1 |h 2 |H 3 |v 4 + |V 1 |h 2 |H 3 |h 4 + |H 1 |v 2 |H 3 |h 4 ) = |W 4 (B3)

Measurement and quantum state tomography
If the photon is encoded either in polarization DOF or path DOF, we first perform measurement on polarization DOF and then the path DOF. As illustrated in Fig. 6, the measurement basis of qubit on polarization DOF α|H + β|V is determined by HWP at θ 1 and QWP at θ 2 . The measurement basis of qubit on path DOF γ|h + δ|v is determined by HWP at θ 3 and QWP at θ 4 . Finally, a PBS is applied before the photon arrives detector. With this setting, the measurement on basis (α|H + β|V ) ⊗ (γ|h + δ|v ) is achieved. The specific calculations are shown in Eq. (B4). If the photon is only encoded in polarization DOF, the measurement on basis α|H + β|V is implemented by a QWP, HWP and PBS. With this experimental setting, we can perform measurement on arbitrary basis. The experimental results of S ψ i and d ψ expt are shown in Fig. 7a and b, respectively. Moreover, we reconstruct experimentally generated states ρ GHZ4 expt , ρ C4 expt , ρ W4 expt , ρ GHZ3 expt and ρ W3 expt by quantum state tomography [48]. The results are shown in Fig. 8, from which we calculate the fidelity F ψ = tr(ρ ψ expt |ψ ψ|). We observe that F GHZ3 = 0.9643 ± 0.0003, F W3 = 0.9589 ± 0.0005, F GHZ4 = 0.9571 ± 0.0003, F C4 = 0.9497 ± 0.0002 and F W4 = 0.915 ± 0.001. Also, the relative entropy of coherence C RE (ρ ψ expt ) of ρ ψ expt can be calculated by C RE (ρ) = S VN (ρ d ) − S VN (ρ), according to which we obtain C RE (ρ GHZ3 expt ) = 0.875 ± 0.002, C RE (ρ W3 expt ) = 1.479 ± 0.004, C RE (ρ GHZ4 expt ) = 0.906 ± 0.002, C RE (ρ C4 expt ) = 1.806 ± 0.002 and C RE (ρ W4 expt ) = 1.964 ± 0.003. channel, which causes that there is only 5% difference in state fidelity (shown in Appendix B) but 15% difference in normalized distance (shown in Fig. 3c).