Optimal Distributed quantum sensing using Gaussian states

We find and investigate the optimal scheme of quantum distributed Gaussian sensing for estimation of the average of independent phase shifts. We show that the ultimate sensitivity is achievable by using an entangled symmetric Gaussian state, which can be generated using a single-mode squeezed vacuum state, a beam-splitter network, and homodyne detection on each output mode in the absence of photon loss. Interestingly, the maximal entanglement of a symmetric Gaussian state is not optimal although the presence of entanglement is advantageous as compared to the case using a product symmetric Gaussian state. It is also demonstrated that when loss occurs, homodyne detection and other types of Gaussian measurements compete for better sensitivity, depending on the amount of loss and properties of a probe state. None of them provide the ultimate sensitivity, indicating that non-Gaussian measurements are required for optimality in lossy cases. Our general results obtained through a full-analytical investigation will offer important perspectives to the future theoretical and experimental study for quantum distributed Gaussian sensing.


I. INTRODUCTION
Quantum resources are known to be useful for further enhancing the precision and the sensitivity of estimation of various physical quantities beyond the standard quantum limit [1][2][3][4][5]. A number of studies on single-parameter estimation have been performed over the past few decades [6], but much attention has begun to be paid to estimation of multiparameters in recent years [7]. Quantum-enhanced sensitivity in simultaneous estimation of multiple phases has been investigated to explain the role of quantum entanglement and identify optimal and realistic setups saturating the ultimate theoretical sensitivity [8][9][10][11][12]. The advantage of exploiting quantum entanglement becomes more significant when sensing takes place in different locations and the parameter of interest is a global feature of the network, e.g., the average of distributed independent phases [13][14][15][16][17][18]. Such distributed sensing is related to applications such as global clock synchronization [19] and phase imaging [8]. These inspire the use of more practical quantum resources that are feasible in a well-controlled manner with current technology, e.g., Gaussian systems [20]. Very recently, the sensitivities of distributed quantum sensing with Gaussian states were studied under specific conditions [16,17]. The ultimate sensitivity and feasible optimal schemes, however, are not yet found and studied in the class of Gaussian metrology [21][22][23].
In this paper, we investigate the ultimate sensitivity for the average phase estimation in distributed quantum sensing * h.jeong37@gmail.com with Gaussian states, where the phases are encoded onto a multimode Gaussian probe state, as described in Fig. 1. We find an optimal probe state and measurement setup that achieve the ultimate sensitivity, which are shown to be experimentally feasible with current technology. Interestingly, we demonstrate that the optimal symmetric Gaussian probe state is not a maximally entangled state. For practical relevance, we further analyze the effect of loss, the entanglement-enhanced gain, and other Gaussian measurements in various conditions.
We begin with a brief introduction to the formalism describing Gaussian states and multiparameter estimation. Gaussian states are defined as states whose Wigner functions are Gaussian distributions and thus characterized by the first moment vector d i = Tr[ρQ i ] and the covariance matrix Here, a quadrature operator vector of a M-mode continuous variable quantum system is defined asQ = (x 1 ,p 1 , ...,x M ,p M ) T , satisfying the canonical commutation relation,

II. DISTRIBUTED SENSING
Consider estimation of M-parameter φ = (φ 1 , φ 2 , ..., φ M ) T based on measurement outcomes x, obtained with a conditional probability p(x|φ). The multiparameter Cramér-Rao inequality states that the M × M estimation error matrix The conditional probability p(x|φ) = Tr[ρ φˆ x ] is given by a positive operator-valued measureˆ x for a given parameterencoded stateρ φ . The quantum Cramér-Rao inequality sets a lower bound for the error of an unbiased estimator, i.e., Schematic of distributed sensing under investigation. A multimode probe stateρ probe generated from the first beam splitter network (BSN) for a given product state input ⊗ M i=1ρ i undergoes the individual phase shifts on each mode. The parameter-imprinted stateρ φ is fed into the second BSN, followed by measurement. The measurement outcomes are used in post-processing to estimate the Fisher information matrix (QFIM), withL i being a symmetric logarithmic derivative operator associated with ith parameter φ i [25]. When a linear combination of φ i 's, φ * = w T φ = M i=1 w i φ i with the weight vector w, is of particular interest, the estimation error is bounded as [26] Here, F −1 and H −1 are understood as the inverse on their support if the matrices are singular. Throughout this paper, we assume the normalization M i=1 |w i | = 1 for simplicity.

A. Quantum Fisher information matrix
Consider a distributed phase sensor in which a product Gaussian input state ⊗ M i=1ρ i is injected into a beam splitter network (BSN), preparing a probe stateρ probe , the multiphase information is encoded ontoρ probe by a unitary operation jâ j ), and the output stateρ φ is measured after the second BSN, as depicted in Fig. 1. Note that configuration of the first BSN enables one to generate any probe Gaussian states [22,27]. We also implicitly assume a strong reference beam to define the phases, accessible in each mode for measurement [28]. Here, we aim to investigate the sensitivity of Gaussian states for estimation of the parameter φ * . When the probe stateρ probe after the first BSN is a pure Gaussian state characterized by ( , d ), the elements of the QFIM are written as [29][30][31][32][33][34] where A (i, j) denotes the 2 × 2 submatrix in the ith row and jth column of the M × M block matrix A, and similar for the vector d (i) . The derivation of the QFIM of Eq. (2) is provided in Appendix A. The convexity of QFIM makes it sufficient to consider only pure probe states to find an optimal state maximizing the QFIM [14], but one can find the analytical form of the QFIM for general Gaussian states [29][30][31][32][33][34]. The quantum Cramér-Rao bound in Eq. (1) can be saturated since the generators of parameters commute [12].
FIG. 2. Estimation errors when probing the phases with product states. The top curve represents the standard quantum limit, 2 φ * SQL , whereas the other curves show the error 2 φ * OPGS when M = 1, 3, 100. The error 2 φ * OPGS increases with M for a fixedN, but is always below the error 2 φ * SQL , which approaches 1/8N as M → ∞.

B. Optimal product Gaussian state
Let us first consider the case where the probe state is a product state and thus the QFIM is evidently a diagonal block matrix. Without loss of generality, we assume that the block matrix of the covariance matrix for ith mode is (i,i) = diag(e 2r i , e −2r i )/2, simplifying the estimation error of φ * to be 2 When probing with a product coherent state, the error bound becomes , and the best strategy for a given total average photon numberN is to distribute the energyN over the modes according to the weight |w i |, i.e., The estimation error is thus where the lower bound defines the standard quantum limit. When w i = 1/M, i.e., φ * is the average phase, 2 φ * SQL = 1/4Mn, wheren ≡N/M represents an equal average number of photons hitting each phase shifter.
Among all product Gaussian states, the best strategy under the energy constraintN is to prepare the probe state in a product squeezed vacuum state with 8N 2 i . Thus, particularly when w i = 1/M, in which φ * is the average phase, the estimation error becomes where we have set r i = r for all i andN = M sinh 2 r. Note that the Heisenberg scaling withn orN is achieved. We refer to the above product squeezed vacuum state as the optimal product Gaussian state (OPGS) throughout this paper. The error 2 φ * OPGS grows with the number of modes M, over which the probe state is distributed for a givenN, as shown in Fig. 2. When an equal energy can be used in all the modes, i.e., for a fixedn, the error 2 φ * OPGS decreases with M, which is obvious since the total energy being used increases by M times. It can be easily shown that the estimation error 2 φ * OPGS can be achieved by performing homodyne detection on each mode without the second BSN [35].

C. Optimal entangled Gaussian state
We now turn to the case when the first BSN is configured to create mode correlation for an injected product input state. In order to find the ultimate sensitivity in distributed sensing using Gaussian states and an optimal probe state, one can further develop the inequality of Eq. (1) as . From now on, let us focus on the estimation of the average phase, i.e., w i = 1/M. Using a series of inequalities, we show that the error for the average phase estimation is given by (see Appendix B for the detail) We note that the ultimate error 2 φ * OEGS scales withN −2 orn −2 , and is smaller than the error 2 φ * OPGS . A similar scaling has been discussed in Ref. [14] but with different quantification of the resource.
We show that the ultimate error 2 φ * OEGS can be achieved by using symmetric Gaussian probe states with zero displacement. The covariance matrix of pure symmetric Gaussian probe states can be written as a M × M partitioned matrix probe with submatrices (i,i) probe = diag(γ 1 , γ 2 ) for all i and [30,[36][37][38]. Since the states are assumed to be pure, the components obey the relations [36][37][38]. The QFIM for symmetric Gaussian states is evidently a symmetric matrix with H ii = H 11 for all i and H i j = H 12 for all i = j. Finally, using Eqs. (1) and (2), the estimation error is reduced to where H 11 = 2(γ 2 1 + γ 2 2 ) − 1 and H 12 = 2( 2 1 + 2 2 ) (see Appendix C for details). It is clear that the correlation quantified by 1 and 2 , H 12 , plays an important role, but the sensitivity is eventually determined by an interplay with the term H 11 that is not independent of H 12 for a given energy. After minimizing the lower bound in Eq. (4) under the energy constraintN = M(γ 1 + γ 2 − 1)/2 (see Appendix D for the detail), we recover the ultimate error 2 φ * OEGS when γ 1,2 = 1/2 + 1,2 and 1,2 = [N ± N (N + 1)]/M, leading to H 11 = 4N (2N + M + 1)/M 2 and H 12 = 4N (2N + 1)/M 2 . Therefore, the ultimate estimation error 2 φ * OEGS can be achieved by the optimal symmetric Gaussian state, which we call the optimal entangled Gaussian state (OEGS) throughout this paper. Most importantly, in contrast to the error 2 φ * OPGS , the error 2 φ * OEGS is independent of the number of modes M for a fixed energyN and scales with M −2 for a fixed n, evidently resulting from exploiting entanglement. Thus, the mode entanglement enables one to prevent the estimation error from growing with M. (b) However, the OPGS is optimal for the simultaneous phase estimation.

D. Role of entanglement
One might wonder whether the OEGS is the maximally entangled Gaussian state, for which the entropy of the reduced state is maximized. We now demonstrate that it is not the case. The entropy of the single-mode reduced state having a diagonal covariance matrix γ is given by S(γ ) =n T ln(1 + 1/n T ) + ln(n T + 1) [21], wheren T = √ γ 1 γ 2 − 1/2 is the average thermal photon number of the reduced single-mode state. The entropy S(γ ) increases with the entanglement of the total system under investigation, where pure symmetric Gaussian states are only considered [39]. Interestingly, the OEGS achieving the ultimate sensitivity does not have the maximal entropy, as shown in Fig. 3(a). This is surprising and in contrast to other cases, where maximally entangled states have shown to lead to the optimal sensitivity, e.g., the GHZ state of qubits exhibiting the maximal entropy of the reduced state [13]. In our scenario, the state often referred to as the continuous variable GHZ-type state having the maximal reduced entropy [40][41][42] exhibits worse sensitivity than the OEGS. A similar result has been reported for estimation of unitarily generated parameters in Ref. [10].
It is worth comparing with the error of simultaneous phase estimation, 2 where the bound 2 φ * OPGS can be achieved by the OPGS. When using the OEGS, however, the estimation error is given by It is clear that the error 2 φ OEGS is larger than the error 2 φ OPGS . More generally, any entangled symmetric Gaussian states exhibit worse sensitivity than the OPGS, as shown in Fig. 3(b).

A. Physical implementation of the optimal scheme
We have shown above that the ultimate estimation error 2 φ * OEGS is achieved by the OEGS. Generation of the latter is experimentally feasible with current technology as we provide here. Suppose that a product state of a p-squeezed vacuum and (M − 1) vacua is injected into the first BSN, configured asÛ BSN and θ j = arccos(M − j + 1) −1/2 . Consequently, one can show that the output state of the BSN is the OEGS. Notice that different configurations of BSN can be employed to generate the OEGS [16].
We demonstrate here that homodyne detection on each mode is sufficient to achieve the ultimate error 2 φ * OEGS without using the second BSN. The resultant probability distribution of homodyne detection follows a Gaussian distribution with the zero first-moment vector and the M × M covariance matrix HD with diagonal components [ HD ] ii = γ 1 cos 2 ϕ i + γ 2 sin 2 ϕ i and off-diagonal components [ HD ] i j = 1 cos ϕ i cos ϕ j + 2 sin ϕ i sin ϕ j , where ϕ i = φ i − θ HD,i with homodyne angles θ HD,i on ith mode. The error is thus given by 2 It can be easily shown that the lower bound is equal to 2 φ * OEGS when ϕ i = ϕ opt ≡ π/2 − cot −1 [2 N (N + 1)]/2 for all i. Such optimal phase setting can be made by adjusting the homodyne angles θ HD,i = φ i − ϕ opt .

B. Effects of loss
From a practical perspective, we analyze the effect of photon loss on the sensitivity. When loss is assumed to occur in each mode with an equal η, the covariance matrix of the probe state is transformed as probe → η probe + (1 − η)1 2M /2, i.e., γ 1,2 → ηγ 1,2 + (1 − η)/2 and 1,2 → η 1,2 [21,30]. Consequently, the theoretical optimal error bounds 2 φ * OPGS and respectively. When homodyne detection is performed, the resulting error bounds are, respectively, given as for which the homodyne angles have been appropriately chosen. One may also seek other type of Gaussian measurement that could outperform the case yielding 2 φ * OEGS,HD (η) in the presence of loss. We exemplify the latter by performing an appropriate general-dyne detection on the first output mode and heterodyne detection on the other output modes of the second BSN that is set to realizeÛ −1 BSN . The associated error FIG. 4. (a) Comparison among the estimation errors 2 φ * OEGS (η) (red curve), 2 φ * OEGS,HD (η) (green curve), and 2 φ * OEGS,GD (η) (orange curve) with loss η forN = 10. The estimation error of the optimal scheme 2 φ * OEGS (η) achieves an improvement by an order of magnitude compared to homodyne detection and general Gaussian detection when loss is significant. Note that there exist two crossing points ( ) between 2 φ * OEGS,HD (η) and 2 φ * OEGS,GD (η) as η increases. (b) The ratio of 2 φ * OEGS,HD (η) to 2 φ * OEGS,GD (η) as a function of N and η. The boundary, represented by a solid line, is given bȳ , at which homodyne detection and general Gaussian detection after the second BSN (Û −1 BSN ) yields the same sensitivity. General Gaussian detection scheme becomes significant only whenN > 2(1 + √ 2), and exhibits the most advantage over the homodyne detection at η = 0.5. bound when probing with the OEGS is given as whose derivation and detailed setup are provided in Appendix E. Figure 4(a) reveals that the error 2 φ * OEGS,HD (η) is competitive with 2 φ * OEGS,GD (η) depending on η, and none of them attain the ultimate error 2 φ * OEGS (η) when η < 1. Comparable behaviors between 2 φ * OEGS,HD (η) and 2 φ * OEGS,GD (η) are elaborated in terms ofN and η in Fig. 4(b), identifying the regions in which one prevails over the other. It shows that homodyne detection is advantageous whenN > (1 + √ 2)/2η(1 − η). Interestingly, the error bound 2 φ * OEGS,HD (η) is exactly the same as that of a singlemode phase estimation using a squeezed thermal state [43]. One could further reduce the error by having displacement as in Ref. [16], or seek for non-Gaussian measurements to achieve the ultimate error 2 φ * OEGS (η) in lossy cases [32,43]. The enhancement of sensitivity by entanglement can be quantified by the relative error ratio R opt = 2 φ * OPGS (η)/ 2 φ * OEGS (η) for the case that an optimal measurement is assumed, and the error ratio R HD = 2 φ * OPGS,HD (η)/ 2 φ * OEGS,HD (η) for the case that homodyne detection is performed. Figure 5(a) shows that the R opt slightly decreases with a moderate loss η and monotonically increases withn, while the R HD drastically drops with η and exhibits the optimum atn = 1/2 √ Mη(1 − η), where the relative enhancement is maximal, when η < 1. The behaviors of R opt and R HD with increasing M are presented in Fig. 5(b) for n = 6. Remarkably, both R opt and R HD are always greater than unity in all cases with any η, stressing the usefulness of entanglement in Gaussian distributed sensing against loss.

V. DISCUSSION
We have investigated the ultimate sensitivity of the average phase estimation in distributed quantum sensing using The quantum enhancement offered by the optimal scheme is robust against photon loss. Overall, loss is obviously always detrimental for a givenn and M, i.e., the error ratios decrease with η. Lines connecting dots in panel (b) are to guide the eyes.
Gaussian states. The ultimate sensitivity has been shown to be achievable by the OEGS possessing partial entanglement between the modes and by performing homodyne detection on each mode in the absence of loss. When photon loss occurs, homodyne detection ceases to be optimal, but non-Gaussian measurement would be required for achieving the ultimate sensitivity. Alternatively, a slightly better sensitivity can be obtained by conducting other type of Gaussian measurement on the output modes of the second BSN that implements the inverse transformation of the first BSN. Although the sensitivity decreases with loss in all the cases considered in this work, we have revealed that using the OEGS is always advantageous for average phase estimation as compared to the case using unentangled symmetric Gaussian states. While we have focused on identification of the ultimate sensitivity and the optimal setup for the average phase estimation in this work, finding those for estimation of other linear combinations of phases would also be an interesting future study. Another interesting open question is to explain the enhancement of sensitivity in a more intuitive manner such as using the multiparameter squeezing parameter [18] other than entanglement, which we leave for future study. Finally, finding the optimal measurement achieving the best sensitivity in the presence of photon-loss is also an important remaining task as in the recent study, where the optimal non-Gaussian measurement in a special case of M = 1 is theoretically found [43]. The experimental implementation of the optimal measurement needs to be devised.
It is worthwhile to discuss our results in relation to recent results in distributed sensing. First of all, a recent experiment successfully showed an enhancement by entanglement in distributed Gaussian quantum sensing [16]. The theory behind the experiment in Ref. [16] assumed that the phase shifts of interest were extremely small and the estimation error was quantified by the linear error propagation analysis from homodyne detection. However, our work identifies the ultimate estimation error in distributed Gaussian sensing by proposing the optimal Gaussian probe and it can be applied to phase shifts of arbitrary degrees. Thus, the experimental results could be understood better and interpreted from a broader perspective of distributed Gaussian sensing. In addition, the optimal entangled Gaussian state has been proven to be optimal for distributed quantum sensing of field-quadrature displacement [44][45][46]. In this section, we derive the quantum Fisher information matrix (QFIM) for distributed sensing using isothermal Gaussian states. When a phase-encoded state is the isothermal Gaussian quantum states characterized by [ (φ), d(φ)] with a isothermal photon numbern, the QFIM is given by [32]

ACKNOWLEDGMENTS
are the covariance matrix and the first moment vector of the quantum state after the unitary operation encoding φ corresponding to the symplectic matrix S(φ), respectively. In the distributed phase sensor, the symplectic transformation corresponds to Note that symplectic transformation S is defined as ones that preserve the canonical commutation relation, S T 2M S = 2M , corresponding to a Gaussian unitary operationÛ applied to density matrices by the relationÛ †QÛ = SQ.
The first term in Eq. (A1) can be simplified as = Tr[P i probe P j probe + P i probe 2M probe 2M P j where (i, j) probe = P i probe P j . Here, we have set φ = 0 without loss of generality since the QFIM is independent of φ under unitary transformation, and we have used which is the projection onto the ith mode, The second term in Eq. (A1) is Using the fact thatĜ * is invariant under any passive transformation, one can assume that â † iâ iâ † jâ j − â † iâ i â † jâ j = 0 for i = j without loss of generality and get which shows the variance of 2Ĝ * is the sum of the photon number variance in all the modes. Since a squeezed vacuum state exhibits the maximum photon number variance among Gaussian states, which is (cosh 4r − 1)/4, Under the constraint for the total mean photon number of the stateN, one can prove that the upper bound of 2Ĝ * in Eq. (B1) is given by 2N (N + 1), i.e.,

APPENDIX C: PROPERTIES OF THE QFIM FOR SYMMETRIC GAUSSIAN STATES
Let us consider the QFIM having diagonal elements H 11 and off-diagonal elements H 12 , which is then written as represents the standard basis. By introducing

APPENDIX D: MINIMIZATION OF THE ESTIMATION ERROR WHEN PROBING WITH SYMMETRIC GAUSSIAN STATES
For pure symmetric Gaussian states, the elements of the covariance matrix satisfy where 0 r cosh −1 (2N/M + 1)/2. Under the above constraint, our task boils down to finding parameters that maximize whose maximum value can be shown to be 8N (N + 1), in general.
One can easily check that if we use which is exactly the same as the lower bound for singlemode phase estimation using a squeezed thermal probe state, as shown in Ref. [43]. Note that preparing the squeezed thermal state input without a photon-loss channel is equivalent to using the optimal entangled Gaussian state with a photon-loss channel after adjusting appropriate parameters when the photon-loss rates are equal to each other [47]. Let us find the implementation of the Gaussian measurement corresponding to M . Noticing that mixing a p-squeezed state and (M − 1) vacua by the first beam splitter network (BSN) in the main text generates the optimal entangled state,ˆ 0 can be represented bŷ 0 =Û BSN |r, 0, ..., 0 r, 0, ..., 0|Û † BSN , where |r r| represents a p-squeezed state with a squeezing parameter r and |0 0| is a vacuum state. Since a BSN transforms a displacement operator into another displacement operator and a single-mode Gaussian measurement of ζ =D(ζ)ˆ 0D † (ζ)/π can be implemented by the generaldyne measurement [48], the Gaussian measurement can be performed by general-dyne measurement on M modes after the second BSN that processes the reverse of the first BSN generating the optimal entangled state. Especially whenˆ 0 is a vacuum, the general-dyne measurement reduces to a heterodyne measurement.