0em1.6em.

# Quantum Fisher information matrix and multiparameter estimation

Jing Liu1, Haidong Yuan2, Xiao-Ming Lu3, Xiaoguang Wang4 1MOE Key Laboratory
of Fundamental Physical Quantities Measurement & Hubei Key Laboratory of
Gravitation and Quantum Physics, PGMF and School of Physics, Huazhong
University of Science and Technology, Wuhan 430074, P. R. China  2Department
of Mechanical and Automation Engineering, The Chinese University of Hong Kong,
Shatin, Hong Kong  3Department of Physics, Hangzhou Dianzi University,
Hangzhou 310018, China  4Zhejiang Institute of Modern Physics, Department of
Physics, Zhejiang University, Hangzhou 310027, China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Quantum Fisher information matrix (QFIM) is a core concept in theoretical
quantum metrology due to the significant importance of quantum Cramér-Rao
bound in quantum parameter estimation. However, studies in recent years have
revealed wide connections between QFIM and other aspects of quantum mechanics,
including quantum thermodynamics, quantum phase transition, entanglement
witness, quantum speed limit and non-Markovianity. These connections indicate
that QFIM is more than a concept in quantum metrology, but rather a
fundamental quantity in quantum mechanics. In this paper, we summarize the
properties and existing calculation techniques of QFIM for various cases, and
review the development of QFIM in some aspects of quantum mechanics apart from
quantum metrology. On the other hand, as the main application of QFIM, the
second part of this paper reviews the quantum multiparameter Cramér-Rao bound,
its attainability condition and the associated optimal measurements. Moreover,
recent developments in a few typical scenarios of quantum multiparameter
estimation and the quantum advantages are also thoroughly discussed in this
part.

###### Contents

  1. 1 Introduction
  2. 2 Quantum Fisher information matrix
    1. 2.1 Definition
    2. 2.2 Parameterization processes
    3. 2.3 Calculating QFIM
      1. 2.3.1 General methods
      2. 2.3.2 Pure states
      3. 2.3.3 Few-qubit states
      4. 2.3.4 Unitary processes
      5. 2.3.5 Gaussian states
    4. 2.4 QFIM and geometry of quantum mechanics
      1. 2.4.1 Fubini-Study metric
      2. 2.4.2 Fidelity and Bures metric
      3. 2.4.3 Quantum geometric tensor
    5. 2.5 QFIM and thermodynamics
    6. 2.6 QFIM in quantum dynamics
      1. 2.6.1 Quantum speed limit
      2. 2.6.2 Non-Markovianity
  3. 3 Quantum multiparameter estimation
    1. 3.1 Quantum multiparameter Cramér-Rao bound
      1. 3.1.1 Introduction
      2. 3.1.2 Attainability
      3. 3.1.3 Optimal measurements
    2. 3.2 Phase estimation in the Mach-Zehnder interferometer
      1. 3.2.1 Double-phase estimation
      2. 3.2.2 Multi-phase estimation
    3. 3.3 Waveform estimation
    4. 3.4 Control-enhanced multiparameter estimation
    5. 3.5 Estimation of a magnetic field
    6. 3.6 Other applications and alternative mathematical tools
  4. 4 Conclusion and outlook
  5. A Derivation of traditional form of QFIM
  6. B Derivation of QFIM for arbitrary-rank density matrices
  7. C QFIM in Bloch representation
  8. D One-qubit basis-independent expression of QFIM
  9. E Derivation of SLD operator for Gaussian states
    1. E.1 SLD operator for multimode Gaussian states
    2. E.2 SLD operator for single-mode Gaussian state
  10. F QFIM and Bures metric
  11. G Relation between QFIM and cross-correlation functions
  12. H Derivation of quantum multiparameter Cramér-Rao bound
  13. I Construction of optimal measurement for pure states
  14. J Gradient in GRAPE for Hamiltonian estimation
    1. J.1 Gradient of CFIM
    2. J.2 Gradient of QFIM

##  1 Introduction

After decades of rapid development, quantum mechanics has now gone deep into
almost every corner of modern science, not only as a fundamental theory, but
also as a technology. The technology originated from quantum mechanics is
usually referred to as quantum technology, which is aiming at developing brand
new technologies or improving current existing technologies with the
association of quantum resources, quantum phenomena or quantum effects. Some
aspects of quantum technology, such as quantum communications, quantum
computation, quantum cryptography and quantum metrology, have shown great
power in both theory and laboratory to lead the next industrial revolution.
Among these aspects, quantum metrology is the most promising one that can step
into practice in the near future.

Quantum metrology focuses on making high precision measurements of given
parameters using quantum systems and quantum resources. Generally, a complete
quantum metrological process contains four steps: (1) preparation of the probe
state; (2) parameterzation; (3) measurement and (4) classical estimation, as
shown in figure 1. The last step has been well studies in classical
statistics, hence, the major concern of quantum metrology is the first three
steps.

Figure 1: Schematic of a complete quantum metrological process, which contains
four steps: (1) preparation of the probe state; (2) parameterzation; (3)
measurement; (4) classical estimation.

Quantum parameter estimation is the theory for quantum metrology, and quantum
Cramér-Rao bound is the most well-studied mathematical tool for quantum
parameter estimation [1, 2]. In quantum Cramér-Rao bound, the quantum Fisher
information (QFI) and quantum Fisher information matrix (QFIM) are the key
quantities representing the precision limit for single parameter and
multiparameter estimations. In recent years, several outstanding reviews on
quantum metrology and quantum parameter estimation have been provided from
different perspectives and at different time, including the ones given by
Giovannetti et al. on the quantum-enhanced measurement [3] and the advances in
quantum metrology [4], the ones given by Paris [5] and Toth et al. [6] on the
QFI and its applications in quantum metrology, the one by Braun et al. on the
quantum enhanced metrology without entanglement [7], the ones by Pezzè et al.
[8] and Huang et al. [9] on quantum metrology with cold atoms, the one by
Degen et al. on general quantum sensing [10], the one by Pirandola et al. on
the photonic quantum sensing [11], the ones by Dowling on quantum optical
metrology with high-N00N state [12] and Dowling and Seshadreesan on
theoretical and experimental optical technologies in quantum metrology,
sensing and imaging [13], the one by Demkowicz-Dobrzański et al. on the
quantum limits in optical interferometry [14], the one by Sidhu and Kok on
quantum parameter estimation from a geometric perspective [15], and the one by
Szczykulska et al. on simultaneous multiparameter estimation [16]. Petz et al.
also wrote a thorough technical introduction on QFI [17].

Apart from quantum metrology, the QFI also connects to other aspects of
quantum physics, such as quantum phase transition [18, 19, 20] and
entanglement witness [21, 22]. The widespread application of QFI may be due to
its connection to the Fubini-study metric, a Kähler metric in the complex
projective Hilbert space. This relation gives the QFI a strong geometric
meaning and makes it a fundamental quantity in quantum physics. Similarly, the
QFIM shares this connection since the diagonal entries of QFIM simply gives
the QFI. Moreover, the QFIM also connects to other fundamental quantity like
the quantum geometric tensor [23]. Thus, besides the role in multiparameter
estimation, the QFIM should also be treated as a fundamental quantity in
quantum mechanics.

In recent years, the calculation techniques of QFIM have seen a rapid
development in various scenarios and models. However, there lack papers that
thoroughly summarize these techniques in a structured manner for the
community. Therefore, this paper not only reviews the recent developments of
quantum multiparameter estimation, but also provides comprehensive techniques
on the calculation of QFIM in a variety of scenarios. For this purpose, this
paper is presented in a way similar to a textbook with many technical details
given in the appendices, which could help the readers to follow and better
understand the corresponding results.

##  2 Quantum Fisher information matrix

###  2.1 Definition

Consider a vector of parameters
$\vec{x}=(x_{0},x_{1},...,x_{a},...)^{\mathrm{T}}$ with $x_{a}$ the $a$th
parameter. $\vec{x}$ is encoded in the density matrix $\rho=\rho(\vec{x})$. In
the entire paper we denote the QFIM as $\mathcal{F}$, and an entry of
$\mathcal{F}$ is defined as [1, 2]

$\mathcal{F}_{ab}:=\frac{1}{2}\mathrm{Tr}\left(\rho\\{L_{a},L_{b}\\}\right),$
(1)

where $\\{\cdot,\cdot\\}$ represents the anti-commutation and $L_{a}$
($L_{b}$) is the symmetric logarithmic derivative (SLD) for the parameter
$x_{a}$ ($x_{b}$), which is determined by the equation 111In the entire paper
the notation $\partial_{a}$ ($\partial_{t}$) is used as an abbreviation of
$\frac{\partial}{\partial x_{a}}$ $\left(\frac{\partial}{\partial t}\right)$.

$\partial_{a}\rho=\frac{1}{2}\left(\rho L_{a}+L_{a}\rho\right).$ (2)

The SLD operator is a Hermitian operator and the expected value
$\mathrm{Tr}(\rho L_{a})=0$. Utilizing the equation above, $\mathcal{F}_{ab}$
can also be expressed by [24]

$\mathcal{F}_{ab}=\mathrm{Tr}\left(L_{b}\partial_{a}\rho\right)=-\mathrm{Tr}(\rho\partial_{a}L_{b}).$
(3)

Based on equation (1), the diagonal entry of QFIM is

$\mathcal{F}_{aa}=\mathrm{Tr}\left(\rho L^{2}_{a}\right),$ (4)

which is exactly the QFI for parameter $x_{a}$.

The definition of Fisher information matrix originated from classical
statistics. For a probability distribution $\\{p(y|\vec{x})\\}$ where
$p(y|\vec{x})$ is the conditional probability for the outcome result $y$, an
entry of Fisher information matrix is defined as

$\mathcal{I}_{ab}:=\int\frac{\left[\partial_{a}p(y|\vec{x})\right]\\!\left[\partial_{b}p(y|\vec{x})\right]}{p(y|\vec{x})}~{}\mathrm{d}y.$
(5)

For discrete outcome results, it becomes
$\mathcal{I}_{ab}:=\sum_{y}\frac{[\partial_{a}p(y|\vec{x})][\partial_{b}p(y|\vec{x})]}{p(y|\vec{x})}$.
With the development of quantum metrology, the Fisher information matrix
concerning classical probability distribution is usually referred to as
_classical Fisher information matrix_ (CFIM), with the diagonal entry referred
to as  _classical Fisher information_ (CFI). In quantum mechanics, it is well
known that the choice of measurement will affect the obtained probability
distribution, and thus result in different CFIM. This fact indicates the CFIM
is actually a function of measurement. However, while the QFI is always
attained by optimizing over the measurements [30], i.e.,
$\mathcal{F}_{aa}=\max_{\\{\Pi_{y}\\}}\mathcal{I}_{aa}(\rho,\\{\Pi_{y}\\})$,
where $\\{\Pi_{y}\\}$ represents a positive-operator valued measure (POVM), in
general there may not be any measurement that can attain the QFIM.

The QFIM based on SLD is not the only quantum version of CFIM. Another well-
used ones are based on the right and left logarithmic derivatives [2, 25],
defined by $\partial_{a}\rho=\rho R_{a}$ and
$\partial_{a}\rho=R^{\dagger}_{a}\rho$, with the corresponding QFIM
$\mathcal{F}_{ab}=\mathrm{Tr}(\rho R_{a}R^{\dagger}_{b})$. Different with the
one based on SLD, which are real symmetric, the QFIM based on right and left
logarithmic derivatives are complex and Hermitian. All versions of QFIMs
belong to a family of Riemannian monotone metrics established by Petz [26, 27]
in 1996, which will be further discussed in section 2.4. All the QFIMs can
provide quantum versions of Cramér-Rao bound, yet with different
achievability. For instance, for the D-invariant models only the one based on
right logarithmic derivative provides an achievable bound [28]. The quantum
Cramér-Rao bound will be further discussed in section 3. For pure states,
Fujiwara and Nagaoka [29] also extended the SLD to a family via
$\partial_{a}\rho=\frac{1}{2}(\rho L_{a}+L^{\dagger}_{a}\rho)$, in which
$L_{a}$ is not necessarily to be Hermitian, and when it is, it reduces to the
SLD. An useful example here is the anti-symmetric logarithmic derivative
$L^{\dagger}_{a}=-L_{a}$. This paper focuses on the QFIM based on the SLD,
thus, the QFIM in the following only refers to the QFIM based on SLD without
causing any confusion.

The properties of QFI have been well organized by G. Tóth  _et al._ in
reference [6]. Similarly, the QFIM also has some powerful properties that have
been widely applied in practice. Here we organize these properties as below.

###### Proposition 2.1

Properties and useful formulas of the QFIM.

  * •

$\mathcal{F}$ is real symmetric, i.e.,
$\mathcal{F}_{ab}=\mathcal{F}_{ba}\in\mathbb{R}$ 111$\mathbb{R}$ represents
the set of real numbers..

  * •

$\mathcal{F}$ is positive semi-definite, i.e., $\mathcal{F}\geq 0$. If
$\mathcal{F}>0$, then $[\mathcal{F}^{-1}]_{aa}\geq 1/\mathcal{F}_{aa}$ for any
$a$.

  * •

$\mathcal{F}(\rho)=\mathcal{F}(U\rho U^{\dagger})$ for a $\vec{x}$-independent
unitary operation $U$.

  * •

If $\rho=\bigotimes_{i}\rho_{i}(\vec{x})$, then
$\mathcal{F}(\rho)=\sum_{i}\mathcal{F}(\rho_{i})$.

  * •

If $\rho=\bigoplus_{i}\mu_{i}\rho_{i}(\vec{x})$ with $\mu_{i}$ a
$\vec{x}$-independent weight, then
$\mathcal{F}(\rho)=\sum_{i}\mu_{i}\mathcal{F}(\rho_{i})$.

  * •

Convexity: $\mathcal{F}(p\rho_{1}+(1-p)\rho_{2})\leq
p\mathcal{F}(\rho_{1})+(1-p)\mathcal{F}(\rho_{2})$ for $p\in[0,1]$.

  * •

$\mathcal{F}$ is monotonic under completely positive and trace preserving map
$\Phi$, i.e., $\mathcal{F}(\Phi(\rho))\leq\mathcal{F}(\rho)$ [27].

  * •

If $\vec{y}$ is function of $\vec{x}$, then the QFIMs with respect to
$\vec{y}$ and $\vec{x}$ satisfy
$\mathcal{F}(\rho(\vec{x}))=J^{\mathrm{T}}\mathcal{F}(\rho(\vec{y}))J$, with
$J$ the Jacobian matrix, i.e., $J_{ij}=\partial y_{i}/\partial x_{j}$.

###  2.2 Parameterization processes

Generally, the parameters are encoded into the probe state via a parameter-
dependent dynamics. According to the types of dynamics, there exist three
types of parameterization processes: Hamiltonian parameterization, channel
parameterization and hybrid parameterization, as shown in figure 2. In the
Hamiltonian parameterization, $\vec{x}$ is encoded in the probe state
$\rho_{0}$ through the Hamiltonian $H_{\vec{x}}$. The dynamics is then
governed by the Schrödinger equation

$\partial_{t}\rho=-i[H_{\vec{x}},\rho],$ (6)

and the parameterized state can be written as

$\rho=e^{-iH_{\vec{x}}t}\rho_{0}e^{iH_{\vec{x}}t}.$ (7)

Thus, the Hamiltonian parameterization is a unitary process. In some other
scenarios the parameters are encoded via the interaction with another system,
which means the probe system here has to be treated as an open system and the
dynamics is governed by the master equation. This is the channel
parameterization. The dynamics for the channel parameterization is

$\partial_{t}\rho=-i[H,\rho]+\mathcal{L}_{\vec{x}}(\rho),$ (8)

where $\mathcal{L}_{\vec{x}}(\rho)$ represents the decay term dependent on
$\vec{x}$. A well-used form of $\mathcal{L}_{\vec{x}}$ is the Lindblad form

$\mathcal{L}(\rho)=\sum_{j}\gamma_{j}\left(\Gamma_{j}\rho\Gamma^{\dagger}_{j}-\frac{1}{2}\left\\{\Gamma^{\dagger}_{j}\Gamma_{j},\rho\right\\}\right),$
(9)

where $\Gamma_{j}$ is the $j$th Lindblad operator and $\gamma_{j}$ is the
$j$th decay rate. All the decay rates are unknown parameters to be estimated.
The third type is the hybrid parameterization, in which both the Hamiltonian
parameters and decay rates in equation (8) are unknown and need to be
estimated.

Figure 2: The schematic of multiparameter parameterization processes. (a)
Hamiltonian parameterization (b) Channel parameterization (c) Hybrid
parameterization.

###  2.3 Calculating QFIM

In this section we review the techniques in the calculation of QFIM and some
analytic results for specific cases.

####  2.3.1 General methods

The traditional derivation of QFIM usually assumes the rank of the density
matrix is full, i.e., all the eigenvalues of the density matrix are positive.
Specifically if we write
$\rho=\sum^{\dim(\rho)-1}_{i=0}\lambda_{i}|\lambda_{i}\rangle\langle\lambda_{i}|$,
with $\lambda_{i}$ and $|\lambda_{i}\rangle$ the eigenvalue and the
corresponding eigenstate, it is usually assumed that $\lambda_{i}>0$ for all
$0\leq i\leq\dim(\rho)-1$. Under this assumption the QFIM can be obtained as
follows.

###### Theorem 2.1

The entry of QFIM for a full-rank density matrix with the spectral
decomposition
$\rho=\sum^{d-1}_{i=0}\lambda_{i}|\lambda_{i}\rangle\langle\lambda_{i}|$ can
be written as

$\mathcal{F}_{ab}=\sum^{d-1}_{i,j=0}\frac{2\mathrm{Re}(\langle\lambda_{i}|\partial_{a}\rho|\lambda_{j}\rangle\langle\lambda_{j}|\partial_{b}\rho|\lambda_{i}\rangle)}{\lambda_{i}+\lambda_{j}},$
(10)

where $d:=\dim(\rho)$ is the dimension of the density matrix.

One can easily see that if the density matrix is not of full rank, there can
be divergent terms in the above equation. To extend it to the general density
matrices which may not have full rank, we can manually remove the divergent
terms as

$\mathcal{F}_{ab}=\sum^{d-1}_{i,j=0,\lambda_{i}+\lambda_{j}\neq
0}\frac{2\mathrm{Re}(\langle\lambda_{i}|\partial_{a}\rho|\lambda_{j}\rangle\langle\lambda_{j}|\partial_{b}\rho|\lambda_{i}\rangle)}{\lambda_{i}+\lambda_{j}}.$
(11)

By substituting the spectral decomposition of $\rho$ into the equation above,
it can be rewritten as [5]

$\mathcal{F}_{ab}=\\!\\!\\!\sum^{d-1}_{i=0}\\!\frac{(\partial_{a}\lambda_{i})(\partial_{b}\lambda_{i})}{\lambda_{i}}+\\!\\!\\!\sum_{i\neq
j,\lambda_{i}+\lambda_{j}\neq
0}\\!\\!\\!\\!\frac{2(\lambda_{i}-\lambda_{j})^{2}}{\lambda_{i}+\\!\lambda_{j}}\mathrm{Re}\\!\left(\langle\lambda_{i}|\partial_{a}\lambda_{j}\rangle\langle\partial_{b}\lambda_{j}|\lambda_{i}\rangle\right)\\!.$
(12)

Recently, it has been rigorously proved that the QFIM for a finite dimensional
density matrix can be expressed with the support of the density matrix [31].
The support of a density matrix, denoted by $\mathcal{S}$, is defined as
$\mathcal{S}:=\\{\lambda_{i}\in\\{\lambda_{i}\\}|\lambda_{i}\neq 0\\}$
($\\{\lambda_{i}\\}$ is the full set of $\rho$’s eigenvalues), and the
spectral decomposition can then be modified as
$\rho=\sum_{\lambda_{i}\in\mathcal{S}}\lambda_{i}|\lambda_{i}\rangle\langle\lambda_{i}|$.
The QFIM can then be calculated via the following theorem.

###### Theorem 2.2

Given the spectral decomposition of a density matrix,
$\rho=\sum_{\lambda_{i}\in\mathcal{S}}|\lambda_{i}\rangle\langle\lambda_{i}|$
where $\mathcal{S}=\\{\lambda_{i}\in\\{\lambda_{i}\\}|\lambda_{i}\neq 0\\}$ is
the support, an entry of QFIM can be calculated as [31]

$\displaystyle\mathcal{F}_{ab}$ $\displaystyle=$
$\displaystyle\sum_{\lambda_{i}\in\mathcal{S}}\frac{(\partial_{a}\lambda_{i})(\partial_{b}\lambda_{i})}{\lambda_{i}}+\sum_{\lambda_{i}\in\mathcal{S}}4\lambda_{i}\mathrm{Re}\left(\langle\partial_{a}\lambda_{i}|\partial_{b}\lambda_{i}\rangle\right)$
(13)
$\displaystyle-\sum_{\lambda_{i},\lambda_{j}\in\mathcal{S}}\frac{8\lambda_{i}\lambda_{j}}{\lambda_{i}+\lambda_{j}}\mathrm{Re}(\langle\partial_{a}\lambda_{i}|\lambda_{j}\rangle\langle\lambda_{j}|\partial_{b}\lambda_{i}\rangle).$

The detailed derivation of this equation can be found in B. It is a general
expression of QFIM for a finite-dimensional density matrix of arbitrary rank.
Due to the relation between the QFIM and QFI, one can easily obtain the
following corollary.

###### Corollary 2.2.1

Given the spectral decomposition of a density matrix,
$\rho=\sum_{\lambda_{i}\in\mathcal{S}}|\lambda_{i}\rangle\langle\lambda_{i}|$,
the QFI for the parameter $x_{a}$ can be calculated as [32, 33, 34, 36, 35]

$\mathcal{F}_{aa}=\sum_{\lambda_{i}\in\mathcal{S}}\frac{(\partial_{a}\lambda_{i})^{2}}{\lambda_{i}}+\sum_{\lambda_{i}\in\mathcal{S}}4\lambda_{i}\langle\partial_{a}\lambda_{i}|\partial_{a}\lambda_{i}\rangle-\\!\\!\sum_{\lambda_{i},\lambda_{j}\in\mathcal{S}}\frac{8\lambda_{i}\lambda_{j}}{\lambda_{i}+\lambda_{j}}|\langle\partial_{a}\lambda_{i}|\lambda_{j}\rangle|^{2}\\!\\!.$
(14)

The first term in equations (12) and (13) can be viewed as the counterpart of
the classical Fisher information as it only contains the derivatives of the
eigenvalues which can be regarded as the counterpart of the probability
distribution. The other terms are purely quantum [5, 36]. The derivatives of
the eigenstates reflect the local structure of the eigenspace on $\vec{x}$.
The effect of this local structure on QFIM can be easily observed via
equations (12) and (13).

The SLD operator is important since it is not only related to the calculation
of QFIM, but also contains the information of the optimal measurements and the
attainability of the quantum Cramér-Rao bound, which will be further discussed
in sections 3.1.2 and 3.1.3. In terms of the eigen-space of $\rho$, the
entries of the SLD operator can be obtained as follows 222The derivation is in
B.

$\langle\lambda_{i}|L_{a}|\lambda_{j}\rangle=\delta_{ij}\frac{\partial_{a}\lambda_{i}}{\lambda_{i}}+\frac{2(\lambda_{j}-\lambda_{i})}{\lambda_{i}+\lambda_{j}}\langle\lambda_{i}|\partial_{a}\lambda_{j}\rangle;$
(15)

for $\lambda_{i}\in\mathcal{S}$ and $\lambda_{j}\not\in\mathcal{S}$,
$\langle\lambda_{i}|L_{a}|\lambda_{j}\rangle=-2\langle\lambda_{i}|\partial_{a}\lambda_{j}\rangle$;
and for $\lambda_{i},\lambda_{j}\not\in\mathcal{S}$,
$\langle\lambda_{i}|L_{a}|\lambda_{j}\rangle$ can take arbitrary values.
Fujiwara and Nagaoka [29, 37] first proved that this randomness does not
affect the value of QFI and all forms of SLD provide the same QFI. As a matter
of fact, this conclusion can be extended to the QFIM for any quantum state
[31, 38], i.e., the entries that can take arbitrary values do not affect the
value of QFIM. Hence, if we focus on the calculation of QFIM we can just set
them zeros. However, this randomness plays a role in the search of optimal
measurement, which will be further discussed in section 3.1.3.

In control theory, equation (2) is also referred to as the Lyapunov equation
and the solution can be obtained as [5]

$L_{a}=2\int^{\infty}_{0}e^{-\rho s}\left(\partial_{a}\rho\right)e^{\rho
s}\mathrm{d}s,$ (16)

which is independent of the representation of $\rho$. This can also be written
in an expanded form [38]

$L_{a}=-2\lim_{s\rightarrow\infty}\sum_{n=0}^{\infty}\frac{\left(-s\right)^{n+1}}{(n+1)!}\mathcal{R}^{n}_{\rho}(\partial_{a}\rho),$
(17)

here $\mathcal{R}_{\rho}(\cdot):=\\{\rho,\cdot\\}$ denotes the anti-
commutator. Using the fact that
$\mathcal{R}^{n}_{\rho}(\partial_{a}\rho)=\sum^{n}_{m=0}{n\choose
m}\rho^{m}\left(\partial_{a}\rho\right)\rho^{n-m}$, where ${n\choose
m}=\frac{n!}{m!(n-m)!}$, equation (17) can be rewritten as

$L_{a}=-2\lim_{s\rightarrow\infty}\sum_{n=0}^{\infty}\sum^{n}_{m=0}\frac{\left(-s\right)^{n+1}}{(n+1)!}{n\choose
m}\rho^{m}\left(\partial_{a}\rho\right)\rho^{n-m}.$ (18)

This form of SLD can be easy to calculate if
$\rho^{m}(\partial_{a}\rho)\rho^{n-m}$ is only non-zero for limited number of
terms or has some recursive patterns.

Recently, Safránek [39] provided another method to compute the QFIM utilizing
the density matrix in Liouville space. In Liouville space, the density matrix
is a vector containing all the entries of the density matrix in Hilbert space.
Denote $\mathrm{vec}(A)$ as the _column_ vector of $A$ in Liouville space and
$\mathrm{vec}(A)^{\dagger}$ as the conjugate transpose of $\mathrm{vec}(A)$.
The entry of $\mathrm{vec}(A)$ is $[\mathrm{vec}(A)]_{id+j}=A_{ij}$
($i,j\in[0,d-1]$). The QFIM can be calculated as follows.

###### Theorem 2.3

For a full-rank density matrix, the QFIM can be expressed by [39]

$\mathcal{F}_{ab}=2\mathrm{vec}(\partial_{a}\rho)^{\dagger}\left(\rho\otimes{\leavevmode\hbox{\small
1\kern-4.3pt\normalsize 1}}+{\leavevmode\hbox{\small 1\kern-4.3pt\normalsize
1}}\otimes\rho^{*}\right)^{-1}\mathrm{vec}(\partial_{b}\rho),$ (19)

where $\rho^{*}$ is the conjugate of $\rho$, and the SLD operator in Liouville
space, denoted by $\mathrm{vec}(L_{a})$, reads

$\mathrm{vec}(L_{a})=2\left(\rho\otimes{\leavevmode\hbox{\small
1\kern-4.3pt\normalsize 1}}+{\leavevmode\hbox{\small 1\kern-4.3pt\normalsize
1}}\otimes\rho^{*}\right)^{-1}\mathrm{vec}(\partial_{a}\rho).$ (20)

This theorem can be proved by using the facts that
$\mathrm{vec}(AB\leavevmode\hbox{\small 1\kern-4.3pt\normalsize
1})=(A\otimes\leavevmode\hbox{\small 1\kern-4.3pt\normalsize
1})\mathrm{vec}(B)=(\leavevmode\hbox{\small 1\kern-4.3pt\normalsize 1}\otimes
B^{\mathrm{T}})\mathrm{vec}(A)$ ($B^{\mathrm{T}}$ is the transpose of $B$)
[40, 41, 42] and
$\mathrm{Tr}(A^{\dagger}B)=\mathrm{vec}(A)^{\dagger}\mathrm{vec}(B)$.

Bloch representation is another well-used tool in quantum information theory.
For a $d$-dimensional density matrix, it can be expressed by

$\rho=\frac{1}{d}\left(\leavevmode\hbox{\small 1\kern-4.3pt\normalsize
1}+\sqrt{\frac{d(d-1)}{2}}\vec{r}\cdot\vec{\kappa}\right),$ (21)

where $\vec{r}=(r_{1},r_{2}...,r_{m},...)^{\mathrm{T}}$ is the Bloch vector
($|\vec{r}|^{2}\leq 1$) and $\vec{\kappa}$ is a $(d^{2}-1)$-dimensional vector
of $\mathfrak{su}(d)$ generator satisfying $\mathrm{Tr}(\kappa_{i})=0$. The
anti-commutation relation for them is
$\\{\kappa_{i},\kappa_{j}\\}=\frac{4}{d}\delta_{ij}\leavevmode\hbox{\small
1\kern-4.3pt\normalsize 1}+\sum^{d^{2}-1}_{m=1}\mu_{ijm}\kappa_{m}$, and the
commutation relation is
$\left[\kappa_{i},\kappa_{j}\right]=i\sum^{d^{2}-1}_{m=1}\epsilon_{ijm}\kappa_{m}$,
where $\mu_{ijm}$ and $\epsilon_{ijm}$ are the symmetric and antisymmetric
structure constants. Watanabe et al. recently [43, 44, 45] provided the
formula of QFIM for a general Bloch vector by considering the Bloch vector
itself as the parameters to be estimated. Here we extend their result to a
general case as the theorem below.

###### Theorem 2.4

In the Bloch representation of a $d$-dimensional density matrix, the QFIM can
be expressed by

$\mathcal{F}_{ab}=(\partial_{b}\vec{r})^{\mathrm{T}}\left(\frac{d}{2(d-1)}G-\vec{r}\,\vec{r}^{\,\mathrm{T}}\right)^{-1}\partial_{a}\vec{r},$
(22)

where $G$ is a real symmetric matrix with the entry

$G_{ij}=\frac{1}{2}\mathrm{Tr}(\rho\\{\kappa_{i},\kappa_{j}\\})=\frac{2}{d}\delta_{ij}+\sqrt{\frac{d-1}{2d}}\sum_{m}\mu_{ijm}r_{m}.$
(23)

The most well-used scenario of this theorem is single-qubit systems, in which
$\rho=(\leavevmode\hbox{\small 1\kern-4.3pt\normalsize
1}+\vec{r}\cdot\vec{\sigma})/2$ with
$\vec{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ the vector of Pauli
matrices. For a single-qubit system, we have the following corollary.

###### Corollary 2.4.1

For a single-qubit mixed state, the QFIM in Bloch representation can be
expressed by

$\mathcal{F}_{ab}=(\partial_{a}\vec{r})\cdot(\partial_{b}\vec{r})+\frac{(\vec{r}\cdot\partial_{a}\vec{r})(\vec{r}\cdot\partial_{b}\vec{r})}{1-|\vec{r}|^{2}},$
(24)

where $|\vec{r}|$ is the norm of $\vec{r}$. For a single-qubit pure state,
$\mathcal{F}_{ab}=(\partial_{a}\vec{r})\cdot(\partial_{b}\vec{r})$.

The diagonal entry of equation (24) is exactly the one given by reference
[47]. The proofs of the theorem and corollary are provided in C.

####  2.3.2 Pure states

A pure state satisfies $\rho=\rho^{2}$, i.e., the purity
$\mathrm{Tr}(\rho^{2})$ equals 1. For a pure state $|\psi\rangle$, the
dimension of the support is 1, which means only one eigenvalue is non-zero (it
has to be 1 since $\mathrm{Tr\rho=1}$), with which the corresponding
eigenstate is $|\psi\rangle$. For pure states, the QFIM can be obtained as
follows.

###### Theorem 2.5

The entries of the QFIM for a pure parameterized state
$|\psi\rangle:=|\psi(\vec{x})\rangle$ can be obtained as [1, 2]

$\mathcal{F}_{ab}=4\mathrm{Re}(\langle\partial_{a}\psi|\partial_{b}\psi\rangle-\langle\partial_{a}\psi|\psi\rangle\langle\psi|\partial_{b}\psi\rangle).$
(25)

The QFI for the parameter $x_{a}$ is just the diagonal element of the QFIM,
which is given by

$\mathcal{F}_{aa}=4(\langle\partial_{a}\psi|\partial_{a}\psi\rangle-|\langle\partial_{a}\psi|\psi\rangle|^{2}),$
(26)

and the SLD operator corresponding to $x_{a}$ is
$L_{a}=2\left(|\psi\rangle\langle\partial_{a}\psi|+|\partial_{a}\psi\rangle\langle\psi|\right)$.

The SLD formula is obtained from the fact $\rho^{2}=\rho$ for a pure state,
then $\partial_{a}\rho=\rho\partial_{a}\rho+(\partial_{a}\rho)\rho$. Compared
this equation to the definition equation, it can be seen that
$L=2\partial_{a}\rho$. A simple example is
$|\psi\rangle=e^{-i\sum_{j}H_{j}x_{j}t}|\psi_{0}\rangle$ with
$[H_{a},H_{b}]=0$ for any $a$ and $b$, here $|\psi_{0}\rangle$ denotes the
initial probe state. In this case, the QFIM reads

$\mathcal{F}_{ab}=4t^{2}\mathrm{cov}_{|\psi_{0}\rangle}(H_{a},H_{b}),$ (27)

where $\mathrm{cov}_{|\varphi\rangle}(A,B)$ denotes the covariance between $A$
and $B$ on $|\varphi\rangle$, i.e.,

$\mathrm{cov}_{|\varphi\rangle}(A,B):=\frac{1}{2}\langle\varphi|\\{A,B\\}|\varphi\rangle-\langle\varphi|A|\varphi\rangle\langle\varphi|B|\varphi\rangle.$
(28)

A more general case where $H_{a}$ and $H_{b}$ do not commute will be discussed
in section 2.3.4.

####  2.3.3 Few-qubit states

The simplest few-qubit system is the single-qubit system. A single-qubit pure
state can always be written as $\cos\theta|0\rangle+\sin\theta
e^{i\phi}|1\rangle$ ($\\{|0\rangle,|1\rangle\\}$ is the basis), i.e., it only
has two degrees of freedom, which means only two independent parameters
($\vec{x}=(x_{0},x_{1})^{\mathrm{T}}$) can be encoded in a single-qubit pure
state. Assume $\theta$, $\phi$ are the parameters to be estimated, the QFIM
can then be obtained via equation (25) as

$\mathcal{F}_{\theta\theta}=4,~{}\mathcal{F}_{\phi\phi}=\sin^{2}(2\theta),~{}\mathcal{F}_{\theta\phi}=0.$
(29)

If the unknown parameters are not $\theta,\phi$, but functions of
$\theta,\phi$, the QFIM can be obtained from formula above with the assistance
of Jacobian matrix.

For a single-qubit mixed state, when the number of encoded parameters is
larger than three, the determinant of QFIM would be zero, indicating that
these parameters cannot be simultaneously estimated. This is due to the fact
that there only exist three degrees of freedom in a single-qubit mixed state,
thus, only three or fewer independent parameters can be encoded into the
density matrix $\rho$. However, more parameters may be encoded if they are not
independent. Since $\rho$ here only has two eigenvalues $\lambda_{0}$ and
$\lambda_{1}$, equation (13) then reduces to

$\mathcal{F}_{ab}=\frac{(\partial_{a}\lambda_{0})(\partial_{b}\lambda_{0})}{\lambda_{0}(1-\lambda_{0})}+4(1-2\lambda_{0})^{2}\mathrm{Re}\left(\langle\partial_{a}\lambda_{0}|\lambda_{1}\rangle\langle\lambda_{1}|\partial_{b}\lambda_{0}\rangle\right).$
(30)

In the case of single qubit, equation (30) can also be written in a basis-
independent formula [46] below.

###### Theorem 2.6

The basis-independent expression of QFIM for a single-qubit mixed state $\rho$
is of the following form

$\mathcal{F}_{ab}=\mathrm{Tr}\left[(\partial_{a}\rho)(\partial_{b}\rho)\right]+\frac{1}{\det(\rho)}\mathrm{Tr}\left[\rho(\partial_{a}\rho)\rho(\partial_{b}\rho)\right],$
(31)

where $\det(\rho)$ is the determinant of $\rho$. For a single-qubit pure
state, $\mathcal{F}_{ab}=2\mathrm{Tr}[(\partial_{a}\rho)(\partial_{b}\rho)]$.

Equation (31) is the reduced form of the one given in reference [46]. The
advantage of the basis-independent formula is that the diagonalization of the
density matrix is avoided. Now we show an example for single-qubit. Consider a
spin in a magnetic field which is in the $z$-axis and suffers from dephasing
noise also in the $z$-axis. The dynamics of this spin can then be expressed by

$\partial_{t}\rho=-i[B\sigma_{z},\rho]+\frac{\gamma}{2}(\sigma_{z}\rho\sigma_{z}-\rho),$
(32)

where $\sigma_{z}$ is a Pauli matrix. $B$ is the amplitude of the field. Take
$B$ and $\gamma$ as the parameters to be estimated. The analytical solution
for $\rho(t)$ is

$\rho(t)=\left(\begin{array}[]{cc}\rho_{00}(0)&\rho_{01}(0)e^{-i2Bt-\gamma
t}\\\ \rho_{10}(0)e^{i2Bt-\gamma t}&\rho_{11}(0)\end{array}\right).$ (33)

The derivatives of $\rho(t)$ on both $B$ and $\gamma$ are simple in this
basis. Therefore, the QFIM can be directly calculated from equation (31),
which is a diagonal matrix ($\mathcal{F}_{B\gamma}=0$) with the diagonal
entries

$\displaystyle\mathcal{F}_{BB}$ $\displaystyle=$ $\displaystyle
16|\rho_{01}(0)|^{2}e^{-2\gamma t}t^{2},$ (34)
$\displaystyle\mathcal{F}_{\gamma\gamma}$ $\displaystyle=$
$\displaystyle\frac{4\rho_{00}(0)\rho_{11}(0)|\rho_{01}(0)|^{2}t^{2}}{\rho_{00}(0)\rho_{11}(0)e^{2\gamma
t}-|\rho_{01}(0)|^{2}}.$ (35)

For a general two-qubit state, the calculation of QFIM requires the
diagonalization of a 4 by 4 density matrix, which is difficult to solve
analytically. However, some special two-qubit states, such as the X state, can
be diagonalized analytically. An X state has the form (in the computational
basis $\\{|00\rangle,|01\rangle,|10\rangle,|11\rangle\\}$) of

$\rho=\left(\begin{array}[]{cccc}\rho_{00}&0&0&\rho_{03}\\\
0&\rho_{11}&\rho_{12}&0\\\ 0&\rho_{21}&\rho_{22}&0\\\
\rho_{30}&0&0&\rho_{33}\end{array}\right).$ (36)

By changing the basis into
$\\{|00\rangle,|11\rangle,|01\rangle,|10\rangle\\}$, this state can be
rewritten in the block diagonal form as $\rho=\rho^{(0)}\oplus\rho^{(1)}$,
where $\oplus$ represents the direct sum and

$\rho^{(0)}=\left(\begin{array}[]{cc}\rho_{00}&\rho_{03}\\\
\rho_{30}&\rho_{33}\\\
\end{array}\right),\quad\rho^{(1)}=\left(\begin{array}[]{cc}\rho_{11}&\rho_{12}\\\
\rho_{21}&\rho_{22}\\\ \end{array}\right).$ (37)

Note that $\rho^{(0)}$ and $\rho^{(1)}$ are not density matrices as their
trace is not normalized. The QFIM for this block diagonal state can be written
as $\mathcal{F}_{ab}=\mathcal{F}^{(0)}_{ab}+\mathcal{F}^{(1)}_{ab}$ [36],
where $\mathcal{F}^{(0)}_{ab}$ ($\mathcal{F}^{(1)}_{ab}$) is the QFIM for
$\rho^{(0)}$ ($\rho^{(1)}$). The eigenvalues of $\rho^{(i)}$ are
$\lambda^{(i)}_{\pm}=\frac{1}{2}\Big{(}\mathrm{Tr}\rho^{(i)}\pm\sqrt{\mathrm{Tr}^{2}\rho^{(i)}-4\det\rho^{(i)}}\Big{)}$
and corresponding eigenstates are

$|\lambda^{(i)}_{\pm}\rangle=\mathcal{N}^{(i)}_{\pm}\left(\frac{1}{2\mathrm{Tr}(\rho^{(i)}\sigma_{+})}\left[\mathrm{Tr}\left(\rho^{(i)}\sigma_{z}\right)\pm\sqrt{\mathrm{Tr}^{2}\rho^{(i)}-4\det\rho^{(i)}}\right]\\!\\!,1\right)^{\mathrm{T}},$
(38)

for non-diagonal $\rho^{(i)}$ with $\mathcal{N}^{(i)}_{\pm}$ ($i=0,1$) the
normalization coefficient. Here the specific form of $\sigma_{z}$ and
$\sigma_{+}$ are

$\sigma_{z}=\left(\begin{array}[]{cc}1&0\\\ 0&-1\\\
\end{array}\right),\quad\sigma_{+}=\left(\begin{array}[]{cc}0&1\\\ 0&0\\\
\end{array}\right).$ (39)

Based on above information, $\mathcal{F}^{(i)}_{ab}$ can be specifically
written as

$\displaystyle\mathcal{F}^{(i)}_{ab}$ $\displaystyle=$
$\displaystyle\sum_{k=\pm}\frac{(\partial_{a}\lambda^{(i)}_{k})(\partial_{b}\lambda^{(i)}_{k})}{\lambda^{(i)}_{k}}+\lambda^{(i)}_{k}\mathcal{F}_{ab}\left(|\lambda^{(i)}_{k}\rangle\right)$
(40)
$\displaystyle-\frac{16\det\rho^{(i)}}{\mathrm{Tr}\rho^{(i)}}\mathrm{Re}(\langle\partial_{a}\lambda^{(i)}_{+}|\lambda^{(i)}_{-}\rangle\langle\lambda^{(i)}_{-}|\partial_{b}\lambda^{(i)}_{+}\rangle),$

where $\mathcal{F}_{ab}(|\lambda^{(i)}_{k}\rangle)$ is the QFIM entry for the
state $|\lambda^{(i)}_{k}\rangle$. For diagonal $\rho^{(i)}$,
$|\lambda^{(i)}_{\pm}\rangle$ is just $(0,1)^{\mathrm{T}}$ and only the
classical contribution term remains in above equation.

####  2.3.4 Unitary processes

Unitary processes are the most fundamental dynamics in quantum mechanics since
it can be naturally obtained via the Schrödinger equation. For a
$\vec{x}$-dependent unitary process $U=U(\vec{x})$, the parameterized state
$\rho$ can be written as $\rho=U\rho_{0}U^{\dagger}$, where $\rho_{0}$ is the
initial probe state which is $\vec{x}$-independent. For such a process, the
QFIM can be calculated via the following theorem.

###### Theorem 2.7

For a unitary parametrization process $U$, the entry of QFIM can be obtained
as [48]

$\displaystyle\mathcal{F}_{ab}$ $\displaystyle=$
$\displaystyle\sum_{\eta_{i}\in\mathcal{S}}4\eta_{i}\mathrm{cov}_{|\eta_{i}\rangle}(\mathcal{H}_{a},\mathcal{H}_{b})$
(41) $\displaystyle-\sum_{\eta_{i},\eta_{j}\in\mathcal{S},i\neq
j}\frac{8\eta_{i}\eta_{j}}{\eta_{i}+\eta_{j}}\mathrm{Re}\left(\langle\eta_{i}|\mathcal{H}_{a}|\eta_{j}\rangle\langle\eta_{j}|\mathcal{H}_{b}|\eta_{i}\rangle\right),$

where $\eta_{i}$ and $|\eta_{i}\rangle$ are $i$th eigenvalue and eigenstate of
the initial probe state $\rho_{0}$.
$\mathrm{cov}_{|\eta_{i}\rangle}(\mathcal{H}_{a},\mathcal{H}_{b})$ is defined
in equation (28). The operator $\mathcal{H}_{a}$ is defined as [49, 50]

$\mathcal{H}_{a}:=i\left(\partial_{a}U^{\dagger}\right)U=-iU^{\dagger}\left(\partial_{a}U\right).$
(42)

$\mathcal{H}_{a}$ is a Hermitian operator for any parameter $x_{a}$ due to
above definition.

For the unitary processes, the parameterized state will remain pure for a pure
probe state. The QFIM for this case is given as follows.

###### Corollary 2.7.1

For a unitary process $U$ with a pure probe state $|\psi_{0}\rangle$, the
entry of QFIM is in the form

$\mathcal{F}_{ab}=4\mathrm{cov}_{|\psi_{0}\rangle}(\mathcal{H}_{a},\mathcal{H}_{b}),$
(43)

where $\mathrm{cov}_{|\psi_{0}\rangle}(\mathcal{H}_{a},\mathcal{H}_{b})$ is
defined by equation (28) and the QFI for $x_{a}$ can then be obtained as
$\mathcal{F}_{aa}=4\mathrm{var}_{|\psi_{0}\rangle}(\mathcal{H}_{a})$. Here
$\mathrm{var}_{|\psi_{0}\rangle}(\mathcal{H}_{a}):=\mathrm{cov}_{|\psi_{0}\rangle}(\mathcal{H}_{a},\mathcal{H}_{a})$
is the variance of $\mathcal{H}_{a}$ on $|\psi_{0}\rangle$.

For a single-qubit mixed state $\rho_{0}$ under a unitary process, the QFIM
can be written as

$\mathcal{F}_{ab}=4\left[2\mathrm{Tr}(\rho^{2}_{0})-1\right]\mathrm{cov}_{|\eta_{0}\rangle}(\mathcal{H}_{a},\mathcal{H}_{b})$
(44)

with $|\eta_{0}\rangle$ an eigenstate of $\rho_{0}$. This equation is
equivalent to

$\mathcal{F}_{ab}=4\left[2\mathrm{Tr}(\rho^{2}_{0})-1\right]\mathrm{Re}\left(\langle\eta_{0}|\mathcal{H}_{a}|\eta_{1}\rangle\langle\eta_{1}|\mathcal{H}_{b}|\eta_{0}\rangle\right).$
(45)

The diagonal entry reads
$\mathcal{F}_{aa}=4\left[2\mathrm{Tr}(\rho^{2}_{0})-1\right]|\langle\eta_{0}|\mathcal{H}_{a}|\eta_{1}\rangle|^{2}.$
Recall that Theorem 2.6 provides the basis-independent formula for single-
qubit mixed state, which leads to the next corollary.

###### Corollary 2.7.2

For a single-qubit mixed state $\rho_{0}$ under a unitary process, the basis-
independent formula of QFIM is

$\displaystyle\mathcal{F}_{ab}$ $\displaystyle=$
$\displaystyle\mathrm{Tr}(\rho^{2}_{0}\\{\mathcal{H}_{a},\mathcal{H}_{b}\\})-2\mathrm{Tr}(\rho_{0}\mathcal{H}_{a}\rho_{0}\mathcal{H}_{b})$
(46)
$\displaystyle+\frac{1}{\det\rho_{0}}\left[\mathrm{Tr}\\!\left(\rho_{0}\mathcal{H}_{a}\rho_{0}\left\\{\rho_{0}^{2},\mathcal{H}_{b}\right\\}\right)-2\mathrm{Tr}(\rho^{2}_{0}\mathcal{H}_{a}\rho^{2}_{0}\mathcal{H}_{b})\right].$

The diagonal entry reads

$\displaystyle\mathcal{F}_{aa}$ $\displaystyle=$ $\displaystyle
2\mathrm{Tr}(\rho^{2}_{0}\mathcal{H}^{2}_{a})-2\mathrm{Tr}[(\rho_{0}\mathcal{H}_{a})^{2}]$
(47)
$\displaystyle+\frac{2}{\det\rho_{0}}\left[\mathrm{Tr}\\!\left(\rho_{0}\mathcal{H}_{a}\rho^{3}_{0}\mathcal{H}_{a}\right)-\mathrm{Tr}[(\rho^{2}_{0}\mathcal{H}_{a})^{2}]\right].$

Under the unitary process, the QFIM for pure probe states, as given in
equation (1), can be rewritten as

$\mathcal{F}_{ab}=\frac{1}{2}\langle\psi_{0}|\left\\{L_{a,\mathrm{eff}},L_{b,\mathrm{eff}}\right\\}|\psi_{0}\rangle,$
(48)

where $L_{a,\mathrm{eff}}:=U^{\dagger}L_{a}U$ can be treated as an effective
SLD operator, which leads to the following theorem.

###### Theorem 2.8

Given a unitary process, $U$, with a pure probe state, $|\psi_{0}\rangle$, the
effective SLD operator $L_{a,\mathrm{eff}}$ can be obtained as

$L_{a,\mathrm{eff}}=i2\left[\mathcal{H}_{a},|\psi_{0}\rangle\langle\psi_{0}|\right].$
(49)

From equation (41), all the information of the parameters is involved in the
operator set $\\{\mathcal{H}_{a}\\}$, which might benefit the analytical
optimization of the probe state in some scenarios. Generally, the unitary
operator can be written as $\exp(-itH)$ where $H=H(\vec{x})$ is the
Hamiltonian for the parametrization. $\mathcal{H}_{a}$ can then be calculated
as

$\mathcal{H}_{a}=-\int^{t}_{0}e^{isH}\left(\partial_{a}H\right)e^{-isH}\mathrm{d}s,$
(50)

where the technique
$\partial_{x}e^{A}=\int^{1}_{0}e^{sA}\partial_{x}Ae^{(1-s)A}ds$ ($A$ is an
operator) is applied. Denote $H^{\times}(\cdot):=[H,\cdot]$, the expression
above can be rewritten in an expanded form [48]

$\mathcal{H}_{a}=-\sum^{\infty}_{n=0}\frac{t^{n+1}}{(n+1)!}\left(iH^{\times}\right)^{n}\partial_{a}H.$
(51)

In some scenarios, the recursive commutations in expression above display
certain patterns, which can lead to analytic expressions for the $\mathcal{H}$
operator. The simplest example is $H=\sum_{a}x_{a}H_{a}$, with all $H_{a}$
commute with each other. In this case $\mathcal{H}_{a}=-tH_{a}$ since only the
zeroth order term in equation (51) is nonzero. Another example is the
interaction of a collective spin system with a magnetic field with the
Hamiltonian $H=-BJ_{\vec{n}_{0}}$, where $B$ is the amplitude of the external
magnetic field, $J_{\vec{n}_{0}}=\vec{n}_{0}\cdot\vec{J}$ with
$\vec{n}_{0}=(\cos\theta,0,\sin\theta)$ and $\vec{J}=(J_{x},J_{y},J_{z})$.
$\theta$ is the angle between the field and the collective spin.
$J_{i}=\sum_{k}\sigma^{(k)}_{i}/2$ for $i=x,y,z$ is the collective spin
operator. $\sigma^{(k)}_{i}$ is the Pauli matrix for $k$th spin. In this case,
the $\mathcal{H}$ operator for $\theta$ can be analytically calculated via
equation (51), which is [48]

$\mathcal{H}_{\theta}=-2\sin\left(\frac{1}{2}Bt\right)J_{\vec{n}_{1}},$ (52)

where $J_{n_{1}}=\vec{n}_{1}\cdot\vec{J}$ with
$\vec{n}_{1}=(\cos(Bt/2)\sin\theta,\sin(Bt/2),-\cos(Bt/2)\cos\theta)$.

Recently, Sidhu and Kok [51, 52] use this $\mathcal{H}$-representation to
study the spatial deformations, epecially the grid deformations of classical
and quantum light emitters. By calculating and analyzing the QFIM, they showed
that the higher average mode occupancies of the classical states performs
better in estimating the deformation when compared with single photon
emitters.

An alternative operator that can be used to characterize the precision limit
of unitary process is [50, 53, 54]

$\mathcal{K}_{a}:=i\left(\partial_{a}U\right)U^{\dagger}=-iU\left(\partial_{a}U^{\dagger}\right).$
(53)

As a matter of fact, this operator is the infinitesimal generator of $U$ of
parameter $x_{a}$. Assume $\vec{x}$ is shifted by $dx_{a}$ along the direction
of $x_{a}$ and other parameters are kept unchanged. Then $U(\vec{x}+dx_{a})$
can be expanded as $U(\vec{x}+dx_{a})\partial_{a}U(\vec{x})$. The density
matrix $\rho_{\vec{x}+dx_{a}}$ can then be approximately calculated as
$\rho_{\vec{x}+dx_{a}}=e^{-i\mathcal{K}_{a}dx_{a}}\rho
e^{i\mathcal{K}_{a}dx_{a}}$ [53], which indicates that $\mathcal{K}_{a}$ is
the generator of $U$ along parameter $x_{a}$. The relation between
$\mathcal{H}_{a}$ and $\mathcal{K}_{a}$ can be easily obtained as

$\mathcal{K}_{a}=-U\mathcal{H}_{a}U^{\dagger}.$ (54)

With this relation, the QFIM can be easily rewritten with $\mathcal{K}_{a}$ as

$\displaystyle\mathcal{F}_{ab}$ $\displaystyle=$
$\displaystyle\sum_{\lambda_{i}\in\mathcal{S}}4\lambda_{i}\mathrm{cov}_{|\lambda_{i}\rangle}(\mathcal{K}_{a},\mathcal{K}_{b})$
(55) $\displaystyle-\sum_{\lambda_{i},\lambda_{j}\in\mathcal{S},i\neq
j}\frac{8\lambda_{i}\lambda_{j}}{\lambda_{i}+\lambda_{j}}\mathrm{Re}\left(\langle\lambda_{i}|\mathcal{K}_{a}|\lambda_{j}\rangle\langle\lambda_{j}|\mathcal{K}_{b}|\lambda_{i}\rangle\right),$

where $|\lambda_{i}\rangle=U|\eta_{i}\rangle$ is the $i$th eigenstate of the
parameterized state $\rho$. And
$\mathrm{cov}_{|\lambda_{i}\rangle}(\mathcal{K}_{a},\mathcal{K}_{b})$ is
defined by equation (28). The difference between the calculation of QFIM with
$\\{\mathcal{K}_{a}\\}$ and $\\{\mathcal{H}_{a}\\}$ is that the expectation is
taken with the eigenstate of the probe state $\rho_{0}$ for the use of
$\\{\mathcal{H}_{a}\\}$ but with the parameterized state $\rho$ for
$\\{\mathcal{K}_{a}\\}$. For a pure probe state $|\psi_{0}\rangle$, the
expression above reduces to
$\mathcal{F}_{ab}=4\mathrm{cov}_{|\psi\rangle}(\mathcal{K}_{a},\mathcal{K}_{b})$
with $|\psi\rangle=U|\psi_{0}\rangle$. Similarly, for a mixed state of single
qubit, the QFIM reads
$\mathcal{F}_{ab}=4\left(2\mathrm{Tr}\rho^{2}_{0}-1\right)\mathrm{cov}_{|\lambda_{0}\rangle}(\mathcal{K}_{a},\mathcal{K}_{b})$.

####  2.3.5 Gaussian states

Gaussian state is a widely-used quantum state in quantum physics, particularly
in quantum optics, quantum metrology and continuous variable quantum
information processes. Consider a $m$-mode bosonic system with $a_{i}$
($a^{\dagger}_{i}$) as the annihilation (creation) operator for the $i$th
mode. The quadrature operators are [55, 56]
$\hat{q}_{i}:=\frac{1}{\sqrt{2}}(a_{i}+a^{\dagger}_{i})$ and
$\hat{p}_{i}:=\frac{1}{i\sqrt{2}}(a_{i}-a^{\dagger}_{i})$, which satisfy the
commutation relation $\left[\hat{q}_{i},\hat{p}_{j}\right]=i\delta_{ij}$
($\hbar=1$). A vector of quadrature operators,
$\vec{R}=(\hat{q}_{1},\hat{p}_{1},...,\hat{q}_{m},\hat{p}_{m})^{\mathrm{T}}$
satisfies

$\left[R_{i},R_{j}\right]=i\Omega_{ij}$ (56)

for any $i$ and $j$ where $\Omega$ is the symplectic matrix defined as
$\Omega:=i\sigma^{\oplus m}_{y}$ with $\oplus$ denote the direct sum. Now we
introduce the covariance matrix $C(\vec{R})$ with the entries defined as
$C_{ij}:=\mathrm{cov}_{\rho}(R_{i},R_{j})=\frac{1}{2}\mathrm{Tr}(\rho\\{R_{i},R_{j}\\})-\mathrm{Tr}(\rho
R_{i})\mathrm{Tr}(\rho R_{j})$. $C$ satisfies the uncertainty relation
$C+\frac{i}{2}\Omega\geq 0$ [57, 58]. According to the Williamson’s theorem,
the covariance matrix can be diagonalized utilizing a symplectic matrix $S$
[58, 59], i.e.,

$C=SC_{\mathrm{d}}S^{\mathrm{T}},$ (57)

where $C_{\mathrm{d}}=\bigoplus^{m}_{k=1}c_{k}\leavevmode\hbox{\small
1\kern-4.3pt\normalsize 1}_{2}$ with $c_{k}$ the $k$th symplectic eigenvalue.
$S$ is a $2m$-dimensional real matrix which satisfies $S\Omega
S^{\mathrm{T}}=\Omega$.

A very useful quantity for Gaussian states is the characteristic function

$\chi(\vec{s})=\mathrm{Tr}\\!\left(\rho~{}e^{i\vec{R}^{\,\mathrm{T}}\Omega\vec{s}}\right),$
(58)

where $\vec{s}$ is a $2m$-dimensional real vector. Another powerful function
is the Wigner function, which can be obtained by taking the Fourier transform
of the characteristic function

$W(\vec{R})=\frac{1}{(2\pi)^{2m}}\int_{\mathbb{R}^{2m}}\\!e^{-i\vec{R}^{\,\mathrm{T}}\Omega\vec{s}}~{}\chi(\vec{s})~{}\mathrm{d}^{2m}\vec{s}.$
(59)

Considering the scenario with first and second moments, a state is a Gaussian
state if $\chi(\vec{s})$ and $W(\vec{R})$ are Gaussian, i.e., [58, 55, 56, 60,
61]

$\displaystyle\chi(\vec{s})$ $\displaystyle=$ $\displaystyle
e^{-\frac{1}{2}\vec{s}^{\,\mathrm{T}}\Omega
C\Omega^{\mathrm{T}}\vec{s}-i(\Omega\langle\vec{R}\rangle)^{\mathrm{T}}\vec{s}},$
(60) $\displaystyle W(\vec{R})$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{m}\sqrt{\det
C}}e^{-\frac{1}{2}\left(\vec{R}-\langle\vec{R}\rangle\right)^{\\!\\!\mathrm{T}}C^{-1}\left(\vec{R}-\langle\vec{R}\rangle\right)},$
(61)

where $\langle\vec{R}\rangle_{j}=\mathrm{Tr}(R_{j}\rho)$ is the first moment.
A pure state is Gaussian if and only if its Wigner function is non-negative
[58].

The study of QFIM for Gaussian states started from the research of QFI. The
expression of QFI was first given in 2013 by Monras for the multi-mode
case[62] and Pinel et al. for the single-mode case [63]. In 2018, Nichols et
al. [65] and Šafránek [66] provided the expression of QFIM for multi-mode
Gaussian states independently, which was obtained based on the calculation of
SLD [7]. The SLD operator for Gaussian states has been given in references
[62, 65, 67], and we organize the corresponding results in the following
theorem.

###### Theorem 2.9

For a continuous variable bosonic $m$-mode Gaussian state with the
displacement vector (first moment) $\langle\vec{R}\rangle$ and the covariance
matrix (second moment) $C$, the SLD operator is [62, 65, 67, 66, 64, 69]

$L_{a}=L^{(0)}_{a}\leavevmode\hbox{\small 1\kern-4.3pt\normalsize
1}_{2m}+\vec{L}^{(1),\mathrm{T}}_{a}\vec{R}+\vec{R}^{\,\mathrm{T}}G_{a}\vec{R},$
(62)

where $\leavevmode\hbox{\small 1\kern-4.3pt\normalsize 1}_{2m}$ is the
$2m$-dimensional identity matrix and the coefficients read

$\displaystyle G_{a}$ $\displaystyle=$
$\displaystyle\sum^{m}_{j,k=1}\sum^{3}_{l=0}\frac{g^{(jk)}_{l}}{4c_{j}c_{k}+(-1)^{l+1}}\left(S^{\mathrm{T}}\right)^{-1}A^{(jk)}_{l}S^{-1},$
(63) $\displaystyle\vec{L}^{(1)}_{a}$ $\displaystyle=$ $\displaystyle
C^{-1}(\partial_{a}{\langle\vec{R}\rangle})-2G_{a}\langle\vec{R}\rangle,$ (64)
$\displaystyle L^{(0)}_{a}$ $\displaystyle=$
$\displaystyle\langle\vec{R}\rangle^{\mathrm{T}}G_{a}\langle\vec{R}\rangle-(\partial_{a}\langle\vec{R}\rangle)^{\mathrm{T}}C^{-1}\langle\vec{R}\rangle-\mathrm{Tr}(G_{a}C).$
(65)

Here

$A^{(jk)}_{l}=\frac{1}{\sqrt{2}}i\sigma^{(jk)}_{y},~{}\frac{1}{\sqrt{2}}\sigma^{(jk)}_{z},~{}\frac{1}{\sqrt{2}}\leavevmode\hbox{\small
1\kern-4.3pt\normalsize 1}^{(jk)}_{2},~{}\frac{1}{\sqrt{2}}\sigma^{(jk)}_{x}$
(66)

for $l=0,1,2,3$ and
$g^{(jk)}_{l}=\mathrm{Tr}[S^{-1}(\partial_{a}C)(S^{\mathrm{T}})^{-1}A^{(jk)}_{l}]$.
$\sigma^{(jk)}_{i}$ is a $2m$-dimensional matrix with all the entries zero
expect a $2\times 2$ block, shown as below

$\sigma^{(jk)}_{i}=\left(\begin{array}[]{ccccc}&1\mathrm{st}&\cdots&k\mathrm{th}&\cdots\\\
1\mathrm{st}&0_{2\times 2}&0_{2\times 2}&0_{2\times 2}&0_{2\times 2}\\\
\vdots&0_{2\times 2}&\vdots&\vdots&\vdots\\\ j\mathrm{th}&0_{2\times
2}&\cdots&\sigma_{i}&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\vdots\end{array}\right),$ (67)

where $0_{2\times 2}$ represents a 2 by 2 block with zero entries.
$\leavevmode\hbox{\small 1\kern-4.3pt\normalsize 1}^{(jk)}_{2}$ is similar to
$\sigma^{(jk)}_{i}$ but replace the block $\sigma_{i}$ with
$\leavevmode\hbox{\small 1\kern-4.3pt\normalsize 1}_{2}$ 333The derivation of
this theorem is in E.1..

Being aware of the expression of SLD operator given in Theorem 2.9, the QFIM
can be calculated via equation (1). Here we show the result explicitly in
following theorem.

###### Theorem 2.10

For a continuous variable bosonic $m$-mode Gaussian state with the
displacement vector (first moment) $\langle\vec{R}\rangle$ and the covariance
matrix (second moment) $C$, the entry of QFIM can be expressed by [64, 65, 66]

$\mathcal{F}_{ab}=\mathrm{Tr}\left(G_{a}\partial_{b}C\right)+(\partial_{a}\langle\vec{R}\rangle^{\mathrm{T}})C^{-1}\partial_{b}\langle\vec{R}\rangle,$
(68)

and the QFI for an $m$-mode Gaussian state with respect to $x_{a}$ can be
immediately obtained as [62]

$\mathcal{F}_{aa}=\mathrm{Tr}\left(G_{a}\partial_{a}C\right)+(\partial_{a}\langle\vec{R}\rangle^{\mathrm{T}})C^{-1}\partial_{a}\langle\vec{R}\rangle.$
(69)

The expression of right logarithmic derivative for a general Gaussian state
and the corresponding QFIM was provided by Gao and Lee [64] in 2014, which is
an appropriate tool for the estimation of complex numbers [25], such as the
number $\alpha$ of a coherent state $|\alpha\rangle$. The simplest case is a
single-mode Gaussian state. For such a state, $G_{a}$ can be calculated as
following.

###### Corollary 2.10.1

For a single-mode Gaussian state $G_{a}$ can be expressed as 444The derivation
of this corollary is in E.2.

$G_{a}=\frac{4c^{2}-1}{4c^{2}+1}\Omega(\partial_{a}J)\Omega,$ (70)

where $c=\sqrt{\det C}$ is the symplectic eigenvalue of $C$ and

$J=\frac{1}{4c^{2}-1}C.$ (71)

For pure states, $\det C$ is a constant, $G_{a}$ then reduces to

$G_{a}=\frac{1}{4c^{2}+1}\Omega(\partial_{a}C)\Omega.$ (72)

From this $G_{a}$, $\vec{L}^{(1)}_{a}$ and $L^{(0)}_{a}$ can be further
obtained, which can be used to obtain the SLD operator via equation (62) and
the QFIM via equation (68).

Another widely used method to obtain the QFI for Gaussian states is through
the fidelity (see section 2.4.2 for the relation between fidelity and QFIM).
The QFI for pure Gaussian states is studied in reference [68]. The QFI for
single-mode Gaussian states has been obtained through the fidelity by Pinel et
al. in 2013 [63], and for two-mode Gaussian states by Šafránek et al. in 2015
[69] and Marian et al. [70] in 2016, based on the expressions of the fidelity
given by Scutaru [71] and Marian et al. in 2012 [72]. The expressions of the
QFI and the fidelity for multi-mode Gaussian states are given by Monras [62],
Safranek et al. [69] and Banchi et al. [73], and reproduced by Oh et al. [74]
with a Hermitian operator related to the optimal measurement of the fidelity.

There are other approaches, such as the exponential state [76], Husimi Q
function [77], that can obtain the QFI and the QFIM for some specific types of
Gaussian states. Besides, a general method to find the optimal probe states to
optimize the QFIM of Gaussian unitary channels is also provided by Šafránek
and Fuentes in 2016 [78], and Matsubara et al. [75] in 2019. Matsubara et al.
performed the optimization of the QFI for Gaussian states in a passive linear
optical circuit. For a fixed total photon number, the optimal Gaussian state
is proved to be a single-mode squeezed vacuum state and the optimal
measurement is a homodyne measurement.

###  2.4 QFIM and geometry of quantum mechanics

####  2.4.1 Fubini-Study metric

In quantum mechanics, the pure states is a normalized vector because of the
basic axiom that the norm square of its amplitude represents the probability.
The pure states thus can be represented as rays in the projective Hilbert
space, on which Fubini-Study metric is a Kähler metric. The squared
infinitesimal distance here is usually expressed as [79]

$\mathrm{d}s^{2}=\frac{\langle\mathrm{d}\psi|\mathrm{d}\psi\rangle}{\langle\psi|\psi\rangle}-\frac{\langle\mathrm{d}\psi|\psi\rangle\langle\psi|\mathrm{d}\psi\rangle}{\langle\psi|\psi\rangle^{2}}.$
(73)

As $\langle\psi|\psi\rangle=1$ and
$|\mathrm{d}\psi\rangle=\sum_{\mu}|\partial_{x_{\mu}}\psi\rangle\mathrm{d}x_{\mu}$,
$\mathrm{d}s^{2}$ can be expressed as

$\mathrm{d}s^{2}=\sum_{\mu\nu}\frac{1}{4}\mathcal{F}_{\mu\nu}\mathrm{d}x_{\mu}\mathrm{d}x_{\nu},$
(74)

here $\mathcal{F}_{\mu\nu}$ is the $\mu\nu$ element of the QFIM. This means
the Fubini-Study metric is a quarter of the QFIM for pure states. This is the
intrinsic reason why the QFIM can depict the precision limit. Intuitively, the
precision limit is just a matter of distinguishability. The best precision
means the maximum distinguishability, which is naturally related to the
distance between the states. The counterpart of Fubini-study metric for mixed
states is the Bures metric, a well-known metric in quantum information and
closely related to the quantum fidelity, which will be discussed below.

####  2.4.2 Fidelity and Bures metric

In quantum information, the fidelity $f(\rho_{1},\rho_{2})$ quantifies the
similarity between two quantum states $\rho_{1}$ and $\rho_{2}$, which is
defined as [80]

$f(\rho_{1},\rho_{2}):=\mathrm{Tr}\sqrt{\sqrt{\rho_{1}}\rho_{2}\sqrt{\rho_{1}}}.$
(75)

Here $f\in[0,1]$ and $f=1$ only when $\rho_{1}=\rho_{2}$. Although the
fidelity itself is not a distance measure, it can be used to construct the
Bures distance, denoted as $D_{\mathrm{B}}$, as [80]

$D^{2}_{\mathrm{B}}(\rho_{1},\rho_{2})=2-2f(\rho_{1},\rho_{2}).$ (76)

The relationship between the fidelity and the QFIM has been well studied in
the literature [30, 81, 82, 83, 31, 84]. In the case that the rank of
$\rho(\vec{x})$ is unchanged with the varying of $\vec{x}$, the QFIM is
related to the infinitestmal Bures distance in the same way as the QFIM
related to the Fubini-study metric 555The derivation is given in F.

$D^{2}_{\mathrm{B}}(\rho(\vec{x}),\rho({\vec{x}+\mathrm{d}\vec{x}}))=\frac{1}{4}\sum_{\mu\nu}\mathcal{F}_{\mu\nu}dx_{\mu}dx_{\nu}.$
(77)

In recent years it has been found that the fidelity susceptibility, the
leading order (the second order) of the fidelity, can be used as an indicator
of the quantum phase transitions [18]. Because of this deep connection between
the Bures metric and the QFIM, it is not surprising that the QFIM can be used
in a similar way. On the other hand, the enhancement of QFIM at the critical
point indicates that the precision limit of the parameter can be improved near
the phase transition, as shown in references [85, 86].

In the case that the rank of $\rho(\vec{x})$ does not equal to that of
$\rho(\vec{x}+\mathrm{d}\vec{x})$, Šafránek recently showed [87] that the QFIM
does not exactly equal to the fidelity susceptibility. Later, Seveso et al.
further suggested [88] that the quantum Cramér-Rao bound may also fail at
those points.

Besides the Fubini-Study metric and the Bures metric, the QFIM is also closely
connected to the Riemannian metric due to the fact that the state space of a
quantum system is actually a Riemannian manifold. In more concrete terms, the
QFIM belongs to a family of contractive Riemannian metric [26, 24, 89, 90],
associated with which the infinitesimal distance in state space is
$\mathrm{d}s^{2}=\sum_{\mu}g_{\mu\nu}\mathrm{d}x_{\mu}\mathrm{d}x_{\nu}$ with
$g_{\mu\nu}$ as the contractive Riemannian metric. In the eigenbasis of the
density matrix $\rho$, $g_{\mu\nu}$ takes the form as [91, 92, 93]

$g_{\mu\nu}=\frac{1}{4}\sum_{i}\frac{\langle\lambda_{i}|\mathrm{d}\rho|\lambda_{i}\rangle^{2}}{\lambda_{i}}+\frac{1}{2}\sum_{i<j}\frac{|\langle\lambda_{i}|\mathrm{d}\rho|\lambda_{j}\rangle|^{2}}{\lambda_{j}h(\lambda_{i}/\lambda_{j})},$
(78)

where $h(\cdot)$ is the Morozova-Čencov function, which is an operator
monotone (for any positive semi-definite operators), self inverse
($xh(1/x)=1/h(x)$) and normalized ($h(1)=1$) real function. When
$h(x)=(1+x)/2$, the metric above reduces to the QFIM (based on the SLD). The
QFIMs based on right and left logarithmic derivatives can also be obtained by
taking $h(x)=x$ and $h(x)=1$. The Wigner-Yanase information metric can be
obtained from it by taking $h(x)=\frac{1}{4}(\sqrt{x}+1)^{2}$.

####  2.4.3 Quantum geometric tensor

The quantum geometric tensor originates from a complex metric in the
projective Hilbert space, and is a powerful tool in quantum information
science that unifies the QFIM and the Berry connection. For a pure state
$|\psi\rangle=|\psi(\vec{x})\rangle$, the quantum geometric tensor $Q$ is
defined as [95, 23]

$Q_{\mu\nu}=\langle\partial_{\mu}\psi|\partial_{\nu}\psi\rangle-\langle\partial_{\mu}\psi|\psi\rangle\langle\psi|\partial_{\nu}\psi\rangle.$
(79)

Recall the expression of QFIM for pure states, given in equation (25), the
real part of $Q_{\mu\nu}$ is actually the QFIM up to a constant factor, i.e.,

$\mathrm{Re}(Q_{\mu\nu})=\frac{1}{4}\mathcal{F}_{\mu\nu}.$ (80)

In the mean time, due to the fact that

$(\langle\partial_{\mu}\psi|\psi\rangle\langle\psi|\partial_{\nu}\psi\rangle)^{*}=\langle\psi|\partial_{\mu}\psi\rangle\langle\partial_{\nu}\psi|\psi\rangle=\langle\partial_{\mu}\psi|\psi\rangle\langle\psi|\partial_{\nu}\psi\rangle,$
(81)

i.e.,
$\langle\partial_{\mu}\psi|\psi\rangle\langle\psi|\partial_{\nu}\psi\rangle$
is real, the imaginary part of $Q_{\mu\nu}$ then reads

$\mathrm{Im}(Q_{\mu\nu})=\mathrm{Im}(\langle\partial_{\mu}\psi|\partial_{\nu}\psi\rangle)=-\frac{1}{2}\left(\partial_{\mu}\mathcal{A}_{\nu}-\partial_{\nu}\mathcal{A}_{\mu}\right),$
(82)

where $\mathcal{A}_{\mu}:=i\langle\psi|\partial_{\mu}\psi\rangle$ is the Berry
connection [96] and
$\Upsilon_{\mu\nu}:=\partial_{\mu}\mathcal{A}_{\nu}-\partial_{\nu}\mathcal{A}_{\mu}$
is the Berry curvature. The geometric phase can then be obtained as [97]

$\phi=\ointop\mathcal{A}_{\mu}\mathrm{d}x_{\mu},$ (83)

where the integral is taken over a closed trajectory in the parameter space.

Recently, Guo et al. [98] connected the QFIM and the Berry curvature via the
Robertson uncertainty relation. Specifically, for a unitary process with two
parameters,
$\Upsilon_{\mu\nu}=i\langle\psi_{0}|[\mathcal{H}_{\mu},\mathcal{H}_{\nu}]|\psi_{0}\rangle$
with $\mathcal{H}_{\mu}$ defined in equation (42) and $|\psi_{0}\rangle$ the
probe state, the determinants of the QFIM and the Berry curvature should
satisfy

$\det\mathcal{F}+4\det\Upsilon\geq 0.$ (84)

###  2.5 QFIM and thermodynamics

The density matrix of a quantum thermal state is

$\rho=\frac{1}{Z}e^{-\beta H},$ (85)

where $Z=\mathrm{Tr}(e^{-\beta H})$ is the partition function and
$\beta=1/(k_{\mathrm{B}}T)$. $k_{\mathrm{B}}$ is the Boltzmann constant and
$T$ is the temperature. For such state we have
$\partial_{T}\rho=\frac{1}{T^{2}}\left(\langle H\rangle-H\right)\rho$, where
we have set $k_{\mathrm{B}}=1$. If we take the temperature as the unknown
parameter, the SLD, which is the solution to
$\partial_{T}\rho=\frac{1}{2}(\rho L_{T}+L_{T}\rho)$, can then be obtained as

$L_{T}=\frac{1}{T^{2}}\left(\langle H\rangle-H\right),$ (86)

which commutes with $\rho$. The QFI for the temperature hence reads

$\mathcal{F}_{TT}=\frac{1}{T^{4}}(\langle H^{2}\rangle-\langle H\rangle^{2}),$
(87)

i.e., $\mathcal{F}_{TT}$ is proportional to the fluctuation of the
Hamiltonian. Compared to the specific heat $C_{v}=\frac{\partial_{T}\langle
H\rangle}{\partial T}=\frac{1}{T^{2}}(\langle H^{2}\rangle-\langle
H\rangle^{2})$, we have [18, 99, 100, 82, 101]

$\mathcal{F}_{TT}=\frac{1}{T^{2}}C_{v},$ (88)

i.e., for a quantum thermal state, the QFI for the temperature is proportional
to the specific heat of this system. For a system of which the Hamiltonian has
no interaction terms, the relation above still holds for its subsystems [102].

The correlation function is an important concept in quantum physics and
condensed matter physics due to the wide applications of the linear response
theory. It is well known that the static susceptibility between two
observables $A$ and $B$, which represents the influence of $\langle
A\rangle$’s perturbation on $\langle B\rangle$ under the thermal equilibrium,
is proportional to the canonical correlation
$\frac{1}{\beta}\int^{\beta}_{0}\frac{1}{Z}\mathrm{Tr}(e^{-\beta
H}e^{sH}Ae^{-sH}B)\mathrm{d}s$ [103], which can be further written into
$\int^{1}_{0}\mathrm{Tr}(\rho^{s}A\rho^{1-s}B)\mathrm{d}s$. Denote
$\int^{1}_{0}\rho^{s}\tilde{L}_{a}\rho^{1-s}\mathrm{d}s=\partial_{a}\rho$ and
replace $A$, $B$ with $\tilde{L}_{a}$ and $\tilde{L}_{b}$, the canonical
correlation reduces to the so-called Bogoliubov–Kubo–Mori Fisher information
matrix $\int^{1}_{0}\rho^{s}\tilde{L}_{a}\rho^{1-s}\tilde{L}_{b}\mathrm{d}s$
[104, 105, 27, 106, 107]. However, this relation does not suggest how to
connect the linear response function with the QFIM based on SLD. In 2007, You
et al. [18, 108] first studied the connection between the fidelity
susceptibility and the correlation function, and then in 2016, Hauke et al.
[109] extended this connection between the QFI and the symmetric and
asymmetric correlation functions to the thermal states. Here we use their
methods to establish the relation between the QFIM and the cross-correlation
functions.

Consider a thermal state corresponding to the Hamiltonian
$H=\sum_{a}x_{a}O_{a}$, where $O_{a}$ is a Hermitian generator for $x_{a}$ and
$[O_{a},O_{b}]=0$ for any $a$ and $b$, the QFIM can be expressed as 666The
details of the derivation can be found in G.

$\mathcal{F}_{ab}=\frac{4}{\pi}\int^{\infty}_{-\infty}\tanh^{2}\left(\frac{\omega}{2T}\right)\mathrm{Re}(S_{ab}(\omega))\mathrm{d}\omega,$
(89)

or equivalently,

$\mathcal{F}_{ab}=\frac{4}{\pi}\int^{\infty}_{-\infty}\tanh\left(\frac{\omega}{2T}\right)\mathrm{Im}(\chi_{ab}(\omega))\mathrm{d}\omega.$
(90)

Here $S_{ab}(\omega)$ is the symmetric cross-correlation spectrum defined as

$S_{ab}(\omega)=\int^{\infty}_{-\infty}\frac{1}{2}\langle\\{Q_{a}(t),O_{b}\\}\rangle
e^{i\omega t}\mathrm{d}t,$ (91)

where $\langle\cdot\rangle=\mathrm{Tr}(\rho\cdot)$ and
$O_{a(b)}(t)=e^{iHt}O_{a(b)}e^{-iHt}$. Its real part can also be written as

$\mathrm{Re}(S_{ab}(\omega))=\int^{\infty}_{-\infty}\frac{1}{2}\langle
Q_{a}(t)O_{b}+O_{b}(t)O_{a}\rangle e^{i\omega t}\mathrm{d}t.$ (92)

$\chi_{ab}$ is the asymmetric cross-correlation spectrum defined as

$\chi_{ab}(\omega)=\int^{\infty}_{-\infty}\frac{i}{2}\langle[O_{a}(t),O_{b}]\rangle
e^{i\omega t}\mathrm{d}t.$ (93)

Because of equations (89) and (90), and the fact that $S_{ab}(\omega)$ and
$\chi_{ab}(\omega)$ can be directly measured in the experiments [110, 111,
112, 113, 114], $\mathcal{F}_{ab}$ becomes measurable in this case, which
breaks the previous understanding that QFI is not observable since the
fidelity is not observable. Furthermore, due to the fact that the QFI is a
witness for multipartite entanglement [22], and a large QFI can imply Bell
correlations [115], equations (89) and (90) provide an experimentally-friendly
way to witness the quantum correlations in the thermal systems. As a matter of
fact, Shitara and Ueda [94] showed that the relations in equations (89) and
(90) can be further extended to the family of metric described in equation
(78) by utilizing the generalized fluctuation-dissipation theorem.

###  2.6 QFIM in quantum dynamics

Quantum dynamics is not only a fundamental topic in quantum mechanics, but
also widely connected to various topics in quantum information and quantum
technology. Due to some excellent mathematical properties, the QFIM becomes a
good candidate for the characterization of certain behaviors and phenomena in
quantum dynamics. In the following we show the roles of QFIM in quantum speed
limit and the characterization of non-Markovianity.

####  2.6.1 Quantum speed limit

Quantum speed limit aims at obtaining the smallest evolution time for quantum
processes [49, 6, 116, 117, 118, 119, 120, 121, 93]. It is closely related to
the geometry of quantum states since the dynamical trajectory with the minimum
evolution time is actually the geodesic in the state space, which indicates
that the QFIM should be capable to quantify the speed limit. As a matter of
fact, the QFI and the QFIM have been used to bound the quantum speed limit in
recent studies [49, 6, 116, 117]. For a unitary evolution, $U=\exp(-iHt)$, to
steer a state away from the initial position with a Bures angle
$D_{\mathrm{B}}=\arccos(f)$ ($f$ is the fidelity defined in equation (75)),
the evolution time $t$ needs to satisfy [6, 49, 122]

$t\geq\frac{2D_{\mathrm{B}}}{\sqrt{\mathcal{F}_{tt}}},$ (94)

where $\mathcal{F}_{tt}$ is the QFI for the time $t$. For a more general case,
that the Hamiltonian is time-dependent, Taddei et al. [49] provided an
implicit bound based on the QFI,

$\sqrt{\frac{1}{2}\frac{\mathrm{d}^{2}\mathcal{D}(f)}{\mathrm{d}f^{2}}\left(\frac{\mathrm{d}\mathcal{D}}{\mathrm{d}f}\right)^{-3}}\Bigg{|}_{f\rightarrow
1}\mathcal{D}(f(\rho(0),\rho(t)))\leq\int^{t}_{0}\sqrt{\frac{\mathcal{F}_{\tau\tau}}{4}}\mathrm{d}\tau,$
(95)

where $\mathcal{D}$ is any metric on the space of quantum states via the
fidelity $f$. In 2017, Beau and del Campo [123] discussed the nonlinear
metrology of many-body open systems and established the relation between the
QFI for coupling constants and the quantum speed limit, which indicates that
the quantum speed limit directly determines the amplitude of the estimation
error in such cases.

Recently, Pires et al. [93] established an infinite family of quantum speed
limits based on a contractive Riemannian metric discussed in section 2.4.2. In
the case that $\vec{x}$ is time-dependent, i.e., $\vec{x}=\vec{x}(t)$, the
geodesic distance $D(\rho_{0},\rho_{t})$ gives a lower bound of general
trajectory,

$D(\rho_{0},\rho_{t})\leq\int^{t}_{0}\left(\frac{\mathrm{d}s}{\mathrm{d}t}\right)\mathrm{d}t=\int^{t}_{0}\mathrm{d}t\sqrt{\sum_{\mu\nu}g_{\mu\nu}\frac{\mathrm{d}x_{\mu}}{dt}\frac{\mathrm{d}x_{\nu}}{dt}},$
(96)

where $g_{\mu\nu}$ is defined in equation (78). Taking the maximum Morozova-
Čencov function, the above inequality leads to the quantum speed limit with
time-dependent parameters given in reference [49].

####  2.6.2 Non-Markovianity

Non-Markovianity is an emerging concept in open quantum systems. Many
different quantification of the non-Markovianity based on monotonic quantities
under the completely positive and trace-preserving maps have been proposed
[124, 125, 126]. The QFI can also be used to characterize the non-Markovianity
since it also satisfies the monotonicity [6]. For the master equation

$\partial_{t}\rho=-i[H_{\vec{x}},\rho]+\sum_{j}\gamma_{j}(t)\left(\Gamma_{j}\rho\Gamma^{\dagger}_{j}-\frac{1}{2}\left\\{\Gamma^{\dagger}_{j}\Gamma_{j},\rho\right\\}\right),$
(97)

the quantum Fisher information flow

$\partial_{t}\mathcal{F}_{aa}=-\sum_{j}\gamma_{j}(t)\mathrm{Tr}\left(\rho[L_{a},\Gamma_{j}]^{\dagger}[L_{a},\Gamma_{j}]\right)$
(98)

given by Lu et al. [127] in 2010 is a valid witness for non-Markovianity.
Later in 2015, Song et al. [128] utilized the maximum eigenvalue of average
QFIM flow to construct a quantitative measure of non-Markovianity. The average
QFIM flow is the time derivative of average QFIM
$\bar{\mathcal{F}}=\int\mathcal{F}\mathrm{d}\vec{x}$. Denote
$\lambda_{\max}(t)$ as the maximum eigenvalue of
$\partial_{t}\bar{\mathcal{F}}$ at time $t$, then the non-Markovianity can be
alternatively defined as

$\mathcal{N}:=\int_{\lambda_{\max}>0}\lambda_{\max}\mathrm{d}t.$ (99)

One may notice that this is not the only way to define non-Makovianity with
the QFIM, similar constructions would also be qualified measures for non-
Markovianity.

##  3 Quantum multiparameter estimation

###  3.1 Quantum multiparameter Cramér-Rao bound

####  3.1.1 Introduction

The main application of the QFIM is in the quantum multiparameter estimation,
which has shown very different properties and behaviors compared to its
single-parameter counterpart [16]. The quantum multiparameter Cramér-Rao
bound, also known as Helstrom bound, is one of the most widely used asymptotic
bound in quantum metrology [2, 1].

###### Theorem 3.1

For a density matrix $\rho$ in which a vector of unknown parameters
$\vec{x}=(x_{0},x_{1},...,x_{m},...)^{\mathrm{T}}$ is encoded, the covariance
matrix $\mathrm{cov}(\hat{\vec{x}},\\{\Pi_{y}\\})$ of an unbiased estimator
$\hat{\vec{x}}$ under a set of POVM, $\\{\Pi_{y}\\}$, satisfies the following
inequality 777The derivation of this theorem is given in H.

$\mathrm{cov}(\hat{\vec{x}},\\{\Pi_{y}\\})\geq\frac{1}{n}\mathcal{I}^{-1}(\\{\Pi_{y}\\})\geq\frac{1}{n}\mathcal{F}^{-1},$
(100)

where $\mathcal{I}(\\{\Pi_{y}\\})$ is the CFIM, $\mathcal{F}$ is the QFIM and
$n$ is the repetition of the experiment.

The second inequality is called the quantum multiparameter Cramér-Rao bound.
In the derivation, we assume the QFIM can be inverted, which is reasonable
since a singular QFIM usually means not all the unknown parameters are
independent and the parameters cannot be estimated simultaneously. In such
cases one should first identify the set of parameters that are independent,
then calculate the corresponding QFIM for those parameters.

For cases where the number of unknown parameters is large, it may be difficult
or even meaningless to know the error of every parameter, and the total
variance or the average variance is a more appropriate macroscopic quantity to
study. Recall that the $a$th diagonal entry of the covariance matrix is
actually the variance of the parameter $x_{a}$. Thus, the bound for the total
variance can be immediately obtained as following.

###### Corollary 3.1.1

Denote $\mathrm{var}(\hat{x}_{a},\\{\Pi_{y}\\})$ as the variance of $x_{a}$,
then the total variance $\sum_{a}\mathrm{var}(x_{a},\\{\Pi_{y}\\})$ is bounded
by the trace of $\mathcal{F}^{-1}$, i.e.,

$\sum_{a}\mathrm{var}(\hat{x}_{a},\\{\Pi_{y}\\})\geq\frac{1}{n}\mathrm{Tr}\left(\mathcal{I}^{-1}(\\{\Pi_{y}\\})\right)\geq\frac{1}{n}\mathrm{Tr}\left(\mathcal{F}^{-1}\right).$
(101)

The inverse of QFIM sometimes is difficult to obtain analytically and one may
need a lower bound of $\mathrm{Tr}(\mathcal{F}^{-1})$ to roughly evaluate the
precision limit. Being aware of the property of QFIM (given in section 2.1)
that $[\mathcal{F}^{-1}]_{aa}\geq 1/\mathcal{F}_{aa}$, one can easily obtain
the following corollary.

###### Corollary 3.1.2

The total variance is bounded as

$\sum_{a}\mathrm{var}(\hat{x}_{a},\\{\Pi_{y}\\})\geq\frac{1}{n}\mathrm{Tr}\left(\mathcal{F}^{-1}\right)\geq\sum_{a}\frac{1}{n\mathcal{F}_{aa}}.$
(102)

The second inequality can only be attained when $\mathcal{F}$ is diagonal.
Similarly,

$\sum_{a}\mathrm{var}(\hat{x}_{a},\\{\Pi_{y}\\})\geq\frac{1}{n}\mathrm{Tr}\left(\mathcal{I}^{-1}(\\{\Pi_{y}\\})\right)\geq\sum_{a}\frac{1}{n\mathcal{I}_{aa}(\\{\Pi_{y}\\})}.$
(103)

The simplest example for the multi-parameter estimation is the case with two
parameters. In this case, $\mathcal{F}^{-1}$ can be calculated analytically as

$\mathcal{F}^{-1}=\frac{1}{\det(\mathcal{F})}\left(\begin{array}[]{cc}\mathcal{F}_{bb}&-\mathcal{F}_{ab}\\\
-\mathcal{F}_{ab}&\mathcal{F}_{aa}\\\ \end{array}\right).$ (104)

Here $\det(\cdot)$ denotes the determinant. With this equation, the corollary
above can reduce to the following form.

###### Corollary 3.1.3

For two-parameter quantum estimation, corollary 3.1.1 reduces to

$\sum_{a}\mathrm{var}(\hat{x}_{a},\\{\Pi_{y}\\})\geq\frac{1}{nI_{\mathrm{eff}}(\\{\Pi_{y}\\})}\geq\frac{1}{nF_{\mathrm{eff}}},$
(105)

where
$I_{\mathrm{eff}}(\\{\Pi_{y}\\})=\det(\mathcal{I})/\mathrm{Tr}(\mathcal{I})$
and $F_{\mathrm{eff}}=\det(\mathcal{F})/\mathrm{Tr}(\mathcal{F})$ can be
treated as effective classical and quantum Fisher information.

####  3.1.2 Attainability

Attainability is a crucial problem in parameter estimation. An unattainable
bound usually means the given precision limit is too optimistic to be realized
in physics. In classical statistical estimations, the classical Cramér-Rao
bound can be attained by the maximum likelihood estimator in the asymptotic
limit, i.e.,
$\lim_{n\rightarrow\infty}n\mathrm{cov}(\hat{\vec{x}}_{\mathrm{m}}(n))=\mathcal{I}^{-1}$,
where $\hat{\vec{x}}_{\mathrm{m}}(n)$ is the local maximum likelihood
estimator and a function of repetition or sample number, which is unbiased in
the asymptotic limit. Because of this, the parameter estimation based on
Cramér-Rao bound is an asymptotic theory and requires infinite samples or
repetition of the experiment. For a finite sample case, the maximum likelihood
estimator is not unbiased and the Cramér-Rao bound may not be attainable as
well. Therefore, the true attainability in quantum parameter estimation should
also be considered in the sense of asymptotic limit since the attainability of
quantum Cramér-Rao bound usually requires the QFIM equals the CFIM first. The
general study of quantum parameter estimation from the asymptotic aspect is
not easy, and the recent progress can be found in reference [90] and
references therein. Here in the following, the attainability majorly refers to
that if the QFIM coincides with the CFIM in theory. Besides, another thing
that needs to emphasize is that the maximum likelihood estimator is optimal in
a local sense [129], i.e., the estimated value is very close to the true
value, and a locally unbiased estimator only attains the bound locally but not
globally in the parameter space, thus, the attainability and optimal
measurement discussed below are referred to the local ones.

For the single-parameter quantum estimations, the quantum Cramér-Rao bound can
be attained with a theoretical optimal measurement. However, for multi-
parameter quantum estimation, different parameters may have different optimal
measurements, and these optimal measurements may not commute with each other.
Thus there may not be a common measurement that is optimal for the estimation
for all the unknown parameters. The quantum Cramér-Rao bound for the
estimation of multiple parameters is then not necessary attainable, which is a
major obstacle for the utilization of this bound in many years. In 2002,
Matsumoto [130] first provided the necessary and sufficient condition for pure
states. After this, its generalization to mixed states was discussed in
several specific scenarios [131, 132, 133, 28] and rigorously proved via the
Holevo bound firstly with the theory of local asymptotic normality [134] and
then the direct minimization of one term in Holevo bound [135]. We first show
this condition in the following theorem.

###### Theorem 3.2

The necessary and sufficient condition for the saturation of the quantum
multiparameter Cramér-Rao bound is

$\mathrm{Tr}\left(\rho[L_{a},L_{b}]\right)=0,~{}~{}\forall~{}a,b.$ (106)

For a pure parameterized state $|\psi\rangle:=|\psi(\vec{x})\rangle$, this
condition reduces to

$\langle\psi|[L_{a},L_{b}]|\psi\rangle=0,~{}~{}\forall~{}a,b,$ (107)

which is equivalent to the form

$\mathrm{Im}(\langle\partial_{a}\psi|\partial_{b}\psi\rangle)=0.$ (108)

When this condition is satisfied, the Holevo bound is also attained and
equivalent to the Cramér-Rao bound [134, 135]. Recall that the Berry curvature
introduced in section 2.4.3 is of the form

$\displaystyle\Upsilon_{ab}$ $\displaystyle=$ $\displaystyle
i\partial_{a}(\langle\psi|\partial_{b}\psi\rangle)-i\partial_{b}(\langle\psi|\partial_{a}\psi\rangle)$
(109) $\displaystyle=$
$\displaystyle-2\mathrm{Im}(\langle\partial_{a}\psi|\partial_{b}\psi\rangle).$

Hence, the above condition can also be expressed as following.

###### Corollary 3.2.1

The multi-parameter quantum Cramér-Rao bound for a pure parameterized state
can be saturated if and only if

$\Upsilon=0,$ (110)

i.e., the matrix of Berry curvature is a null matrix.

For a unitary process $U$ with a pure probe state $|\psi_{0}\rangle$, this
condition can be expressed with the operator $\mathcal{H}_{a}$ and
$\mathcal{H}_{b}$, as shown in the following corollary [48].

###### Corollary 3.2.2

For a unitary process $U$ with a pure probe state $|\psi_{0}\rangle$, the
necessary and sufficient condition for the attainability of quantum
multiparameter Cramér-Rao bound is

$\langle\psi_{0}|\left[\mathcal{H}_{a},\mathcal{H}_{b}\right]|\psi_{0}\rangle=0,~{}~{}\forall~{}a,b.$
(111)

Here $\mathcal{H}_{a}$ was introduced in section 2.3.4.

####  3.1.3 Optimal measurements

The satisfaction of attainability condition theoretically guarantees the
existence of some CFIM that can reach the QFIM. However, it still requires an
optimal measurement. The search of practical optimal measurements is always a
core mission in quantum metrology, and it is for the best that the optimal
measurement is independent of the parameter to be estimated. The most well
studied measurement strategies nowadays include the individual measurement,
adaptive measurement and collective measurement, as shown in figure 3. The
individual measurement refers to the measurement on a single copy (figure
3(a)) or local systems (black lines in figure 3(d)), and can be easily extend
the sequential scenario (figure 3(c)), which is the most common scheme for
controlled quantum metrology. The collective measurement, or joint
measurement, is the one performed simultaneously on multi-copies or on the
global system (orange lines in figure 3(d)) in parallel schemes. A typical
example for collective measurement is the Bell measurement. The adaptive
measurement (figure 3(b)) usually uses some known tunable operations to adjust
the outcome. A well-studied case is the optical Mach-Zehnder interferometer
with a tunable path in one arm. The Mach-Zehnder interferometer will be
thoroughly introduced in the next section.

Figure 3: Schematics for basic measurement schemes in quantum metrology,
including individual measurment, adaptive measurement and collective
measurement.

For the single parameter case, a possible optimal measurement can be
constructed with the eigenstates of the SLD operator. Denote
$\\{|l_{i}\rangle\langle l_{i}|\\}$ as the set of eigenstates of $L_{a}$, if
we choose the set of POVM as the projections onto these eigenstates, then the
probability for the $i$th measurement result is $\langle
l_{i}|\rho|l_{i}\rangle$. In the case where $|l_{i}\rangle$ is independent of
$x_{a}$, the CFI then reads

$\mathcal{I}_{aa}=\sum_{i}\frac{\langle
l_{i}|\partial_{a}\rho|l_{i}\rangle^{2}}{\langle l_{i}|\rho|l_{i}\rangle}.$
(112)

Due to the equation $2\partial_{a}\rho=\rho L_{a}+L_{a}\rho$, the equation
above reduces to

$\mathcal{I}_{aa}=\sum_{i}l^{2}_{i}\langle
l_{i}|\rho|l_{i}\rangle=\mathrm{Tr}(\rho L^{2}_{a})=\mathcal{F}_{aa},$ (113)

which means the POVM $\\{|l_{i}\rangle\langle l_{i}|\\}$ is the optimal
measurement to attain the QFI. However, if the eigenstates of the SLD are
dependent on $x_{a}$, it is no longer the optimal measurement. In the case
with a high prior information, the CFI with respect to
$\\{|l_{i}(\hat{x}_{a})\rangle\langle l_{i}(\hat{x}_{a})|\\}$ ($\hat{x}_{a}$
is the estimated value of $x_{a}$) may be very close to the QFI. In practice,
this measurement has to be used adaptively. Once we obtain a new estimated
value $\hat{x}_{a}$ via the measurement, we need to update the measurement
with the new estimated value and then perform the next round of measurement.
For a non-full rank parameterized density matrix, the SLD operator is not
unique, as discussed in section 2.3.1, which means the optimal measurement
constructed via the eigenbasis of SLD operator is not unique. Thus, finding a
realizable and simple optimal measurement is always the core mission in
quantum metrology. Update to date, only known states in single-parameter
estimation that own parameter-independent optimal measurement is the so-called
quantum exponential family [27, 90], which is of the form

$\rho=e^{\frac{1}{2}\left(\int^{x}_{0}c(x^{\prime})\mathrm{d}x^{\prime}O-\int^{x}_{0}x^{\prime}c(x^{\prime})\mathrm{d}x^{\prime}\right)}\rho_{0}e^{\frac{1}{2}(\int^{x}_{0}c(x^{\prime})\mathrm{d}x^{\prime}O-\int^{x}_{0}x^{\prime}c(x^{\prime})\mathrm{d}x^{\prime})},$
(114)

where $c(x)$ is a function of the unknown parameter $x$, $\rho_{0}$ is a
parameter-independent density matrix, and $O$ is an unbiased observable of
$x$, i.e., $\langle O\rangle=x$. For this family of states, the SLD is
$L_{x}=c(x)(O-x)$ and the optimal measurement is the eigenstates of $O$.

For multiparameter estimation, the SLD operators for different parameters may
not share the same eigenbasis, which means
$\\{|l_{i}(\hat{x}_{a})\rangle\langle l_{i}(\hat{x}_{a})|\\}$ is no longer an
optimal choice for the estimation of all unknown parameters, even with the
adaptive strategy. Currently, most of the studies in multiparameter estimation
focus on the construction of the optimal measurements for a pure parameterized
state $|\psi\rangle$. In 2013, Humphreys et al. [151] proposed a method to
construct the optimal measurement, a complete set of projectors containing the
operator
$|\psi_{\vec{x}_{\mathrm{true}}}\rangle\langle\psi_{\vec{x}_{\mathrm{true}}}|$
($\vec{x}_{\mathrm{true}}$ is the true value of $\vec{x}$). Here
$|\psi_{\vec{x}_{\mathrm{true}}}\rangle$ equals the value of $|\psi\rangle$ by
taking $\vec{x}=\vec{x}_{\mathrm{true}}$ ($\vec{x}_{\mathrm{true}}$). All the
other projectors can be constructed via the Gram-Schmidt process. In practice,
since the true value is unknown, the measurement has to be performed
adaptively with the estimated values $\hat{\vec{x}}$, similar as the single
parameter case. Recently, Pezzè et al. [136] provided the specific conditions
this set of projectors should satisfy to be optimal, which is organized in the
following three theorems.

###### Theorem 3.3

Consider a parameterized pure state $|\psi\rangle$.
$|\psi_{\vec{x}_{\mathrm{true}}}\rangle:=|\psi(\vec{x}=\vec{x}_{\mathrm{true}})\rangle$
with $\vec{x}_{\mathrm{true}}$ the true value of $\vec{x}$. The set of
projectors $\\{|m_{k}\rangle\langle
m_{k}|,|m_{0}\rangle=|\psi_{\vec{x}_{\mathrm{true}}}\rangle\\}$ is an optimal
measurement to let the CFIM reach QFIM if and only if [136]

$\lim_{\vec{x}\rightarrow\vec{x}_{\mathrm{true}}}\frac{\mathrm{Im}(\langle\partial_{a}\psi|m_{k}\rangle\langle
m_{k}|\psi\rangle)}{|\langle m_{k}|\psi\rangle|}=0,~{}~{}\forall
x_{a}~{}\mathrm{and}~{}k\neq 0,$ (115)

which is equivalent to

$\mathrm{Im}(\langle\partial_{a}\psi|m_{k}\rangle\langle
m_{k}|\partial_{b}\psi\rangle)=0,~{}\forall x_{a},x_{b}~{}\mathrm{and}~{}k\neq
0.$ (116)

The proof is given in I. This theorem shows that if the quantum Cramér-Rao
bound can be saturated then it is always possible to construct the optimal
measurement with the projection onto the probe state itself at the true value
and a suitable choice of vectors on the orthogonal subspace [136].

###### Theorem 3.4

For a parameterized state $|\psi\rangle$, the set of projectors
$\\{|m_{k}\rangle\langle m_{k}|,\langle\psi|m_{k}\rangle\neq 0~{}\forall k\\}$
is an optimal measurement to let the CFIM reach QFIM if and only if [136]

$\mathrm{Im}(\langle\partial_{a}\psi|m_{k}\rangle\langle
m_{k}|\psi\rangle)=|\langle\psi|m_{k}\rangle|^{2}\mathrm{Im}(\langle\partial_{a}\psi|\psi\rangle),~{}~{}\forall
k,x_{a}.$ (117)

For the most general case that some projectors are vertical to $|\psi\rangle$
and some not, we have following theorem.

###### Theorem 3.5

For a parameterized pure state $|\psi\rangle$, assume a set of projectors
$\\{|m_{k}\rangle\langle m_{k}|\\}$ include two subsets
$A=\\{|m_{k}\rangle\langle m_{k}|,\langle\psi|m_{k}\rangle=0~{}\forall k\\}$
and $B=\\{|m_{k}\rangle\langle m_{k}|,\langle\psi|m_{k}\rangle\neq 0~{}\forall
k\\}$, i.e., $\\{|m_{k}\rangle\langle m_{k}|\\}=A\cup B$, then it is an
optimal measurement to let the CFIM to reach the QFIM if and nly if [136]
equation (115) is fulfilled for all the projectors in set $A$ and (117) is
fulfilled for all the projectors in set $B$.

Apart from the QFIM, the CFIM is also bounded by a measurement-dependent
matrix with the $ab$th entry $\sum_{k}\mathrm{Re}[\mathrm{Tr}(\rho
L_{a}\Pi_{k}L_{b})]$ [137], where $\\{\Pi_{k}\\}$ is a set of POVM. Recently,
Yang et al. [137] provided the attainable conditions for the CFIM to attain
this bound by generalizing the approach in reference [30].

###  3.2 Phase estimation in the Mach-Zehnder interferometer

Mach-Zehnder interferometer (MZI) is one of the most important model in
quantum technology. It was first proposed in 1890th as an optical
interferometer, and its quantum description was given in 1986 [138]. With the
development of quantum mechanics, now it is not only a model for optical
interferometer, but can also be realized via other systems like spin systems
and cold atoms. The recently developed gravitational wave detector GEO 600 can
also be mapped as a MZI in the absence of noise [139]. It is a little bit more
complicated when the noise is involved, for which a valid bound has been
provided by Branford et al. [140] and is attainable by a frequency-dependent
homodyne detection. Phase estimation in MZI is the earliest case showing
quantum advantages in metrology. In 1981, Caves [141] showed that there exists
an unused port in the MZI due to quantum mechanics and the vacuum fluctuation
in that port actually affects the phase precision and limit it to the standard
quantum limit (also known as shot-noise limit), which is the ultimate limit
for a classical apparatus. He continued to point out that injecting a squeezed
state in the unused port can lead to a high phase precision beating the
standard quantum limit. This pioneer work proved that quantum technologies can
be powerful in the field of precision measure, which was experimentally
confirmed in 1987 [142, 143]. Since then, quantum metrology has been seeing a
rapid development and grown into a major topic in quantum technology.

Figure 4: (a) Schematic for a standard Mach-Zehnder interferometer, which
generally consists of two beam splitters and a phase shift; (b) Schematic for
a double phase Mach-Zehnder interferometer in which the phase shifts are put
in all two paths. Here we emphasize that though the schematic is given in an
optical setup, the MZI model can also be realized via other systems like spin
systems and cold atoms.

MZI is a two-path interferometer, which generally consists of two beam
splitters and a phase shift between them, as shown in Fig. 4(a). In theory,
the beam splitters and phase shift are unitary evolutions. A 50:50 beam
splitter can be theoretically expressed by $B=\exp(-i\frac{\pi}{2}J_{x})$
where$J_{x}=\frac{1}{2}(\hat{a}^{\dagger}\hat{b}+\hat{b}^{\dagger}\hat{a})$ is
a Schwinger operator with $\hat{a}$, $\hat{b}$ ($\hat{a}^{\dagger}$,
$\hat{b}^{\dagger}$) the annihilation (creation) operators for two paths.
Other Schwinger operators are
$J_{y}=\frac{1}{2i}(\hat{a}^{\dagger}\hat{b}-\hat{b}^{\dagger}\hat{a})$ and
$J_{z}=\frac{1}{2}(\hat{a}^{\dagger}\hat{a}-\hat{b}^{\dagger}\hat{b})$. The
Schwinger operators satisfy the commutation
$[J_{i},J_{j}]=i\sum_{k}\epsilon_{ijk}J_{k}$ with $\epsilon_{ijk}$ the Levi-
Civita symbol. The second beam splitter usually takes the form as
$B^{\dagger}$. For the standard MZI, the phase shift can be modeled as
$P=\exp(i\theta J_{z})$, which means the entire operation of MZI is
$BPB^{\dagger}=\exp(-i\theta J_{y})$, which is a SU(2) rotation, thus, this
type of MZI is also referred as SU(2) interferometer. If the input state is a
pure state $|\psi_{0}\rangle$, the QFI for $\theta$ is just the variance of
$J_{y}$ with respect to $|\psi_{0}\rangle$, i.e.,
$\mathrm{var}_{|\psi_{0}\rangle}(J_{y})$.

####  3.2.1 Double-phase estimation

A double-phase MZI consists of two beam splitters and a phase shift in each
path, as shown in figure 4(b). In this setup, the operator of the two phase
shifts is
$P(\phi_{a},\phi_{b})=\exp(i(\phi_{a}\hat{a}^{\dagger}\hat{a}+\phi_{b}\hat{b}^{\dagger}\hat{b}))$.
According to Proposition 2.1, the QFIM for the phases is not affected by the
second beam splitter $B^{\dagger}$ since it is independent of the phases.
Thus, the second beam splitter can be neglected for the calculation of QFIM.
Now take $\phi_{\mathrm{tot}}=\phi_{a}+\phi_{b}$ and
$\phi_{\mathrm{d}}=\phi_{a}-\phi_{b}$ as the parameters to be estimated. For a
separable input state $|\alpha\rangle\otimes|\chi\rangle$ with
$|\alpha\rangle$ as a coherent state, the entries of QFIM reads [144]

$\displaystyle\mathcal{F}_{\phi_{\mathrm{tot}},\phi_{\mathrm{tot}}}$
$\displaystyle=$
$\displaystyle|\alpha|^{2}+\mathrm{var}_{|\chi\rangle}(b^{\dagger}b),$ (118)
$\displaystyle\mathcal{F}_{\phi_{\mathrm{d}},\phi_{\mathrm{d}}}$
$\displaystyle=$ $\displaystyle
2|\alpha|^{2}\mathrm{cov}_{|\chi\rangle}(b,b^{\dagger})-2\mathrm{Re}\left(\alpha^{2}\mathrm{var}_{|\chi\rangle}(b^{\dagger})\right)+\langle
b^{\dagger}b\rangle,$ (119)
$\displaystyle\mathcal{F}_{\phi_{\mathrm{tot}},\phi_{\mathrm{d}}}$
$\displaystyle=$ $\displaystyle-i\alpha^{*}\langle
b\rangle-i\mathrm{Im}(\alpha^{*}(\langle b^{\dagger}b^{2}\rangle-\langle
b^{\dagger}b\rangle\langle b\rangle)).$ (120)

Focusing on the maximization of
$\mathcal{F}_{\phi_{\mathrm{d}},\phi_{\mathrm{d}}}$ subject to a constraint of
fixed mean photon number on $b$ mode, Lang and Caves [144] proved that the
squeezed vacuum state is the optimal choice for $|\chi\rangle$ which leads to
the maximum sensitivity of $\phi_{\mathrm{d}}$.

Since Caves already pointed out that the vacuum fluctuation will affect the
phase sensitivity [141], it is then interesting to ask the question that how
bad the phase sensitivity could be if one input port keeps vacuum? Recently
Takeoka et al. [145] considered this question and gave the answer by proving a
no-go theorem stating that in the double phase MZI, if one input port is the
vacuum state, the sensitivity can never be better than the standard quantum
limit regardless of the choice of quantum state in the other port and the
detection scheme. However, this theorem does not hold for a single phase shift
scenario in figure 4(a). Experimentally, Polino et al. [146] recently
demonstrated quantum-enhanced double-phase estimation with a photonic chip.

Besides the two-phase estimation, another practical two-parameter scenario in
interometry is the joint measurement of phase and phase diffusion. In 2014,
Vidrighin et al. [131] provide a trade-off bound on the statistical variances
for the joint estimation of phase and phase diffusion. Later in 2015, Altorio
et al. [147] addressed the usefulness of weak measurements in this case and in
2018 Hu et al. [148] discussed the SU(1,1) interferometry in the presence of
phase diffusion. In the same year, Roccia et al. [149] experimentally showed
that some collective measurement, like Bell measurement, can benefit the joint
estimation of phase and phase diffusion.

####  3.2.2 Multi-phase estimation

Figure 5: Schematics of (a) simultaneous multi-phase estimation and (b)
independent estimation of multiple phases. $0$-mode light is the reference
mode. In the independent estimation, the phase in the $i$th mode is estimated
via the MZI consisting of $0$-mode and $i$-mode lights.

Apart from double phase estimation, multi-phase estimation is also an
important scenario in multiparameter estimation. The multi-phase estimation is
usually considered in the multi-phase interferometer shown in figure 5(a).
Another recent review on this topic is reference [16]. In this case, the total
variance of all phases is the major concern, which is bounded by
$\frac{1}{n}\mathrm{Tr}(\mathcal{F}^{-1})$ according to Corollary 3.1.1. For
multiple phases, the probe state undergoes the parameterization, which can be
represented by $U=\exp(i\sum_{j}x_{j}H_{j})$ where $H_{j}$ is the generator of
the $j$th mode. In the optical scenario, this operation can be chosen as
$\exp(i\sum_{j}x_{j}a^{\dagger}_{j}a_{j})$ with $a^{\dagger}_{j}$ ($a_{j}$) as
the creation (annihilation) operator for the $j$th mode. For a separable state
$\sum_{k}p_{k}\bigotimes^{d}_{i=0}\rho^{(i)}_{k}$ where $\rho^{(i)}_{k}$ is a
state of the $i$th mode, and $p_{k}$ is the weight with $p_{k}>0$ and
$\sum_{k}p_{k}=1$, the QFI satisfies $\mathcal{F}_{jj}\leq
d(h_{j,\mathrm{max}}-h_{j,\mathrm{min}})^{2}$ [150] with $h_{j,\mathrm{max}}$
($h_{j,\mathrm{min}}$) as the maximum (minimum) eigenvalue of $H_{j}$. From
Corollary 3.1.2, the total variance is then bounded by [150]

$\sum_{j}\mathrm{var}(\hat{x}_{j},\\{\Pi_{y}\\})\geq\frac{1}{d}\sum^{d}_{j=1}\frac{1}{(h_{j,\mathrm{max}}-h_{j,\mathrm{min}})^{2}}.$
(121)

This bound indicates that the entanglement is crucial in the multi-phase
estimation to beat the standard quantum limit.

N00N state is a well known entangled state in quantum metrology which can
saturate the Heisenberg limit. In 2013, Humphreys et al. [151] discussed a
generalized $(d+1)$-mode N00N state in the form
$c_{0}|N_{0}\rangle+c_{1}\sum_{i}|N_{i}\rangle$ where
$|N_{i}\rangle=|0...0N0...0\rangle$ is the state in which only the $i$th mode
corresponds to a Fock state and all other modes are left vacuum. $c_{0}$ and
$c_{1}$ are real coefficients satisfying $c^{2}_{0}+dc^{2}_{1}=1$. The 0-mode
is the reference mode and the parametrization process is
$\exp(i\sum^{d}_{j=1}x_{j}a^{\dagger}_{j}a_{j})$. The generation of this state
has been proposed in reference [152]. Since this process is unitary, the QFIM
can be calculated via Corollary 2.7.1 with
$\mathcal{H}_{j}=-a^{\dagger}_{j}a_{j}$. It is then easy to see that the entry
of QFIM reads [151] $\mathcal{F}_{ij}=4N^{2}(\delta_{ij}c^{2}_{1}-c^{4}_{1})$,
which further gives

$\min_{c_{1}}\mathrm{Tr}(\mathcal{F}^{-1})=\frac{(1+\sqrt{d})^{2}d}{4N^{2}},$
(122)

where the optimal $c_{1}$ is $1/\sqrt{d+\sqrt{d}}$.

Apart from the scheme in figure 5(a) with the simultaneous estimation, the
multi-phase estimation can also be performed by estimating the phases
independently, using the $i$th mode and the reference mode, as shown in figure
5(b). In this scheme, the phases are estimated one by one with the N00N state,
which provides the total precision limit $d^{3}/N^{2}$ [151]. Thus, the
simultaneous estimation scheme in Fig. 5(a) shows a $\mathcal{O}(d)$ advantage
compared to the independent scheme in figure 5(b). However, the performance of
the simultaneous estimation may be strongly affected by the noise [153, 154]
and the $\mathcal{O}(d)$ advantage may even disappear under the photon loss
noise [155].

In the single phase estimation, it is known that the entangled coherent state
$\mathcal{N}(|0\alpha\rangle+|\alpha 0\rangle)$ ($|\alpha\rangle$ is a
coherent state and $\mathcal{N}$ is the normalization) can provide a better
precision limit than the N00N state [156]. Hence, it is reasonable to think
the generalization of entangled coherent state may also outperform the
generalized N00N state. This result was theoretically confirmed in reference
[157]. For a generalized entangled coherent state written in the form
$c_{0}|\alpha_{0}\rangle+c_{1}\sum^{d}_{j=1}|\alpha_{j}\rangle$ with
$|\alpha_{j}\rangle=|0...0\alpha 0...0\rangle$, the QFIM is
$\mathcal{F}_{ij}=4|c_{1}|^{2}|\alpha|^{2}[\delta_{ij}(1+|\alpha|^{2})-|c_{1}|^{2}|\alpha|^{2}]$
[157]. For most values of $d$ and $|\alpha|$, the minimum
$\mathrm{Tr}(\mathcal{F}^{-1})$ can be written as

$\min_{c_{1}}\mathrm{Tr}(\mathcal{F}^{-1})=\frac{(1+\sqrt{d})^{2}d}{4(1+|\alpha|^{2})^{2}},$
(123)

which is smaller than the counterpart of the generalized N00N state,
indicating the performance of the generalized entangled coherent state is
better than the generalized N00N state. For a large $|\alpha|$, the total
particle number $\sim|\alpha|^{2}$ and the performances of both states
basically coincide. With respect to the independent scheme, the generalized
entangled coherent also shows a $\mathcal{O}(d)$ advantage in the absence of
noise.

In 2017, Zhang et al. [158] considered a general balanced state
$c\sum_{j}|\psi_{j}\rangle$ with $|\psi_{j}\rangle=|0...0\psi 0...0\rangle$.
$c$ is the normalization coefficient. Four specific balanced states: N00N
state, entangled coherent state, entangled squeezed vacuum state and entangled
squeezed coherent state, were calculated and they found that the entangled
squeezed vacuum state shows the best performance, and the balanced type of
this state outperforms the unbalanced one in some cases.

A linear network is another common structure for multi-phase estimation [159,
160]. Different with the structure in figure 5(a), a network requires the
phase in each arm can be detected independently, meaning that each arm needs a
reference beam. The calculation of Cramér-Rao bound in this case showed [160]
that the linear network for miltiparameter metrology behaves classically even
though endowed with well-distributed quantum resources. It can only achieve
the Heisenberg limit when the input photons are concentrated in a small number
of input modes. Moreover, the performance of a mode-separable state
$\mathcal{N}(|N\rangle+v|0\rangle)$ may also shows a high theoretical
precision limit in this case if $v\propto\sqrt{d}$.

Recently, Gessner [161] considered a general case for multi-phase estimation
with multimode interferometers. The general Hamiltonian is
$H=\sum^{N}_{j=1}h^{(j)}_{k}$ with $h^{(j)}_{k}$ a local Hamiltonian for the
$i$th particle in the $k$th mode. For the particle-separable states, the
maximum QFIM is a diagonal matrix with the $i$th diagonal entry $\langle
n_{i}\rangle$ being the $i$th average particle number ($n_{i}$ is the $i$th
particle number operator). This bound could be treated as the shot-noise limit
for the multi-phase estimation. For the mode-separable states, the maximum
QFIM is also diagonal, with the $i$th entry $\langle n^{2}_{i}\rangle$, which
means the maximum QFIM for mode-separable state is larger than the particle-
separate counterpart. Taking into account both the particle and mode
entanglement, the maximum QFIM is in the form
$\mathcal{F}_{ij}=\mathrm{sgn}(\langle n_{i}\rangle)\mathrm{sgn}(\langle
n_{j}\rangle)\sqrt{\langle n^{2}_{i}\rangle\langle n^{2}_{j}\rangle}$, which
gives the Heisenberg limit in this case.

The comparison among the performances of different strategies usually requires
the same amount of resources, like the same sensing time, same particle number
and so on. Generally speaking, the particle number here usually refers to the
average particle number, which is also used to define the standard quantum
limit and Heisenberg limit in quantum optical metrology. Although the
Heisenberg limit is usually treated as a scaling behavior, people still
attempt to redefine it as an ultimate bound given via quantum mechanics by
using the expectation of the square of number operator to replace the square
of average particle number [162], i.e., the variance should be involved in the
Heisenberg limit [162, 163] and be treated as a resource.

###  3.3 Waveform estimation

In many practical problems, like the detection of gravitational waves [140] or
the force detection [164], what needs to be estimated is not a parameter, but
a time-varing signal. The QFIM also plays an important role in the estimation
of such signals, also known as waveform estimation [165]. By discretizing time
into small enough intervals, the estimation of a time-continuous signal $x(t)$
becomes the estimation of multiparameters $x_{j}$’s. The prior information of
the waveform has to be taken into account in the estimation problem, e.g.,
restricting the signal to a finite bandwidth, in order to make the estimation
error well-defined. The estimation-error covariance matrix $\Sigma$ is then
given by

$\Sigma_{jk}=\int[\hat{x}_{j}(y)-x_{j}][\hat{x}_{k}(y)-x_{k}]\,\mathrm{d}y\mathrm{d}\vec{x},$
(124)

where $\mathrm{d}\vec{x}=\prod_{j}\mathrm{d}x_{j}$ and $y$ denotes the
measurement outcome. Tsang proved the most general form of Bayesian quantum
Cramér-Rao bound

$\Sigma\geq\left(\mathcal{F}^{\mathrm{(C)}}+\mathcal{F}^{\mathrm{(Q)}}\right)^{-1},$
(125)

where the classical part

$\mathcal{F}^{\mathrm{(C)}}_{ab}=\int\left[\partial_{a}\ln
p(\vec{x})\right]\left[\partial_{b}\ln
p(\vec{x})\right]p(\vec{x})\,\mathrm{d}\vec{x}$ (126)

depends only on the prior information about the vector parameter $\vec{x}$ and
the quantum part [29, 165]

$\mathcal{F}^{\mathrm{(Q)}}_{ab}=\int\mathrm{Re}\left[\mathrm{Tr}(\tilde{L}_{a}^{\dagger}\tilde{L}_{b}\rho)\right]p(\vec{x})\,\mathrm{d}\vec{x}$
(127)

depends on the parametric family $\rho$ of density operators with
$\tilde{L}_{a}$ being determined via

$\partial_{a}\rho=\frac{1}{2}\left(\tilde{L}_{a}\rho+\rho\tilde{L}_{a}^{\dagger}\right).$
(128)

Note that $\tilde{L}_{a}$ is an extended version of SLD and is not necessary
to be Hermitian. When all $\tilde{L}_{a}$’s are Hermitian,
$\mathcal{F}^{\mathrm{(Q)}}$ is the average of the QFIM over the prior
distribution of $\vec{x}$. $\tilde{L}_{a}$ can also be anti-Hermiaition [29].
In this case, the corresponding bound is equivalent to the SLD one for pure
states but a potentially looser one for mixed states. For a unitary evolution,
the entire operator $U$ can be discretized into $U=U_{m}U_{m-1}\cdots
U_{1}U_{0}$, where $U_{a}=\exp(-iH(x_{a})\delta t)$. Denote
$h_{a}=U^{\dagger}_{0}U^{\dagger}_{1}\cdots
U^{\dagger}_{a}(\partial_{a}H)U_{a}\cdots U_{1}U_{0}$ and $\Delta
h_{a}=h_{a}-\mathrm{Tr}(\rho_{0}h_{a})$ with $\rho_{0}$ the probe state, the
quantum part in equation (127) reads [165]

$\mathcal{F}^{(Q)}_{ab}=2(\delta t)^{2}\int\mathrm{Tr}\left(\left\\{\Delta
h_{a},\Delta h_{b}\right\\}\rho_{0}\right)p(\vec{x})\mathrm{d}\vec{x}.$ (129)

Taking the continuous-time limit, it can be rewritten into [165]

$\mathcal{F}^{(Q)}(t,t^{\prime})=2\int\mathrm{Tr}\left(\left\\{\Delta
h(t),\Delta h(t^{\prime})\right\\}\rho_{0}\right)p(\vec{x})\mathrm{d}\vec{x}.$
(130)

Together with $\mathcal{F}^{(C)}(t,t^{\prime})$, the fundamental quantum limit
to waveform estimation based on QFIM is established [165]. The waveform
estimation can also be solved with other tools, like the Bell-Ziv-Zakai lower
bounds [166].

###  3.4 Control-enhanced multiparameter estimation

The dynamics of many artificial quantum systems, like the Nitrogen-vacancy
center, trapped ion and superconducting circuits, can be precisely altered by
control. Hence, control provides another freedom for the enhancement of the
precision limit in these apparatuses. It is already known that quantum control
can help to improve the QFI to the Heisenberg scaling with the absence of
noise in some scenarios [168, 169, 170]. However, like other resources, this
improvement could be sensitive to the noise, and the performance of optimal
control may strongly depend on the type of noise [171, 172]. In general, the
master equation for a noisy quantum system is described by

$\partial_{t}\rho=\mathcal{E}_{\vec{x}}\rho.$ (131)

Here $\mathcal{E}_{\vec{x}}$ is a $\vec{x}$-dependent superoperator. For the
Hamiltonian estimation under control, the dynamics is

$\mathcal{E}_{\vec{x}}\rho=-i[H_{0}(\vec{x})+H_{\mathrm{c}},\rho]+\mathcal{L}\rho,$
(132)

where $H_{\mathrm{c}}=\sum^{p}_{k=1}V_{k}(t)H_{k}$ is the control Hamiltonian
with $H_{k}$ the $k$th control and $V_{k}(t)$ the corresponding time-dependent
control amplitude. For a general Hamiltonian, the optimal control can only be
tackled via numerical methods. One choice for this problem is the Gradient
ascent pulse engineering algorithm.

Figure 6: Flow chart of GRAPE algorithm for controlled single-parameter and
multiparameter estimation. Reprinted figure with permission from [171, 172].
Copyright (2019) by the American Physical Society.

Gradient ascent pulse engineering (GRAPE) algorithm is a gradient-based
algorithm, which was originally developed to search the optimal control for
the design of nuclear magnetic resonance pulse sequences [173], and now is
extended to the scenario of quantum parameter estimation [171, 172]. As a
gradient-based method, GRAPE requires an objective function and analytical
expression of the corresponding gradient. For quantum single-parameter
estimation, the QFI is a natural choice for the objective function. Of course,
in some specific systems the measurement methods might be very limited, which
indicates the CFI is a proper objective function in such cases. For quantum
multiparameter estimation, especially those with a large parameter number, it
is difficult or even impossible to take into account the error of every
parameter, and thus the total variance
$\sum_{a}\mathrm{var}(\hat{x}_{a},\\{\Pi_{y}\\})$ which represents the average
error of all parameters, is then a good index to show system precision.
According to Corollary 3.1.1, $\mathrm{Tr}(\mathcal{F}^{-1})$ and
$\mathrm{Tr}(\mathcal{I}^{-1})$ (for fixed measurement) could be proper
objective functions for GRAPE if we are only concern with the total variance.
For two-parameter estimation, $\mathrm{Tr}(\mathcal{F}^{-1})$ and
$\mathrm{Tr}(\mathcal{I}^{-1})$ reduce to the effective QFI $F_{\mathrm{eff}}$
and CFI $I_{\mathrm{eff}}$ according to Corollary 3.1.3. However, as the
parameter number increases, the inverse matrix of QFIM and CFIM are difficult
to obtain analytically, and consequently superseded objective functions are
required. $(\sum_{a}\mathcal{I}_{aa})^{-1}$ or
$(\sum_{a}\mathcal{F}_{aa})^{-1}$, based on corollary 3.1.2, is then a
possible superseded objective function for fixed measurement. However, the use
of $(\sum_{a}\mathcal{F}_{aa})^{-1}$ should be very cautious since
$\mathrm{Tr}(\mathcal{F}^{-1})$ cannot always be achievable. The specific
expressions of gradients for these objective functions in Hamiltonian
estimation are given in J.

The flow of the algorithm, shown in figure 6, is formulated as follows [171,
172]:

  1. 1.

Guess a set of initial values for $V_{k}(j)$ [$V_{k}(j)$ is the $k$th control
at the $j$th time step].

  2. 2.

Evolve the dynamics with the controls.

  3. 3.

Calculate the objective function. For single-parameter estimation, the
objective function can be chosen as QFI $\mathcal{F}_{aa}$ or CFI
$\mathcal{I}_{aa}$ (for mixed measurement). For two-parameter estimation, it
can be chosen as $I_{\mathrm{eff}}$ or $F_{\mathrm{eff}}$. For large parameter
number, it can be chosen as $f_{0}=\frac{1}{\sum_{a}\mathcal{I}^{-1}_{aa}}$ or
$\frac{1}{\sum_{a}\mathcal{F}^{-1}_{aa}}$.

  4. 4.

Calculate the gradient.

  5. 5.

Update $V_{k}(j)$ to $V_{k}(j)+\epsilon\cdot\mathrm{gradient}$ ($\epsilon$ is
a small quantity) for all $j$ simultaneously.

  6. 6.

Go back to step (ii) until the objective function converges.

Figure 7: (a) In GRAPE, the entire evolution time $T$ is cut into $m$ parts
with time interval $\Delta t$, i.e., $m\Delta t=T$. $V_{k}(t)$ within the
$j$th time interval is denoted as $V_{k}(j)$ and is assumed to be a constant.
All the $V_{k}(j)$ are simultaneously updated in each iteration. (b) Krotov’s
method updates one $V_{k}(j)$ in each iteration.

In step (v), all $V_{k}(j)$ are updated simultaneously in each time of
iteration [173], as shown in figure 7(a), which could improve the speed of
convergence in some cases. However, this parallel update method does not
promise the monotonicity of convergence, and the choice of initial guess of
$V_{k}(j)$ is then important for a convergent result. The specific codes for
this algorithm can be found in reference [175].

Krotov’s method is another gradient-based method in quantum control, which
promises the monotonicity of convergence during the iteration [174]. Different
from GRAPE, only one $V_{k}(j)$ is updated in each iteration in Krotov’s
method, as shown in figure 7(b). This method can also be extended to quantum
parameter estimation with the aim of searching the optimal control for a high
precision limit. It is known that the gradient-based methods can only harvest
the local extremals. Thus, it is also useful to involve gradient-free methods,
including Monte Carlo method and Particle Swarm Optimization, into the quantum
parameter estimation, or apply a hybrid method combining gradient-based and
gradient-free methods as discussed in reference [176].

###  3.5 Estimation of a magnetic field

Measurement of the magnetic fields is an important application of quantum
metrology, as it promises better performances than the classical counterparts.
Various physical systems have been used as quantum magnetometers, including
but not limited to nitrogen-vacancy centers [177, 178], optomechanical systems
[179], and cold atoms [180].

A magnetic field can be represented as a vector,
$\vec{B}=B(\cos\theta\cos\phi,\cos\theta\sin\phi,\sin\phi)$, in the spherical
coordinates, where $B$ is the amplitude, and $\theta,\phi$ define the
direction of the field. Therefore, the estimation of a magentic field is in
general a three-parameter estimation problem. The estimation of the amplitude
with known angles is the most widely-studied case in quantum metrology, both
in theory and experiments. For many quantum systems, it is related to the
estimation of the strength of system-field coupling. In the case where one
angle, for example $\phi$ is known, the estimation of the field becomes a two-
parameter estimation problem. The simplest quantum detector for this case is a
single-spin system [53, 48]. An example of the interaction Hamiltonian is
$H=-B\vec{n}_{0}\cdot\vec{\sigma}$ with $n_{0}=(\cos\theta,0,\sin\theta)$ and
$\vec{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$. The QFIM can then be
obtained by making use of the expressions of $\mathcal{H}$ operators [53, 48]

$\mathcal{H}_{B}=t\vec{n}_{0}\cdot\vec{\sigma},~{}~{}\mathcal{H}_{\theta}=-\frac{1}{2}\sin(Bt)\vec{n}_{1}\cdot\vec{\sigma},$
(133)

where $n_{1}=\left(\cos(Bt)\sin\theta,\sin(Bt),-\cos(Bt)\cos\theta\right)$.
The entries of QFIM then read

$\displaystyle\mathcal{F}_{\theta\theta}$ $\displaystyle=$
$\displaystyle\sin^{2}(Bt)[1-(\vec{n}_{1}\cdot\vec{r}_{\mathrm{in}})^{2}],$
(134) $\displaystyle\mathcal{F}_{BB}$ $\displaystyle=$ $\displaystyle
4t^{2}[1-(\vec{n}_{0}\cdot\vec{r}_{\mathrm{in}})^{2}],$ (135)
$\displaystyle\mathcal{F}_{B\theta}$ $\displaystyle=$ $\displaystyle
2t\sin(Bt)(\vec{n}_{0}\cdot\vec{r}_{\mathrm{in}})(\vec{n}_{1}\cdot\vec{r}_{\mathrm{in}}),$
(136)

where $\vec{r}_{\mathrm{in}}$ is the Bloch vector of the probe state. The
maximal values of $\mathcal{F}_{\theta\theta}$ and $\mathcal{F}_{BB}$ can be
attained when $\vec{r}_{\mathrm{in}}$ is vertical to both $\vec{n}_{0}$ and
$\vec{n}_{1}$. However, as a two-parameter estimation problem, we have to
check the value of
$\langle\psi_{0}|[\mathcal{H}_{B},\mathcal{H}_{\theta}]|\psi_{0}\rangle$ due
to Corollary 3.2.2. Specifically, for a pure probe state it reads
$\langle\psi_{0}|[\mathcal{H}_{B},\mathcal{H}_{\theta}]|\psi_{0}\rangle=-\sin(Bt)(\vec{n}_{0}\times\vec{n}_{1})\cdot\vec{\sigma}$.
For the time $t=n\pi/B$ ($n=1,2,3...$), the expression above vanishes.
However, another consequence of this condition is that
$\mathcal{F}_{\theta\theta}$ also vanishes. Thus, single-qubit probe may not
be an ideal magnetometer when at least one angle is unknown. One possible
candidate is a collective spin system, in which the $\mathcal{H}$ operator and
QFIM are provided in reference [181]. Another simple and practical-friendly
candidate is a two-qubit system. One qubit is the probe and the other one is
an ancilla, which does not interact with the field. Consider the Hamiltonian
$H=-\vec{B}\cdot\vec{\sigma}$, the maximal QFIM in this case is (in the basis
$\\{B,\theta,\phi\\}$) [84]

$\max\mathcal{F}=4\left(\begin{array}[]{ccc}t^{2}&0&0\\\ 0&\sin^{2}(Bt)&0\\\
0&0&\sin^{2}(Bt)\cos^{2}\theta\end{array}\right),$ (137)

which can be attained by any maximally entangled state and the Bell
measurement as the optimal measurement. With the assistance of an ancilla, all
three parameters $B,\theta,\phi$ of the field can be simultaneously estimated.
However, $\mathcal{F}_{\theta\theta}$ and $\mathcal{F}_{\phi\phi}$ are only
proportional to $\sin^{2}(Bt)$, indicating that unlike the estimation of the
amplitude, a longer evolution time does not always lead to a better precision
for the estimation of $\theta$ or $\phi$.

In 2016, Baumgratz and Datta [182] provided a framework for the estimation of
a multidimensional field with noncommuting unitary generators. Analogous to
the optical multi-phase estimation, simultaneous estimation shows a better
performance than separate and individual estimation.

To improve the performance of the probe system, additional control could be
employed. For unitary evolution, reference [84] showed that the performance
can be significantly enhanced by inserting the anti-evolution operator as the
control. Specifically if after each evolution of a period of $\delta t$, we
insert a control which reverses the evolution as $U_{c}=U^{\dagger}(\delta
t)=e^{iH(\vec{B})\delta t}$, the QFIM can reach

$\mathcal{F}=4N^{2}\left(\begin{array}[]{ccc}(\delta t)^{2}&0&0\\\
0&\sin^{2}(B\delta t)&0\\\ 0&0&\sin^{2}(B\delta
t)\cos^{2}\theta\end{array}\right),$ (138)

where $N$ is the number of the injected control pulses, and $N\delta t=t$. In
practice the control needs to be applied adaptively as
$U_{c}=e^{iH(\hat{\vec{B}})\delta t}$ with $\hat{\vec{B}}$ as the estimated
value obtained from previous measurement results. In the controlled scheme,
the number of control pulses becomes a resource for the estimation of all
three parameters. For the amplitude, the precision limit is the same as non-
controlled scheme. But for the angles, the precision limit is significantly
improved, especially when $N$ is large. When $\delta t\rightarrow 0$, the QFIM
becomes

$\mathcal{F}=4\left(\begin{array}[]{ccc}t^{2}&0&0\\\ 0&B^{2}t^{2}&0\\\
0&0&B^{2}t^{2}\cos^{2}\theta\end{array}\right),$ (139)

which reaches the highest precision for the estimation of $B$, $\theta$ and
$\phi$ simultaneously. Recently, Hou et al. [183] demonstrated this scheme up
to eight controls in an optical platform and demonstrate a precision near the
Heisenberg limit.

Taking into account the effect of noise, the optimal control for the
estimation of the magnetic field is then hard to obtain analytically. The
aforementioned numerical methods like GRAPE or Krotov’s method could be used
in this case to find the optimal control. The performance of the controlled
scheme relies on the specific form of noise [171, 172], which could be
different for different magnetometers.

In practice, the magnetic field may not be a constant field in space. In this
case, the gradient of the field also needs to be estimated. The corresponding
theory for the estimation of gradients based on Cramér-Rao bound has been
established in references [184, 185].

###  3.6 Other applications and alternative mathematical tools

Apart from the aforementioned cases, quantum multiparameter estimation is also
studied in other scenarios, including the spinor systems [187], unitary
photonic systems [188], networked quantum sensors [189], and the quantum
thermometry [190]. In the case of noise estimation, the simultaneous
estimation of loss parameters has been studied in reference [191].

In 2016, Tsang et al. [192] used quantum multiparameter estimation for a
quantum theory of superresolution for two incoherent optical point sources.
With weak-source approximation, the density operator for the optical field in
the imaging plane can be expressed as
$\rho=(1-\epsilon)|\mathrm{vac}\rangle\langle\mathrm{vac}|+\frac{\epsilon}{2}(|\psi_{1}\rangle\langle\psi_{1}|+|\psi_{2}\rangle\langle\psi_{2}|)$,
with $\epsilon$ the average photon number in a temporal mode,
$|{\mathrm{vac}}\rangle$ the vacuum state,
$|\psi_{j}\rangle=\int\mathrm{d}y\,\psi_{x_{j}}(y)|y\rangle$ the quantum state
of an arrival photon from the point source located at $x_{j}$ of the object
plane, $\psi_{x_{j}}(y)$ the wave function in the image plane, and $|y\rangle$
the photon image-plane position eigenstate. Taking the position coordinates to
be one-dimensional, the parameters to be estimated are $x_{1}$ and $x_{2}$.
The performance of the underlying quantum measurement can be assessed by the
CFIM and its fundamental quantum limit can be revealed by the QFIM. In this
case, the QFIM is analytically calculated in the 4-dimensional Hilbert
subspace spanned by $|\psi_{1}\rangle$, $|\psi_{2}\rangle$,
$|\partial_{1}\psi_{1}\rangle$, and $|\partial_{2}\psi_{2}\rangle$. For
convenience, the two parameters $x_{1}$ and $x_{2}$ can be transformed to the
centriod $\theta_{1}=(x_{1}+x_{2})/2$ and the separation
$\theta_{2}=x_{2}-x_{1}$. Assuming the imaging system is spatially invariant
such that the point-spread function of the form
$\psi_{x_{j}}(y)=\psi(y-x_{j})$ with real-valued $\psi(y)$, the QFIM with
respect to $\theta_{1}$ and $\theta_{2}$ is then given by

$\mathcal{F}=\left(\begin{array}[]{cc}4\epsilon(\Delta k^{2}-\gamma^{2})&0\\\
0&\epsilon\Delta k^{2}\end{array}\right),$ (140)

where $\Delta k^{2}=\int_{-\infty}^{\infty}[\partial\psi(y)/\partial
y]^{2}\,\mathrm{d}y$ and
$\gamma=\int_{-\infty}^{\infty}\psi(y-\theta_{2})\partial\psi(y)/\partial
y\,\mathrm{d}y$ [192]. With this result, it has been shown that direct imaging
performs poorly for estimating small separations and other elaborate methods
like spatial-mode demultiplexing and image inversion interferometry can be
used to achieve the quantum limit [192, 193].

The design of enhanced schemes for quantum parameter estimation would
inevitably face optimization processes, including the optimization of probe
states, the parameterization trajectories, the measurement and so on. If
control is involved, one would also need to harvest the optimal control.
Therefore, the stochastic optimization methods, including convex optimization,
Monte Carlo method and machine learning, could be used in quantum parameter
estimation.

Machine learning is one of the most promising and cutting-edge methods
nowadays. In recent years, it has been applied to condensed matter physics for
phase transitions [194, 195], design of a magneto-optical trap [196] and other
aspects [197]. With respect to quantum parameter estimation, in 2010,
Hentschel and Sanders [198] applied the particle swarm optimization algorithm
[199] in the adaptive phase estimation to determine the best estimation
strategy, which was experimentally realized by Lumino et al. [200] in 2018.
Meanwhile, in 2017, Greplova et al. [201] proposed to use neural networks to
estimate the rates of coherent and incoherent processes in quantum systems
with continuous measurement records. In the case of designing controlled
schemes, similar to gradient-based methods, machine learning could also help
use to find the optimal control protocal in order to achieve the best
precision limit [186]. In particular, for quantum multiparameter estimation,
there still lacks of systematic research on the role of machine learning. For
instance, we are not assured if there exists different performance compared to
single-parameter estimation.

As a mathematical tool, quantum Cramér-Rao bound is not the only one for
quantum multiparameter estimation. Bayesian approach [204, 205, 206, 207,
208], including Ziv-Zakai family [167, 202, 166] and Weiss-Weinstein family
[165, 209], and some other tools [209, 210, 211] like Holevo bound [2, 28,
203], have also shown their validity in various multiparameter scenarios.
These tools beyond the Cramér-Rao bound will be thoroughly reviewed by M. Guţă
et al. [212] in the same special issue and hence we do not discuss them in
this paper. Please see reference [212] for further reading of this topic. In
some scenarios in quantum parameter estimations, there may exist some
parameters that not of interest but affect the precision of estimating other
parameters of interest, which is usually referred to as the nuisance
parameters. Suzuki et al. [213] has thoroughly reviewed the quantum state
estimation with nuisance parameters in the same special issue. Please see
reference [213] for further reading.

##  4 Conclusion and outlook

Quantum metrology has been recognized as one of the most promising quantum
technologies and has seen a rapid development in the past few decades. Quantum
Cramér-Rao bound is the most studied mathematical tool for quantum parameter
estimation, and as the core quantity of quantum Cramér-Rao bound, the QFIM has
drawn plenty of attentions. It is now well-known that the QFIM is not only a
quantity to quantify the precision limit, but also closely connected to many
different subjects in quantum physics. This makes it an important and
fundamental concept in quantum mechanics.

In recent years, many quantum multiparameter estimation schemes have been
proposed and discussed with regard to various quantum systems, and some of
them have shown theoretical advances compared to single-parameter schemes.
However, there are still many open problems, such as the design of optimal
measurement, especially simple and practical measurement which are independent
of the unknown parameters, as well as the robustness of precision limit
against the noises and imperfect controls. We believe that some of the issues
will be solved in the near future. Furthermore, quantum multiparameter
metrology will then go deeply into the field of applied science and become a
basic technology for other aspects of sciences.

The authors thank Zibo Miao, Christiane Koch, Rafał Demkowicz-Dobrzański,
Yanming Che, Dominic Branford, Dominik Šafránek, Jasminder Sidhu, Wenchao Ge,
Jesus Rubio, Wojciech Rzadkowski, Marco Genoni, Haggai Landa and two anonymous
referees for for helpful discussions and suggestions. J.L. particularly thanks
Yanyan Shao and Jinfeng Qin for the assistance on drawing some of the
schematics. J.L. acknowledges the support from National Natural Science
Foundation of China (Grant No. 11805073) and the startup grant of HUST. H.Y.
acknowledges the support from the Research Grants Council of Hong Kong (GRF
No. 14207717). X.M.L. acknowledges the support from the National Natural
Science Foundation of China (Grant Nos. 61871162 and 11805048), and Zhejiang
Provincial Natural Science Foundation of China (Grant No. LY18A050003). X.W.
acknowledges the support from National Key Research and Development Program of
China (Grant Nos. 2017YFA0304202 and 2017YFA0205700), the National Natural
Science Foundation of China (Grant Nos. 11935012 and 11875231), and the
Fundamental Research Funds for the Central Universities (Grant No.
2017FZA3005).

##  Appendix A Derivation of traditional form of QFIM

The traditional calculation of QFIM usually assume a full rank density matrix,
of which the spectral decompostion is in the form

$\rho=\sum^{\dim\rho-1}_{i=0}\lambda_{i}|\lambda_{i}\rangle\langle\lambda_{i}|,$
(1)

where $\lambda_{i}$ and $|\lambda_{i}\rangle$ are $i$th eigenvalue and
eigenstate of $\rho$. $\dim\rho$ is the dimension of $\rho$.
$|\lambda_{i}\rangle$ satisfies
$\sum^{\dim\rho-1}_{i=0}|\lambda_{i}\rangle\langle\lambda_{i}|=\leavevmode\hbox{\small
1\kern-4.3pt\normalsize 1}$ with 11 the identity matrix. Substituting equation
(1) into the equation of SLD ($\partial_{x_{a}}$ is the abbreviation of
$\partial/\partial_{x_{a}}$)

$\partial_{x_{a}}\rho=\frac{1}{2}\left(\rho L_{x_{a}}+L_{x_{a}}\rho\right),$
(2)

and taking $\langle\lambda_{i}|\cdot|\lambda_{j}\rangle$ on both sides of
above equation, one can obtain

$\langle\lambda_{i}|L_{x_{a}}|\lambda_{j}\rangle=\frac{2\langle\lambda_{i}|\partial_{x_{a}}\rho|\lambda_{j}\rangle}{\lambda_{i}+\lambda_{j}}.$
(3)

Next, utlizing equation (1), the QFIM

$\mathcal{F}_{ab}=\frac{1}{2}\mathrm{Tr}\left(\rho\\{L_{x_{a}},L_{x_{b}}\\}\right)$
(4)

can be rewritten as

$\mathcal{F}_{ab}=\frac{1}{2}\sum^{\mathrm{dim}\rho-1}_{i=0}\lambda_{i}\left(\langle\lambda_{i}|L_{x_{a}}L_{x_{b}}|\lambda_{i}\rangle+\langle\lambda_{i}|L_{x_{b}}L_{x_{a}}|\lambda_{i}\rangle\right).$
(5)

Inserting the equation $\leavevmode\hbox{\small 1\kern-4.3pt\normalsize
1}=\sum^{\dim\rho-1}_{i=0}|\lambda_{i}\rangle\langle\lambda_{i}|$ into above
equation, one can obtain

$\mathcal{F}_{ab}=\sum^{\dim\rho-1}_{i,j=0}\lambda_{i}\mathrm{Re}(\langle\lambda_{i}|L_{x_{a}}|\lambda_{j}\rangle\langle\lambda_{j}|L_{x_{b}}|\lambda_{i}\rangle).$
(6)

Substituting equation (3) into this equation, one has

$\mathcal{F}_{ab}=\sum^{\dim\rho-1}_{i,j=0}4\lambda_{i}\frac{\mathrm{Re}(\langle\lambda_{i}|\partial_{x_{a}}\rho|\lambda_{j}\rangle\langle\lambda_{j}|\partial_{x_{b}}\rho|\lambda_{i}\rangle)}{(\lambda_{i}+\lambda_{j})^{2}}.$
(7)

Exchange subscripts $i$ and $j$, the traditional form of QFIM is obtained as
below

$\mathcal{F}_{ab}=\sum^{\dim\rho-1}_{i,j=0}\frac{2\mathrm{Re}(\langle\lambda_{i}|\partial_{x_{a}}\rho|\lambda_{j}\rangle\langle\lambda_{j}|\partial_{x_{b}}\rho|\lambda_{i}\rangle)}{\lambda_{i}+\lambda_{j}}.$
(8)

##  Appendix B Derivation of QFIM for arbitrary-rank density matrices

In this appendix we show the detailed derivation of QFIM for arbitrary-rank
density matrices. Here the spectral decompostion of $\rho$ is

$\rho=\sum^{N-1}_{i=0}\lambda_{i}|\lambda_{i}\rangle\langle\lambda_{i}|,$ (1)

where $\lambda_{i}$ and $|\lambda_{i}\rangle$ are $i$th eigenvalue and
eigenstate of $\rho$. Notice $N$ here is the dimension of $\rho$’s support.
For a full-rank density matrix, $N$ equals to $\dim\rho$, the dimension of
$\rho$, and for a non-full rank density matrix, $N<\dim\rho$, and
$\sum^{N-1}_{i=0}|\lambda_{i}\rangle\langle\lambda_{i}|\neq\leavevmode\hbox{\small
1\kern-4.3pt\normalsize 1}$ (11 is the identity matrix and
$\leavevmode\hbox{\small 1\kern-4.3pt\normalsize
1}=\sum^{\dim\rho-1}_{i=0}|\lambda_{i}\rangle\langle\lambda_{i}|$).
Furthermore, $\lambda_{i}\neq 0$ for $i\in[0,N-1]$. With these notations,
$\mathcal{F}_{ab}$ can be expressed by

$\displaystyle\mathcal{F}_{ab}$ $\displaystyle=$
$\displaystyle\sum^{N-1}_{i=0}\frac{(\partial_{x_{a}}\lambda_{i})(\partial_{x_{b}}\lambda_{i})}{\lambda_{i}}+\sum^{N-1}_{i=0}4\lambda_{i}\mathrm{Re}\left(\langle\partial_{x_{a}}\lambda_{i}|\partial_{x_{b}}\lambda_{i}\rangle\right)$
(2)
$\displaystyle-\sum^{N-1}_{i,j=0}\frac{8\lambda_{i}\lambda_{j}}{\lambda_{i}+\lambda_{j}}\mathrm{Re}(\langle\partial_{x_{a}}\lambda_{i}|\lambda_{j}\rangle\langle\lambda_{j}|\partial_{x_{b}}\lambda_{i}\rangle).$

Followings are the detailed proof of this equation.

Based on the equation of SLD ($\partial_{x_{a}}$ is the abbreviation of
$\partial/\partial x_{a}$)