# Multi-parameter quantum estimation of single- and two-mode pure Gaussian
states

Gabriele Bressanini  Blackett Laboratory, Imperial College London, London SW7
2AZ, United Kingdom  Marco G. Genoni  Dipartimento di Fisica Aldo Pontremoli,
Universitá degli Studi di Milano, I-20133 Milano, Italy  M. S. Kim  Blackett
Laboratory, Imperial College London, London SW7 2AZ, United Kingdom  Matteo G.
A. Paris  Dipartimento di Fisica Aldo Pontremoli, Universitá degli Studi di
Milano, I-20133 Milano, Italy

###### Abstract

We discuss the ultimate precision bounds on the multiparameter estimation of
single- and two-mode pure Gaussian states. By leveraging on previous
approaches that focused on the estimation of a complex displacement only, we
derive the Holevo Cramér-Rao bound (HCRB) for both displacement and squeezing
parameter characterizing single and two-mode squeezed states. In the single-
mode scenario, we obtain an analytical bound and find that it degrades
monotonically as the squeezing increases. Furthermore, we prove that
heterodyne detection is nearly optimal in the large squeezing limit, but in
general the optimal measurement must include non-Gaussian resources. On the
other hand, in the two-mode setting, the HCRB improves as the squeezing
parameter grows and we show that it can be attained using double-homodyne
detection.

##  I Introduction

Quantum devices of interest for quantum technology are complex systems
characterized by many parameters. Their use in effective protocols requires
the precise characterization the underlying physical systems. This process,
vital from both a theoretical and technological standpoint, constitutes the
current challenge in quantum sensing and metrology. Quantum estimation theory
provides the tools to establish the ultimate limits on the precision of
parameter estimation in the quantum domain, and aims to identify potential
advantages with respect to classical protocols by leveraging quantum
resources, including entanglement and squeezing [1, 2, 3, 4, 5, 6, 7, 8, 9].
Multiparameter quantum metrology [10, 11, 12, 13] has received much attention
in the last years, ranging from the joint estimation of unitary parameters
[14, 15, 16, 17, 18, 19, 20, 21, 22], of unitary and loss parameters [23, 24,
25, 26, 27, 28], and for both spatial and time superresolution imaging [29,
30, 31, 32, 33, 34, 35]. From the theoretical point of view, the derivations
of the ultimate bounds on the estimation precision relies on the seminal works
by Helstrom [36] and Holevo [37]; by inspecting these derivations it is
immediate clear how in the quantum realm the multiparameter bounds are not a
trivial generalization of the single-parameter ones, as it is indeed the case
in the classical scenario. In fact, the potential non-commutativity of
quantum-mechanical observables may lead to the incompatibility of the optimal
measurements corresponding to each single parameter. For this reason, the
standard quantum Cramér-Rao bound based on the symmetric logarithmic
derivative (SLD) operators is in general not tight, and the departure from
this bound has been indicated as a signature of quantumness of the
corresponding quantum statistical model [38, 39, 40]. The Holevo Cramér-Rao
bound (HCRB) [37] is a tighter bound which may differ by no more than a factor
of two from the SLD Cramér-Rao bound [38, 41] and it is in principle
attainable by performing a collective measurement on an asymptotically large
number of copies of the quantum state encoding the parameters. If one
restricts to separable measurements on single copies, a tighter bound than the
HCRB exists which is referred to as the Nagaoka-Hayashi bound (NHB) [42, 43,
21]. Nevertheless, it was shown that the HCRB is achievable at the single-copy
level, and thus equal to the NHB, for quantum statistical models corresponding
to pure states [44] and for displacement estimation tasks with Gaussian probe
states [37]. However, its evaluation relies on a non-trivial function
minimization, which can be recast as a semidefinite program [45]. Numerical
results have been presented in the context of error-corrected multiparameter
quantum metrology [46] and for 3D magnetometry [47]. Analytical closed
formulas for the HCRB are hard to obtain and, to the best of our knowledge,
they have been derived for two-parameter estimation with pure states [44], for
generic qubit systems [48], for displacement estimation with Gaussian states
[16, 15], for light polarization [49] and for three-parameter rotations to
two-qubit states [21].

In this work, we leverage on the methods introduced in Ref. [16] and extend
the evaluation of the HCRB to the estimation of both displacement and
squeezing, for single- and two-mode displaced squeezed states. In the single-
mode case, we provide an analytical expression of the HCRB for the three-
parameter quantum statistical model. Furthermore, we prove that the optimal
measurement scheme must include non-Gaussian resources. We also identify the
best-performing general-dyne measurement and observe how that becomes nearly
optimal in the large squeezing regime. In the two-mode example, we resort to
numerical methods for the evaluation of the HCRB and prove that the double
homodyne measurement scheme is optimal in all the regimes considered, allowing
us to infer the analytical form of the bound.

The manuscript is organized as follows. In Sec. II we give a brief
introduction to the Gaussian quantum state formalism for bosonic quantum
systems, while in Sec. III we introduce the basic ingredients of
multiparameter quantum estimation theory. In Sections IV and V we present our
main results, i.e. we evaluate the HCRB for single- and two-mode displaced
squeezed vacuum states respectively, along with the analysis of the
performance of Gaussian general-dyne measurements schemes. We conclude the
manuscript in Sec. VI with some final remarks and outlooks.

##  II Gaussian quantum states

In this section we review key aspects of the Gaussian formalism, including
Gaussian states and Gaussian measurements. We refer to the following
references for a more detailed introduction [50, 51, 52]. Let us consider a
continuous variables (CV) quantum systems made up of $d$ bosonic modes
described by annihilation operators $\hat{a}_{j}$ with $j=1,\dots,d$, that
satisfy the standard bosonic commutation relations
$[\hat{a}_{j},\hat{a}^{\dagger}_{k}]=\delta_{jk}\hat{\mathcal{I}}$. Here,
$\delta_{jk}$ denotes the Kronecker delta and $\hat{\mathcal{I}}$ is the
identity operator. For each mode, one can then introduce the quadrature
operators

$\hat{q}_{j}=\frac{\hat{a}+\hat{a}^{\dagger}}{\sqrt{2}}\,,\quad\hat{p}_{j}=\frac{\hat{a}-\hat{a}^{\dagger}}{i\sqrt{2}}$
(1)

that satisfy the canonical commutation relations (CCR)
$[\hat{q}_{j},\hat{p}_{k}]=i\delta_{jk}\hat{\mathcal{I}}$, where we have set
$\hbar=1$. It is convenient to arrange these operators in a vector

$\bm{\hat{r}}=(\hat{q}_{1},\hat{p}_{1},\dots,\hat{q}_{d},\hat{p}_{d})^{\intercal}\,,$
(2)

and the CCR can now be expressed compactly as

$[\bm{\hat{r}},\bm{\hat{r}}^{\intercal}]=i\Omega\,,$ (3)

where the symplectic matrix $\Omega$ is given by

$\Omega=\bigoplus_{j=1}^{d}\Omega_{1}\,,\quad\quad\Omega_{1}=\begin{pmatrix}0&1\\\
-1&0\end{pmatrix}\,.$ (4)

A generic Gaussian state $\hat{\rho}_{G}$ is fully characterized by its vector
of first moments $\bm{\overline{r}}$ and its covariance matrix $\bm{\sigma}$,
defined respectively as

$\bm{\overline{r}}=\Tr[\hat{\rho}_{G}\bm{\hat{r}}]\,,$ (5)
$\bm{\sigma}=\Tr[\hat{\rho}_{G}\\{\bm{\hat{r}}-\bm{\overline{r}},(\bm{\hat{r}}-\bm{\overline{r}})^{\intercal}\\}]\,.$
(6)

Note that the conventions we use are such that the covariance matrix of a
single-mode coherent state is the identity matrix, i.e.
$\bm{\sigma}=\mathbb{I}_{2}$. One can show that a generic $d$-mode Gaussian
state can always be written as

$\hat{\rho}_{G}=\hat{D}^{\dagger}_{\bm{\overline{r}}}\hat{U}^{\dagger}(\bigotimes_{j}\hat{\nu}_{j}^{th}(\overline{n}_{j}))\hat{U}\hat{D}_{\bm{\overline{r}}}\,,$
(7)

where
$\hat{D}_{\bm{\overline{r}}}=e^{i\bm{\overline{r}}^{T}\Omega\bm{\hat{r}}}$ is
the displacement operator, $\hat{U}$ is a unitary operator generated by a
purely quadratic Hamiltonian in the mode operators, and
$\hat{\nu}_{j}^{th}(\overline{n}_{j})$ is the thermal state of a single-mode
free bosonic field. The latter reads

$\hat{\nu}_{j}^{th}(\overline{n}_{j})=\frac{1}{1+\overline{n}_{j}}\left(\frac{\overline{n}_{j}}{1+\overline{n}_{j}}\right)^{\hat{a}_{j}^{\dagger}\hat{a}_{j}}\,,$
(8)

where $\overline{n}_{j}$ is the mean number of bosons. It is useful to also
introduce the following parametrization of the displacement operator for the
single-mode case $d=1$, namely

$\hat{D}(\alpha)=e^{\alpha\hat{a}^{\dagger}-\alpha^{*}\hat{a}}=\hat{D}_{-\overline{\bm{r}}},\quad\overline{\bm{r}}=\sqrt{2}\left(\Re{\alpha},\Im{\alpha}\right)^{\intercal}\,.$
(9)

By virtue of the Euler decomposition of symplectic matrices, we can express a
generic single-mode Gaussian state as a displaced squeezed thermal state, i.e.

$\hat{\rho}_{G}=\hat{D}(\alpha)\hat{S}(\xi)\hat{\nu}^{th}(\overline{n})\hat{S}^{\dagger}(\xi)\hat{D}^{\dagger}(\alpha)\,,$
(10)

where $\hat{S}(\xi)=e^{\frac{1}{2}(\xi\hat{a}^{\dagger
2}-\xi^{*}\hat{a}^{2})}$ is the single-mode squeezing operator and
$\xi=re^{i\chi}\in\mathbb{C}$ is the (complex) squeezing parameter.

Next, we introduce general-dyne measurements, a general class of Gaussian
measurements whose corresponding positive-operator-valued measure (POVM)
originates from a generalization of the well known over-completeness relation
of coherent states [50, 51], namely

$\frac{1}{(2\pi)^{d}}\int
d^{2d}\bm{r}_{m}\,\hat{D}_{\bm{r}_{m}}^{\dagger}\,\hat{\rho}_{m}\,\hat{D}_{\bm{r}_{m}}=\hat{\mathcal{I}}\,,$
(11)

where $\hat{\rho}_{m}$ is a generic $d-$mode Gaussian state with null vector
of first moments and covariance matrix $\bm{\sigma}_{m}$. A general-dyne
detection is said to be ideal if $\hat{\rho}_{m}$ is a pure state. The POVM
associated with this measurement thus reads

$\hat{\Pi}_{\bm{r}_{m}}=\frac{\hat{D}_{\bm{r}_{m}}^{\dagger}\,\hat{\rho}_{m}\,\hat{D}_{\bm{r}_{m}}}{(2\pi)^{d}}\,.$
(12)

Hence, $\bm{\sigma}_{m}$ defines the specific measurement scheme, while
$\bm{r}_{m}$ represents the measurement outcome. The homodyne and heterodyne
detections are easily recovered in this formalism. For example, in a single-
mode scenario, the covariance matrix associated with the $\hat{q}$ quadrature
measurement is

$\bm{\sigma}_{m}=\lim_{s\rightarrow-\infty}\matrixquantity(e^{2s}&0\\\
0&e^{-2s})=\lim_{z\rightarrow 0}\matrixquantity(z&0\\\ 0&\frac{1}{z}),$ (13)

while the heterodyne detection is retrieved with the substitution
$\bm{\sigma}_{m}=\mathbb{I}_{2}$ (the POVM elements are projectors on coherent
states). Lastly, the probability distribution of outcomes of a general-dyne
detection on a Gaussian input state $\hat{\rho}_{G}$ characterized by vector
of first moments $\overline{\bm{r}}$ and covariance matrix $\bm{\sigma}$
corresponds to a multi-variate Gaussian distribution centered in
$\overline{\bm{r}}$ and with covariance matrix
$\bm{\Sigma}=(\bm{\sigma}+\bm{\sigma}_{m})/2$.

##  III Multi-parameter Quantum Estimation Theory

Let us consider a quantum statistical model $\hat{\rho}_{\bm{\theta}}$, i.e. a
family of quantum states labelled by a vector of $d$ real parameters
$\bm{\theta}=(\theta_{1},\dots,\theta_{d})^{\intercal}$. We want to estimate
the value of $\bm{\theta}$ from the outcomes of $M$ independent measurements
$\bm{x}=(x_{1},\dots,x_{M})^{\intercal}$ described by the POVM
$\hat{\Pi}_{x}$, using a suitable unbiased estimator
$\tilde{\bm{\theta}}(\bm{x})$. The accuracy of the latter may be evaluated in
terms of its covariance matrix, defined as

$\bm{V}(\tilde{\bm{\theta}})=\int
d\bm{x}\,p(\bm{x}|\bm{\theta})(\tilde{\bm{\theta}}(\bm{x})-\bm{\theta})(\tilde{\bm{\theta}}(\bm{x})-\bm{\theta})^{\intercal}\,.$
(14)

Here,
$p(\bm{x}|\bm{\theta})=\Pi_{j=1}^{M}p(x_{j}|\bm{\theta})=\Pi_{j=1}^{M}\Tr[\Pi_{x}\rho_{\bm{\theta}}]$,
where we have implicitly assumed that the $M$ measurements are independent of
each other. The covariance matrix satisfies the Cramér-Rao bound (CRB)

$\bm{V}(\tilde{\bm{\theta}})\geq\frac{1}{M}\bm{F}^{-1}\,,$ (15)

where $\bm{F}$ is the Fisher Information (FI) matrix, whose elements are
defined as

$\bm{F}_{\mu\nu}=\int
d\bm{x}\,p(\bm{x}|\bm{\theta})\left(\partial_{\mu}\log{p(\bm{x}|\bm{\theta})}\right)\left(\partial_{\nu}\log{p(\bm{x}|\bm{\theta})}\right)\,.$
(16)

Throughout this work we will use the notation
$\partial_{\mu}\equiv\partialderivative{\theta_{\mu}}$. In particular, for a
Gaussian probability distribution $p({\bm{x}}|\bm{\theta})$ centered in
$\overline{\bm{x}}$ and with covariance matrix $\bm{\Sigma}$, one has

$\bm{F}_{\mu\nu}=(\partial_{\mu}\overline{\bm{x}}^{\intercal})\bm{\Sigma}^{-1}(\partial_{\nu}\overline{\bm{x}})+\frac{1}{2}\hbox{Tr}\left[\bm{\Sigma}^{-1}(\partial_{\mu}\bm{\Sigma})\bm{\Sigma}^{-1}(\partial_{\nu}\bm{\Sigma})\right]\,.$
(17)

The CRB is attainable, at least asymptotically, by choosing a suitable
efficient estimator. Furthermore, quantum mechanics allows us to find tighter
precision bounds, which only depend on the statistical model. Let us introduce
the symmetric logarithmic derivative (SLD) operators, implicitly defined by
the following Lyapunov equation

$\partial_{\mu}\hat{\rho}_{\bm{\theta}}=\frac{\hat{L}_{\mu}\hat{\rho}_{\bm{\theta}}+\hat{\rho}_{\bm{\theta}}\hat{L}_{\mu}}{2}\,.$
(18)

For pure statistical models
$\hat{\rho}_{\bm{\theta}}=\outerproduct{\psi_{\bm{\theta}}}{\psi_{\bm{\theta}}}$
one can easily solve the equation above and show that the SLD assumes the
following simple form:

$\hat{L}_{\mu}=2\partial_{\mu}\hat{\rho}_{\bm{\theta}}=2(\ket{\psi_{\bm{\theta}}}\\!\\!\bra{\partial_{\mu}\psi_{\bm{\theta}}}+\ket{\partial_{\mu}\psi_{\bm{\theta}}}\\!\\!\bra{\psi_{\bm{\theta}}}).$
(19)

We then use the SLDs to define the Quantum Fisher Information (QFI) matrix

$\bm{Q}_{\mu\nu}=\Tr[\hat{\rho}_{\bm{\theta}}\frac{\hat{L}_{\mu}\hat{L}_{\nu}+\hat{L}_{\nu}\hat{L}_{\mu}}{2}]\,,$
(20)

which is in turn used to derive the matrix quantum Cramér-Rao bound [36]

$\bm{V}(\tilde{\bm{\theta}})\geq\frac{1}{M}\bm{Q}^{-1}\,.$ (21)

In the single-parameter scenario, the Cramér-Rao inequality Eq. (21) is a
scalar bound that can be attained by projective measurement of the (Hermitian)
SLD operator. On the other hand, in the multi-parametric setting the quantum
Cramér-Rao bound is in general not tight, because of the possible non-
commutativity of the SLD operators. As it is usually more convenient to deal
with scalar bounds, we use a real, positive, $d\times d$ weight matrix
$\bm{W}$ to introduce the following scalar inequality, which we refer to as
the SLD Cramér-Rao bound (SLD-CRB):

$\Tr[\bm{W}\bm{V}]\geq\Tr[\bm{W}\bm{Q}^{-1}]\equiv
C^{S}(\bm{\theta},\bm{W})\,.$ (22)

Different $\bm{W}$ matrices are used to weigh the uncertainties of the
parameters differently. In this paper we will set $\bm{W}=\mathbb{I}_{d}$, so
that the inequality above bounds $\Tr[\bm{V}]$, i.e. the sum of the variances
of each parameter’s estimate. Like the corresponding matrix bound, the SLD-CRB
in Eq. (22) is also in general not attainable. We can also define the most
informative bound as a minimization of the classical bound over all possible
quantum measurements, i.e.

$C^{\textrm{MI}}(\bm{\theta})=\min_{\textrm{POVM}}\Tr[\bm{F}^{-1}]\,,$ (23)

which is in general larger than the SLD-CRB. However, there exists a tighter
bound, known as the HCRB, such that the following chain of inequalities holds

$\Tr[\bm{V}]\geq C^{\textrm{MI}}(\bm{\theta})\geq C^{H}(\bm{\theta})\geq
C^{S}(\bm{\theta})\,.$ (24)

The HCRB $C^{H}(\bm{\theta})$ is defined via the following minimization [37]

$C^{H}(\bm{\theta})=\min_{\hat{\bm{X}}}h_{\bm{\theta}}\left[\hat{\bm{X}}\right]\,,$
(25)

where $\hat{\bm{X}}=(\hat{X}_{1},\dotsc,\hat{X}_{d})$ is a vector of Hermitian
operators satisfying the locally unbiased conditions

$\displaystyle\Tr[\hat{\rho}_{\bm{\theta}}\hat{X}_{j}]$ $\displaystyle=0\,,$
(26) $\displaystyle\Tr[(\partial_{j}{\hat{\rho}_{\bm{\theta}}})\hat{X}_{k}]$
$\displaystyle=\delta_{jk}\,.$ (27)

Finally, the function to minimize reads

$h_{\bm{\theta}}[\hat{\bm{X}}]=\Tr[\Re{\bm{Z}_{\bm{\theta}}[\hat{\bm{X}}]}]+\norm{\Im{\bm{Z}_{\bm{\theta}}[\hat{\bm{X}}]}}_{1}\,,$
(28)

where $\norm{A}_{1}\equiv\Tr[\sqrt{A^{\dagger}A}]$ and
$\bm{Z}_{\bm{\theta}}[\hat{\bm{X}}]$ is a matrix of operators defined as

$\bm{Z}_{\bm{\theta}}[\hat{\bm{X}}]_{jk}=\Tr[\hat{\rho}_{\bm{\theta}}\hat{X}_{j}\hat{X}_{k}]\,.$
(29)

The HCRB is typically regarded as the most fundamental scalar bound in multi-
parameter quantum estimation theory, as it can be shown to be attainable by
performing a collective measurement on an asymptotically large number of
copies of the quantum state $\hat{\rho}_{\bm{\theta}}$ encoding the
parameters. Nevertheless, it was shown that the HCRB is actually attainable at
the single copy level for displacement estimation tasks of Gaussian states and
for quantum statistical models encoded in pure quantum states. The quantity
$C^{H}(\bm{\theta})$ may also be bounded [38, 41] as follows

$C^{S}(\bm{\theta})\leq C^{H}(\bm{\theta})\leq(1+R)C^{S}(\bm{\theta})\leq
2C^{S}(\bm{\theta})\,,$ (30)

where

$\displaystyle R=\norm{i\bm{Q}^{-1}(\bm{\theta})\bm{D}(\bm{\theta})}_{\infty}$
(31)

has been referred to as the asymptotic incompatibility (AI) of the
corresponding quantum statistical model [39, 40]. In the formula above,
$\norm{\bm{A}}_{\infty}$ denotes the largest eigenvalue in modulus of
$\bm{A}$, and $\bm{D}(\bm{\theta})$ is known as the Uhlmann curvature, whose
matrix elements are defined as

$\bm{D}_{\mu\nu}=-\frac{i}{2}\Tr[\hat{\rho}_{\bm{\theta}}[\hat{L}_{\mu},\hat{L}_{\nu}]]\,.$
(32)

The AI measure provides a bound for the true normalized gap between the HCRB
and the SLD-CRB bounds, as follows

$\displaystyle\frac{C^{H}(\bm{\theta})-C^{S}(\bm{\theta})}{C^{S}(\bm{\theta})}\leq
R\,.$ (33)

One can also prove [38] that $0\leq R\leq 1$, hence the SLD-CRB gives an
estimate of the HCRB up to a factor two. As it is also clear from Eq. (30) and
Eq. (32), the HCRB and the SLD-CRB coincide whenever the Uhlmann curvature
matrix elements are all equal to zero, that is whenever the average values of
the commutators of the SLD operators on the quantum statistical model
$\hat{\rho}_{\bm{\theta}}$ are equal to zero.

In the following Sections, we apply the tools of quantum estimation theory to
fully characterize general classes of single-mode and two-mode displaced
squeezed vacuum states, which correspond to three-parameter quantum
statistical models. The choice of focusing on pure statistical models
originates from the two following observations. First, from an experimental
point of view, the description of optical frequencies radiation modes at
temperatures accessible in a lab with the vacuum state is a very good
approximation. Furthermore, we recall that the HCRB is defined as a
minimization problem over a set of peculiar Hermitian matrices, however this
difficulty is eased by the particular structure of the quantum statistical
models considered and we are thus able to compute the bound analytically, or
at least to greatly simplify its numerical evaluation.

##  IV Multi-parameter quantum estimation of single mode pure Gaussian state

Let us consider a quantum statistical model given by single-mode displaced
squeezed vacuum states

$\ket{\psi_{\bm{\theta}}}=\hat{D}(\alpha)\hat{S}(r)\ket{0}\,,$ (34)

where $\alpha$ is a complex displacement and $r\geq 0$ is the (real) squeezing
parameter. As already mentioned, Eq. (34) represents the most general single-
mode pure Gaussian state (up to a phase shift), corresponding to Eq. (10) with
$\overline{n}=0$. We want to assess the ultimate bound to the joint estimation
precision of the three real parameters that characterize the statistical
model, namely

$\bm{\theta}=(\Re{\alpha},\Im{\alpha},r)^{\intercal}=(\theta_{1},\theta_{2},\theta_{3})^{\intercal}\,.$
(35)

In the following, we compute the SLD-CRB, the HCRB and discuss the
experimental attainability of the latter. In order to simplify the notation we
introduce the following Hilbert space basis

$\ket{e_{n}}=\hat{D}(\alpha)\hat{S}(r)\ket{n}\,,$ (36)

where $\ket{n}$ is the $n-$boson Fock state. Hence, the quantum statistical
model simply reads $\ket{\psi_{\bm{\theta}}}=\ket{e_{0}}$. In Appendix A, we
outline the detailed calculation of the SLDs and the QFI matrix. We find the
latter to be

$\bm{Q}=\matrixquantity(4e^{-2r}&0&0\\\ 0&4e^{2r}&0\\\ 0&0&2)\,.$ (37)

As expected due to symmetry reasons, $\bm{Q}$ does not depend on the
displacement $\alpha$. Additionally, the diagonal matrix elements of reveal
that, at the single parameter estimation level, squeezing does not affect the
estimation of $r$ itself, however it plays a role in the estimation of the
displacement. We can understand this intuitively: squeezing decreases the
variance of the $\hat{p}$ quadrature and increases that of the $\hat{q}$
quadrature. Consequently, even before computing the CRB, we may conclude that
squeezing is not a useful resource for the simultaneous estimation of the
three parameters: as $r$ grows, the estimation error of $\theta_{1}$
decreases, eventually reaching zero, while the variance of $\theta_{2}$ grows
indefinitely. This simple argument tells us that $r$ is not a useful resource
already in a displacement-only estimation scenario. We can now compute the
SLD-CRB, namely

$\Tr[\bm{V}^{-1}]\geq\Tr[\bm{Q}^{-1}]=C^{S}=\frac{1+\cosh{(2r)}}{2}\,.$ (38)

We also compute the quantumness $R$ of our statistical model. In particular,
the Uhlmann curvature $\bm{D}$ reads

$\bm{D}=\matrixquantity(0&4&0\\\ -4&0&0\\\ 0&0&0)\,,$ (39)

hence, via Eq. (31), we obtain $R=1$, i.e. the maximum degree of
incompatibility between the parameters. We remark that the same result would
be obtained for the estimation of the displacement parameters only, and the
maximum incompatibility can be ascribed to the incompatibility of the
$\hat{q}$ and $\hat{p}$ quadrature operators.

We are now ready to tackle the calculation of the HCRB. The derivatives of the
quantum statistical model with respect to the three parameters can be found in
the Appendix A (see Eqs. (84), (85) and (86) ). We recall that the HCRB is
defined as a minimization of the function $h_{\bm{\theta}}[\hat{\bm{X}}]$ in
Eq. (28) over all possible vectors of Hermitian operators $\hat{\bm{X}}$
satisfying the locally unbiased conditions Eq.(26) and Eq.(27). In our case
$\hat{\bm{X}}=\\{\hat{X}_{1},\hat{X}_{2},\hat{X}_{3}\\}$, and the constraint
imposed by Eq. (26) implies

$\bra{e_{0}}\hat{X}_{1}\ket{e_{0}}=\bra{e_{0}}\hat{X}_{2}\ket{e_{0}}=\bra{e_{0}}\hat{X}_{3}\ket{e_{0}}=0\,.$
(40)

Since our quantum statistical model is pure, Eq. (27) can be expressed as

$\delta_{jk}=\Tr[(\partial_{j}{\hat{\rho}_{\bm{\theta}}})\hat{X}_{k}]=\bra{\partial_{j}\psi_{\bm{\theta}}}\hat{X}_{k}\ket{\psi_{\bm{\theta}}}+\bra{\psi_{\bm{\theta}}}\hat{X}_{k}\ket{\partial_{j}\psi_{\bm{\theta}}}\,.$
(41)

The latter imposes three additional conditions for each of our Hermitian
operators. As an example, we find that $\hat{X}_{1}$ must satisfy

$\displaystyle 1$
$\displaystyle=\Tr[(\partial_{1}{\hat{\rho}_{\bm{\theta}}})\hat{X}_{1}]\Rightarrow\Re{\bra{e_{0}}\hat{X}_{1}\ket{e_{1}}}=\frac{e^{r}}{2}\,,$
(42) $\displaystyle 0$
$\displaystyle=\Tr[(\partial_{2}{\hat{\rho}_{\bm{\theta}}})\hat{X}_{1}]\Rightarrow\Im{\bra{e_{0}}\hat{X}_{1}\ket{e_{1}}}=0\,,$
(43) $\displaystyle 0$
$\displaystyle=\Tr[(\partial_{3}{\hat{\rho}_{\bm{\theta}}})\hat{X}_{1}]\Rightarrow\Re{\bra{e_{0}}\hat{X}_{1}\ket{e_{2}}}=0\,.$
(44)

These in turn imply

${\bra{e_{0}}\hat{X}_{1}\ket{e_{1}}}=\frac{e^{r}}{2}\,,\quad{\bra{e_{0}}\hat{X}_{1}\ket{e_{2}}}=i\beta\,,$
(45)

where $\beta\in\mathbb{R}$. Analogous calculations yield the following
conditions on the matrix elements of $\hat{X}_{2}$ and $\hat{X}_{3}$:

$\displaystyle{\bra{e_{0}}\hat{X}_{2}\ket{e_{1}}}$
$\displaystyle=-\frac{ie^{-r}}{2}\,,\quad$
$\displaystyle{\bra{e_{0}}\hat{X}_{2}\ket{e_{2}}}$ $\displaystyle=i\gamma\,,$
(46) $\displaystyle{\bra{e_{0}}\hat{X}_{3}\ket{e_{1}}}$
$\displaystyle=0\,,\quad$ $\displaystyle{\bra{e_{0}}\hat{X}_{3}\ket{e_{2}}}$
$\displaystyle=\frac{1}{\sqrt{2}}+i\delta\,,$ (47)

where $\gamma$ and $\delta$ are free real parameters. Clearly, the hermiticity
of $\hat{X}_{j}$ imposes additional constraints on the matrix elements of
these operators. Now comes the great simplification of dealing with pure
states statistical models: in Ref. [44] it is shown that we can set to zero
all the matrix elements of $\hat{X}_{j}$ not involved in the aforementioned
constraints. This is equivalent to saying that it is sufficient to study these
Hermitian operators in the Hilbert subspace spanned by the partial derivatives
of the statistical model. After some algebra, we obtain the matrix

$\bm{Z}_{\bm{\theta}}[\hat{\bm{X}}]=\matrixquantity(\frac{e^{2r}}{4}+\beta^{2}&\frac{i}{4}+\beta\gamma&\frac{i\beta}{\sqrt{2}}+\beta\delta\\\
-\frac{i}{4}+\beta\gamma&\frac{e^{-2r}}{4}+\gamma^{2}&\frac{i\gamma}{\sqrt{2}}+\gamma\delta\\\
-\frac{i\beta}{\sqrt{2}}&-\frac{i\gamma}{\sqrt{2}}+\gamma\delta&\frac{1}{2}+\delta^{2})\,.$
(48)

We can then use the definition in Eq. (28) to compute the function
$h_{\bm{\theta}}[\hat{\bm{X}}]$, namely

$h_{\bm{\theta}}[\hat{\bm{X}}]=\frac{\cosh{(2r)}+1}{2}+\beta^{2}+\gamma^{2}+\delta^{2}+\frac{1}{2}\sqrt{1+8\beta^{2}+8\gamma^{2}}\,,$
(49)

and perform a minimization over the real parameters $\beta,\gamma$ and
$\delta$. The optimization is trivial and yields the HCRB, i.e.

$C^{H}=\frac{\cosh{(2r)}+2}{2}\,.$ (50)

As we can see, optimal performance is achieved for $r=0$ and we also notice
that $C^{H}>C^{S}\,\,\forall r$, thus confirming that the SLD-CRB is in
general not tight. The latter coincides with $C^{H}$ only in the asymptotic
limit $r\rightarrow\infty$. This also means that in this case the AI measure
$R$ overestimates the true gap between $C^{H}$ and $C^{S}$. As mentioned
above, it has been proven in Ref. [44] that the HCRB is achievable via single-
copy measurements for pure quantum statistical models. However, we do not know
in general what measurement saturates said bound. In Ref. [14] it was shown
that if we are interested in estimating the complex displacement $\alpha$
only, then the optimal measurement corresponds to heterodyne detection. As we
now want to evaluate the the ultimate bound including the squeezing parameter,
we study the more general class of general-dyne detections, and compare their
performances to the HCRB. We remind the reader that the covariance matrix
$\bm{\sigma}_{m}$ which characterizes the specific general-dyne measurement
scheme is that of a generic Gaussian state with zero displacement. As we are
looking for the optimal measurement, we also restrict ourselves to ideal
(noiseless) measurements, thus corresponding to pure Gaussian states. This
means that the measurement is fully described by a covariance matrix of the
form

$\bm{\sigma}_{m}=\bm{S}\bm{S}^{\intercal}\,,$ (51)

where $\bm{S}$ is a symplectic matrix corresponding to a unitary evolution
induced by a purely quadratic Hamiltonian in the quadrature operators. The
Euler decomposition theorem states that any symplectic matrix can be expressed
in terms of symplectic matrices corresponding to phase rotations, beam
splitters and local squeezing only. However, since we are dealing with single-
mode statistical models, we do not need to include beam splitters and it is
easy to convince ourselves that phase rotations do not play a role (we have
also verified this intuition by numerical means), hence we may discard those
as well. As a result, $\bm{S}$ is just the usual symplectic matrix that
corresponds to single-mode squeezing, hence $\bm{\sigma}_{m}$ reads

$\bm{\sigma}_{m}=\matrixquantity(z&0\\\ 0&z^{-1})\,,$ (52)

with $z>0$. The covariance matrix of the statistical model is

$\bm{\sigma}=\matrixquantity(e^{2r}&0\\\ 0&e^{-2r})\,,$ (53)

and the conditional probability distribution of outcomes has Gaussian form,
centered in $\overline{\bm{x}}=(\theta_{1},\theta_{2})^{\intercal}$ and with
covariance matrix given by $\bm{\Sigma}=(\bm{\sigma}+\bm{\sigma}_{m})/2$. The
FI matrix can then be obtained via Eq. (17), yielding the following non-zero
matrix elements:

$\displaystyle\bm{F}_{\textrm{gen},11}$ $\displaystyle=\frac{4}{e^{2r}+z}\,,$
(54) $\displaystyle\bm{F}_{\textrm{gen},22}$
$\displaystyle=\frac{4}{e^{-2r}+z^{-1}}\,,$ (55)
$\displaystyle\bm{F}_{\textrm{gen},33}$
$\displaystyle=\frac{2(e^{4r}+z^{2})}{(e^{2r}+z)^{2}}\,.$ (56)

The bound on the estimation precision is thus given by

$\Tr[\bm{F}_{\textrm{gen}}^{-1}]=\frac{1}{4}\left(2+e^{-2r}+z+\frac{1}{z}+e^{2r}\left(1+\frac{4z}{e^{4r}+z^{2}}\right)\right)\,.$
(57)

For $z=1$ we obtain the heterodyne detection precision bound

$\Tr[\bm{F}_{\textrm{het}}^{-1}]=\frac{2\cosh^{4}{r}}{\cosh{(2r)}}\,,$ (58)

which coincides with the HCRB just in the $r\rightarrow\infty$ limit, thus
making the heterodyne detection nearly optimal in the large squeezing regime.
However, we also point out that as the squeezing parameter grows, the overall
estimation precision gets worse. We want to investigate whether Eq. (57) is
equal to the HCRB Eq. (50) for some values of $z$ and $r$. This is equivalent
to finding the zeros of the positive definite function

$f(z,r)=\Tr[\bm{F}_{\textrm{gen}}^{-1}]-C^{H}\,.$ (59)

By resorting to numerical methods, we discover that the equation $f(z,r)=0$
does not have any real solutions. Therefore, we conclude that the optimal
measurement must involve some non-Gaussian features, such as photon counting,
as it is often the case in several other single-parameter estimation problems
with Gaussian quantum statistical models [53]. Despite this, we may still ask
ourselves what is the best performance we can achieve with Gaussian resources
only. This translates to finding the value of $z$ that minimizes $f(z,r)$ at
fixed $r$. We call this value $z_{opt}(r)$ and display its behavior in
Fig.(1).

Figure 1:  The implicit function $z_{opt}(r)$, obtained from numerically
minimizing Eq. (59) at fixed value of the squeezing parameter $r$. We remind
that $z_{opt}$ determines the optimal measurement - among general-dyne
detection schemes - for our three-parameter estimation problem and that
$z_{opt}=1$ corresponds to heterodyne measurement.

We notice that, in the absence of squeezing, the heterodyne detection scheme
is the best one amongst general-dyne measurements, i.e. $z_{opt}(r=0)=1$. As
$r$ increases, at first $z_{opt}$ decreases until it reaches $z\sim 0.6$, and
then grows monotonically to asymptotically reach $1$ again. We also point out
that, as $r$ grows, the minimum of $f(z,r)$, i.e. $f(z_{opt}(r),r)\equiv
f_{opt}(r)$, decreases monotonically and approaches zero in the
$r\rightarrow\infty$ limit, consistently with our previous findings, as shown
in Fig. (2).

Figure 2: The function $f_{opt}(r)=f(z_{opt}(r),r)$, representing the
difference between the optimal estimation precision obtained with a POVM
belonging to the general-dyne class and the HCRB. The optimization is carried
out at fixed $r$. We notice that the function is monotonically decreasing in
$r$, hence the estimation performance achievable by the optimal general-dyne
matches the HCRB only for large values of squeezing.

##  V Multiparameter quantum estimation of displaced two-mode squeezed vacuum
states

We now consider the two-mode quantum statistical model described by displaced
two-mode squeezed vacuum states, where the displacement operator to one mode
only:

$\ket{\psi_{\bm{\theta}}}=\left(\hat{D}(\alpha)\otimes\hat{\mathcal{I}}\right)\hat{S}^{(2)}(r)\ket{0}\,,$
(60)

where $\hat{D}(\alpha)$ is a displacement operator acting on the first mode,
$\hat{S}^{(2)}(r)=e^{r(\hat{a}^{\dagger}\hat{b}^{\dagger}-\hat{a}\hat{b})}$ is
the two-mode squeezing operator with a real squeezing parameter $r$, and
$\hat{a}$ and $\hat{b}$ denote the annihilation operators of the first and
second mode, respectively. The HCRB for the estimation of the displacement
parameters only have been derived analytically in Ref. [16], and via semi-
definite programming in the case of mixed-state quantum statistical model in
Ref. [15]. Here we are interested in jointly estimating both the displacement
and squeezing, hence the vector of parameters $\bm{\theta}$ is still given by
Eq. (35). It is possible to manipulate the statistical model in order to make
the calculations easier and physical interpretation of the problem clearer. In
fact, the two-mode state Eq. (60) can actually be transformed into a tensor
product of single-mode states by means of a balanced beam splitter. The
justification to perform this unitary transformation lies in the fact that we
can always think of the latter as part of the detection scheme. We remind the
reader that the unitary evolution associated with a balanced beam splitter,
i.e. a beam splitter with trasmissivity $T=1/2$, is given by

$\hat{U}_{\textrm{BS}}=e^{\frac{\pi}{4}(\hat{a}^{\dagger}\hat{b}-\hat{a}\hat{b}^{\dagger})}\,.$
(61)

The evolution of the bosonic operators under the action of Eq. (61) is given
by

$\hat{U}^{\dagger}_{\textrm{BS}}\,\hat{a}\,\hat{U}_{\textrm{BS}}=\frac{\hat{a}+\hat{b}}{\sqrt{2}},\quad\hat{U}^{\dagger}_{\textrm{BS}}\,\hat{b}\,\hat{U}_{\textrm{BS}}=\frac{\hat{b}-\hat{a}}{\sqrt{2}}\,.$
(62)

Hence, the evolution of the state in Eq. (60) under $\hat{U}_{\textrm{BS}}$
reads

$\hat{U}_{\textrm{BS}}\ket{\psi_{\bm{\theta}}}=\hat{D}(\tfrac{\alpha}{\sqrt{2}})\hat{S}(r)\ket{0}\otimes\hat{D}(\tfrac{-\alpha}{\sqrt{2}})\hat{S}(-r)\ket{0}\,.$
(63)

As the ultimate bound on the estimation precision cannot depend on the
application of unitary operation on the original statistical model
$\ket{\psi_{\bm{\theta}}}$, we will now consider Eq. (63) to be our new
statistical model and, with a slight abuse of notation, we will still denote
it with $\ket{\psi_{\bm{\theta}}}$. We also introduce the following basis for
the two single-mode Hilbert spaces

$\displaystyle\ket{e_{n}}$
$\displaystyle=\hat{D}(\tfrac{\alpha}{\sqrt{2}})\hat{S}(r)\ket{n}\,,$ (64)
$\displaystyle\ket{f_{n}}$
$\displaystyle=\hat{D}(\tfrac{-\alpha}{\sqrt{2}})\hat{S}(-r)\ket{n}\,.$ (65)

Hence the statistical model simply reads
$\ket{\psi_{\bm{\theta}}}=\ket{e_{0}}\otimes\ket{f_{0}}\equiv\ket{e_{0}f_{0}}$.
We can now compute the SLD-CRB, the HCRB, and discuss their experimental
attainability. In Appendix B we outline the detailed calculation of the SLDs
and the QFI matrix. After some calculations, one finds

$\bm{Q}=\matrixquantity(4\cosh{(2r)}&0&0\\\ 0&4\cosh{(2r)}&0\\\ 0&0&4)\,.$
(66)

As expected from symmetry arguments, similarly to the single-mode case, the
QFI matrix does not depend on the displacement $\alpha$. The SLD-CRB is then
given by

$C^{S}=\Tr[\bm{Q}^{-1}]=\frac{1}{4}+\frac{1}{2\cosh{(2r)}}\,.$ (67)

Contrary to the single-mode scenario, in this setting squeezing is a useful
resource for our estimation protocol, i.e. the bigger the squeezing parameter,
the better the achievable overall estimation precision. Note that this is
already true when we are interested in the estimation of the displacement
only, as can be seen by the fact that $\bm{Q}_{33}$ does not depend on $r$. We
can understand this behaviour intuitively by giving a closer look at our
statistical model Eq. (63). Fist of all, we notice that both modes carry
complete information about the displacement $\alpha$, encoded in the unitary
operators $\hat{D}(\frac{\pm\alpha}{\sqrt{2}})$. We also know from the
previous section that single-mode squeezing allows for the estimation of the
real (imaginary) part of $\alpha$ with arbitrary precision, in the
$r\rightarrow-\infty$ ($r\rightarrow+\infty$) limit, at the cost of ignoring
the imaginary (real) part altogether. However, Eq. (63) is factorized into two
single-mode states which are squeezed in orthogonal directions in the phase
space. We can thus exploit structure this to optimally estimate
$\theta_{1}=\Re{\alpha}$ and $\theta_{2}=\Im{\alpha}$ simultaneously by
performing proper quadrature measurements on the two modes. The Uhlmann
curvature $\bm{D}$ reads

$\bm{D}=\matrixquantity(0&4&0\\\ -4&0&0\\\ 0&0&0)\,,$ (68)

which is equal to that of the single-mode scenario. Therefore the quantumness
of the statistical model is $R=1$, i.e. the maximum degree of incompatibility
between the parameters. We suspect that this is due to the fact that we are
dealing with a system at zero temperature, as seen in previous sections.

Following the techniques introduced in the previous section for the single-
mode case, we know that in order to compute the HCRB it is sufficient to
restrict the study of the Hermitian operators
$\hat{\bm{X}}=\\{\hat{X}_{1},\hat{X}_{2},\hat{X}_{3}\\}$ satisfying the
locally unbiased conditions Eq. (26) and Eq. (41), to the Hilbert subspace
spanned by the partial derivatives of the statistical model. We compute these
in Appendix B, i.e. Equations (95-97), hence it is clear that the Hilbert
subspace of interest is spanned by

$\begin{split}\\{\ket{e_{0}f_{0}}&,\ket{e_{0}f_{1}},\ket{e_{1}f_{0}},\ket{e_{0}f_{2}},\ket{e_{2}f_{0}}\\}\equiv\\\
&\\{\ket{\lambda_{1}},\ket{\lambda_{2}},\ket{\lambda_{3}},\ket{\lambda_{4}},\ket{\lambda_{5}}\\}\,.\end{split}$
(69)

With this new notation the statistical model reads
$\ket{\psi_{\bm{\theta}}}=\ket{\lambda_{1}}$, and its partial derivatives are
given by

$\displaystyle\ket{\partial_{1}\psi_{\bm{\theta}}}$
$\displaystyle=\frac{e^{-r}}{\sqrt{2}}\ket{\lambda_{3}}-i\theta_{2}\ket{\lambda_{1}}-\frac{e^{r}}{\sqrt{2}}\ket{\lambda_{2}}\,,$
(70) $\displaystyle\ket{\partial_{2}\psi_{\bm{\theta}}}$
$\displaystyle=\frac{ie^{r}}{\sqrt{2}}\ket{\lambda_{3}}+i\theta_{1}\ket{\lambda_{1}}-\frac{ie^{-r}}{\sqrt{2}}\ket{\lambda_{2}}\,,$
(71) $\displaystyle\ket{\partial_{3}\psi_{\bm{\theta}}}$
$\displaystyle=\frac{1}{\sqrt{2}}\ket{\lambda_{5}}-\frac{1}{\sqrt{2}}\ket{\lambda_{4}}\,.$
(72)

The constraint imposed by Eq. (26) implies

$\bra{\lambda_{1}}\hat{X}_{1}\ket{\lambda_{1}}=\bra{\lambda_{1}}\hat{X}_{2}\ket{\lambda_{1}}=\bra{\lambda_{1}}\hat{X}_{3}\ket{\lambda_{1}}=0\,.$
(73)

In Appendix C we display the additional constraints on each Hermitian operator
$\hat{X}_{k}$ imposed by Eq. (41). We end up with the vector of operators
$\hat{\bm{X}}$ parameterized by fifteen free real variables. One can then work
out the three by three Hermitian matrix $\bm{Z}_{\bm{\theta}}[\hat{\bm{X}}]$
and use the latter to compute the function $h_{\bm{\theta}}[\hat{\bm{X}}]$,
whose minimization yields the HCRB. We do not report the analytic expression
of $h_{\bm{\theta}}[\hat{\bm{X}}]$ for brevity and also because the
optimization is now non-trivial and must be carried out with numerical
methods. Before displaying the results of such minimization, we want to
further discuss how to exploit the fact that the statistical model can be
expressed - by means of a balanced beam-splitter - as a tensor product of two
single-mode states squeezed in orthogonal directions in the phase space. As we
have hinted before, this suggests that performing quadrature measurements on
the two modes would allow us to optimally estimate the displacement. In fact,
it has been shown in Ref. [15] that a double-homodyne measurement achieves
optimality whenever the quantum state is entangled, if we are interested in
the estimation of the displacement $\alpha$ only. Hence, it seems reasonable
to choose this detection scheme as the starting point of our analysis for the
three-parameter estimation case. The covariance matrix $\bm{\sigma}$ of the
statistical model Eq. (63) is simply the direct sum of the covariance matrices
associated with single-mode squeezed vacuum states, namely

$\bm{\sigma}=\text{diag}(e^{2r},e^{-2r},e^{-2r},e^{2r}),$ (74)

while the statistical model’s vector of first moments reads

$\overline{\bm{r}}=(\theta_{1},\theta_{2},-\theta_{1},-\theta_{2})^{\intercal}\,.$
(75)

The double-homodyne detection scheme consists in performing a $\hat{p}$
quadrature measurement on the first mode, and a $\hat{q}$ quadrature
measurement on the second one (we should remind ourselves that these
measurements are performed after a balanced beam-splitter, if we consider the
state in Eq. (60)) as our initial state). Consequently, the covariance matrix
$\bm{\sigma}_{m}$ that characterizes this POVM reads

$\bm{\sigma}_{m}=\lim_{z\rightarrow 0}\text{diag}(z^{-1},z,z,z^{-1})\,.$ (76)

Hence, the corresponding conditional probability distribution of outcomes is a
Gaussian centered in $\overline{\bm{x}}=\overline{\bm{r}}$, with covariance
matrix given by $\bm{\Sigma}=(\bm{\sigma}+\bm{\sigma}_{m})/2$. By taking the
limit $z\to 0$ after the matrix inversion, one can prove that

$(\bm{\sigma}+\bm{\sigma}_{m})^{-1}=\text{diag}(0,e^{2r},e^{2r},0)\,.$ (77)

The FI matrix can then be evaluated using Eq. (17), yielding

$\bm{F}=\matrixquantity(2e^{2r}&0&0\\\ 0&2e^{2r}&0\\\ 0&0&4)\,.$ (78)

The corresponding estimation precision of the double-homodyne detection scheme
thus reads

$\Tr[\bm{F}^{-1}]=\frac{1}{4}+e^{-2r}\,.$ (79)

It turns out that this bound coincides with the HCRB we obtained from
numerical optimization, as it can be seen in Fig. (3).

Figure 3: The dashed orange line denotes the SLD-CRB Eq. (67), while the solid
blue line denotes the estimation precision achieved with the double-homodyne
detection scheme Eq. (79). The dots are the result of the numerical
optimization that yields the HCRB. The plot shows that $C^{H}$ coincides with
the estimation precision obtained with the double-homodyne detection scheme,
thus making the latter the optimal measurement for the joint estimation
problem at study.

Hence, we may infer that the expression of the HCRB is given by

$C^{H}=\frac{1}{4}+e^{-2r}$ (80)

and that the double-homodyne detection scheme represents the optimal POVM for
our estimation problem, for every value $r$. Also in this case, the AI measure
$R$ largely overestimates the true gap between the HCRB and the SLD-CRB bound.
We recall that, given the statistical model Eq. (63), the double-homodyne
represents the optimal measurement for the estimation of the displacement
only. Here we have shown that from the measurement outcomes of that same POVM
one can _optimally_ estimate the squeezing parameter as well, without
incurring in any additional cost.

##  VI Conclusions

We have analyzed in details three-parameter estimation problems involving pure
Gaussian states, focusing on the joint estimation of complex displacement and
real squeezing, encoded unitarily in single- and two-mode pure states. We have
explicitly evaluated the Holevo-Cramér-Rao bound (HCRB) and explored its
dependence on the degree of squeezing.

Our findings reveal that for single-mode states, the HCRB deteriorates as
squeezing increases. In scenarios of substantial squeezing, heterodyne
detection approaches near-optimality. However, our analysis reveals that, more
broadly, achieving the optimal measurement requires non-Gaussian measurements.
Conversely, in the two-mode scenario, we discover that squeezing—or
equivalently, entanglement—serves to enhance the HCRB. Remarkably, this
enhanced bound is attainable through double-homodyne detection across all
levels of squeezing, showcasing the intrinsic value of entanglement in quantum
metrology.

Our results reaffirms the role of entanglement as a resource in quantum
metrology, setting the stage for advanced quantum sensing technologies. By
clarifying the relationship between squeezing, measurement optimality, and
quantum state estimation, we pave the way to further improve metrological
precision with realistic and reliable detection schemes.

##  VII Acknowledgments

MGG acknowledges useful discussions with Francesco Albarelli. GB is part of
the AppQInfo MSCA ITN which received funding from the European Union’s Horizon
2020 research and innovation programme under the Marie Sklodowska-Curie grant
agreement No 956071. MGG has been partially supported by MUR through Project
No. PRIN22-2022KB2JJM-CONTRABASS. The work of MGAP has been done under the
auspices of GNFM-INdAM and has been partially supported by MUR through Project
No. PRIN22-2022T25TR3-RISQUE.

## References

  * [1] Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone.  Quantum-enhanced measurements: Beating the standard quantum limit.  Science, 306(5700):1330–1336, 2004.
  * [2] Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone.  Quantum Metrology.  Phys. Rev. Lett., 96(1):010401, jan 2006.
  * [3] Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone.  Advances in quantum metrology.  Nat. Photonics, 5(4):222–229, apr 2011.
  * [4] Carlton M. Caves.  Quantum-mechanical noise in an interferometer.  Phys. Rev. D, 23:1693–1708, Apr 1981.
  * [5] Matteo G. A. Paris.  Quantum estimation for quantum technology.  Int. J. Quantum Inf., 07(supp01):125–137, jan 2009.
  * [6] Rafał Demkowicz-Dobrzański, Marcin Jarzyna, and Jan Kołodyński.  Quantum Limits in Optical Interferometry.  In Emil Wolf, editor, Prog. Opt. Vol. 60, chapter 4, pages 345–435. Elsevier, Amsterdam, 2015.
  * [7] C. L. Degen, F Reinhard, and Paola Cappellaro.  Quantum sensing.  Rev. Mod. Phys., 89(3):035002, jul 2017.
  * [8] Stefano Pirandola, Bhaskar Roy Bardhan, Tobias Gehring, Christian Weedbrook, and Seth Lloyd.  Advances in photonic quantum sensing.  Nat. Photonics, 12(12):724–733, dec 2018.
  * [9] Marco Barbieri.  Optical quantum metrology.  PRX Quantum, 3:010202, Jan 2022.
  * [10] F. Albarelli, M. Barbieri, M.G. Genoni, and I. Gianani.  A perspective on multiparameter quantum metrology: From theoretical tools to applications in quantum imaging.  Physics Letters A, 384(12):126311, 2020.
  * [11] Magdalena Szczykulska, Tillmann Baumgratz, and Animesh Datta.  Multi-parameter quantum metrology.  Adv. Phys. X, 1(4):621–639, jul 2016.
  * [12] Jing Liu, Haidong Yuan, Xiao-Ming Lu, and Xiaoguang Wang.  Quantum fisher information matrix and multiparameter estimation.  Journal of Physics A: Mathematical and Theoretical, 53(2):023001, dec 2019.
  * [13] Rafał Demkowicz-Dobrzański, Wojciech Górecki, and Mădălin Guţă.  Multi-parameter estimation beyond quantum fisher information.  Journal of Physics A: Mathematical and Theoretical, 53(36):363001, aug 2020.
  * [14] Marco G. Genoni, Matteo G. A. Paris, Gerardo Adesso, Hyunchul Nha, Peter L. Knight, and M. S. Kim.  Optimal estimation of joint parameters in phase space.  Phys. Rev. A, 87(1):012107, jan 2013.
  * [15] Mark Bradshaw, Ping Koy Lam, and Syed M. Assad.  Ultimate precision of joint quadrature parameter estimation with a Gaussian probe.  Phys. Rev. A, 97(1):012106, jan 2018.
  * [16] Mark Bradshaw, Syed M. Assad, and Ping Koy Lam.  A tight Cramér–Rao bound for joint parameter estimation with a pure two-mode squeezed probe.  Phys. Lett. A, 381(32):2598–2607, aug 2017.
  * [17] Peter C. Humphreys, Marco Barbieri, Animesh Datta, and Ian A. Walmsley.  Quantum Enhanced Multiple Phase Estimation.  Phys. Rev. Lett., 111(7):070403, 2013.
  * [18] Christos N. Gagatsos, Dominic Branford, and Animesh Datta.  Gaussian systems for quantum-enhanced multiple phase estimation.  Phys. Rev. A, 94(4):042342, oct 2016.
  * [19] Paul A. Knott, Timothy J. Proctor, Anthony J. Hayes, Jason F. Ralph, Pieter Kok, and Jacob A. Dunningham.  Local versus global strategies in multiparameter estimation.  Phys. Rev. A, 94(6):062312, dec 2016.
  * [20] Luca Pezzè, Mario A. Ciampini, Nicolò Spagnolo, Peter C. Humphreys, Animesh Datta, Ian A. Walmsley, Marco Barbieri, Fabio Sciarrino, and Augusto Smerzi.  Optimal Measurements for Simultaneous Quantum Estimation of Multiple Phases.  Phys. Rev. Lett., 119(13):130504, sep 2017.
  * [21] Lorcán O. Conlon, Jun Suzuki, Ping Koy Lam, and Syed M. Assad.  Efficient computation of the Nagaoka–Hayashi bound for multiparameter estimation with separable measurements.  npj Quantum Inf, 7(1):1–8, July 2021.  Bandiera_abtest: a Cc_license_type: cc_by Cg_type: Nature Research Journals Number: 1 Primary_atype: Research Publisher: Nature Publishing Group Subject_term: Quantum information;Quantum metrology Subject_term_id: quantum-information;quantum-metrology.
  * [22] Lorcán O. Conlon, Ping Koy Lam, and Syed M. Assad.  Multiparameter estimation with two-qubit probes in noisy channels.  Entropy, 25(8), 2023.
  * [23] G. Mauro D’Ariano, Matteo G. A. Paris, and Massimiliano F. Sacchi.  Parameter estimation in quantum optics.  Phys. Rev. A, 62:023815, Jul 2000.
  * [24] Sergey I. Knysh and Gabriel A. Durkin.  Estimation of phase and diffusion: Combining quantum statistics and classical noise, 2013.
  * [25] Mihai D. Vidrighin, Gaia Donati, Marco G. Genoni, Xian-Min Jin, W. Steven Kolthammer, M. S. Kim, Animesh Datta, Marco Barbieri, and Ian A. Walmsley.  Joint estimation of phase and phase diffusion for quantum metrology.  Nature Communications, 5(1):3532, 2014.
  * [26] Matteo Altorio, Marco G. Genoni, Mihai D. Vidrighin, Fabrizia Somma, and Marco Barbieri.  Weak measurements and the joint estimation of phase and phase diffusion.  Phys. Rev. A, 92:032114, Sep 2015.
  * [27] Magdalena Szczykulska, Tillmann Baumgratz, and Animesh Datta.  Reaching for the quantum limits in the simultaneous estimation of phase and phase diffusion.  Quantum Science and Technology, 2(4):044004, aug 2017.
  * [28] Emanuele Roccia, Ilaria Gianani, Luca Mancino, Marco Sbroscia, Fabrizia Somma, Marco G Genoni, and Marco Barbieri.  Entangling measurements for multiparameter estimation with two qubits.  Quantum Science and Technology, 3(1):01LT01, oct 2017.
  * [29] Mankei Tsang, Ranjith Nair, and Xiao-Ming Lu.  Quantum theory of superresolution for two incoherent optical point sources.  Phys. Rev. X, 6:031033, Aug 2016.
  * [30] Andrzej Chrostowski, Rafał Demkowicz-Dobrzański, Marcin Jarzyna, and Konrad Banaszek.  On super-resolution imaging as a multiparameter estimation problem.  International Journal of Quantum Information, 15(08):1740005, 2017.
  * [31] J. Rehacek, M. Paúr, B. Stoklasa, Z. Hradil, and L. L. Sánchez-Soto.  Optimal measurements for resolution beyond the rayleigh limit.  Opt. Lett., 42(2):231–234, Jan 2017.
  * [32] J. Řeháček, Z. Hradil, D. Koutný, J. Grover, A. Krzic, and L. L. Sánchez-Soto.  Optimal measurements for quantum spatial superresolution.  Phys. Rev. A, 98:012103, Jul 2018.
  * [33] Carmine Napoli, Samanta Piano, Richard Leach, Gerardo Adesso, and Tommaso Tufarelli.  Towards superresolution surface metrology: Quantum estimation of angular and axial separations.  Phys. Rev. Lett., 122:140505, Apr 2019.
  * [34] Lukas J. Fiderer, Tommaso Tufarelli, Samanta Piano, and Gerardo Adesso.  General expressions for the quantum fisher information matrix with applications to discrete quantum imaging.  PRX Quantum, 2:020308, Apr 2021.
  * [35] Vahid Ansari, Benjamin Brecht, Jano Gil-Lopez, John M. Donohue, Jaroslav Řeháček, Zden ěk Hradil, Luis L. Sánchez-Soto, and Christine Silberhorn.  Achieving the ultimate quantum timing resolution.  PRX Quantum, 2:010301, Jan 2021.
  * [36] Carl W. Helstrom.  Quantum detection and estimation theory.  Academic Press, New York, 1976.
  * [37] Alexander S. Holevo.  Probabilistic and Statistical Aspects of Quantum Theory.  Edizioni della Normale, Pisa, 2nd edition, 2011.
  * [38] Angelo Carollo, Bernardo Spagnolo, Alexander A Dubkov, and Davide Valenti.  On quantumness in multi-parameter quantum estimation.  J. Stat. Mech. Theory Exp., 2019(9):094010, sep 2019.
  * [39] Sholeh Razavian, Matteo G. A. Paris, and Marco G. Genoni.  On the quantumness of multiparameter estimation problems for qubit systems.  Entropy, 22(11), 2020.
  * [40] Alessandro Candeloro, Matteo G A Paris, and Marco G Genoni.  On the properties of the asymptotic incompatibility measure in multiparameter quantum estimation.  Journal of Physics A: Mathematical and Theoretical, 54(48):485301, nov 2021.
  * [41] Mankei Tsang, Francesco Albarelli, and Animesh Datta.  Quantum semiparametric estimation.  Phys. Rev. X, 10:031023, Jul 2020.
  * [42] Hiroshi Nagaoka.  A New Approach to Cramér-Rao Bounds for Quantum State Estimation, pages 100–112.  2005\.
  * [43] Masahito Hayashi.  A Linear Programming Approach to Attainable Cramésr-Rao Type Bounds, pages 150–161.  2005\.
  * [44] K Matsumoto.  A new approach to the cramér-rao-type bound of the pure-state model.  Journal of Physics A: Mathematical and General, 35(13):3111–3123, mar 2002.
  * [45] Francesco Albarelli, Jamie F. Friel, and Animesh Datta.  Evaluating the holevo cramér-rao bound for multiparameter quantum metrology.  Phys. Rev. Lett., 123:200503, Nov 2019.
  * [46] Wojciech Górecki, Sisi Zhou, Liang Jiang, and Rafał Demkowicz-Dobrzański.  Optimal probes and error-correction schemes in multi-parameter quantum metrology.  Quantum, 4:288, July 2020.
  * [47] Jamie Friel, Pantita Palittapongarnpim, Francesco Albarelli, and Animesh Datta.  Attainability of the holevo-cramér-rao bound for two-qubit 3d magnetometry, 2020.
  * [48] Jun Suzuki.  Explicit formula for the Holevo bound for two-parameter qubit-state estimation problem.  J. Math. Phys., 57(4):042201, apr 2016.
  * [49] Marcin Jarzyna.  Quantum limits to polarization measurement of classical light, 2021.
  * [50] Marco G. Genoni, Ludovico Lami, and Alessio Serafini.  Conditional and unconditional Gaussian quantum dynamics.  Contemp. Phys., 57(3):331–349, July 2016.
  * [51] Alessio Serafini.  Quantum Continuous Variables : A Primer of Theoretical Methods.  CRC Press, Boca Raton, 2017.
  * [52] A. Ferraro, S. Olivares, and M. G. A. Paris.  Gaussian states in Quantum Information.  Bibliopolis, Napoli, 2005.
  * [53] Changhun Oh, Changhyoup Lee, Carsten Rockstuhl, Hyunseok Jeong, Jaewan Kim, Hyunchul Nha, and Su-Yong Lee.  Optimal gaussian measurements for phase estimation in single-mode gaussian metrology.  npj Quantum Information, 5(1):10, 2019.

##  Appendix A QFI matrix evaluation for single-mode Gaussian quantum state

In this appendix we derive the QFI matrix for the single-mode quantum
statistical model Eq. (34). We remind the reader that for pure statistical
models the SLDs can be expressed as in Eq. (19). Hence, in order to compute
$\hat{L}_{\mu}$ we just need to differentiate the statistical mode with
respect to the parameter to estimated $\theta_{\mu}$. To this end we first
compute the derivatives of the displacement operator and the squeezing
operator with respect to the parameters to estimate, namely

$\displaystyle\partial_{1}\hat{D}(\alpha)$
$\displaystyle=(\hat{a}^{\dagger}-\hat{a}+i\theta_{2})\hat{D}(\alpha)\,,$ (81)
$\displaystyle\partial_{2}\hat{D}(\alpha)$
$\displaystyle=i(\hat{a}^{\dagger}+\hat{a}-\theta_{1})\hat{D}(\alpha)\,,$ (82)
$\displaystyle\partial_{3}\hat{S}(r)$ $\displaystyle=\frac{\hat{a}^{\dagger
2}-\hat{a}^{2}}{2}\hat{S}(r)\,.$ (83)

Exploiting these and the identities
$\hat{D}^{\dagger}(\alpha)\hat{a}\hat{D}(\alpha)=\hat{a}+\alpha$ and
$\hat{S}^{\dagger}(r)\hat{a}\hat{S}(r)=\mu\hat{a}+\nu\hat{a}^{\dagger}$, where
$\mu=\cosh{r}$ and $\nu=\sinh{r}$, we finally obtain

$\displaystyle\ket{\partial_{1}\psi_{\bm{\theta}}}$
$\displaystyle=-i\theta_{2}\ket{e_{0}}+e^{-r}\ket{e_{1}}\,,$ (84)
$\displaystyle\ket{\partial_{2}\psi_{\bm{\theta}}}$
$\displaystyle=i\theta_{1}\ket{e_{0}}+ie^{r}\ket{e_{1}}\,,$ (85)
$\displaystyle\ket{\partial_{3}\psi_{\bm{\theta}}}$
$\displaystyle=\frac{1}{\sqrt{2}}\ket{e_{2}}\,.$ (86)

Substituting these expressions into Eq. (19) yields the SLDs, i.e.

$\displaystyle\hat{L}_{1}$
$\displaystyle=2e^{-r}\left(\ket{e_{1}}\\!\\!\bra{e_{0}}+\ket{e_{0}}\\!\\!\bra{e_{1}}\right)\,,$
(87) $\displaystyle\hat{L}_{2}$
$\displaystyle=2ie^{r}\left(\ket{e_{1}}\\!\\!\bra{e_{0}}-\ket{e_{0}}\\!\\!\bra{e_{1}}\right)\,,$
(88) $\displaystyle\hat{L}_{3}$
$\displaystyle=\sqrt{2}\left(\ket{e_{2}}\\!\\!\bra{e_{0}}+\ket{e_{0}}\\!\\!\bra{e_{2}}\right)\,.$
(89)

We can then evaluate the QFI matrix, whose expression for our pure statistical
model $\ket{\psi_{\bm{\theta}}}=\ket{e_{0}}$ reads

$\bm{Q}_{\mu\nu}=\frac{1}{2}\bra{e_{0}}\hat{L}_{\mu}\hat{L}_{\nu}+\hat{L}_{\nu}\hat{L}_{\mu}\ket{e_{0}}\,.$
(90)

After some calculations we obtain the following diagonal matrix

$\bm{Q}=\matrixquantity(4e^{-2r}&0&0\\\ 0&4e^{2r}&0\\\ 0&0&2)\,.$ (91)

##  Appendix B QFI matrix evaluation for single-mode Gaussian quantum state

In this appendix we derive the QFI matrix for the two-mode quantum statistical
model in Eq. (63). Let us first work out the derivatives of
$\hat{D}(\tfrac{\pm\alpha}{\sqrt{2}})$ and $\hat{S}(\pm r)$ with respect to
the parameters to be estimated $\theta_{\mu}$, namely

$\displaystyle\partial_{1}\hat{D}(\tfrac{\pm\alpha}{\sqrt{2}})$
$\displaystyle=\left(\frac{\pm\hat{a}^{\dagger}\mp\hat{a}}{\sqrt{2}}+\frac{i\theta_{2}}{2}\right)\hat{D}(\tfrac{\pm\alpha}{\sqrt{2}})\,,$
(92) $\displaystyle\partial_{2}\hat{D}(\tfrac{\pm\alpha}{\sqrt{2}})$
$\displaystyle=i\left(\pm\frac{\hat{a}+\hat{a}^{\dagger}}{\sqrt{2}}-\frac{\theta_{1}}{2}\right)\hat{D}(\tfrac{\pm\alpha}{\sqrt{2}})\,,$
(93) $\displaystyle\partial_{3}\hat{S}(\pm r)$
$\displaystyle=\pm\left(\frac{\hat{a}^{\dagger
2}-\hat{a}^{2}}{2}\right)\hat{S}(\pm r)\,.$ (94)

Using these expressions, together with the previously introduced
transformation properties of the mode operators under displacement and local
squeezing unitaries, we cam compute the derivative of our statistical model.
After cumbersome calculations one obtains

$\displaystyle\ket{\partial_{1}\psi_{\bm{\theta}}}$
$\displaystyle=\frac{e^{-r}}{\sqrt{2}}\ket{e_{1}f_{0}}-i\theta_{2}\ket{e_{0}f_{0}}-\frac{e^{r}}{\sqrt{2}}\ket{e_{0}f_{1}}\,,$
(95) $\displaystyle\ket{\partial_{2}\psi_{\bm{\theta}}}$
$\displaystyle=\frac{ie^{r}}{\sqrt{2}}\ket{e_{1}f_{0}}+i\theta_{1}\ket{e_{0}f_{0}}-\frac{ie^{-r}}{\sqrt{2}}\ket{e_{0}f_{1}}\,,$
(96) $\displaystyle\ket{\partial_{3}\psi_{\bm{\theta}}}$
$\displaystyle=\frac{1}{\sqrt{2}}\ket{e_{2}f_{0}}-\frac{1}{\sqrt{2}}\ket{e_{0}f_{2}}\,.$
(97)

Substituting these expressions into Eq. (19) yields the SLDs

$\displaystyle\hat{L}_{1}$
$\displaystyle=\sqrt{2}\left(e^{-r}\ket{e_{1}f_{0}}\\!\\!\bra{e_{0}f_{0}}-e^{r}\ket{e_{0}f_{1}}\\!\\!\bra{e_{0}f_{0}}+h.c.\right)\,,$
(98) $\displaystyle\hat{L}_{2}$
$\displaystyle=i\sqrt{2}\left(e^{r}\ket{e_{1}f_{0}}\\!\\!\bra{e_{0}f_{0}}-e^{-r}\ket{e_{0}f_{1}}\\!\\!\bra{e_{0}f_{0}}-h.c.\right)\,,$
(99) $\displaystyle\hat{L}_{3}$
$\displaystyle=\sqrt{2}\left(\ket{e_{2}f_{0}}\\!\\!\bra{e_{0}f_{0}}-\ket{e_{0}f_{2}}\\!\\!\bra{e_{0}f_{0}}+h.c.\right)\,.$
(100)

We can then evaluate the QFI matrix $\bm{Q}$, whose matrix elements are given
by

$\bm{Q}_{\mu\nu}=\frac{1}{2}\bra{e_{0}f_{0}}\hat{L}_{\mu}\hat{L}_{\nu}+\hat{L}_{\nu}\hat{L}_{\mu}\ket{e_{0}f_{0}}\,.$
(101)

After some calculations we obtain the following diagonal matrix

$\bm{Q}=\matrixquantity(4\cosh{(2r)}&0&0\\\ 0&4\cosh{(2r)}&0\\\ 0&0&4)\,.$
(102)

##  Appendix C HCRB evaluation for two-mode displaced squeezed vacuum states

In this appendix we outline the computation of the HCRB for the two-mode
quantum statistical model Eq. (63). The locally unbiased condition Eq. (41)
imposes three additional constraints for each Hermitian operator
$\hat{X}_{\mu}$. In particular, one finds that $\hat{X}_{1}$ must satisfy

$\displaystyle\sqrt{2}\left[e^{-r}\Re{\bra{\lambda_{1}}\hat{X}_{1}\ket{\lambda_{3}}}-e^{r}\Re{\bra{\lambda_{1}}\hat{X}_{1}\ket{\lambda_{2}}}\right]=1\,,$
$\displaystyle-\sqrt{2}\left[e^{r}\Im{\bra{\lambda_{1}}\hat{X}_{1}\ket{\lambda_{3}}}-e^{-r}\Im{\bra{\lambda_{1}}\hat{X}_{1}\ket{\lambda_{2}}}\right]=0\,,$
$\displaystyle\sqrt{2}\left[\Re{\bra{\lambda_{1}}\hat{X}_{1}\ket{\lambda_{5}}}-\Re{\bra{\lambda_{1}}\hat{X}_{1}\ket{\lambda_{4}}}\right]=0\,.$

These conditions in turn imply that

$\displaystyle\bra{\lambda_{1}}\hat{X}_{1}\ket{\lambda_{2}}$
$\displaystyle=x_{1}+ie^{2r}x_{2}\,,$
$\displaystyle\bra{\lambda_{1}}\hat{X}_{1}\ket{\lambda_{3}}$
$\displaystyle=\frac{e^{r}}{\sqrt{2}}+e^{2r}x_{1}+ix_{2}\,,$
$\displaystyle\bra{\lambda_{1}}\hat{X}_{1}\ket{\lambda_{4}}$
$\displaystyle=x_{3}+ix_{4}\,,$
$\displaystyle\bra{\lambda_{1}}\hat{X}_{1}\ket{\lambda_{5}}$
$\displaystyle=x_{3}+ix_{5}\,,$

where $\\{x_{1},x_{2},x_{3},x_{4},x_{5}\\}$ are five free real parameters.
Analogously, the additional constraints on $\hat{X}_{2}$ read

$\displaystyle\sqrt{2}\left[e^{-r}\Re{\bra{\lambda_{1}}\hat{X}_{2}\ket{\lambda_{3}}}-e^{r}\Re{\bra{\lambda_{1}}\hat{X}_{2}\ket{\lambda_{2}}}\right]=0\,,$
$\displaystyle-\sqrt{2}\left[e^{r}\Im{\bra{\lambda_{1}}\hat{X}_{2}\ket{\lambda_{3}}}-e^{-r}\Im{\bra{\lambda_{1}}\hat{X}_{2}\ket{\lambda_{2}}}\right]=1\,,$
$\displaystyle\sqrt{2}\left[\Re{\bra{\lambda_{1}}\hat{X}_{2}\ket{\lambda_{5}}}-\Re{\bra{\lambda_{1}}\hat{X}_{2}\ket{\lambda_{4}}}\right]=0\,.$

These are equivalent to

$\displaystyle\bra{\lambda_{1}}\hat{X}_{2}\ket{\lambda_{2}}$
$\displaystyle=x_{6}+i\left(\frac{e^{r}}{\sqrt{2}}+e^{2r}x_{6}\right)\,,$
$\displaystyle\bra{\lambda_{1}}\hat{X}_{2}\ket{\lambda_{3}}$
$\displaystyle=e^{2r}x_{6}+ix_{7}\,,$
$\displaystyle\bra{\lambda_{1}}\hat{X}_{2}\ket{\lambda_{4}}$
$\displaystyle=x_{8}+ix_{9}\,,$
$\displaystyle\bra{\lambda_{1}}\hat{X}_{2}\ket{\lambda_{5}}$
$\displaystyle=x_{8}+ix_{10}\,,$

with $\\{x_{6},x_{7},x_{8},x_{9},x_{10}\\}$ five other free real parameters.
Finally, the constraints on $\hat{X}_{3}$ are

$\displaystyle\sqrt{2}\left[e^{-r}\Re{\bra{\lambda_{1}}\hat{X}_{3}\ket{\lambda_{3}}}-e^{r}\Re{\bra{\lambda_{1}}\hat{X}_{3}\ket{\lambda_{2}}}\right]=0\,,$
$\displaystyle-\sqrt{2}\left[e^{r}\Im{\bra{\lambda_{1}}\hat{X}_{3}\ket{\lambda_{3}}}-e^{-r}\Im{\bra{\lambda_{1}}\hat{X}_{3}\ket{\lambda_{2}}}\right]=0\,,$
$\displaystyle\sqrt{2}\left[\Re{\bra{\lambda_{1}}\hat{X}_{3}\ket{\lambda_{5}}}-\Re{\bra{\lambda_{1}}\hat{X}_{3}\ket{\lambda_{4}}}\right]=1\,.$

Hence, the matrix elements of $\hat{X}_{3}$ must satisfy

$\displaystyle\bra{\lambda_{1}}\hat{X}_{3}\ket{\lambda_{2}}$
$\displaystyle=x_{11}+ie^{2r}x_{12}\,,$
$\displaystyle\bra{\lambda_{1}}\hat{X}_{3}\ket{\lambda_{3}}$
$\displaystyle=e^{2r}x_{11}+ix_{12}\,,$
$\displaystyle\bra{\lambda_{1}}\hat{X}_{3}\ket{\lambda_{4}}$
$\displaystyle=x_{13}+ix_{14}\,,$
$\displaystyle\bra{\lambda_{1}}\hat{X}_{3}\ket{\lambda_{5}}$
$\displaystyle=x_{13}+\frac{1}{\sqrt{2}}+ix_{15}\,,$

where $\\{x_{11},x_{12},x_{13},x_{14},x_{15}\\}$ are free real parameters. The
hermiticity of these operators imposes additional constraints on the matrix
elements and, as explained Section IV of the main text, we can set to zero all
other matrix elements of $\hat{X}_{\mu}$ not involved in the above-mentioned
constraints. We therefore have a total of fifteen free real variables.