Quantum parameter estimation using multi-mode Gaussian states

One of the main aims of quantum metrology is to find the ultimate precision bound on the estimation of a physical parameter encoded in a quantum state. Of special interest are parameters that cannot be measured directly, since they do not correspond to observables of the system. However, they can be estimated by finding an appropriate measurement strategy. The estimation also involves choosing an estimator $\hat{\epsilon }$ which maps the set of the measurement results onto the set of possible parameters. The ultimate precision limit is given by the quantum Cramér–Rao bound [1, 2] which gives a lower bound on the mean squared error of any locally unbiased estimator $\hat{\epsilon }$. The local unbiasedness means that in the limit where the number of measurements goes to infinity, the value of the estimator converges to the real value of the parameter. The bound is given by the number of measurements taken on the identical copies of the state $\rho (\epsilon )$ and a quantity $H(\epsilon )$ called the quantum Fisher information. Higher precision is achieved by increasing the number of measurements and maximizing the quantum Fisher information. Calculating the quantum Fisher information thus gives us an idea of how well we can estimate the parameter when only a fixed amount of measurements are available. This technique has been applied, for example, in large interferometers like VIRGO [3] and LIGO [4] assigned to measure gravitational waves, or a current proposal [5] of measuring gravitation waves using phonons in Bose–Einstein condensates, magnetometers [6, 7], and gravimeters [8].

Calculating the quantum Fisher information is not always an easy task. Although a general formula for the quantum Fisher information exists, it is written in a terms of the density matrix [2]. On the other hand, many applications use a special kind of a continuous-variable systems called Gaussian states, for which the description using density matrices seems particularly ineffective. Gaussian states can be conveniently described in terms of the first and the second moments of the so-called quadrature operators. This description is usually called a phase-space or the covariance matrix formalism [9].

Despite of the importance and practical usage of the quantum Fisher information, the theory for estimation using Gaussian states in the phase-space formalism is far from complete, and only partial results are known. The first leap in deriving general formulae has been taken by Pinel et al [10], who found a formula for pure states, i.e., for the states which are pure at point and remain pure even if the slightly changes. The same year Marian and Marian found the formula for the fidelity between one-mode and two-mode Gaussian states [11], which allowed for the derivation of the general formula for the one-mode state [12]. Also, Spedalieri et al found a formula for the fidelity between one pure and one mixed Gaussian state [13], from which one can derive a slightly more general formula for pure states, i.e., for the states which are pure at the point but the small change in introduces impurity. A different path has been followed by Monras [14], who connected the quantum Fisher information to the solution of the so-called Stein-equation. Using this approach, he derived the quantum Fisher information for a generalization of pure states called iso-thermal states, and a general formula for any multi-mode Gaussian state in terms of an infinite series. Using the previous result, Jiang derived a formula [15] for Gaussian states in exponential form and simplified a known formula for pure states. Quite recently, Gao and Lee derived an exact formula [16] for the quantum Fisher information for the multi-mode Gaussian states in terms of the inverse of certain tensor products, elegantly generalizing the previous results, however with some possible drawbacks, especially in the necessity of inverting relatively large matrices.

In this article, we first introduce a phase-space description of the Gaussian states. Then we derive an exact formula for two-mode Gaussian states in analogy with [12]. We simplify the result of Monras for multi-mode Gaussian states and we show that the infinite sum involved converges as a geometrical series. However, the series may not be easy to evaluate, and for that reason we calculate the error when only a finite number of terms of the sum are taken into account. On the other hand, in the case when the Williamson decomposition of the covariance matrix is known, the infinite series can be evaluated. This gives a general and exact formula for the quantum Fisher information of any multi-mode Gaussian state in terms of its decomposition. Finally, we use the derived formula to find the optimal Gaussian states for estimating the squeezing parameter of a single mode squeezing channel, and we demonstrate that this estimation can be exponentially enhanced by the initial squeezing when followed by an appropriate rotation. In appendix A we connect the so-called real and the complex form of the covariance matrix, in appendix B we study the case of pure states and the regularization procedure which allows us to use our results for generally mixed states, even for states where some or all modes remain pure. Appendices C and D contain detailed proofs for some results and in appendix E we show a full derivation of our single mode squeezing channel example. In appendix F we provide a table of frequently used notation that are constant throughout the paper and appear repeatedly.

1.1. Phase-space formalism of Bosonic modes and the Williamson decomposition

In this section we recapitulate the phase-space description of a Bosonic system, which will be particularly useful for the continuous parameter states known as Gaussian states. First, let us consider a Bosonic system with the set of annihilation and creation operators $\{{\hat{a}}_{n},{\hat{a}}_{n}^{\dagger }\}$. We collect them into a vector of operators $\hat{{\boldsymbol{A}}}:= {({\hat{a}}_{1},{\hat{a}}_{2},...,{\hat{a}}_{1}^{\dagger },{\hat{a}}_{2}^{\dagger },...)}^{T}$. Now we can write the commutation relation between these operators in a compact form

$$\begin{eqnarray}[{\hat{{\boldsymbol{A}}}}_{i},{\hat{{\boldsymbol{A}}}}_{j}^{\dagger }]={K}_{{ij}}\mathrm{id}\quad \Rightarrow \quad K=\left[\begin{array}{cc}I & 0\\ 0 & -I\end{array}\right],\end{eqnarray} \tag{ 1 }$$

where $\mathrm{id}$ denotes the identity element of an algebra and I is the identity matrix. Note that ${K}^{-1}={K}^{\dagger }=K$ and ${K}^{2}=I$.

One way of representing a state in quantum mechanics is by using a density matrix $\hat{\rho }$, however, for Bosonic systems an alternative and completely equivalent description exists, which is particularly useful in a description of Gaussian states. Given a state $\hat{\rho }$ we define the symmetric characteristic function as

$$\begin{eqnarray}&&\chi ({\boldsymbol{\xi }})=\mathrm{tr}[\hat{\rho }\;\hat{D}({\boldsymbol{\xi }})],\end{eqnarray} \tag{ 2 }$$

where $\hat{D}({\boldsymbol{\xi }})={{\rm{e}}}^{{\hat{{\boldsymbol{A}}}}^{\dagger }K{\boldsymbol{\xi }}}$ is the Weyl displacement operator with the variable of the form ${\boldsymbol{\xi }}={\boldsymbol{\gamma }}\oplus \bar{{\boldsymbol{\gamma }}}$. Gaussian states are those whose characteristic function is, by definition, of Gaussian form, i.e.

$$\begin{eqnarray}&&\chi ({\boldsymbol{\xi }})={{\rm{e}}}^{-\frac{1}{4}{{\boldsymbol{\xi }}}^{\dagger }\sigma {\boldsymbol{\xi }}-{\rm{i}}\;{{\boldsymbol{d}}}^{\dagger }K{\boldsymbol{\xi }}}.\end{eqnarray} \tag{ 3 }$$

In analogy with classical probability theory, Gaussian states are completely described by the first and the second statistical moments ${\boldsymbol{d}}$ and σ, where vector ${\boldsymbol{d}}$ and the positive-definite Hermitian matrix σ are defined as

$$\begin{eqnarray}&&{{\boldsymbol{d}}}_{i}=\mathrm{tr}\left[\rho {\hat{{\boldsymbol{A}}}}_{i}\right],\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&{\sigma }_{{ij}}=\mathrm{tr}\left[\hat{\rho }\;\{\Delta {\hat{{\boldsymbol{A}}}}_{i},\Delta {\hat{{\boldsymbol{A}}}}_{j}^{\dagger }\}\right],\end{eqnarray} \tag{ 4b }$$

where $\Delta \hat{{\boldsymbol{A}}}:= \hat{{\boldsymbol{A}}}-{\boldsymbol{d}}$, and $\{\cdot ,\cdot \}$ denotes the anti-commutator. We call ${\boldsymbol{d}}$ the displacement vector and σ the covariance matrix. Note that other authors use different conventions. We choose the covariance matrix defined by the anti-commutator of annihilation and creation operators which is known as the 'complex form', while some authors define it using the correlations between position and momenta operators. In our convention the vacuum is represented by the identity matrix I, i.e., the variance of the quadrature operators ${\hat{x}}_{i}$ and ${\hat{p}}_{i}$ are $\mathrm{var}({\hat{x}}_{i})\;=\;\mathrm{var}({\hat{p}}_{i})=+1$ (some authors define the vacuum variances as $+1/2$). This is of course only a definition and does not affect any physical interpretation of the results.

According to Williamson theorem [17–19] any positive-definite $2N\times 2N$ Hermitian matrix σ can be diagonalized using symplectic matrices, i.e. $\sigma ={{SDS}}^{\dagger }$, where S is an element of a complex representation [19] of the real symplectic group $\mathrm{Sp}(2N,{\mathbb{R}})$, i.e., S is an element of a group isomorphic to the $\mathrm{Sp}(2N,{\mathbb{R}})$, and D is a diagonal matrix. S and D take the form

$$\begin{eqnarray}S=\left[\begin{array}{cc}\alpha & \beta \\ \bar{\beta } & \bar{\alpha }\end{array}\right],\quad D=\left[\begin{array}{cc}L & 0\\ 0 & L\end{array}\right],\end{eqnarray} \tag{ 5 }$$

where S additionally satisfies the relation ${{SKS}}^{\dagger }=K$ with K defined by equation (1), and $L={\bf{diag}}({\lambda }_{1},...,{\lambda }_{N})$ is a diagonal matrix consisting of the so-called symplectic eigenvalues of a covariance matrix σ. This result will be used through this article.

The symplectic eigenvalues can be found by solving the usual eigenvalue problem for the matrix $K\sigma$. Eigenvalues of $K\sigma$ always appear in pairs. If ${\lambda }_{i}$ is an eigenvalue of $K\sigma$, then also $-{\lambda }_{i}$ is an eigenvalue of the same operator. The symplectic spectrum is then defined as a collection of the positive eigenvalues of $K\sigma$. In other words, ${\lambda }_{i}$ is a symplectic eigenvalue of σ if and only if it is positive and $\pm {\lambda }_{i}$ are the eigenvalues of the operator $K\sigma$. Symplectic eigenvalues are always greater or equal to one and are related to the purity of the Gaussian state. The state is pure if and only if for all i, ${\lambda }_{i}=1$, and the larger the symplectic eigenvalues are, the more mixed the state is. Knowing this, we can say symplectic eigenvalues are analogous to the eigenvalues of the density matrix ρ in the density-matrix formalism. On the other hand, symplectic matrices S usually represent some form of a squeezing or an entangling operation and are analogous to the unitary operators in the density matrix formalism. Given the special form in equation (5) of the symplectic matrices and the relation ${{SKS}}^{\dagger }=K$, one can easily prove [20] that the complex form of the symplectic matrices forms a subgroup of the more general pseudo-unitary group $U(N,N)=\{S\in {GL}(2N,{\mathbb{C}})| {{SKS}}^{\dagger }=K\}$.

For more details about the complex and the real form of the covariance matrix see appendix A, for a more detailed analysis of Gaussian states see [9, 21].

In this section we derive an exact expression for the quantum Fisher information for any two-mode Gaussian state. There are numerous ways to compute this quantity, however, for the purpose of this section we adopt the definition via the Bures distance [22]. The Bures distance is a measure of distinguishability between two quantum states ${\rho }_{\mathrm{1,2}}$ and is defined through the Uhlmann fidelity [23] ${\mathcal{F}}({\rho }_{1},{\rho }_{2}):= {\left(\mathrm{tr}\sqrt{\sqrt{{\rho }_{1}}\;{\rho }_{2}\;\sqrt{{\rho }_{1}}}\right)}^{2}$ as ${d}_{B}^{2}=2\left(1-\sqrt{{\mathcal{F}}({\rho }_{1},{\rho }_{2})}\right)$. The quantum Fisher information which measures how well we can distinguish two neighbouring states ${\rho }_{\epsilon }$ and ${\rho }_{\epsilon +d\epsilon }$ is defined as a limit [24]

$$\begin{eqnarray}&&H(\epsilon )=8\underset{d\epsilon \to 0}{\mathrm{lim}}\displaystyle \frac{1-\sqrt{{\mathcal{F}}({\rho }_{\epsilon },{\rho }_{\epsilon +d\epsilon })}}{d{\epsilon }^{2}}.\end{eqnarray} \tag{ 6 }$$

The problem of finding the quantum Fisher information thus reduces to expanding the fidelity around the point . As stated before, for Gaussian states the density matrix can be represented by a couple of the first and the second moments, ${\rho }_{1}\equiv ({{\boldsymbol{d}}}_{1},{\sigma }_{1})$, ${\rho }_{2}\equiv ({{\boldsymbol{d}}}_{2},{\sigma }_{2})$. In the case of a two-mode Gaussian state the fidelity can be written as [11]

$$\begin{eqnarray}&&{\mathcal{F}}({\rho }_{1},{\rho }_{2})=\displaystyle \frac{4{{\rm{e}}}^{-\delta {{\boldsymbol{d}}}^{\dagger }{\left({\sigma }_{1}+{\sigma }_{2}\right)}^{-1}\delta {\boldsymbol{d}}}}{\left(\sqrt{\Gamma }+\sqrt{\Lambda }\right)-\sqrt{{\left(\sqrt{\Gamma }+\sqrt{\Lambda }\right)}^{2}-\Delta }},\end{eqnarray} \tag{ 7 }$$

where $\delta {\boldsymbol{d}}={{\boldsymbol{d}}}_{1}-{{\boldsymbol{d}}}_{2}$ is the relative displacement and $\Delta ,\Gamma ,\Lambda$ denotes three determinants defined as

$$\begin{eqnarray}&&\Delta =| {\sigma }_{1}+{\sigma }_{2}| ,\end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray}&&\Gamma =| I+K{\sigma }_{1}K{\sigma }_{2}| ,\end{eqnarray} \tag{ 8b }$$

$$\begin{eqnarray}&&\Lambda =| {\sigma }_{1}+K| | {\sigma }_{2}+K| ,\end{eqnarray} \tag{ 8c }$$

with $K=I\oplus -I$ already introduced in the previous section. We note that our definitions differs from [11] by the factor of 2, which results in the factor of 4 in equation (7).

Let us denote an expansion of an arbitrary matrix around point up to second order in $d\epsilon$ as

$$\begin{eqnarray}&&M={M}_{0}+{M}_{1}d\epsilon +{M}_{2}d{\epsilon }^{2}+O(d{\epsilon }^{3}),\end{eqnarray} \tag{ 9 }$$

where $M=M(\epsilon +d\epsilon )$, ${M}_{0}=M(\epsilon )$, ${M}_{1}=\dot{M}(\epsilon )$, ${M}_{2}=\displaystyle \frac{1}{2}\ddot{M}(\epsilon )$, where dot denotes the derivative with respect to . Using this notation and the definition (6) we can write $H(\epsilon )=-4{{\mathcal{F}}}_{2}(\epsilon )$, where ${\mathcal{F}}$ from now on denotes the fidelity between the two close states ${\rho }_{\epsilon }$ and ${\rho }_{\epsilon +d\epsilon }$. The problem is that if we try to use equations (8) to expand equation (6) directly to find the expression of ${{\mathcal{F}}}_{2}$, we arrive at a complicated expression that depends on the second derivatives of σ. However, expressions of the quantum Fisher information in [2, 14] do not depend on second derivatives. This is true most of the time but there are places where the second derivatives appear. This fact is usually ignored in the literature and we will explain those scenarios later. First, to find an expression that depends only on first derivatives we use the Williamson theorem $\sigma ={{SDS}}^{\dagger }$ to rewrite equations (8) as

$$\begin{eqnarray}&&\Delta =| {D}_{0}+{{PDP}}^{\dagger }| ,\end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray}&&\Gamma =| I+{D}_{0}{{PDP}}^{\dagger }| ,\end{eqnarray} \tag{ 10b }$$

$$\begin{eqnarray}&&\Lambda =| {D}_{0}+K| | D+K| ,\end{eqnarray} \tag{ 10c }$$

where the matrix P is defined as $P={S}_{0}^{-1}S$. Following equation (9), P has expansion $P=I+{P}_{1}d\epsilon +{P}_{2}d{\epsilon }^{2}+O(d{\epsilon }^{3})$, where ${P}_{1}={S}_{0}^{-1}{S}_{1}=S{(\epsilon )}^{-1}\dot{S}(\epsilon )$ and ${P}_{2}={S}_{0}^{-1}{S}_{2}$. A useful property of P is that it is symplectic, i.e., ${{PKP}}^{\dagger }=K$, giving us conditions on the first and second derivatives,

$$\begin{eqnarray}&&{P}_{1}^{\dagger }=-{{KP}}_{1}K,\end{eqnarray} \tag{ 11a }$$

$$\begin{eqnarray}&&{P}_{2}^{\dagger }=-{{KP}}_{2}K+{{KP}}_{1}^{2}\;K.\end{eqnarray} \tag{ 11b }$$

Using the expansion of the determinant,

$$\begin{eqnarray}| {M}_{0}+{M}_{1}d\epsilon +{M}_{2}d{\epsilon }^{2}+O(d{\epsilon }^{3})| & = & | {M}_{0}| \left(1+\left(\mathrm{tr}{M}_{0}^{-1}{M}_{1}\right)d\epsilon \right.\\ & & +\left.\displaystyle \frac{1}{2}\left(2\mathrm{tr}{M}_{0}^{-1}{M}_{2}+{(\mathrm{tr}{M}_{0}^{-1}{M}_{1})}^{2}-\mathrm{tr}[{({M}_{0}^{-1}{M}_{1})}^{2}]\right)d{\epsilon }^{2}\right)\\ & & +O(d{\epsilon }^{3}),\end{eqnarray} \tag{ 12 }$$

which holds for an invertible matrix M₀, equations (11), and the cyclic property of the trace we remove the dependence on second derivatives P₂ and obtain

$$\begin{eqnarray}\Delta & = & | {D}_{0}| \left(16+8\mathrm{tr}[{D}_{0}^{-1}{D}_{1}]d\epsilon +\left(4\mathrm{tr}[{P}_{1}^{2}]-4\mathrm{tr}[{D}_{0}^{-1}{{KP}}_{1}{D}_{0}{{KP}}_{1}]+8\mathrm{tr}[{D}_{0}^{-1}{D}_{2}]\right.\right.\\ & & +\left.\left.2{(\mathrm{tr}[{D}_{0}^{-1}{D}_{1}])}^{2}-2\mathrm{tr}[{({D}_{0}^{-1}{D}_{1})}^{2}]\right)d{\epsilon }^{2}\right)+O(d{\epsilon }^{3}),\end{eqnarray} \tag{ 13a }$$

$$\begin{eqnarray}\Gamma & = & | C| \left(1+\mathrm{tr}[{C}^{-1}{D}_{0}{D}_{1}]d\epsilon +\left(\mathrm{tr}[{C}^{-1}{P}_{1}^{2}]-\mathrm{tr}[{({C}^{-1}{P}_{1})}^{2}]-\mathrm{tr}[({C}^{-1}{D}_{0}{{KP}}_{1}{)}^{2}]\right.\right.\\ & & +\left.\left.\mathrm{tr}[{C}^{-1}{D}_{0}{D}_{2}]+\displaystyle \frac{1}{2}{(\mathrm{tr}[{C}^{-1}{D}_{0}{D}_{1}])}^{2}-\displaystyle \frac{1}{2}\mathrm{tr}[{({C}^{-1}{D}_{0}{D}_{1})}^{2}]\right)d{\epsilon }^{2}\right)+O(d{\epsilon }^{3}),\end{eqnarray} \tag{ 13b }$$

$$\begin{eqnarray}\Lambda & = & {| E| }^{2}\left(1+\mathrm{tr}[{E}^{-1}{D}_{1}]d\epsilon +\left(\mathrm{tr}[{E}^{-1}{D}_{2}]\right.\right.+\left.\left.\displaystyle \frac{1}{2}{(\mathrm{tr}[{E}^{-1}{D}_{1}])}^{2}-\displaystyle \frac{1}{2}\mathrm{tr}[{({E}^{-1}{D}_{1})}^{2}]\right)d{\epsilon }^{2}\right)+O(d{\epsilon }^{3}),\end{eqnarray} \tag{ 13c }$$

where we have denoted $C=I+{D}_{0}^{2}$ and $E={D}_{0}+K$. In the above we assumed that there exist a Taylor expansion around point in this particular form, which is not true for Λ in the case where at least one of the symplectic eigenvalues of σ is equal to one. We address this subtle issue in appendix B. Assuming all symplectic eigenvalues are larger than one, we insert expressions (13) into the Uhlmann fidelity (7). We derive that the zeroth order sums to 1, the first order vanishes, and the second order provides the quantum Fisher information

$$\begin{eqnarray}&&H(\epsilon )=\displaystyle \frac{{\Delta }_{2}-8({(\sqrt{\Gamma })}_{2}+{(\sqrt{\Lambda })}_{2})}{2(\sqrt{{\Gamma }_{0}}+\sqrt{{\Lambda }_{0}}-4)}+2{{\boldsymbol{d}}}_{1}^{\dagger }{\sigma }_{0}^{-1}{{\boldsymbol{d}}}_{1}.\end{eqnarray} \tag{ 14 }$$

The denominator in the above expression actually has a compact expression in terms of either a) the determinant of the state or b) the symplectic eigenvalues, $\sqrt{{\Gamma }_{0}}+\sqrt{{\Lambda }_{0}}-4=2(| \sigma | -1)=2(| {D}_{0}| -1)$, which helps with computations considerably. To derive an expression which depends only on the first derivatives, we insert equations (13) into the above formula. The terms proportional to D₂ vanish giving us an alternative expression for the quantum Fisher information,

$$\begin{eqnarray}H(\epsilon ) & = & \displaystyle \frac{1}{| {D}_{0}| -1}\left(| {D}_{0}| \left(\mathrm{tr}[{P}_{1}^{2}]-\mathrm{tr}[{D}_{0}^{-1}{{KP}}_{1}{D}_{0}{{KP}}_{1}]\right)\right.\\ & & +\left.\sqrt{| C| }\left(\mathrm{tr}[{({C}^{-1}{P}_{1})}^{2}]+\mathrm{tr}[{({C}^{-1}{D}_{0}{{KP}}_{1})}^{2}]-\mathrm{tr}[{C}^{-1}{P}_{1}^{2}]\right)\right)\\ & & +\displaystyle \frac{1}{2}\mathrm{tr}\left[{({D}_{0}+K)}^{-1}{D}_{0}^{-1}{D}_{1}^{2}\right]+2{{\boldsymbol{d}}}_{1}^{\dagger }{\sigma }_{0}^{-1}{{\boldsymbol{d}}}_{1},\end{eqnarray} \tag{ 15 }$$

where ${P}_{1}={S}_{0}^{-1}{S}_{1}=S{(\epsilon )}^{-1}\dot{S}(\epsilon )$, ${D}_{0}=D(\epsilon )$, ${D}_{1}=\dot{D}(\epsilon )$, $C=I+{D}_{0}^{2}$, and ${{\boldsymbol{d}}}_{1}=\dot{{\boldsymbol{d}}}(\epsilon )$. The above formula is useful when we know the initial Williamson decomposition. For example when we are trying to estimate squeezing in the case of two-mode squeezed thermal state. However, in general finding the Williamson decomposition is not an easy task. That is why we find an alternative expression only in terms of the covariance matrix σ, displacement ${\boldsymbol{d}}$ and symplectic eigenvalues. For convenience we use the dot-notation, where dot denotes the derivative with respect to . For convenience we also denote $A:= K\sigma (\epsilon )$. The quantum Fisher information for a two-mode Gaussian state is given by

$$\begin{eqnarray}H(\epsilon ) & = & \displaystyle \frac{1}{2(| A| -1)}\left(| A| \mathrm{tr}\left[{\left({A}^{-1}\dot{A}\right)}^{2}\right]+\sqrt{| I+{A}^{2}| }\mathrm{tr}\left[{\left({(I+{A}^{2})}^{-1}\dot{A}\right)}^{2}\right]\right.\\ & & +\left.4\left({\lambda }_{1}^{2}-{\lambda }_{2}^{2}\right)\left(-\displaystyle \frac{{\dot{{\lambda }_{1}}}^{2}}{{\lambda }_{1}^{4}-1}+\displaystyle \frac{{\dot{{\lambda }_{2}}}^{2}}{{\lambda }_{2}^{4}-1}\right)\right)+2{\dot{{\boldsymbol{d}}}}^{\dagger }{\sigma }^{-1}\dot{{\boldsymbol{d}}},\end{eqnarray} \tag{ 16 }$$

where the symplectic eigenvalues of σ can be calculated as

$$\begin{eqnarray}&&{\lambda }_{\mathrm{1,2}}=\displaystyle \frac{1}{2}\sqrt{\mathrm{tr}[{A}^{2}]\pm \sqrt{{(\mathrm{tr}[{A}^{2}])}^{2}-16| A| }},\end{eqnarray} \tag{ 17 }$$

and $| \cdot |$ denotes the determinant. Using the Williamson theorem one can prove that the above formula reduces to equation (15). Also, strictly speaking, the above formula is defined only for covariance matrices with both symplectic eigenvalues larger than one. However, we can use a regularization procedure which allows us to use this formula in any case. This consists of multiplying the original covariance matrix σ by a mixedness parameter $\nu \gt 1$, using the formula (16) to calculate the quantum Fisher information for the state $\nu \sigma$, performing the limit $\nu \;\longrightarrow \;1$, and adding the second derivative $\ddot{{\lambda }_{i}}(\epsilon )$ for every symplectic eigenvalue ${\lambda }_{i}(\epsilon )$ which is equal to one. We need to add these second derivatives because by performing the limit $\nu \to 1$ we set the problematic terms $\displaystyle \frac{{\dot{{\lambda }_{i}}}^{2}}{{\lambda }_{i}^{4}-1}$ of equation (16) to zero. However, in cases when the symplectic eigenvalue ${\lambda }_{i}$ is equal to one, such terms have a non-zero contribution which needs to be accounted for. Altogether, we have

$$\begin{eqnarray}&&H(\epsilon )=\underset{\nu \to 1}{\mathrm{lim}}\;H(\nu \sigma (\epsilon ))+\underset{i:{\lambda }_{i}(\epsilon )=1}{\displaystyle \sum }{\ddot{\lambda }}_{i}(\epsilon ).\end{eqnarray} \tag{ 18 }$$

More details about the procedure can be found in appendix B.

Note that in the case where the symplectic eigenvalues do not change with a small variation in , i.e. $\dot{{\lambda }_{1}}=\dot{{\lambda }_{2}}=0$, or $\ddot{{\lambda }_{i}}=0$ for ${\lambda }_{i}=1$, the term depending on the symplectic eigenvalues vanishes. This includes the case where purity does not change or where the parameter of interest was encoded into the initial state by a symplectic transformation. Although all the computations were performed in the complex representation of the covariance matrix, we can easily transform the result to the real representation. This is done by substituting $A={\rm{i}}\Omega {\sigma }_{R}$, ${\dot{{\boldsymbol{d}}}}^{\dagger }{\sigma }^{-1}\dot{{\boldsymbol{d}}}\;\longrightarrow \;{\dot{{\boldsymbol{d}}}}_{R}^{T}{\sigma }_{R}^{-1}{\dot{{\boldsymbol{d}}}}_{R}$, where ${\rm{i}}\Omega$ is the real form equivalent to matrix K, ${\sigma }_{R}$ the real form covariance matrix and ${{\boldsymbol{d}}}_{R}$ the real form displacement. For more details see appendix A.

In the previous section we derived an exact expression for two-mode Gaussian states, however, in recent work [14] a general formula for the quantum Fisher information was derived as a limit of a particular infinite series. Here we simplify the expression for the infinite series. We also find a bound on the remainder of the series when we sum only a finite number of terms. We then go on to simplify an already known formula for pure states and derive an exact expression for the multi-mode quantum Fisher information for the cases when the Williamson decomposition is known.

According to [14], the quantum Fisher information for a general multi-mode case which has all symplectic eigenvalues larger than one can be calculated as

$$\begin{eqnarray}&&H(\epsilon )=\displaystyle \frac{1}{2}\mathrm{tr}\left[\dot{\sigma }Y\right]+2{\dot{{\boldsymbol{d}}}}^{\dagger }{\sigma }^{-1}\dot{{\boldsymbol{d}}},\end{eqnarray} \tag{ 19 }$$

where Y is a solution to the so-called Stein equation [25]

$$\begin{eqnarray}&&\dot{\sigma }=\sigma Y\sigma -{KYK}.\end{eqnarray} \tag{ 20 }$$

If all symplectic eigenvalues of σ are larger than one, the solution is unique and can be written as an infinite series

$$\begin{eqnarray}&&Y=-\underset{n=0}{\overset{\infty }{\displaystyle \sum }}{(K\sigma )}^{-n}\dot{({\sigma }^{-1})}{(\sigma K)}^{-n}.\end{eqnarray} \tag{ 21 }$$

Inserting (21) into (19), using $\dot{({\sigma }^{-1})}=-{\sigma }^{-1}\dot{\sigma }{\sigma }^{-1}$, and the properties of matrix K we find an elegant expression for the multi-mode quantum Fisher information,

$$\begin{eqnarray}&&H(\epsilon )=\displaystyle \frac{1}{2}\underset{n=1}{\overset{M}{\displaystyle \sum }}\mathrm{tr}\left[{({A}^{-n}\dot{A})}^{2}\right]+{R}_{M}+2{\dot{{\boldsymbol{d}}}}^{\dagger }{\sigma }^{-1}\dot{{\boldsymbol{d}}},\end{eqnarray} \tag{ 22 }$$

where we use the notation $A:= K\sigma$ or $A:= {\rm{i}}\Omega {\sigma }_{R}$ for the real form of the covariance matrix and R_M is the remainder of the series $M\;\longrightarrow \;\infty$. As we prove in appendix C, the remainder is bounded $| {R}_{M}| \leqslant \displaystyle \frac{\mathrm{tr}[{(A\dot{A})}^{2}]}{2{\lambda }_{\mathrm{min}}^{2M+2}({\lambda }_{\mathrm{min}}^{2}-1)}$ with ${\lambda }_{\mathrm{min}}:= {\mathrm{min}}_{i}\{{\lambda }_{i}\}$ being the smallest symplectic eigenvalue, i.e., the smallest positive eigenvalue of matrix A. This means that for ${\lambda }_{\mathrm{min}}\gt 1$ we have ${\mathrm{lim}}_{M\longrightarrow \infty }{R}_{M}=0$, the series converges, and we can write

$$\begin{eqnarray}&&H(\epsilon )=\displaystyle \frac{1}{2}\underset{n=1}{\overset{\infty }{\displaystyle \sum }}\mathrm{tr}\left[{({A}^{-n}\dot{A})}^{2}\right]+2{\dot{{\boldsymbol{d}}}}^{\dagger }{\sigma }^{-1}\dot{{\boldsymbol{d}}}.\end{eqnarray} \tag{ 23 }$$

To calculate the quantum Fisher information for the states which have some eigenvalues equal to one we can, once again, use the regularization procedure (18). A small example, which has been already shown in [14], is to consider the class of iso-thermal states given by ${A}^{2}={\nu }^{2}I$, $\nu \gt 1$, which is equivalent to $D(\epsilon )=\nu I$. For such states we can easily evaluate the infinite sum (23) and derive

$$\begin{eqnarray}&&H(\epsilon )=\displaystyle \frac{{\nu }^{2}}{2(1+{\nu }^{2})}\mathrm{tr}\left[{({A}^{-1}\dot{A})}^{2}\right]+2{\dot{{\boldsymbol{d}}}}^{\dagger }{\sigma }^{-1}\dot{{\boldsymbol{d}}},\end{eqnarray} \tag{ 24 }$$

which for $\nu =1$ gives a formula for pure states also derived in [10]. As noted in [15], this formula can be further simplified. Differentiating ${A}^{2}={\nu }^{2}$ and multiplying each side by ${A}^{-1}$ we obtain an anti-commutation relation ${A}^{-1}\dot{A}=-\dot{A}{A}^{-1}$ which together with ${A}^{-2}={\nu }^{-2}I$ gives

$$\begin{eqnarray}&&H(\epsilon )=-\displaystyle \frac{1}{2(1+{\nu }^{2})}\mathrm{tr}\left[{\dot{A}}^{2}\right]+2{\dot{{\boldsymbol{d}}}}^{\dagger }{\sigma }^{-1}\dot{{\boldsymbol{d}}},\end{eqnarray} \tag{ 25 }$$

meaning that for iso-thermal states with fixed displacement we do not need to invert the covariance matrix anymore.

Using $\mathrm{tr}[{({A}^{-n}\dot{A})}^{2}]=2\mathrm{tr}[{({D}_{0}^{-n+1}{K}^{-n+1}{P}_{1})}^{2}]-2\mathrm{tr}[{D}_{0}^{-n+2}{K}^{-n}{P}_{1}{D}_{0}^{-n}{K}^{-n}{P}_{1}]+\mathrm{tr}[{D}_{0}^{-2n}{D}_{1}^{2}]$ we can rewrite formula (23) for the quantum Fisher information in terms of the elements of the symplectic matrices and eigenvalues,

$$\begin{eqnarray}H(\epsilon ) & = & \mathrm{tr}\left[{P}_{1}^{2}\right]-\underset{n=1}{\overset{\infty }{\displaystyle \sum }}\mathrm{tr}\left[\left({D}_{0}^{2}-I\right){\left({D}_{0}^{-n}{K}^{-n}{P}_{1}\right)}^{2}\right]+\displaystyle \frac{1}{2}\mathrm{tr}\left[{({D}_{0}+K)}^{-1}{D}_{0}^{-1}{D}_{1}^{2}\right]+2{{\boldsymbol{d}}}_{1}^{\dagger }{\sigma }_{0}^{-1}{{\boldsymbol{d}}}_{1},\end{eqnarray} \tag{ 26 }$$

where we have used the notation from equation (15). We can see that the diagonal part coincides with the two-mode case indicating the validity of the general multi-mode formula. However, in contrast to equation (23) the infinite series here converges also when symplectic eigenvalues are equal to one. Nevertheless, this does not mean that it is valid to use this formula for such states. To be more specific, plainly inserting ${D}_{0}=I$ does not give the correct formula for pure states—this still needs to be obtained by the regularization procedure. The reason why it does not give the correct result is that in general limits ${\mathrm{lim}}_{M\to \infty }$ and ${\mathrm{lim}}_{\nu \to 1}$ do not commute. On the other hand, the infinite sum in this formula leads to a geometric series that can be evaluated, as shown in appendix D. This allows us to derive a much more elegant, entirely general and exact formula for the quantum Fisher information in terms of elements of the Williamson decomposition $\sigma ={{SDS}}^{\dagger }$ for multi-mode Gausian states. Using the definitions (5) and (11a), we find that matrix P₁ has an elegant structure

$$\begin{eqnarray}{P}_{1}(\epsilon )=\left[\begin{array}{cc}R & Q\\ \bar{Q} & \bar{R}\end{array}\right],\end{eqnarray} \tag{ 27 }$$

where $R={\alpha }^{\dagger }\dot{\alpha }-\bar{{\beta }^{\dagger }\dot{\beta }}$ is a skew-Hermitian and $Q={\alpha }^{\dagger }\dot{\beta }-\bar{{\beta }^{\dagger }\dot{\alpha }}$ a (complex) symmetric matrix. Matrix P₁ is actually an element of the Lie algebra associated with the complex form of the real symplectic group. If the diagonalizing symplectic matrix forms a one-parameter group $S={{\rm{e}}}^{X\epsilon }$, where X is an element of the algebra independent of , then ${P}_{1}=X$. Inserting equation (27) into equation (26) and evaluating the infinite sum we derive the quantum Fisher information for the N-mode Gaussian state

$$\begin{eqnarray}&&H(\epsilon )=\underset{i,j=1}{\overset{N}{\displaystyle \sum }}\displaystyle \frac{{({\lambda }_{i}-{\lambda }_{j})}^{2}}{{\lambda }_{i}{\lambda }_{j}-1}{\left|{R}_{{ij}}\right|}^{2}+\displaystyle \frac{{({\lambda }_{i}+{\lambda }_{j})}^{2}}{{\lambda }_{i}{\lambda }_{j}+1}{\left|{Q}_{{ij}}\right|}^{2}+\underset{i=1}{\overset{N}{\displaystyle \sum }}\displaystyle \frac{{\dot{{\lambda }_{i}}}^{2}}{{\lambda }_{i}^{2}-1}+2{\dot{{\boldsymbol{d}}}}^{\dagger }{\sigma }^{-1}\dot{{\boldsymbol{d}}},\end{eqnarray} \tag{ 28 }$$

where ${\lambda }_{i}$ are the symplectic eigenvalues. Strictly speaking, the formula is not defined for the symplectic eigenvalues equal to one. Under the assumption of differentiability of the covariance matrix, for ${\lambda }_{i}(\epsilon )={\lambda }_{j}(\epsilon )=1$ both $\displaystyle \frac{{({\lambda }_{i}-{\lambda }_{j})}^{2}}{{\lambda }_{i}{\lambda }_{j}-1}(\epsilon )$ and $\displaystyle \frac{{\dot{{\lambda }_{i}}}^{2}}{{\lambda }_{i}^{2}-1}(\epsilon )$ are expressions of type $\displaystyle \frac{0}{0}$. Nevertheless, we can define these problematic points in a way which makes the quantum Fisher information a continuous function, which also corresponds to the knowledge gained from a study of the states close to the points of purity described in appendix B. For ${\lambda }_{i}(\epsilon )=1$ we define $\displaystyle \frac{{\dot{{\lambda }_{i}}}^{2}}{{\lambda }_{i}^{2}-1}(\epsilon ):= \ddot{{\lambda }_{i}}(\epsilon )$, and for ${\lambda }_{i}(\epsilon )={\lambda }_{j}(\epsilon )=1$ we define $\displaystyle \frac{{({\lambda }_{i}-{\lambda }_{j})}^{2}}{{\lambda }_{i}{\lambda }_{j}-1}(\epsilon ):= 0$.

Assuming all symplectic eigenvalues are larger than one, we can define the Hermitian matrix ${\tilde{R}}_{{ij}}:= \displaystyle \frac{{\lambda }_{i}-{\lambda }_{j}}{\sqrt{{\lambda }_{i}{\lambda }_{j}-1}}{R}_{{ij}}$ and the symmetric matrix ${\tilde{Q}}_{{ij}}:= \displaystyle \frac{{\lambda }_{i}+{\lambda }_{j}}{\sqrt{{\lambda }_{i}{\lambda }_{j}+1}}{Q}_{{ij}}$ which allows us to write the quantum Fisher information in a compact form,

$$\begin{eqnarray}&&H(\epsilon )=\mathrm{tr}\left[\tilde{R}{\tilde{R}}^{\dagger }+\tilde{Q}{\tilde{Q}}^{\dagger }\right]+\mathrm{tr}\left[{({L}^{2}-I)}^{-1}{\dot{L}}^{2}\right]+2{\dot{{\boldsymbol{d}}}}^{\dagger }{\sigma }^{-1}\dot{{\boldsymbol{d}}}.\end{eqnarray} \tag{ 29 }$$

Now we can easily interpret each term. The first part corresponds to the change in relative orientation and squeezing, the second to the change in purity and the third to the change in displacement. Added together, they all contribute to the quantum Fisher information and increase the precision in the estimation of a parameter .

Note that formulas (28), (29) respectively, are actually multi-mode generalizations of equation (15). The form (15) takes when we rewrite it in terms of submatrices given by equation (27) is exactly the same as (28) for N = 2. The same holds for the one-mode formula derived in [12] and N = 1, partially validating our general result. For the derivation of the formula see appendix D, for the details why we choose that particular definition of problematic points see appendix B, and for the quantum Fisher information in the real form formalism see appendix A.

To illustrate the use of the general formula (28) we derive the precision with which we can estimate the squeezing parameter of a one-mode squeezing channel. Similar problems have already been studied. For example in [26] the vacuum and coherent probe states were considered, in [27] it was the coherent states and the displaced squeezed vacuum states, which were squeezed in either position or momenta direction. Here we are generalizing the precision bounds in those articles. We also note that a similar work has been done in [28] using the density matrix formalism for a different type of a one-mode squeezing channel ${\tilde{S}}_{r}=\mathrm{exp}(-{\rm{i}}\displaystyle \frac{r}{2}({\hat{a}}^{2}+{\hat{a}}^{\dagger 2}))$. Here, however, we use an entirely general single-mode Gaussian state as an input state.

We choose a one-mode squeezed rotated displaced thermal state as an input state, which is the most general single-mode state [21], with the initial squeezing parameter r, angle of rotation θ, and the initial displacement $\tilde{{d}_{0}}=\left|\tilde{{d}_{0}}\right|{{\rm{e}}}^{{\rm{i}}\phi }$. We feed this state into a one-mode squeezing channel, which encodes the unknown parameter we are trying to estimate, leaving us with the final state,

$$\begin{eqnarray}&&{\rho }_{\epsilon }={S}_{\epsilon }{D}_{\tilde{{d}_{0}}}{R}_{\theta }{S}_{r}{\rho }_{\mathrm{th}}{S}_{r}^{\dagger }{R}_{\theta }^{\dagger }{D}_{\tilde{{d}_{0}}}^{\dagger }{S}_{\epsilon }^{\dagger },\end{eqnarray} \tag{ 30 }$$

where ${S}_{r}=\mathrm{exp}(\displaystyle \frac{r}{2}({\hat{a}}^{2}-{\hat{a}}^{\dagger 2}))$ is the one-mode squeezing operator, ${R}_{\theta }=\mathrm{exp}(-{\rm{i}}\theta {\hat{a}}^{\dagger }\hat{a})$ the rotation operator, and ${D}_{\tilde{{d}_{0}}}=\mathrm{exp}(\tilde{{d}_{0}}{\hat{a}}^{\dagger }-\bar{\tilde{{d}_{0}}}\hat{a})$ the displacement operator. As we show in appendix E, the quantum Fisher information is

$$\begin{eqnarray}H(\epsilon ) & = & \displaystyle \frac{4{\lambda }_{1}^{2}}{{\lambda }_{1}^{2}+1}\left({\mathrm{cosh}}^{4}r+{\mathrm{sinh}}^{4}r-2\mathrm{cos}(4\theta ){\mathrm{cosh}}^{2}r{\mathrm{sinh}}^{2}r\right)\\ & & +\displaystyle \frac{4{\left|\tilde{{d}_{0}}\right|}^{2}}{{\lambda }_{1}}(\mathrm{cosh}(2r)+\mathrm{cos}(2(\theta -\phi ))\mathrm{sinh}(2r)).\end{eqnarray} \tag{ 31 }$$

Now, let us have a few notes on the derived formula. The result does not depend on the parameter we want to estimate, which must be true for every encoding channel which forms a one-parameter unitary group $U(\epsilon )={{\rm{e}}}^{-{\rm{i}}\hat{K}\epsilon }$, where $\hat{K}$ is a Hermitian operator. Also, as pointed out in [29], because the symplectic eigenvalue ${\lambda }_{1}\geqslant 1$ is proportional to temperature, we can immediately see that although the thermality slightly enhances the estimation through squeezing (with the maximal enhancement by the factor of 2), it reduces the estimation from the displacement. More importantly, however, by choosing an appropriate rotation of an input state we can significantly increase the estimation precision. Without loss of generality we assume $r\geqslant 0$. The maximal value of (31) is achieved when we choose to rotate for example by the value $\theta =\displaystyle \frac{\pi }{4}$, and displace in the direction $\phi =\displaystyle \frac{\pi }{4}$,

$$\begin{eqnarray}&&{H}_{\mathrm{max}}(\epsilon )=\displaystyle \frac{4{\lambda }_{1}^{2}}{{\lambda }_{1}^{2}+1}{\mathrm{cosh}}^{2}(2r)+\displaystyle \frac{4}{{\lambda }_{1}}{\left|\tilde{{d}_{0}}\right|}^{2}{{\rm{e}}}^{2r},\end{eqnarray} \tag{ 32 }$$

which signifies an exponential increase in the precision of estimation the unknown parameter when considering an initially squeezed and rotated state. To demonstrate, we derive a formula for a squeezing needed to enhance by k orders of magnitude when assuming zero initial displacement, $r=\mathrm{arcsinh}\sqrt{\frac{{10}^{\frac{k}{2}}-1}{2}}$, which for larger k (or $r\gtrsim 1$) behaves as $r\approx 0.35+0.58k$. The current state-of-the-art [30] achieves the squeezing around r = 1.46, which could hypothetically account for an improvement by a factor of 80 as compared to the case where the initial squeezing is zero. As we show in figure 1, the reason why the amount of distinguishability of the two close states rises is because the final squeezing forces the covariance matrix to turn. Also, note that although the initial squeezing leads to an exponential increase, the initial displacement contributes only quadratically. To conclude, the optimal Gaussian state for the estimation of the parameter of a one-mode squeezing channel is a thermal state infinitely squeezed in the angle of $\theta =\frac{\pi }{4}$ from the direction of the squeezing channel we want to estimate, and infinitely displaced in the direction in which the squeezed state is stretched, $\phi =\frac{\pi }{4}$. For a fixed amount of squeezing and displacement, the optimal temperature which maximizes the quantum Fisher information is given by a solution of $\frac{{\lambda }_{1}^{3}}{{({\lambda }_{1}^{2}+1)}^{2}}=\frac{{\left|\tilde{{d}_{0}}\right|}^{2}{{\rm{e}}}^{2r}}{2{\mathrm{cosh}}^{2}(2r)}$. Also, because ${\lambda }_{1}=1+2{n}_{\mathrm{th}}$, where ${n}_{\mathrm{th}}=\mathrm{tr}[{\rho }_{\mathrm{th}}{\hat{a}}^{\dagger }\hat{a}]$ denotes the mean number of thermal bosons, and the mean total number of bosons in displaced squeezed states is given by

$$\begin{eqnarray}&&n={n}_{\tilde{{d}_{0}}}+{n}_{\mathrm{th}}+(1+2{n}_{\mathrm{th}}){\mathrm{sinh}}^{2}r,\end{eqnarray} \tag{ 33 }$$

where we have denoted ${n}_{\tilde{{d}_{0}}}:= {\left|\tilde{{d}_{0}}\right|}^{2}$ the mean number of bosons coming from displacement, we can rewrite the maximal quantum Fisher information (32) as

$$\begin{eqnarray}{H}_{\mathrm{max}}(\epsilon ) & = & \displaystyle \frac{2{(1+2n-2{n}_{\tilde{{d}_{0}}})}^{2}}{1+2{n}_{\mathrm{th}}(1+{n}_{\mathrm{th}})}+\displaystyle \frac{4{n}_{\tilde{{d}_{0}}}\left(1+2n-2{n}_{\tilde{{d}_{0}}}+2\sqrt{n-{n}_{\tilde{{d}_{0}}}-{n}_{\mathrm{th}}}\sqrt{1+n-{n}_{\tilde{{d}_{0}}}+{n}_{\mathrm{th}}}\right)}{{(1+2{n}_{\mathrm{th}})}^{2}}.\end{eqnarray} \tag{ 34 }$$

For a fixed total number of bosons this function achieves the maximum at ${n}_{\mathrm{th}}={n}_{\tilde{{d}_{0}}}=0$. This means that for a finite amount of available energy, the best protocol is to invest it all into squeezing, which is the same result as in [28, 29]. The optimal probe state is then the $\frac{\pi }{4}$-rotated squeezed vacuum state with the quantum Fisher information $H(\epsilon )=2{(1+2n)}^{2}$, which clearly indicates the Heisenberg scaling.

Zoom In Zoom Out Reset image size

We have derived an exact formula for the quantum Fisher information of an arbitrary two-mode Gaussian state. This has been done using the definition of the infinitesimal Bures distance, the Williamson decomposition of positive-definite matrices and the properties of the real symplectic group. Although the formula is not directly applicable for the states with pure modes, we introduced a regularization procedure which allows us to overcome this problem. Then, using a different approach, we simplified an already known formula for the multi-mode quantum Fisher information in terms of an infinite series. We also estimated the remainder of the series, allowing for an effective numerical calculation. Using the previous result, we showed that for the cases when the Williamson decomposition of the covariance matrix is known, the quantum Fisher information for multi-mode Gaussian state can be computed exactly. The general multi-mode formula is equivalent to the known results for one-mode Gaussian state when setting N = 1 and to the previously mentioned two-mode Gaussian states when N = 2. However, we note that using the requirement of the continuity of the quantum Fisher information and studying the case of the pure states gave us a different definition for the problematic points than is mentioned in [12]. Finally we applied our newly gained formula to study the case of the estimation a squeezing parameter of a one-mode squeezing channel. We showed that a strategy of squeezing and rotating the input state can significantly improve the precision in estimation.

We believe the main achievement of this article is in the usefulness of the derived formulae. It allows for the study of the optimal input states of Gaussian nature [31], it helps predict the ultimate sensitivity of a physical detector's particular implementation [3, 4], or to analyse the effects of temperature on the current gravitational wave detector proposal [5]. It gives a limit in the estimation of time [32] or temperature [33]. Also, since certain objects called the Bogoliubov transformations are isomorphic to the symplectic transformations, the natural application lies wherever these transformations appear. This is for example for quantum field theory in curved spacetime [34, 35] but also Bose–Einstein condensates [36] or scattering problems [37].

Thanks to Carlos Sabín, Mădălin Guţă, Alex Monras, David Edward Bruschi, Jan Kohlrus, and Katarzyna Macieszczak for fruitful discussions. ARL thanks the EPSRC Doctoral prize. IF acknowledges support from EPSRC INSPIRE: EP/M003019/1.

In this section we describe the structure of covariance matrices and displacement. We introduce the real form covariance matrix formalism and find how covariance matrices in the real form transform into its complex form.

The complex form of the covariance matrix and displacement was defined by equations (1) and (4). From the definition we can observe the following structure of the first and second moments:

$$\begin{eqnarray}&&{\boldsymbol{d}}=\left[\begin{array}{c}\tilde{{\boldsymbol{d}}}\\ \bar{\tilde{{\boldsymbol{d}}}}\end{array}\right],\quad \sigma =\left[\begin{array}{cc}X & Y\\ \bar{Y} & \bar{X}\end{array}\right]\end{eqnarray} \tag{ A.1 }$$

and ${\sigma }^{\dagger }=\sigma$, i.e., ${X}^{\dagger }=X$ and ${Y}^{T}=Y$. Using the Williamson theorem [17–19] we can write $\sigma ={{SDS}}^{\dagger }$, where ${{SKS}}^{\dagger }=K$. Although the symplectic matrices S are not necessarily Hermitian, they follow the same structure as σ which is expressed by equation (5).

Construction of the real form of the covariance matrix is analogous to the complex form described in the introduction. It is usually defined with respect to the collection of quadrature operators $\hat{{\boldsymbol{Q}}}:= \hat{{\boldsymbol{x}}}\oplus \hat{{\boldsymbol{p}}}=({\hat{x}}_{1},{\hat{x}}_{2},\ldots ,{\hat{p}}_{1},{\hat{p}}_{2},\ldots )$. The canonical commutation relations of these operators can be conveniently expressed as

$$\begin{eqnarray}[{\hat{Q}}_{i},{\hat{Q}}_{j}]\;=\;+{\rm{i}}\;{\Omega }_{{ij}}\;\mathrm{id}\quad \Rightarrow \quad \Omega =\left[\begin{array}{cc}0 & I\\ -I & 0\end{array}\right].\end{eqnarray} \tag{ A.2 }$$

Properties of Ω are $-{\Omega }^{2}=I$ and ${\Omega }^{T}=-\Omega$, in contrast to the complex form version K. In the real form, the definitions of the first and second moments are

$$\begin{eqnarray}&&{{\boldsymbol{d}}}_{R}=\mathrm{tr}\left[\hat{\rho }\hat{{\boldsymbol{Q}}}\right]=\left[\begin{array}{c}{\boldsymbol{x}}\\ {\boldsymbol{p}}\end{array}\right],\end{eqnarray} \tag{ A.3a }$$

$$\begin{eqnarray}{\sigma }_{R}=\mathrm{tr}\left[\hat{\rho }\{\hat{{\boldsymbol{Q}}},\hat{{\boldsymbol{Q}}}\}\right]=\left[\begin{array}{cc}{X}_{R} & {Y}_{R}\\ {Y}_{R}^{T} & {Z}_{R}\end{array}\right].\end{eqnarray} \tag{ A.3b }$$

The real covariance matrix is symmetric, i.e. ${X}_{R}={X}_{R}^{T}$ and ${Z}_{R}={Z}_{R}^{T}$. The corresponding real symplectic matrices are given by ${\sigma }_{R}={S}_{R}{D}_{R}{S}_{R}^{T}$, where $S\Omega {S}^{T}=\Omega$, which is a defining relation of the real symplectic group $\mathrm{Sp}(2N,{\mathbb{R}})$.

Since the change between real and complex form of the covariance matrix is a simple basis transformation, $\hat{{\boldsymbol{Q}}}\to \hat{{\boldsymbol{A}}}$, we can relate these two using the unitary matrix U,

$$\begin{eqnarray}\hat{{\boldsymbol{A}}}=U\hat{{\boldsymbol{Q}}},\ U=\displaystyle \frac{1}{\sqrt{2}}\left[\begin{array}{cc}I & +{\rm{i}}I\\ I & -{\rm{i}}I\end{array}\right].\end{eqnarray} \tag{ A.4 }$$

The resulting transformation between real and complex covariance matrices and displacement are

$$\begin{eqnarray}&&{\boldsymbol{d}}=U{{\boldsymbol{d}}}_{R},\ \sigma =U{\sigma }_{R}{U}^{\dagger },\end{eqnarray} \tag{ A.5 }$$

and the transformations related to the Williamson decomposition are

$$\begin{eqnarray}&&S={{US}}_{R}{U}^{\dagger },\ D={{UD}}_{R}{U}^{\dagger }={D}_{R},\ K=U{\rm{i}}\Omega {U}^{\dagger }.\end{eqnarray} \tag{ A.6 }$$

We explicitly write the connection between real and complex form of symplectic matrix,

$$\begin{eqnarray}{S}_{R}=\left[\begin{array}{cc}{\alpha }_{R} & {\beta }_{R}\\ {\gamma }_{R} & {\delta }_{R}\end{array}\right]=\left[\begin{array}{cc}\mathrm{Re}\;[\alpha +\beta ] & -\mathrm{Im}\;[\alpha -\beta ]\\ \mathrm{Im}\;[\alpha +\beta ] & \mathrm{Re}\;[\alpha -\beta ]\end{array}\right].\end{eqnarray} \tag{ A.7 }$$

Consequently, $\alpha$ and $\beta$ needed for R and Q from equation (28) can be expressed in the real form symplectic matrix elements as

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{1}{2}({\alpha }_{R}+{\delta }_{R}+{\rm{i}}{\gamma }_{R}-{\rm{i}}{\beta }_{R}),\end{eqnarray} \tag{ A.8a }$$

$$\begin{eqnarray}&&\beta =\displaystyle \frac{1}{2}({\alpha }_{R}-{\delta }_{R}+{\rm{i}}{\gamma }_{R}+{\rm{i}}{\beta }_{R}).\end{eqnarray} \tag{ A.8b }$$

Since all important matrices are related via this unitary transformation and traces and determinants are invariant under such transformations, it is clear that every formula we derived can be easily rewritten in the real form formalism by the formal substitution $\sigma \to {\sigma }_{R}$ and $K\to {\rm{i}}\Omega$. On the other hand, the complex form provides a much more elegant structure and exposes the inner symmetries of symplectic and covariance matrices in more detail. Also, it is much easier to work with K since it is diagonal, unitary, and Hermitian in contrast to non-diagonal skew-Hermitian matrix Ω, providing much more convenient language.

In this section we derive a formula for the QFI for the states around the points of purity, i.e., for the points where ${\lambda }_{i}(\epsilon )=1$, but not necessarily ${\lambda }_{i}(\epsilon +d\epsilon )=1$. We can see that formula (16) is undefined at points ${\lambda }_{i}(\epsilon )=1$. This is the consequence that we assumed that for all $\Delta ,\Lambda ,\Gamma$ the Taylor expansion exist in a form given by equations (13), which is not true when ${\lambda }_{i}(\epsilon )=1$ for some $i\in \{1,2\}$.

For simplicity let us study the case where is such that the both ${\lambda }_{1}(\epsilon )={\lambda }_{2}(\epsilon )=1$ and at the end of the calculation it will be clear how various cases work. For such states we have $D=I+{D}_{1}d\epsilon +{D}_{2}d{\epsilon }^{2}+O(d{\epsilon }^{3})$. From equations (10) we can see that $\Lambda =0$ and $\Delta =\Gamma$ and the Uhlmann fidelity (7) simplifies to

$$\begin{eqnarray}&&{\mathcal{F}}\left({\rho }_{\epsilon },{\rho }_{\epsilon +d\epsilon }\right)=\displaystyle \frac{4}{\sqrt{\Delta }},\end{eqnarray} \tag{ B.1 }$$

where for simplicity we have omitted the part consisting of displacement. This coincides with the general formula [13] for fidelity between one pure and one pure or mixed Gaussian state, which is exactly our case. Now, because ${\lambda }_{i}(\epsilon )=1$ and the symplectic eigenvalue cannot fall below 1, we have that either ${\dot{\lambda }}_{i}(\epsilon )=0$ or the derivative does not exist. If it exists and is, for example, positive, then for $d\epsilon \lt 0$ we have $\lambda (\epsilon +d\epsilon )=1+{\dot{\lambda }}_{i}(\epsilon )d\epsilon +O(d{\epsilon }^{2})\lt 1$ for small enough $d\epsilon$, which cannot be true. From now on we assume that map $\sigma \;:\epsilon \;\longrightarrow \;\sigma (\epsilon )$ is differentiable, thus ${D}_{1}=0$ and $D=I+{D}_{2}d{\epsilon }^{2}+O(d{\epsilon }^{3})$. The simplified formula (B.1) for fidelity gives ${\mathcal{F}}\left({\rho }_{\epsilon },{\rho }_{\epsilon +d\epsilon }\right)=1-\displaystyle \frac{1}{8}\left(\mathrm{tr}\left[{P}_{1}^{2}\right]-\mathrm{tr}\left[{({{KP}}_{1})}^{2}\right]+2\mathrm{tr}[{D}_{2}]\right)d{\epsilon }^{2}+O(d{\epsilon }^{3})$, leading to the QFI,

$$\begin{eqnarray}&&H(\epsilon )=\displaystyle \frac{1}{2}\left(\mathrm{tr}\left[{P}_{1}^{2}\right]-\mathrm{tr}\left[{({{KP}}_{1})}^{2}\right]\right)+\mathrm{tr}[{D}_{2}].\end{eqnarray} \tag{ B.2 }$$

Since ${D}_{2}=\displaystyle \frac{1}{2}\ddot{D}$, we have $\mathrm{tr}[{D}_{2}]={\displaystyle \sum }_{i}{\ddot{\lambda }}_{i}$. We can rewrite the above equation in terms of $A=K\sigma$,

$$\begin{eqnarray}&&H(\epsilon )=\displaystyle \frac{1}{4}\left(2\mathrm{tr}\left[{A}^{-1}\ddot{A}\right]-\mathrm{tr}\left[{({A}^{-1}\dot{A})}^{2}\right]\right).\end{eqnarray} \tag{ B.3 }$$

If the state remains pure with a small variation in , then ${\ddot{\lambda }}_{i}=0$ and the above formula is equivalent to the pure state formula $H(\epsilon )=\displaystyle \frac{1}{4}\mathrm{tr}\left[{({\sigma }^{-1}\dot{\sigma })}^{2}\right]$, which can be also derived from the equation for iso-thermal states (24). We see that although general mixed-state formulae do not depend on the second derivative, formulae (B.2) and (B.3) do. But we can show that we achieve the same result from the mixed state formula when requiring the continuity of the QFI. If we assume that $\sigma \;:\epsilon \;\longrightarrow \;\sigma (\epsilon )\in {C}^{(2)}$, i.e., the second derivative exists everywhere and is continuous, then the QFI is continuous everywhere apart from the points where it is undefined. To make the QFI a continuous function we define the problematic points given by ${\lambda }_{i}(\epsilon )=1$ as limits

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{\lambda }}_{i}^{2}}{{\lambda }_{i}^{2}-1}(\epsilon ):= \underset{d\epsilon \to 0}{\mathrm{lim}}\;\displaystyle \frac{{\dot{\lambda }}_{i}{(\epsilon +d\epsilon )}^{2}}{{\lambda }_{i}{(\epsilon +d\epsilon )}^{2}-1}={\ddot{\lambda }}_{i}(\epsilon ),\end{eqnarray} \tag{ B.4 }$$

since ${\lambda }_{i}$ has expansion ${\lambda }_{i}(\epsilon +d\epsilon )=1+\displaystyle \frac{1}{2}{\ddot{\lambda }}_{i}(\epsilon )d{\epsilon }^{2}+O(d{\epsilon }^{3})$. Now we see that the our definition of problematic values using the second derivative corresponds to the actual value of the QFI given by (B.2). On the other hand, note that authors of [12] choose rather $\displaystyle \frac{{\dot{\lambda }}_{i}{(\epsilon )}^{2}}{{\lambda }_{i}{(\epsilon )}^{2}-1}:= 0$, which would correspond to using an exact formula for states which remain pure (setting $\ddot{\lambda }=0$ in (B.2)), but may lead to discontinuities even for the smooth functions $\sigma (\epsilon )$. Nevertheless, to decide which convention is the right to choose, one should rather examine the validity of Cramér–Rao bound itself for such cases.

Similarly to definition (B.4), for the multi-mode formula (28) and ${\lambda }_{i}(\epsilon )={\lambda }_{j}(\epsilon )=1$ we define

$$\begin{eqnarray}&&\displaystyle \frac{{({\lambda }_{i}-{\lambda }_{j})}^{2}}{{\lambda }_{i}{\lambda }_{j}-1}(\epsilon ):= \underset{d\epsilon \to 0}{\mathrm{lim}}\displaystyle \frac{{({\lambda }_{i}(\epsilon +d\epsilon )-{\lambda }_{j}(\epsilon +d\epsilon ))}^{2}}{{\lambda }_{i}(\epsilon +d\epsilon ){\lambda }_{j}(\epsilon +d\epsilon )-1}=0.\end{eqnarray} \tag{ B.5 }$$

Now we describe the regularization process which allows us to derive the correct value of the quantum Fisher information. In fact, all (15), (16), (23), (26) have problems when there exists at least one symplectic eigenvalue equal to one. The basic idea is that instead of computing $H(\epsilon ):= H(\sigma (\epsilon ))$ directly we calculate the quantum Fisher information for the regularized state $H(\nu \sigma (\epsilon ))$, $\nu \gt 1$ and at the end of the computation we perform the limit $\nu \to 1$. This method must always give a finite result because by doing that we must arrive at an expression similar to equation (28). This becomes clear when we look at the derivation in appendix D. Also, the result it gives is certainly correct for any point where the formula is defined, we just need to make sure that this method gives the same definitions (B.4), (B.5) for the problematic points. For ${\lambda }_{i}(\epsilon )={\lambda }_{j}(\epsilon )=1$ and differentiable σ we have ${\dot{\lambda }}_{i}(\epsilon )={\dot{\lambda }}_{j}(\epsilon )=0$ and

$$\begin{eqnarray}&&\underset{\nu \to 1}{\mathrm{lim}}\;\displaystyle \frac{{(\nu {\dot{\lambda }}_{i}(\epsilon ))}^{2}}{{(\nu {\lambda }_{i}(\epsilon ))}^{2}-1}=\underset{\nu \to 1}{\mathrm{lim}}\;\displaystyle \frac{0}{{\nu }^{2}-1}=0,\end{eqnarray} \tag{ B.6a }$$

$$\begin{eqnarray}&&\underset{\nu \to 1}{\mathrm{lim}}\;\displaystyle \frac{{(\nu {\lambda }_{i}(\epsilon )-\nu {\lambda }_{j}(\epsilon ))}^{2}}{\nu {\lambda }_{i}(\epsilon )\nu {\lambda }_{j}(\epsilon )-1}=\underset{\nu \to 1}{\mathrm{lim}}\;\displaystyle \frac{0}{{\nu }^{2}-1}=0.\end{eqnarray} \tag{ B.6b }$$

Although the second equation gives the same definition as required by equation (B.5), we can see that the first equation does not give the same definition as equation (B.4) because the problematic term here is set to zero by the regularization procedure. As a result, for every such that ${\lambda }_{i}(\epsilon )=1$ we need to add additional ${\ddot{\lambda }}_{i}(\epsilon )$ to the regularized version of the quantum Fisher information. Therefore, we have shown that we can calculate the quantum Fisher information at the points of purity as a limit

$$\begin{eqnarray}&&H(\epsilon )=\underset{\nu \to 1}{\mathrm{lim}}\;H(\nu \sigma (\epsilon ))+\underset{i:{\lambda }_{i}(\epsilon )=1}{\displaystyle \sum }{\ddot{\lambda }}_{i}(\epsilon ).\end{eqnarray} \tag{ B.7 }$$

Here we prove the bound on the remainder of the general multi-mode formula. We consider the Williamson decomposition $\sigma ={{SDS}}^{\dagger }$. An element of the sum equation (22) can be written as

$$\begin{eqnarray}&&{a}_{n}=\mathrm{tr}\left[{A}^{-n}\dot{A}{A}^{-n}\dot{A}\right]=\mathrm{tr}\left[{D}^{-n}{{BD}}^{-n}B\right],\end{eqnarray} \tag{ C.1 }$$

where $B={S}^{\dagger }\dot{A}{({S}^{\dagger })}^{-1}{K}^{-n-1}$. We can derive the inequalities

$$\begin{eqnarray}\left|{a}_{n}\right| & = & \left|\underset{k,l}{\displaystyle \sum }\displaystyle \frac{1}{{\lambda }_{k}^{n}{\lambda }_{l}^{n}}{B}_{{kl}}{B}_{{lk}}\right|\leqslant \left|\underset{k,l}{\displaystyle \sum }\displaystyle \frac{1}{{\lambda }_{k}^{n}{\lambda }_{l}^{n}}\left|{B}_{{kl}}\right|\left|{B}_{{lk}}\right|\right|\leqslant \displaystyle \frac{1}{{\lambda }_{\mathrm{min}}^{2n}}\left|\underset{k,l}{\displaystyle \sum }\left|{B}_{{kl}}\right|\left|{B}_{{lk}}\right|\right|\\ & & \leqslant \displaystyle \frac{1}{{\lambda }_{\mathrm{min}}^{2n}}\left(\underset{k,l}{\displaystyle \sum }{\left|{B}_{{kl}}\right|}^{2}\right)=\displaystyle \frac{1}{{\lambda }_{\mathrm{min}}^{2n}}\left(\mathrm{tr}[{B}^{\dagger }B]\right),\end{eqnarray} \tag{ C.2 }$$

where the last inequality is the Cauchy–Schwarz inequality between B_ij and B_ji considered as vectors with 4N² entries where N is number of modes, ${\lambda }_{\mathrm{min}}:= {\mathrm{min}}_{i}\{{\lambda }_{i}\}$ is the smallest symplectic eigenvalue. Defining the Hermitian matrix $C:= {S}^{\dagger }\dot{A}{KS}$ we have

$$\begin{eqnarray}&&\mathrm{tr}[{(A\dot{A})}^{2}]=\mathrm{tr}[{C}^{\dagger }{DCD}]=\underset{k,l}{\displaystyle \sum }{\left|{C}_{{kl}}\right|}^{2}{\lambda }_{k}{\lambda }_{l}\geqslant {\lambda }_{\mathrm{min}}^{2}\mathrm{tr}[{C}^{\dagger }C]={\lambda }_{\mathrm{min}}^{2}\mathrm{tr}[{B}^{\dagger }B].\end{eqnarray} \tag{ C.3 }$$

Combining (C.2) and (C.3) gives

$$\begin{eqnarray}&&\left|{a}_{n}\right|\leqslant \mathrm{tr}\left[{(A\dot{A})}^{2}\right]{\lambda }_{\mathrm{min}}^{-2n-2}.\end{eqnarray} \tag{ C.4 }$$

For ${\lambda }_{\mathrm{min}}\gt 1$ we can estimate the remainder,

$$\begin{eqnarray}&&\left|{R}_{M}\right|\leqslant \displaystyle \frac{\mathrm{tr}\left[{(A\dot{A})}^{2}\right]}{2}\underset{n=M+1}{\overset{\infty }{\displaystyle \sum }}{\lambda }_{\mathrm{min}}^{-2n-2}=\displaystyle \frac{\mathrm{tr}\left[{(A\dot{A})}^{2}\right]}{2{\lambda }_{\mathrm{min}}^{2M+2}({\lambda }_{\mathrm{min}}^{2}-1)}.\end{eqnarray} \tag{ C.5 }$$

Here we derive the general formula for the QFI given by equation (28). Denoting ${D}_{0{ii}}:= {({D}_{0})}_{{ii}}$, we start by evaluating an infinite sum from the equation (26):

$$\begin{eqnarray}\underset{n=1}{\overset{\infty }{\displaystyle \sum }}\mathrm{tr}\left[\left({D}_{0}^{2}-I\right){\left({D}_{0}^{-n}{K}^{-n}{P}_{1}\right)}^{2}\right] & = & \underset{n=1}{\overset{\infty }{\displaystyle \sum }}\mathrm{tr}\left[\left({D}_{0}^{2}-I\right){\left({D}_{0}^{-2n+1}{{KP}}_{1}\right)}^{2}\right]\\ & & +\mathrm{tr}\left[\left({D}_{0}^{2}-I\right){\left({D}_{0}^{-2n}{P}_{1}\right)}^{2}\right]\\ & = & \underset{i,j}{\displaystyle \sum }\left({D}_{0{ii}}^{2}-1\right)\left(\underset{n=1}{\overset{\infty }{\displaystyle \sum }}{({D}_{0{ii}}{D}_{0{jj}})}^{-2n+1}\right){({{KP}}_{1})}_{{ij}}{({{KP}}_{1})}_{{ji}}\\ & & +\left({D}_{0{ii}}^{2}-1\right)\left(\underset{n=1}{\overset{\infty }{\displaystyle \sum }}{({D}_{0{ii}}{D}_{0{jj}})}^{-2n}\right){({P}_{1})}_{{ij}}{({P}_{1})}_{{ji}}\\ & = & \underset{i,j}{\displaystyle \sum }\displaystyle \frac{{D}_{0{ii}}^{2}-1}{{D}_{0{ii}}^{2}{D}_{0{jj}}^{2}-1}\left({({{KP}}_{1})}_{{ij}}{({{KP}}_{1})}_{{ji}}+{D}_{0{ii}}{D}_{0{jj}}{({P}_{1})}_{{ij}}{({P}_{1})}_{{ji}}\right).\end{eqnarray} \tag{ D.1 }$$

We combine this expression with the first part $\mathrm{tr}[{P}_{1}^{2}]$ in (26):

$$\begin{eqnarray}\mathrm{tr}\left[{P}_{1}^{2}\right]-\underset{n=1}{\overset{\infty }{\displaystyle \sum }}(\cdot ) & = & \underset{i,j}{\displaystyle \sum }\displaystyle \frac{{D}_{0{ii}}^{2}\left({D}_{0{jj}}^{2}-1\right)}{{D}_{0{ii}}^{2}{D}_{0{jj}}^{2}-1}{({P}_{1})}_{{ij}}{({P}_{1})}_{{ji}}-\displaystyle \frac{{D}_{0{ii}}{D}_{0{jj}}\left({D}_{0{ii}}^{2}-1\right)}{{D}_{0{ii}}^{2}{D}_{0{jj}}^{2}-1}{({{KP}}_{1})}_{{ij}}{({{KP}}_{1})}_{{ji}}.\end{eqnarray} \tag{ D.2 }$$

Now we use definitions (5), (27), ${P}_{1}={S}^{-1}\dot{S}={{KS}}^{\dagger }K\dot{S}$, $K=I\oplus -I$ and ${D}_{0}=\mathrm{diag}({\lambda }_{1},...,{\lambda }_{N},{\lambda }_{1},...,{\lambda }_{N})$ to rewrite the expression in terms of the symplectic eigenvalues and sub-matrices R and Q,

$$\begin{eqnarray}\mathrm{tr}\left[{P}_{1}^{2}\right]-\underset{n=1}{\overset{\infty }{\displaystyle \sum }}(\cdot ) & = & \underset{i,j=1}{\overset{N}{\displaystyle \sum }}2\displaystyle \frac{{\lambda }_{i}^{2}({\lambda }_{j}^{2}-1)-{\lambda }_{i}{\lambda }_{j}({\lambda }_{i}^{2}-1)}{{\lambda }_{i}^{2}{\lambda }_{j}^{2}-1}\mathrm{Re}({R}_{{ij}}{R}_{{ji}})\\ & & +2\displaystyle \frac{{\lambda }_{i}^{2}({\lambda }_{j}^{2}-1)+{\lambda }_{i}{\lambda }_{j}({\lambda }_{i}^{2}-1)}{{\lambda }_{i}^{2}{\lambda }_{j}^{2}-1}\mathrm{Re}({\bar{Q}}_{{ij}}{Q}_{{ji}})\\ & = & \underset{i,j=1}{\overset{N}{\displaystyle \sum }}2\displaystyle \frac{{\lambda }_{i}({\lambda }_{j}-{\lambda }_{i})}{{\lambda }_{i}{\lambda }_{j}-1}\mathrm{Re}({R}_{{ij}}{R}_{{ji}})+2\displaystyle \frac{{\lambda }_{i}({\lambda }_{j}+{\lambda }_{i})}{{\lambda }_{i}{\lambda }_{j}+1}\mathrm{Re}({\bar{Q}}_{{ij}}{Q}_{{ji}}).\end{eqnarray} \tag{ D.3 }$$

$\mathrm{Re}({\bar{Q}}_{{ij}}{Q}_{{ji}})$ and $\mathrm{Re}({R}_{{ij}}{R}_{{ji}})$ are both symmetric. This is why we decompose parts consisting of symplectic eigenvalues into its symmetric and antisymmetric parts: $2\displaystyle \frac{{\lambda }_{i}({\lambda }_{j}-{\lambda }_{i})}{{\lambda }_{i}{\lambda }_{j}-1}=-\displaystyle \frac{{({\lambda }_{j}-{\lambda }_{i})}^{2}}{{\lambda }_{i}{\lambda }_{j}-1}+\displaystyle \frac{{\lambda }_{j}^{2}-{\lambda }_{i}^{2}}{{\lambda }_{i}{\lambda }_{j}-1}$, $2\displaystyle \frac{{\lambda }_{i}({\lambda }_{j}+{\lambda }_{i})}{{\lambda }_{i}{\lambda }_{j}+1}=\displaystyle \frac{{({\lambda }_{j}+{\lambda }_{i})}^{2}}{{\lambda }_{i}{\lambda }_{j}+1}+\displaystyle \frac{{\lambda }_{i}^{2}-{\lambda }_{j}^{2}}{{\lambda }_{i}{\lambda }_{j}+1}$ and the antisymmetric part vanishes under the sum. In addition, from identity (11a) we have ${R}^{\dagger }=-R$ and ${Q}^{T}=Q$ giving us $\mathrm{Re}({R}_{{ij}}{R}_{{ji}})=-{\left|{R}_{{ij}}\right|}^{2}$, $\mathrm{Re}({\bar{Q}}_{{ij}}{Q}_{{ji}})={\left|{Q}_{{ij}}\right|}^{2}$. The diagonal part is simply $\displaystyle \frac{1}{2}\mathrm{tr}\left[{({D}_{0}+K)}^{-1}{D}_{0}^{-1}{D}_{1}^{2}\right]={\displaystyle \sum }_{i=1}^{N}\displaystyle \frac{{\dot{{\lambda }_{i}}}^{2}}{{\lambda }_{i}^{2}-1}$. Combining all these yields

$$\begin{eqnarray}H(\epsilon ) & = & \underset{i,j=1}{\overset{N}{\displaystyle \sum }}\displaystyle \frac{{({\lambda }_{i}-{\lambda }_{j})}^{2}}{{\lambda }_{i}{\lambda }_{j}-1}{\left|{R}_{{ij}}\right|}^{2}+\displaystyle \frac{{({\lambda }_{i}+{\lambda }_{j})}^{2}}{{\lambda }_{i}{\lambda }_{j}+1}{\left|{Q}_{{ij}}\right|}^{2}+\underset{i=1}{\overset{N}{\displaystyle \sum }}\displaystyle \frac{{\dot{{\lambda }_{i}}}^{2}}{{\lambda }_{i}^{2}-1}+2{\dot{{\boldsymbol{d}}}}^{\dagger }{\sigma }^{-1}\dot{{\boldsymbol{d}}}.\end{eqnarray} \tag{ D.4 }$$

The definitions of expressions for problematic points ${\lambda }_{i}(\epsilon )={\lambda }_{j}(\epsilon )=1$ come from the assumption of the continuity of the QFI and correspond to the knowledge gained in appendix B.

Here, to illustrate how the general formula (28) works, we calculate the quantum Fisher information for an example of a one-mode Gaussian state. For N = 1 the general formula (28) has a compact form

$$\begin{eqnarray}&&H(\epsilon )=4\displaystyle \frac{{\lambda }_{1}^{2}}{{\lambda }_{1}^{2}+1}{\left|{Q}_{11}\right|}^{2}+\displaystyle \frac{{\dot{{\lambda }_{1}}}^{2}}{{\lambda }_{1}^{2}-1}+2{\dot{{\boldsymbol{d}}}}^{\dagger }{\sigma }^{-1}\dot{{\boldsymbol{d}}},\end{eqnarray} \tag{ E.1 }$$

where ${Q}_{11}=\bar{\alpha }\dot{\beta }-\beta \dot{\bar{\alpha }}$.

We consider a task of estimating a squeezing parameter of a squeezing channel, with a squeezed rotated displaced thermal state as an input state. In the complex form of the covariance matrix formalism, the one-mode squeezing operator with a squeezing parameter r, and the rotation operator via angle θ are

$$\begin{eqnarray}{S}_{r}=\left[\begin{array}{cc}\mathrm{cosh}r & -\mathrm{sinh}r\\ -\mathrm{sinh}r & \mathrm{cosh}r\end{array}\right],\quad {R}_{\theta }=\left[\begin{array}{cc}{{\rm{e}}}^{-{\rm{i}}\theta } & 0\\ 0 & {{\rm{e}}}^{+{\rm{i}}\theta }\end{array}\right].\end{eqnarray} \tag{ E.2 }$$

First, we create an input state by squeezing, rotating, and displacing a thermal state. We obtain a state given by its first and the second moments,

$$\begin{eqnarray}&&{{\boldsymbol{d}}}_{0},\quad {\sigma }_{0}={R}_{\theta }{S}_{r}{{DS}}_{r}^{\dagger }{R}_{\theta }^{\dagger },\end{eqnarray} \tag{ E.3 }$$

where ${{\boldsymbol{d}}}_{0}={(\tilde{{d}_{0}},\bar{\tilde{{d}_{0}}})}^{T}$ is the initial displacement, ${\sigma }_{0}$ the initial covariance matrix, and $D=\mathrm{diag}({\lambda }_{1},{\lambda }_{1})$ is the covariance matrix of a thermal state. For a harmonic oscillator with frequency ω we have ${\lambda }_{1}=\mathrm{coth}(\displaystyle \frac{\omega {\hslash }}{2{kT}})$. Now we feed the prepared state into the channel which we consider to be again a simple squeezing operation with the unknown squeezing parameter we want to estimate. We obtain the final state,

$$\begin{eqnarray}&&{\boldsymbol{d}}={S}_{\epsilon }{{\boldsymbol{d}}}_{0},\quad \sigma ={S}_{\epsilon }{R}_{\theta }{S}_{r}{{DS}}_{r}^{\dagger }{R}_{\theta }^{\dagger }{S}_{\epsilon }^{\dagger }.\end{eqnarray} \tag{ E.4 }$$

It is clear from the construction that the diagonalizing symplectic matrix needed to compute Q₁₁ in formula (E.1) is $S={S}_{\epsilon }{R}_{\theta }{S}_{r}$. Considering definitions (5) and (E.2) we derive

$$\begin{eqnarray}&&\alpha ={{\rm{e}}}^{-{\rm{i}}\theta }\mathrm{cosh}\epsilon \mathrm{cosh}r+{{\rm{e}}}^{+{\rm{i}}\theta }\mathrm{sinh}\epsilon \mathrm{sinh}r,\end{eqnarray} \tag{ E.5a }$$

$$\begin{eqnarray}&&\beta ={{\rm{e}}}^{-{\rm{i}}\theta }\mathrm{cosh}\epsilon \mathrm{sinh}r+{{\rm{e}}}^{+{\rm{i}}\theta }\mathrm{sinh}\epsilon \mathrm{cosh}r.\end{eqnarray} \tag{ E.5b }$$

For the part containing displacement we need to compute the inverse of the covariance matrix, ${\sigma }^{-1}={{KSD}}^{-1}{S}^{\dagger }K$, and the change of displacement, $\dot{{\boldsymbol{d}}}={\dot{S}}_{\epsilon }{{\boldsymbol{d}}}_{0}$. The middle part of expression (E.1) is simply zero, because ${\lambda }_{1}$ is independent of $\epsilon$. In other words, the purity of the system remains the same. Now we have everything prepared to calculate the quantum Fisher information,

$$\begin{eqnarray}H(\epsilon ) & = & \displaystyle \frac{4{\lambda }_{1}^{2}}{{\lambda }_{1}^{2}+1}\left({\mathrm{cosh}}^{4}r+{\mathrm{sinh}}^{4}r-2\mathrm{cos}(4\theta ){\mathrm{cosh}}^{2}r{\mathrm{sinh}}^{2}r\right)\\ & & +\displaystyle \frac{4{\left|\tilde{{d}_{0}}\right|}^{2}}{{\lambda }_{1}}(\mathrm{cosh}(2r)+\mathrm{cos}(2(\theta -\phi ))\mathrm{sinh}(2r)).\end{eqnarray} \tag{ E.6 }$$