Phase space formalism for quantum estimation of Gaussian states

# Phase space formalism for quantum estimation of Gaussian states

Alex Monras Centre for Quantum Technologies, National University of Singapore, 2 Science Drive 3, 117542, Singapore
###### Abstract

We formulate, with full generality, the asymptotic estimation theory for Gaussian states in terms of their first and second moments. By expressing the quantum Fisher information (QFI) and the elusive symmetric logarithmic derivative (SLD) in terms of the state’s moments (and their derivatives) we are able to obtain the noncommutative extension of the well known expression for the Fisher information of a Gaussian probability distribution. Focusing on models with fixed first moments and identical Williamson ’diagonal’ states –which include pure state models–, we obtain their SLD and QFI, and elucidate what features of the Wigner function are fundamentally accessible, and at what rates. In addition, we find the optimal homodyne detection scheme for all such models, and show that for pure state models they attain the fundamental limit.

###### pacs:
03.65.Ta, 03.67.-a, 06.20.Dk, 42.50.St

Estimation theory plays a central role in modern developments of quantum enhanced metrology. From a practical point of view, it allows to assess the ultimate precision limits of given metrological schemes. From a fundamental perspective, it provides a gold standard upon which to asses distinguishability of quantum states. Quantum estimation theory is an old subject Helstrom (1976); Holevo (1982) and it has seen huge developments over the last 20 years Braunstein and Caves (1994); Boixo et al. (2007); Giovannetti et al. (2011). It has played a major role in understanding the fundamental powers and limitations of quantum measurement. The groundbreaking advances in atomic clocks and precision metrology Appel et al. (2009); Ma et al. (2011) suggest that quantum estimation theory will only become more relevant as measurement precision reaches its fundamental limits.

Common to almost all the disciplines of physics where precision metrology can provide significant results, is the fact that they benefit from the simplicity and power of the Gaussian state formalism Weedbrook et al. (2012). The latter has already proven its success and serves as an invaluable tool in describing quantum states of light and atomic ensembles, as well as providing useful insight and intuition.

Despite the great success of both the Gaussian state formalism, and quantum estimation theory, these two have never been successfully merged. Indeed, the nontrivial equations defining central objects in quantum estimation theory often forces to numerical methods Adesso et al. (2009); Dorner et al. (2009); Demkowicz-Dobrzanski et al. (2009) for computing precision bounds and determining optimal measurements, and these difficulties are only aggravated by the infinite-dimensional nature of bosonic systems. Most remarkably, these difficulties are not alleviated by the Gaussian state formalism, but for very particular cases and without explicit harvest of the phase-space structure. Ironically, it is in the seminal book by Holevo Holevo (1982) where one encounters one of the first comprehensive accounts of the phase-space formalism and Gaussian states on the one hand, and the foundations of quantum estimation theory on the other. There, certain essential features of Gaussian states regarding optimal detection were established, but the analysis focused on the Gaussian shift model. This model, despite being extremely relevant, and certainly the first candidate to be studied, is not well suited for modern entanglement-enhanced metrology, where the signal is encoded in the state’s correlations rather than the amplitudes. It is the purpose of the present work to provide a fully general phase-space formulation of the central quantities in quantum estimation theory, namely, the symmetric logarithmic derivative (SLD), and the SLD quantum Fisher information, with focus on general Gaussian states. After presenting this new formulation, we use it to study a wide class of models which include all pure state models. We also address the optimality of Gaussian measurements and derive sufficient conditions under which they are optimal.

Given a quantum model, i.e. a parametrized set of quantum states , the quantum Cramér-Rao bound (QCRB) establishes a lower bound to the variance of any unbiased estimator of parameter ,

 (Δ^θ)2≥IQ(θ)−1 (1)

where is the SLD quantum Fisher information (QFI). The QFI is defined in terms of the symmetric logarithmic derivative (SLD), which is the Hermitian operator that satisfies

 ∂θρθ=12(ρθLθ+Lθρθ), (2)

 IQ(θ)=tr[ρθL2θ]. (3)

Although the QCRB only establishes a bound, for uniparametric models it is asymptotically attainable Gill and Massar (2000); Gill (2001); Paris (2009). Therefore, it establishes the best asymptotic rate at which statistical fluctuations can decrease when measuring parameter  Helstrom (1976); Holevo (1982); Braunstein and Caves (1994). Thus, it is a quantitative measure of distinguishability of a state from its neighbors and as such, is intimately related to the quantum Fidelity and the Bures distance Bengtsson and Życzkowski (2006). On the other hand, the SLD not only has a geometric meaning; it also represents, by construction, an optimal observable Gill (2001); Paris (2009), in the sense that –to leading order– is proportional to the deviation from a reference state , , and it minimizes the statistical fluctuations.

We focus our analysis on systems of bosonic modes, described by Hilbert space . These systems are characterized by canonical variables with canonical commutation relations , where is the symplectic matrix, and is the degenerate symplectic inner product, 111Summation of latin indices and runs from 1 to , except for symbol which labels the modes and thus runs from to . The Weyl (displacement) operators are defined as , so that , and . In addition, we introduce the symmetric product . The symmetric product is not associative. We define it to have precedence over ordinary product, .

Gaussian states are defined as those states having Gaussian characteristic function Let be the space of matrices over . We define the first and second moments as

 di =tr[Riρ] (4a) Γij =2tr[(Ri−di)∘(Rj−dj)ρ], (4b)

where and . With these definitions, a Gaussian characteristic function reads

 χρ(ξ)=exp(iξ⊤¯d−14ξ⊤¯Γξ) (5)

where and .

The main goal of this work is to provide a formulation of Eqs. (2) and (3) in terms of and , and to explore the benefits of such formulation. Let be a Gaussian model with parameter . Any such model is fully described by the first and second moments and of . As can be expected from the structure of Gaussian states, –and we show in Appendix A– the SLD is quadratic in the canonical operators. It has also zero expectation as follows from its definition Eq. (2) and . Taking as ansatz for the expression

 Lθ= ∑ijLij(Ri−di)∘(Rj−dj) +∑ibi(Ri−di)−12tr[LΓ] (6)

in Eq. (2) and taking the characteristic function of both sides one can relate the derivative of w.r.t to , (lhs) to expectations of , and (rhs). This relation in turn, implies a relation between and (lhs) and suitable linear combinations of (rhs), where . Equating terms with equal powers of yields solutions to and in Eq. (6). We report the solutions here and the technical details in Appendix A. Define the linear map as

 DX(Y) =XYX⊤−ωYω⊤. (7)

With this, and are given by

 b =2Γ−1∂d, (8a) L =D−1Γ(∂Γ), (8b)

where the a pseudoinverse (Moore-Penrose inverse) is understood whenever is singular. The map in Eq. (8b) will play a central role in our discussion. Notice that is symmetry preserving. Let . Then satisfies the Stein equation, , where and its spectrum is contained in the unit disk (). Assuming is nonsingular, and the unique solution to the Stein equation is Bhatia (2006)

 Y=−∞∑k=0Fk∂Γ−1F⊤k. (9)

A more detailed analysis of the structure of and its pseudoinverse is given in Appendix A. The following relation will be useful,

 ∂Γ=−DΓ(∂Γ−1)−ω∂Γ−1ω⊤, (10)

as follows from .

As a first observation regarding the structure of the SLD obtained in Eq. (6), setting in Eq. (8b) recovers the Gaussian shift model, already studied extensively in the literature Holevo (1982); Gill and Guta (2012). In this case we obtain , linear in the canonical operators, hence recovering the optimality of homodyne detection.

Beyond the Gaussian shift model, for generic Gaussian models the SLD is at most quadratic in the canonical operators. Defining operators one can write

 Lθ=∑ijLij^Ri∘^Rj+C (11)

where is a scalar. In addition, as follows from Eqs. (8b) and (9) the matrix is the image of under the action of a completely-positive map. Therefore, whenever or is positive semidefinite, is so too, and there is a symplectic transformation such that . Thus, defining and the number operators, one has

 Lθ=n∑k=1αk(Nk−⟨Nk⟩θ). (12)

Therefore, in the case where has definite signature, the SLD reduces to photon counting in suitably defined modes. This result extends and generalizes the findings of Monras and Paris (2007); Monras and Illuminati (2010), where the SLD was shown to have this structure for certain classes of channel estimation problems. Still, despite having a physical interpretation of the operator (Gaussian transformations and photon counting), implementing a measurement of it may still be prohibitive, especially if strong squeezing of the signal is required.

Before considering more practical measurements, let us obtain the fundamental limit to their performance. From Eq. (6) one can readily obtain the expression for the quantum Fisher information for general Gaussian models. It is convenient at this point to endow with an inner product structure, , thus regarding matrices as vectors (kets) , and linear maps thereof as , so that we can write and Eq. (8b) reads . With this inner product, and its (pseudo)inverse are self-adjoint. The QFI reads [see Appendix B for details]

 IQ=12(∂Γ|(Γ⊗Γ−ω⊗ω)−1|∂Γ)+2∂d⊤Γ−1∂d. (13)

This expression is, for the first time, the most general form of the QFI in the Gaussian state formalism, and together with Eqs. (6) and (8) constitutes our main result. It allows to compute precision bounds for a number of situations, and expresses the Fisher information in a form amenable for numerical computation, overcoming the difficulties posed by the infinite-dimensional character of Eq. (2). Notice that putting the dimensions back in and taking the classical limit () yields , thus recovering

 Icl=12tr[∂ΓΓ−1∂ΓΓ−1]+2∂d⊤Γ−1∂d, (14)

the Fisher information of a Gaussian probability distribution centered at with covariance 222The factor of 2 in the second summand is due to the factor 2 in the definition of the covariance matrix, Eq. (4b). Thus, Eq. (13) is the noncommutative generalization of Eq. (14). Indeed, plays an essential role in capturing the geometry and distinguishability properties of Gaussian states, and the term proportional to accounts for the uncertainty due to noncommutativity of the canonical observables.

Eq. (9) gives the unique solution for models with nonsingular , namely, no vacuum modes in the Williamson decomposition. This rules out the important case of pure state models. However, a generic solution to the Stein equation can be given for the class of models with covariance matrix satisfying the relation

 (Γω)2=−ν2, (15)

where is constant and we will treat it as a scalar. We call these models isothermal due to the constant temperature of their Williamson decomposition. These models include all pure state Gaussian models, as well as some mixed models that have recently attracted attention Blandino et al. (2012). It is easy to check that Eq. (15) implies , which combined with Eq. (10) leads to . Despite some technicalities in simplifying when is singular, the resulting expression is also valid for pure state models, as can be shown by a detailed analysis of the singular case (). The QFI for such models is readily obtained

 IQ=12ν21+ν2tr[∂ΓΓ−1∂ΓΓ−1]+2∂d⊤Γ−1∂d. (16)

First, notice that the contribution due to first moments is equal to that of Eq. (14). This is due to the fact that there is always a reference frame for which the derivative is along a set of mutually compatible variables. In addition, setting recovers the known expression for the QFI for pure state models Pinel et al. (2012). Interestingly, the correction factor for thermal models approaches 1 in the large temperature limit (), recovering the Fisher information contained in the Wigner distribution. As noncommutativity of the canonical variables dictates, a faithful sampling of the Wigner distribution is beyond reach except in the high temperature regime, when thermal fluctuations render quantum fluctuations irrelevant. This explains why, in the high temperature limit, approaches the Fisher information of the Wigner distribution. In addition, notice that for pure state models with fixed first moments (), is exactly 1/2 of the Fisher information contained in the state’s Wigner function. This fact deserves further attention. Consider the the model in the Williamson form

 Γθ=SθνS⊤θ. (17)

The essential quantity in Eq. (16) is , which is nothing but expressed in the coordinates for which is diagonal. However, the canonical transformation is only determined up to an orthogonal symplectic transformation. Conveniently, is symmetric and Hamiltonian (), thus there is always a symplectic orthogonal transformation that diagonalizes it. Therefore, defining the symplectic transformation and corresponding canonical coordinates , one has

 TΓT⊤=(ν00ν),T∂ΓT⊤=(νλ00−νλ). (18)

where and we have made the dependency explicit in . This illustrates a characteristic trait of all models of the kind (15), i.e., that there exist canonical coordinates for which variations of the parameter correspond to a collection of single-mode squeezing operations. In addition, it is clear that

 IQ=ν21+ν2tr[λ2]. (19)

Now consider the class of homodyne measurements obtained by measuring quadratures (or ). The outcomes of such measurements are Gaussian distributed with covariance matrix and , hence yielding given by

 I⋆cl=12tr[λ2]. (20)

One immediately sees that for pure states () this is optimal, as follows from Eq. (19), . Hence, no other measurement can perform better. Let us pause for a moment to discuss what we mean by homodyne detection in this general multimode setting. In principle, we regard homodyne detection as the measurement of any set of compatible canonical variables, labelled for a suitably chosen canonical coordinates, . However, the set of variables is only a convenient way to account for the information that the measurement outcomes provide. The real estimator is some linear combination of the outcomes, and thus can be described as a linear combination of the canonical observables , where . Expressing in terms of the original coordinates we get , with . Hence, it remains to be seen that any such linear combination can be implemented by simple linear combinations of the original quadratures or passive Gaussian transformations thereof, and does not depend critically on the simultaneous measurement of incompatible variables or on some impractical active transformations. Parametrizing passive transformations as orthogonal symplectic matrices , we seek such that , where . Writing with blocks we get

 Vα=(cs−sc)(αqαq)=(cαq+sαp−sαq+cαp), (21)

so the condition that is always achievable by e.g., and . Let , then [] shows that any linear combination of the canonical operators is implementable by passive transformations and homodyne detection on the original modes.

Going back to Eqs. (19) and (20), it is clear that for mixed states () there is, potentially, room for improvement, as . More general Gaussian measurements can be implemented by attaching a Gaussian ancilla and performing homodyne detection Giedke and Cirac (2002); Eisert et al. (2002). Adding an ancilla corresponds to replacing and , thus . Therefore, without loss of generality we can focus on homodyne detection. Consider the Fisher information of observables out of a generic set of canonical coordinates , parametrized by a symplectic matrix

 U=(ab⋅⋅),ab⊤=ba⊤. (22)

Outcomes are distributed according to the covariance matrix , where the subindex indicates that the block corresponding to quadratures has to be taken. Likewise, the derivative is of the form . One can check that this model yields Fisher information

 Icl(a,b)=12tr[(ϕa(λ)−ϕb(λ))2] (23)

where , with . The maps are completely positive and is unital. Using the property one can show that is also trace-preserving. Finally, adding to Eq. (23) yields , and using Kadison-Schwartz inequality one gets because is trace-preserving.

This shows that, for models of the form of Eq. (15), with fixed first moments, homodyne detection of a suitable quadrature is always optimal among Gaussian measurements, with Fisher information upper bounded by , and heterodyne detection cannot improve its performance. Moreover, for pure models it is optimal in a wider sense –among all quantum measurements–. The optimality of homodyne detection generalizes and puts in context some earlier results Monras (2006); Pinel et al. (2012).

Further work can be envisaged in various directions. On the more practical side, identifying the most general class of models for which homodyne detection is optimal. On the theoretical side, application of these results to the study of Gaussian channels is a natural way of proceeding. In addition, one expects that some extension of the methods used here may be able to provide insight to paradigmatic non-Gaussian models such as phase diffusion and degaussified states.

We thank H. Cable, M. Hayashi and B.-G. Englert for useful discussions. The Centre for Quantum Technologies is funded by the Singapore Ministry of Education and the National Research Foundation as part of the Research Centres of Excellence programme.

## Appendix A The symmetric logarithmic derivative

We will make use of Einstein’s summation convention. In addition, we will make a distinction between covariant and contravariant indices to bookkeep the transformation rules to which they comply to. In this spirit, the inverse of the symplectic metric is , so that . We give a different symbol to avoid confusion, because componentwise . Also, we define , and with no index refers to .

Taking the trace of Eq. (2) with the Weyl operators we get

 ∂θχρ(ξ)=tr[Lθ∘ρW(ξ)]. (24)

As an ansatz, suppose is at most quadratic in the canonical operators

 Lθ=L(0)+L(1)iRi+L(2)ijRi∘Rj, (25)

so the RHS of Eq. (24) can be written as

 tr[Lθ∘ρW(ξ)]= L(0)+L(1)itr[ρRi∘W(ξ)] +L(2)ijtr[ρ(Ri∘Rj)∘W(ξ)]. (26)

On the other hand, the Gaussian characteristic function allows to write

 ∂θχρ(ξ)=( iξ∂θ¯d−14ξ⊤∂θ¯Γξ)χρ(ξ). (27)

where and in proper covariant-contravariant notation read and .

We begin by developping the RHS. The Weyl operators are defined as

 W(ξ)=exp(iΩ(ξ,R)), (28)

and we can show

 ∂iW(ξ)= iΩijRj∘W(ξ), (29)

or equivalently

 Ri∘W(ξ)=−iωij∂jW(ξ) (30)

It is tempting to define lower-indexed operartors so that one can write . However we will avoid this notation since the metric is antisymmetric and it would be essential to be consistent on which index is contracted with. The second derivative reads

 ∂i ∂jW(ξ)=−ΩikΩjlRk∘(Rl∘W(ξ)). (31)

The symmetric product is not associative . To put Eq. (31) in a manifestly symmetric form we use

 Rk∘(Rl∘W(ξ)) =(Rk∘Rl)∘W(ξ)−14ξkξlW(ξ) (32)

 ∂ijW(ξ)=−ΩikΩjl(Rk∘Rl−14ξkξl)∘W(ξ) (33)

which is manifestly symmetric under . From this we can derive the relation

 (Rk∘Rl)∘W(ξ)= (14ξkξl−ωkiωlj∂ij)W(ξ). (34)

Combining Eqs. (30) and (34) with Eq. (26) we get

 tr[ ρ∘LθW(ξ)] (35) =(L(0)−iL(1)iωik∂k+L(2)ij(14ξiξj−ωikωjl∂kl))χρ(ξ).

Using Eqs. (27) and (35), Eqs. (24) reads

 ( iξ∂¯d−14ξ⊤∂¯Γξ)χρ(ξ) (36) =(L(0)−iL(1)iωik∂k+L(2)ij(14ξiξj−ωikωjl∂kl))χρ(ξ).

It is straightforward to evaluate the derivatives of . These are

 −iL(1)iωik∂kχρ(ξ) =(L(1)idi−i2L(1)iΓikΩklξl)χρ(ξ) (37)

and

 L(2)ij( 14ξiξj−ωikωjl∂kl)χρ(ξ)= (38) = (L(2)ij(didj+12Γij−i2L(2)ij(diΓjk+djΓik)Ωklξl −14L(2)ij(ΓikΓjl−ωikωjl)ΩkrΩlsξrξs)χρ(ξ).

Rewriting Eq. (36), and recalling that is nowhere zero, we get

 iξ∂¯d− 14ξ⊤∂¯Γξ= (39) = L(0)+L(1)idi+L(2)ij(didj+12Γij) −i2(L(1)iΓik+L(2)ij(diΓjk+djΓik))Ωklξl −14L(2)ij(ΓikΓjl−ωikωjl)ΩkrΩlsξrξs,

Equaling the different orders in we get

 L(0)+L(1)d+d⊤L(2)d+12tr[L(2)Γ] =0 (40a) L(1)+(L(2)+L(2)⊤)d =2Γ−1∂d (40b) ΓL(2)Γ−ωL(2)ω⊤ =∂Γ. (40c)

By defining the linear map as

 DX(Y)=XYX⊤−ωYω⊤, (41)

Eq. (40c) can be written as

 DΓ(L(2))=∂Γ, (42)

and a solution exists when lies within the range of . Using to denote the pseudoinverse of and the fact that we rewrite Eqs. (40) as

 L(2)= D−1Γ(∂Γ) (43a) L(1)= 2Γ−1∂d−2D−1Γ(∂Γ)d (43b) L(0)= d⊤D−1Γ(∂Γ)d−2d⊤Γ−1∂d−12tr[D−1Γ(∂Γ)Γ]. (43c)

Putting this back into Eq. (25) we get

 Lθ= D−1Γ(∂Γ)ij(Ri−di)∘(Rj−dj) (44) +2∂diΓ−1ij(Rj−dj) (45) −12tr[D−1Γ(∂Γ)Γ] (46)

This proves Eqs. (6) and (8).

As suggested above may be singular. To analyze the spectrum of it is convenient to consider in Williamson form, where is a (thermal) Gibbs state. In this representation,

 DΓ=(S⊗S)DΓth(S⊗S)⊤ (47)

and the eigenvalue equation reads . Since is diagonal and proportional to the identity in each block corresponding to each mode, the eigenvalue equation decouples into independent equations, labelled by the pair , corresponding to the -th and -th modes in the row and column resp.

 ωXijω⊤=(νiνj−λ)Xij. (48)

The map is an involution, thus its eigenvalues are . Hence, the eigenvalues of are of the form , where are the symplectic eigenvalues of . In addition, matrices and are eigenvectors of corresponding to eigenvalue , and Pauli matrices and correspond to eigenvalue . Hence are eigenvectors of with eigenvalue , and are associated to . Thus, the spectrum of is given by . In particular, one can see that if has symplectic eigenvalues equal to 1 (or is singular), then is singular, and so is . The kernel of consists of sparse matrices populated only in blocks such that , with entries of the form ,

 ker(DΓth)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝0⋯0⋯0⋮⋱⋮⋮⋮⋯α1+βω⋯0⋮⋮⋱⋮0⋯0⋯0⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (49)

Thus is nonsingular only when the symplectic eigenvalues of are strictly greater than 1, and if is singular, its kernel consists of matrices with with nonzero entries in blocks corresponding to vacuum modes in the Williamson decomposition of .

## Appendix B The quantum Fisher information

The QFI is defined as . Using the notation introduced in the text, a map reads and the action is specified as follows

 A⊗B|X)=|AXB⊤). (50)

With this has 4 upper indices , of which the last two are contracted with the argument matrix,

 [DX(Y)]ij =[DX]ijklYkl =(XikXjl−ωikωjl)Ykl =XikYklX⊤lj−ωikYklω⊤lj =[XYX⊤−ωYω⊤]ij. (51)

On the other hand, supposing that is nonsingular, , and

 [D−1X(Y)]ij=[D−1X]ijklYkl, (52)

and

 [D−1X]ijkl[DX]klrs=δriδsj. (53)

We will use the relation Monras and Illuminati (2010)

 tr[ρθ(^Ri∘^Rj)∘(^Rk∘^Rl)]=14[ΓijΓkl+[DΓ]ijk′l′(δkk′δll′+δkl′δlk′)], (54)

where the centered canonical operators have been defined. With this, reads

 IQ= D−1Γ(∂Γ)ijD−1Γ(∂Γ)kltr[ρθ(^Ri∘^Rj)∘(^Rk∘^Rl)] (55a) +(4[Γ−1∂d]i[Γ−1∂d]j−tr[D−1Γ(∂Γ)Γ]D−1Γ(∂Γ)ij)tr[ρθ^Ri∘^Rj]+14tr[D−1Γ(∂Γ)Γ]2 = 14D−1Γ(∂Γ)ijD−1Γ(∂Γ)kl[DΓ]ijk′l′(δkk′δll′+δkl′δlk′)+4[Γ−1∂d]i[Γ−1∂d]jtr[ρθ^Ri∘^Rj] (55b) = 14D−1Γ(∂Γ)ijD−1Γ(∂Γ)kl[DΓ]ijk′l′(δkk′δll′+δkl′δlk′)+2∂d⊤Γ−1∂d. (55c)

Since is symmetric, we can write

 IQ =12D−1Γ(∂Γ)ij[DΓ]ijklD−1Γ(∂Γ)kl +2∂d⊤Γ−1∂d. (56)

Finally, notice that

 [DΓ]ijkl[D−1Γ(∂Γ)]kl=[DΓ∘D−1Γ(∂Γ)]ij=[P(∂Γ)]ij (57)

where is the projector onto the support of . Hence,

 IQ=12tr[D−1Γ(∂Γ)P(∂Γ)]+2∂d⊤Γ−1∂d (58)

To conclude, notice that is self adjoint, i.e., , as a consequence of being symmetric. So is