Quantum Fisher information matrix and multiparameter estimation

TL;DR

This paper reviews the quantum Fisher information matrix's properties, calculation methods, and its broader significance in quantum mechanics, along with recent advances in multiparameter quantum estimation and the associated bounds and measurements.

Contribution

It provides a comprehensive overview of QFIM's properties, calculation techniques, and its role beyond quantum metrology, including recent developments in multiparameter estimation.

Findings

QFIM links to quantum thermodynamics, phase transitions, and entanglement.

The paper discusses conditions for attaining the quantum Cramér-Rao bound.

Recent scenarios show quantum advantages in multiparameter estimation.

Abstract

Quantum Fisher information matrix (QFIM) is a core concept in theoretical quantum metrology due to the significant importance of quantum Cram\'{e}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\'{e}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.

Equations

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

$${\partial }_a\rho =\frac{1}{2}\left(\rho L_a+L_a\rho \right).$$

$${\mathcal{F}}_{ab}=\mathrm{Tr}\left(L_b{\partial }_a\rho \right)=-\mathrm{Tr}(\rho {\partial }_a{L}_b).$$

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

$${\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.$$

$${\partial }_t\rho =-i[H_{\vec{x}},\rho ],$$

$$\rho =e^{-i{H}_{\vec{x}}t}{\rho }_0{e}^{i{H}_{\vec{x}}t}.$$

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

$$\mathcal{L}(\rho )=\sum _j{\gamma }_j\left({\Gamma }_j\rho {\Gamma }_j^{†}-\frac{1}{2}\left\{{\Gamma }_j^{†}{\Gamma }_j,\rho \right\}\right),$$

$${\mathcal{F}}_{ab}=\sum _{i,j=0}^{d-1}\frac{2\mathrm{Re}(⟨{\lambda }_i|{\partial }_a\rho |{\lambda }_j⟩⟨{\lambda }_j|{\partial }_b\rho |{\lambda }_i⟩)}{{\lambda }_i+{\lambda }_j},$$

$${\mathcal{F}}_{ab}=\sum _{i,j=0,{\lambda }_i+{\lambda }_j\ne 0}^{d-1}\frac{2\mathrm{Re}(⟨{\lambda }_i|{\partial }_a\rho |{\lambda }_j⟩⟨{\lambda }_j|{\partial }_b\rho |{\lambda }_i⟩)}{{\lambda }_i+{\lambda }_j}.$$

$${\mathcal{F}}_{ab}=\sum _{i=0}^{d-1}\frac{({\partial }_a{\lambda }_i)({\partial }_b{\lambda }_i)}{{\lambda }_i}+\sum _{i\ne j,{\lambda }_i+{\lambda }_j\ne 0}\frac{2({\lambda }_i-{\lambda }_j{)}^2}{{\lambda }_i+{\lambda }_j}\mathrm{Re}\left(⟨{\lambda }_i|{\partial }_a{\lambda }_j⟩⟨{\partial }_b{\lambda }_j|{\lambda }_i⟩\right).$$

$${\mathcal{F}}_{ab}$$

$${\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⟨{\partial }_a{\lambda }_i|{\partial }_a{\lambda }_i⟩-\sum _{{\lambda }_i,{\lambda }_j\in \mathcal{S}}\frac{8{\lambda }_i{\lambda }_j}{{\lambda }_i+{\lambda }_j}|⟨{\partial }_a{\lambda }_i|{\lambda }_j⟩{|}^2.$$

$$⟨{\lambda }_i|L_a|{\lambda }_j⟩={\delta }_{ij}\frac{{\partial }_a{\lambda }_i}{{\lambda }_i}+\frac{2({\lambda }_j-{\lambda }_i)}{{\lambda }_i+{\lambda }_j}⟨{\lambda }_i|{\partial }_a{\lambda }_j⟩;$$

$$L_a=2{\int }_0^{\infty }{e}^{-\rho s}\left({\partial }_a\rho \right)e^{\rho s}\mathrm{d}s,$$

$$L_a=-2\underset{s\to \infty }{\lim }\sum _{n=0}^{\infty }\frac{{\left(-s\right)}^{n+1}}{(n+1)!}{\mathcal{R}}_{\rho }^n({\partial }_a\rho ),$$

$$L_a=-2\underset{s\to \infty }{\lim }\sum _{n=0}^{\infty }\sum _{m=0}^n\frac{{\left(-s\right)}^{n+1}}{(n+1)!}\left(\genfrac{}{}{0ex}{}{n}{m}\right){\rho }^m\left({\partial }_a\rho \right){\rho }^{n-m}.$$

$${\mathcal{F}}_{ab}=2\mathrm{vec}({\partial }_a\rho {)}^{†}{\left(\rho \otimes \text{leavevmode}11+\text{leavevmode}11\otimes {\rho }^{\ast }\right)}^{-1}\mathrm{vec}({\partial }_b\rho ),$$

$$\mathrm{vec}(L_a)=2{\left(\rho \otimes \text{leavevmode}11+\text{leavevmode}11\otimes {\rho }^{\ast }\right)}^{-1}\mathrm{vec}({\partial }_a\rho ).$$

$$\rho =\frac{1}{d}\left(\text{leavevmode}11+\sqrt{\frac{d(d-1)}{2}}\vec{r}\cdot \vec{\kappa }\right),$$

$${\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},$$

$$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.$$

$${\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},$$

$${\mathcal{F}}_{ab}=4\mathrm{Re}(⟨{\partial }_a\psi |{\partial }_b\psi ⟩-⟨{\partial }_a\psi |\psi ⟩⟨\psi |{\partial }_b\psi ⟩).$$

$${\mathcal{F}}_{aa}=4(⟨{\partial }_a\psi |{\partial }_a\psi ⟩-|⟨{\partial }_a\psi |\psi ⟩{|}^2),$$

$${\mathcal{F}}_{ab}=4t^2{\mathrm{cov}}_{|{\psi }_0⟩}(H_a,H_b),$$

$${\mathrm{cov}}_{|\phi ⟩}(A,B):=\frac{1}{2}⟨\phi |\{A,B\}|\phi ⟩-⟨\phi |A|\phi ⟩⟨\phi |B|\phi ⟩.$$

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

$${\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(⟨{\partial }_a{\lambda }_0|{\lambda }_1⟩⟨{\lambda }_1|{\partial }_b{\lambda }_0⟩\right).$$

Full text

1. 1 Introduction

2. 2 Quantum Fisher information matrix

3. 2.1 Definition

4. 2.2 Parameterization processes

5. 2.3 Calculating QFIM

6. 2.3.1 General methods

7. 2.3.2 Pure states

8. 2.3.3 Few-qubit states

9. 2.3.4 Unitary processes

10. 2.3.5 Gaussian states

11. 2.4 QFIM and geometry of quantum mechanics

12. 2.4.1 Fubini-Study metric

13. 2.4.2 Fidelity and Bures metric

14. 2.4.3 Quantum geometric tensor

15. 2.5 QFIM and thermodynamics

16. 2.6 QFIM in quantum dynamics

17. 2.6.1 Quantum speed limit

18. 2.6.2 Non-Markovianity

19. 3 Quantum multiparameter estimation

20. 3.1 Quantum multiparameter Cramér-Rao bound

21. 3.1.1 Introduction

22. 3.1.2 Attainability

23. 3.1.3 Optimal measurements

24. 3.2 Phase estimation in the Mach-Zehnder interferometer

25. 3.2.1 Double-phase estimation

26. 3.2.2 Multi-phase estimation

27. 3.3 Waveform estimation

28. 3.4 Control-enhanced multiparameter estimation

29. 3.5 Estimation of a magnetic field

30. 3.6 Other applications and alternative mathematical tools

31. 4 Conclusion and outlook

32. A Derivation of traditional form of QFIM

33. B Derivation of QFIM for arbitrary-rank density matrices

34. C QFIM in Bloch representation

35. D One-qubit basis-independent expression of QFIM

36. E Derivation of SLD operator for Gaussian states

37. E.1 SLD operator for multimode Gaussian states

38. E.2 SLD operator for single-mode Gaussian state

39. F QFIM and Bures metric

40. G Relation between QFIM and cross-correlation functions

41. H Derivation of quantum multiparameter Cramér-Rao bound

42. I Construction of optimal measurement for pure states

43. J Gradient in GRAPE for Hamiltonian estimation

44. J.1 Gradient of CFIM

45. 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.

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]

[TABLE]

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)$.

[TABLE]

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]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

and the parameterized state can be written as

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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]

[TABLE]

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]

[TABLE]

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.

[TABLE]

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]

[TABLE]

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

[TABLE]

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

[TABLE]

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]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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]

[TABLE]

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]

[TABLE]

$\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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

The diagonal entry reads

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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]

[TABLE]

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]

[TABLE]

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]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.,

[TABLE]

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

[TABLE]

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

[TABLE]

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]

[TABLE]

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]

[TABLE]

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

[TABLE]

Here

[TABLE]

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

[TABLE]

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]

[TABLE]

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

[TABLE]

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.

[TABLE]

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

[TABLE]

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

[TABLE]

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]

[TABLE]

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

[TABLE]

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]

[TABLE]

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]

[TABLE]

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.

[TABLE]

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]

[TABLE]

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]

[TABLE]

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.,

[TABLE]

In the mean time, due to the fact that

[TABLE]

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

[TABLE]

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]

[TABLE]

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

[TABLE]

2.5 QFIM and thermodynamics

The density matrix of a quantum thermal state is

[TABLE]

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

[TABLE]

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

[TABLE]

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]

[TABLE]

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.

[TABLE]

or equivalently,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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]

[TABLE]

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,

[TABLE]

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,

[TABLE]

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

[TABLE]

the quantum Fisher information flow

[TABLE]

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

[TABLE]

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.

[TABLE]

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.,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

which is equivalent to the form

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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]

[TABLE]

which is equivalent to

[TABLE]

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]

[TABLE]

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.

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]

[TABLE]

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

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]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where the classical part

[TABLE]

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

[TABLE]

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

[TABLE]

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]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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]:

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.

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]

[TABLE]

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

[TABLE]

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]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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}}$)

[TABLE]

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

[TABLE]

Next, utlizing equation (1), the QFIM

[TABLE]

can be rewritten as

[TABLE]

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

[TABLE]

Substituting equation (3) into this equation, one has

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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}$)

[TABLE]

one can easily obtain

[TABLE]

The derivative of $\rho$ reads $\partial_{x_{a}}\rho=\sum^{N-1}_{i=0}\partial_{x_{a}}\lambda_{i}|\lambda_{i}\rangle\langle\lambda_{i}|+\lambda_{i}|\partial_{x_{a}}\lambda_{i}\rangle\langle\lambda_{i}|+\lambda_{i}|\lambda_{i}\rangle\langle\partial_{x_{a}}\lambda_{i}|$, which gives $\langle\lambda_{i}|\partial_{x_{a}}\rho|\lambda_{j}\rangle=\partial_{x_{a}}\lambda_{i}\delta_{ij}+\lambda_{j}\langle\lambda_{j}|\partial_{x_{a}}\lambda_{i}\rangle+\lambda_{i}\langle\partial_{x_{a}}\lambda_{j}|\lambda_{i}\rangle.$ $\delta_{ij}$ is the Kronecker delta function ($\delta_{ij}=1$ for $i=j$ and zero otherwise). Substituting this expression into above equation, $\langle\lambda_{i}|L_{x_{a}}|\lambda_{j}\rangle$ can be calculated as

[TABLE]

With the solution of SLD, we can now further calculate the QFIM. Utlizing equation (1), the QFIM can be rewritten into

[TABLE]

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

[TABLE]

In this equation, it can be seen that the arbitrary part of $\langle\lambda_{i}|L_{x_{a}}|\lambda_{j}\rangle$ does not affect the value of QFIM since it is not involved in above equation, but it provides a freedom for the optimal measurements, which will be further discussed later. Substituting equation (5) into above equation, we have

[TABLE]

where the fact $\langle\lambda_{i}|\partial_{x_{a}}\lambda_{j}\rangle=-\langle\partial_{x_{a}}\lambda_{i}|\lambda_{j}\rangle$ has been applied. The third term in above equation can be further calculated as

[TABLE]

Therefore equation (8) can be rewritten into

[TABLE]

Exchange the subscripts $i$ and $j$ in the second and fourth terms of above equation, it reduces into a symmetric form as shown in equation (2).

Appendix C QFIM in Bloch representation

Bloch representation is a well-used representation of quantum states in quantum mechanics and quantum information theory. For a $d$-dimensional quantum state $\rho$, it can be expressed by

[TABLE]

where $\vec{r}=(r_{1},r_{2}...,r_{k},...)^{\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$ and

[TABLE]

where $\delta_{ij}$ is the Kronecker delta function. $\mu_{ijm}$ and $\epsilon_{ijm}$ are the symmetric and antisymmetric structure constants, which can be calculated via the following equations

[TABLE]

It is easy to find $\mathrm{Tr}(\kappa_{i}\kappa_{j})=2\delta_{ij}$ and $\mu_{ijl}$ keeps still for any order of $ijl$, which will be used in the following.

Now we express the SLD operator with the $\mathfrak{su}(d)$ generators as [43, 44]

[TABLE]

where $z_{a}$ is a number and $\vec{y}_{a}$ is a vector. In the following we only use $z$ and $\vec{y}$ for convenience. From the equation $\partial_{x_{a}}\rho=\frac{1}{2}(\rho L+L\rho)$, one can find

[TABLE]

where $c=\sqrt{(d-1)/(2d)}$. Furthermore,

[TABLE]

The coefficients of 11 in equation (16) from both sides have to be the same, which gives

[TABLE]

which can also be obtained from the equation $\mathrm{Tr}(\rho L_{x_{a}})=0$. Similarly, the coefficients of $\kappa_{i}$ in equation (16) from both sides have to be the same, which gives

[TABLE]

This can also be obtained by multipling $\kappa_{i}$ on both sides of equation (16) and then taking the trace. Now we introduce the matrix $G$ with the entry

[TABLE]

It is easy to see $G$ is real symmetric. Substituting equation (11) into the equation above, it can be further written into

[TABLE]

Hence, the last term in equation (19) can be rewritten into

[TABLE]

Substituting this equation into equation (19), one has $c\partial_{x_{a}}r_{i}=czr_{i}+\frac{1}{2}\sum_{m}G_{im}y_{m}$, which immediately leads to

[TABLE]

Futhermore, one can check that

[TABLE]

Utilizing this equation, equation (23) reduces to

[TABLE]

Therefore,

[TABLE]

Next, we calculate the entry of QFIM and we will use the full notation $z_{a}$ and $\vec{y}_{a}$ instead of $z$, $\vec{y}$. The entry of QFIM reads

[TABLE]

Based on equations (18) and (23), one can obtain

[TABLE]

Hence, the QFIM can be written as

[TABLE]

Utilizing equation (26), the QFIM can be finally written into

[TABLE]

The theorem has been proved.

The simplest case here is single-qubit systems. The corresponding density matrix is $\rho=\frac{1}{2}(\leavevmode\hbox{\small 1\kern-4.3pt\normalsize 1}_{2}+\vec{r}\cdot\vec{\sigma})$ with $\vec{\sigma}=(\sigma_{1},\sigma_{2},\sigma_{3})$ the vector of Pauli matrices and $\leavevmode\hbox{\small 1\kern-4.3pt\normalsize 1}_{2}$ is the 2-dimensional identity matrix. In this case, $G=\leavevmode\hbox{\small 1\kern-4.3pt\normalsize 1}_{3}$ ($\leavevmode\hbox{\small 1\kern-4.3pt\normalsize 1}_{3}$ is the 3-dimensional identity matrix) due to the fact $\{\sigma_{i},\sigma_{j}\}=2\delta_{ij}\leavevmode\hbox{\small 1\kern-4.3pt\normalsize 1}_{2}$. Equation (30) then reduces to

[TABLE]

It can be checked that

[TABLE]

The QFIM then reads

[TABLE]

or equivalently,

[TABLE]

Appendix D One-qubit basis-independent expression of QFIM

This appendix shows the proof of Theorem 2.6. We first prove that the QFIM for one-qubit mixed state can be written as [46]

[TABLE]

Utilizing the spectral decomposition $\rho=\lambda_{0}|\lambda_{0}\rangle\langle\lambda_{0}|+\lambda_{1}|\lambda_{1}\rangle\langle\lambda_{1}|$, the first term reads

[TABLE]

The second term

[TABLE]

Since

[TABLE]

and

[TABLE]

one can finally obtain

[TABLE]

which coincides with the traditional formula of QFIM in Theorem 2.1. Equation (1) is then proved. Furthermore, one may notice that $\langle\lambda_{0}|\partial_{x_{a}}\rho|\lambda_{0}\rangle=\partial_{x_{a}}\lambda_{0}$, $\langle\lambda_{1}|\partial_{x_{a}}\rho|\lambda_{1}\rangle=\partial_{x_{a}}\lambda_{1}$, and $\langle\lambda_{0}|\partial_{x_{a}}\rho|\lambda_{1}\rangle=\langle\partial_{x_{a}}\lambda_{0}|\lambda_{1}\rangle+\langle\lambda_{0}|\partial_{x_{a}}\lambda_{1}\rangle=0$, which gives

[TABLE]

namely, the equality $\mathrm{Tr}\left[(\partial_{x_{a}}\rho)(\partial_{x_{b}}\rho)\right]=\mathrm{Tr}\left(\rho\{\partial_{x_{a}}\rho,\partial_{x_{b}}\rho)\}\right)$ holds for single-qubit mixed states. Thus, equation (1) can further reduce to

[TABLE]

The theorem has been proved.

Appendix E Derivation of SLD operator for Gaussian states

E.1 SLD operator for multimode Gaussian states

The derivation in this appendix is majorly based on the calculation in references [62, 67]. Let us first recall the notations before the derivation. The vector of quadrature operators are defined as $\vec{R}=(\hat{q}_{1},\hat{p}_{1},...,\hat{q}_{m},\hat{p}_{m})^{\mathrm{T}}$. $\Omega$ is a symplectic matrix $\Omega=i\sigma^{\oplus m}_{y}$. $\chi(\vec{s})$ is the characteristic function. Furthermore, denote $d=\langle\vec{R}\rangle$. To keep the calculation neat, we use $\chi$, $\rho$, $L$ instead of $\chi(\vec{s})$, $\rho$, $L_{x_{a}}$ in this appendix. Some other notations are $\langle\cdot\rangle=\mathrm{Tr}(\rho\cdot)$ and $\dot{A}=\partial_{x_{a}}A$.

The characteristic function $\chi=\langle D\rangle$, where $D=e^{i\vec{R}^{\mathrm{T}}\Omega\vec{s}}$ with $\vec{s}\in\mathbb{R}^{2m}$. Substituting $D$ into the equation $\partial_{x_{a}}\rho=\frac{1}{2}(\rho L+L\rho)$ and taking the trace, we obtain

[TABLE]

Next, assume the SLD operator is in the following form

[TABLE]

where $L^{(0)}$ is a real number, $\vec{L}^{(1)}$ is a $2m$-dimensional real vector and $G$ is a $2m$-dimensional real symmetric matrix. These conditions promise the Hermiticity of SLD. With this ansatz, equation (1) can then be rewritten into

[TABLE]

Now we calculate $\langle\{R_{i},D\}\rangle$ in equation (3). Denote $[A,\cdot]=A^{\times}(\cdot)$, one can have

[TABLE]

where we have used the fact $\Omega_{ij}=-\Omega_{ji}$ and $\sum_{j}\Omega_{ij}\Omega_{jk}=-\delta_{ik}$ which comes from the equation $\Omega^{2}=-\leavevmode\hbox{\small 1\kern-4.3pt\normalsize 1}_{2m}$. Substituting the equation

[TABLE]

into the first line of equation (4) and repeat the calculation, one can obtain $\partial_{s_{k}}D=-i\sum_{i}\Omega_{ki}D\left(\frac{1}{2}s_{i}+R_{i}\right)$. Combining this equation with equation (4), $\partial_{s_{k}}D$ can be finally written as

[TABLE]

which further gives $\langle\partial_{s_{k}}D\rangle=-\frac{i}{2}\sum_{i}\Omega_{ki}\langle\{R_{i},D\}\rangle$. Based on this equation, it can be found $\sum_{k}\Omega_{jk}\langle\partial_{s_{k}}D\rangle=-\frac{i}{2}\sum_{ik}\Omega_{jk}\Omega_{ki}\langle\{R_{i},D\}\rangle$. Again since $\sum_{k}\Omega_{jk}\Omega_{ki}=-\delta_{ij}$, it reduces to

[TABLE]

Next, continue to take the derivative on $\partial_{s_{k}}D$, we have

[TABLE]

Due to the Baker-Campbell-Hausdorff formula, $D^{\dagger}R_{i}D=R_{i}+[-i\vec{R}^{\mathrm{T}}\Omega\vec{s},R_{i}]=R_{i}+s_{i}$, the commutation between $R_{i}$ and $D$ can be obtained as $[R_{i},D]=s_{i}D$, which further gives $\{R_{i},\{R_{j},D\}\}=\{\{R_{i},R_{j}\},D\}-s_{i}s_{j}D$. Then equation (8) can be rewritten into

[TABLE]

And one can finally obtain

[TABLE]

With equations (7) and (10), equation (3) can be expressed by

[TABLE]

On the other hand, from the expression of characteristic function

[TABLE]

it can be found that

[TABLE]

and

[TABLE]

Then we have $-i\sum_{ik}L^{(1)}_{i}\Omega_{ik}(\partial_{s_{k}}\chi)=\sum_{i}L^{(1)}_{i}\left(-i\sum_{jk}\Omega_{kj}C_{ij}s_{k}+d_{i}\right)\chi$, and

[TABLE]

Substituting these equations into equation (11), $\partial_{x_{a}}\chi$ can be expressed by

[TABLE]

This equation can be written into a more compact way as below

[TABLE]

Compare this equation with equation (13), it can be found that

[TABLE]

From the second equation above, one can obtain

[TABLE]

Using this equation and the equation (16), it can be found that

[TABLE]

Once we obtain the expression of $G$, $L^{(0)}$ and $\vec{L}^{(1)}$ can be obtained correspondingly, which means we need to solve equation (18), which can be rewritten into following form

[TABLE]

Since $C=SC_{\mathrm{d}}S^{\mathrm{T}}$ and $S\Omega S^{\mathrm{T}}=\Omega$, above equation can be rewritten into

[TABLE]

where $G_{\mathrm{s}}=S^{\mathrm{T}}GS$. Denote the map $\mathcal{E}(G_{\mathrm{s}}):=\Omega G_{\mathrm{s}}\Omega+4C_{\mathrm{d}}G_{\mathrm{s}}C_{\mathrm{d}}$, then $G_{\mathrm{s}}=\mathcal{E}^{-1}(\mathcal{E}(G_{\mathrm{s}}))$. The map $\mathcal{E}(G_{\mathrm{s}})$ can be decomposed via the generators $\{A^{(jk)}_{l}\}$ ($j,k=1,...,m$), where

[TABLE]

for $l=0,1,2,3.$ $\sigma^{(jk)}_{i}$ is a $2m$-dimensional matrix with all the entries zero expect the $2\times 2$ block shown as below

[TABLE]

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}$. $A^{(jk)}_{l}$ satisfies the orthogonal relation $\mathrm{Tr}(A^{(jk)}_{l}A^{(j^{\prime}k^{\prime})}_{l^{\prime}})=\delta_{jj^{\prime}}\delta_{kk^{\prime}}\delta_{ll^{\prime}}$. Recall that $C_{\mathrm{d}}=\bigoplus^{m}_{k=1}c_{k}\leavevmode\hbox{\small 1\kern-4.3pt\normalsize 1}_{2}$, it is easy to check that

[TABLE]

Next decompose $S^{-1}\dot{C}(S^{\mathrm{T}})^{-1}$ with $\{A^{(jk)}_{l}\}$ as

[TABLE]

where $g^{(jk)}_{l}=\mathrm{Tr}[S^{-1}\dot{C}(S^{\mathrm{T}})^{-1}A^{(jk)}_{l}]$. Decomposing $G_{\mathrm{s}}$ as $G_{\mathrm{s}}=\sum_{jkl}\tilde{g}^{(jk)}_{l}A^{(jk)}_{l}$, and substituting it into equation (22), we have

[TABLE]

Utilizing equations (25) and (26), above equation reduces to

[TABLE]

which indicates

[TABLE]

Thus, $G=(S^{\mathrm{T}})^{-1}G_{\mathrm{s}}S^{-1}$ can be solved as below

[TABLE]

In summary, the SLD can be expressed by

[TABLE]

where

[TABLE]

with

[TABLE]

for $l=0,1,2,3$ and $g^{(jk)}_{l}=\mathrm{Tr}[S^{-1}\dot{C}(S^{\mathrm{T}})^{-1}A^{(jk)}_{l}]$. And

[TABLE]

E.2 SLD operator for single-mode Gaussian state

For a single-mode case, a 2-dimensional symplectic $S$ matrix satisfies $\det S=1$. Based on the equation $C=SC_{\mathrm{d}}S^{\mathrm{T}}=cSS^{\mathrm{T}}$ ($c$ is the symplectic value), the following equations can be obtained

[TABLE]

Together with the equation $\det S=S_{00}S_{11}-S_{01}S_{10}=1$, the symplectic value $c$ can be solved as $c=\sqrt{\det C}$. Assume the form

[TABLE]

Substituting above expressions into equation (37), it can be found $\theta,\phi$ need to satisfy

[TABLE]

Since there is only four constrains for these five variables, one of them is free. We take $\theta=\pi/2-\phi$, and the symplectic matrix reduces to

[TABLE]

where $\phi$ satisfies $\sin(2\phi)=\frac{C_{01}}{\sqrt{C_{00}C_{11}}}$ and $\cos(2\phi)=-\frac{c}{\sqrt{C_{00}C_{11}}}$. Based on Theorem 2.9, we obtain

[TABLE]

where $\dot{C}$, $\dot{c}$ are short for $\partial_{x_{a}}C$ and $\partial_{x_{a}}c$. Meanwhile, we have

[TABLE]

These expressions immediately give $G_{x_{a}}$ as below

[TABLE]

Define a matrix

[TABLE]

$[G_{x_{a}}]_{00}$ can be rewritten into

[TABLE]

Similarly, one can obtain

[TABLE]

and

[TABLE]

These expressions indicate

[TABLE]

with

[TABLE]

Appendix F QFIM and Bures metric

The Bures distance between two quantum states $\rho_{1}$ and $\rho_{2}$ is defined as

[TABLE]

where $f(\rho_{1},\rho_{2})=\mathrm{Tr}\sqrt{\sqrt{\rho_{1}}\rho_{2}\sqrt{\rho_{1}}}$ is the quantum fidelity. Now we calcuate the fidelity for two close quantum states $\rho(\vec{x})$ and $\rho(\vec{x}+\mathrm{d}\vec{x})$. The Taylor series of $\rho(\vec{x}+\mathrm{d}\vec{x})$ (up to the second order) reads

[TABLE]

In the following $\rho$ will be used as the abbreviation of $\rho(\vec{x})$. Utilizing the equation above, one can obtain

[TABLE]

Now we assume

[TABLE]

Taking the square of above equation and compare it to equation (3), one can obtain

[TABLE]

where $\partial^{2}\rho$ is short for $\frac{\partial^{2}\rho}{\partial x_{a}\partial x_{b}}$. $\frac{1}{2}\{W_{a},W_{b}\}$ (not $W_{a}W_{b}$) is used to make sure the equation is unchanged when the subscripts $a$ and $b$ exchange. Meawhile, from equation (4), the Bures metric $D_{\mathrm{B}}(\rho(\vec{x}),\rho(\vec{x}+\mathrm{d}\vec{x}))$ (the abbreviation $D_{\mathrm{B}}$ will bu used below) can be calculated as

[TABLE]

As long as we obtain the specific expressions of $W_{a}$ and $Y_{ab}$ from equations (5) and (6), $D_{\mathrm{B}}$ can be obtained immediately. To do that, we utilize the spectral decomposition $\rho=\sum_{i}\lambda_{i}|\lambda_{i}\rangle\langle\lambda_{i}|$. In the basis $\{|\lambda_{i}\rangle\}$, the matrix entries ($[\cdot]_{ij}=\langle\lambda_{i}|\cdot|\lambda_{j}\rangle$) read

[TABLE]

Here $[\partial_{x_{a}}\rho]_{ij}=\delta_{ij}\partial_{x_{a}}\lambda_{i}-(\lambda_{i}-\lambda_{j})\langle\lambda_{i}|\partial_{x_{a}}\lambda_{j}\rangle$, and

[TABLE]

where $\partial^{2}\lambda_{i}$ and $|\partial^{2}\lambda_{i}\rangle$ are short for $\frac{\partial^{2}\lambda_{i}}{\partial x_{a}\partial x_{b}}$ and $\frac{\partial^{2}}{\partial x_{a}\partial x_{b}}|\lambda_{i}\rangle$. Now we denote $\rho$’s dimension as $N$ and $\lambda_{i}\in\mathcal{S}$ for $i=0,1,2...,M-1$. Under the assumption that the support $\mathcal{S}$ is not affected by the values of $\vec{x}$, i.e., the rank of $\rho(\vec{x})$ equals to that of $\rho(\vec{x}+\mathrm{d}\vec{x})$, $\sqrt{\rho}\partial_{x_{a}}\rho\sqrt{\rho}$ and $\sqrt{\rho}\partial^{2}\rho\sqrt{\rho}$ are both block diagonal. Based on equation (5), the $ij$th matrix entry of $W_{a}$ can be calculated as

[TABLE]

for $i,j\in[0,M-1]$ and $[W_{a}]_{ij}=0$ for others. One can observe that $W_{a}$ is a Hermitian matrix, and

[TABLE]

which means there is no first order term in Bures metric. With respect to the second order term, we need to know the value of $[Y_{ab}]_{ii}$. From equation (6), one can obtain

[TABLE]

Due to the fact $\langle\partial^{2}\lambda_{i}|\lambda_{i}\rangle+\langle\lambda_{i}|\partial^{2}\lambda_{i}\rangle=-2\mathrm{Re}(\langle\partial_{x_{a}}\lambda_{i}|\partial_{x_{b}}\lambda_{i}\rangle)$, one can have

[TABLE]

which further gives

[TABLE]

where the fact $\sum_{i}\partial^{2}\lambda_{i}=0$ has been applied. Next, from equation (11) one can obtain

[TABLE]

which means

[TABLE]

Thus, $\mathrm{Tr}Y_{ab}$ can then be expressed by

[TABLE]

With this equation, we finally obtain

[TABLE]

One should notice that this proof shows that this relation is established for density matrices with any rank as long as the rank of $\rho(\vec{x})$ is unchanged with the varying of $\vec{x}$. In the case the rank can change, a thorough discussion can be found in reference [87].

Appendix G Relation between QFIM and cross-correlation functions

This appendix gives the thorough calculation of the relation between QFIM and dynamic susceptibility in reference [109]. Consider the unitary parameterization $U=\exp(i\sum_{a}x_{a}O_{a})$ with a thermal state $\rho=\frac{1}{Z}e^{-\beta H}$. Here $Z=\mathrm{Tr}(e^{-\beta H})$ is the partition function. $O_{a}$ is a Hermitian generator for $x_{a}$. In the following we set $k_{B}=1$ and assume all $O_{a}$ are commutative, i.e., $[O_{a},O_{b}]=0$ for any $a$ and $b$. Denote $O_{a}(t)=e^{iHt}O_{a}e^{-iHt}$, and $\langle\cdot\rangle=\mathrm{Tr}(\rho\cdot)$, the symmetric cross-correlation spectrum in this case reads

[TABLE]

Utilizing the energy basis $\{|E_{i}\rangle\}$ (with $E_{i}$ the $i$th energy), it can be rewritten into

[TABLE]

Further use the equation $\int^{\infty}_{-\infty}e^{i(\omega+E_{i}-E_{j})}\mathrm{d}t=2\pi\delta(\omega+E_{i}-E_{j})$, $S_{ab}(\omega)$ can reduce to

[TABLE]

With this expression, one can find

[TABLE]

Since

[TABLE]

one can obtain

[TABLE]

From the expression of QFIM, it can be found that

[TABLE]

It can also be checked that

[TABLE]

In the mean time, the asymmetric cross-correlation spectrum is in the form

[TABLE]

which directly gives the relation between $\chi_{ab}$ and $S_{ab}$ as below

[TABLE]

namely,

[TABLE]

which is just the fluctuation-dissipation theorem. Using this relation, one can further obtain the result in reference [109] as below

[TABLE]

Appendix H Derivation of quantum multiparameter Cramér-Rao bound

The derivation of quantum multiparameter Cramér-Rao bound is based on the Cauchy-Schwarz inequality below

[TABLE]

which comes from the complete form

[TABLE]

Define $X$ and $Y$ as

[TABLE]

where $f_{m},g_{m}$ are real numbers for any $m$, and $\langle\cdot\rangle=\mathrm{Tr}(\cdot\rho)$. Here $O_{m}$ is an observable defined as

[TABLE]

$\hat{x}_{m}(k)$ is the estimator of $x_{m}$ and is the function of $k$th result. Based on equation (3), it can be calculated that $\mathrm{Tr}(X^{\dagger}X)=\sum_{ml}f_{m}f_{l}\mathrm{Tr}(L_{m}L_{l}\rho)$, which can be symmetrized into $\mathrm{Tr}(X^{\dagger}X)=\frac{1}{2}\sum_{ml}f_{m}f_{l}\mathrm{Tr}(\{L_{m},L_{l}\}\rho)=\sum_{ml}f_{m}f_{l}\mathcal{F}_{ab}$. Define $\vec{f}=(f_{0},f_{1},...,f_{m},...)^{\mathrm{T}}$, $\mathrm{Tr}(X^{\dagger}X)$ can be further rewritten into

[TABLE]

Similarly, define $\vec{g}=(g_{0},g_{1},...,g_{m},...)^{\mathrm{T}}$ and through some algebra, $\mathrm{Tr}(Y^{\dagger}Y)$ can be calculated as

[TABLE]

where $C$ is the covariance matrix for $\{O_{m}\}$, and is defined as $C_{ml}=\frac{1}{2}\langle\{O_{m},O_{l}\}\rangle-\langle O_{m}\rangle\langle O_{l}\rangle$. Furthermore, one can also obtain

[TABLE]

where the entry of $B$ is defined as $B_{ml}=\frac{1}{2}\mathrm{Tr}(\rho\{L_{m},\Pi_{l}\})=\frac{1}{2}\mathrm{Tr}(\{\rho,L_{m}\}O_{l})$. Utilizing $\partial_{x_{a}}\rho=\frac{1}{2}\{\rho,L_{m}\}$, $B_{ml}$ reduces to $B_{ml}=\mathrm{Tr}(\Pi_{l}\partial_{x_{m}}\rho)$. Since we consider unbiased estimators, i.e., $\langle O_{m}\rangle=\sum_{m}\hat{x}_{m}\mathrm{Tr}(\rho\Pi_{k})=x_{m}$, $B_{ml}$ further reduces to $B_{ml}=\partial_{x_{m}}\langle O_{l}\rangle=\delta_{ml}$, with $\delta_{ml}$ the Kronecker delta function. Hence, for unbiased estimators $B$ is actually the identity matrix 11.

Now substituting equations (6), (7) and (8) into the Cauchy-Schwarz inequality (1), one can obtain

[TABLE]

Assuming $\mathcal{F}$ is invertable, i.e., it is positive-definite, and taking $\vec{f}=\mathcal{F}^{-1}\vec{g}$, above inequality reduces to $\vec{g}^{\,\mathrm{T}}\mathcal{F}^{-1}\vec{g}\vec{g}^{\,\mathrm{T}}C\vec{g}\geq\left(\vec{g}^{\,\mathrm{T}}\mathcal{F}^{-1}\vec{g}\right)^{2}$. Since $\mathcal{F}$ is positive-definite, $\mathcal{F}^{-1}$ is also positive-definite, which means $\vec{g}^{\,\mathrm{T}}\mathcal{F}^{-1}\vec{g}$ is a positive number, thus, the above equation can further reduce to $\vec{g}^{\,\mathrm{T}}C\vec{g}\geq\vec{g}^{\,\mathrm{T}}\mathcal{F}^{-1}\vec{g}$, namely,

[TABLE]

Next we discuss the relation between $C$ and $\mathrm{cov}(\hat{\vec{x}},\{\Pi_{k}\})$. Utilizing the defintion of $O_{m}$, $C_{ml}$ can be written as

[TABLE]

and $\mathrm{cov}(\hat{\vec{x}},\{\Pi_{k}\})$ for unbiased estimators reads

[TABLE]

If $\{\Pi_{k}\}$ is a set of projection operators, it satisfies $\Pi_{k}\Pi_{k^{\prime}}=\Pi_{k}\delta_{kk^{\prime}}$. For unbiased estimators, $\sum_{k}\hat{x}_{m(l)}\mathrm{Tr}(\rho\Pi_{k})=x_{m(l)}$, which gives $C_{ml}=\sum_{k}\hat{x}_{m}\hat{x}_{l}\mathrm{Tr}(\rho\Pi_{k})-x_{m}x_{l}$, i.e.,

[TABLE]

If $\{\Pi_{k}\}$ is a set of POVM, one can see

[TABLE]

In the mean time,

[TABLE]

Now define $\mathcal{B}_{k}:=\sqrt{\Pi_{k}}\left(\sum_{m}g_{m}\hat{x}_{m}-\sum_{m}g_{m}O_{m}\right)$. Based on the Cauchy-Schwarz inequality $\mathrm{Tr}(A^{\dagger}A)\geq 0$, which is valid for any operator $A$, one can immediately obtain $\mathrm{Tr}(\rho\mathcal{B}^{\dagger}_{k}\mathcal{B}_{k})=\mathrm{Tr}(\sqrt{\rho}\mathcal{B}^{\dagger}_{k}\mathcal{B}_{k}\sqrt{\rho})\geq 0$, which further gives $\sum_{k}\mathrm{Tr}(\rho\mathcal{B}^{\dagger}_{k}\mathcal{B}_{k})\geq 0$. Through some calculations, the expression of $\sum_{k}\mathrm{Tr}(\rho\mathcal{B}^{\dagger}_{k}\mathcal{B}_{k})$ is in the form

[TABLE]

Now taking the difference of equations (14) and (15), it can be found

[TABLE]

Finally, we obtain the following inequality $\vec{g}^{\,\mathrm{T}}(\mathrm{cov}(\hat{\vec{x}},\{\Pi_{k}\})-C)\vec{g}\geq 0$, namely, for any POVM measurement,

[TABLE]

Based on inequality (10) and the property of quadratic form, we finally obtain $\mathrm{cov}(\hat{\vec{x}},\{\Pi_{k}\})\geq\mathcal{F}^{-1}$. Consider the repetition of experiments (denoted as $n$), above bound needs to add a factor of $1/n$. Hence the quantum multiparameter Cramér-Rao bound can be finally expressed by

[TABLE]

Appendix I Construction of optimal measurement for pure states

Assume the true values of the vector of unknown parameters $\vec{x}$ is $\vec{x}_{\mathrm{true}}$, we now provide the proof that for a pure parameterized state $|\psi\rangle$, a set of projectors containing the state $|\psi_{\vec{x}_{\mathrm{true}}}\rangle:=|\psi(\vec{x}=\vec{x}_{\mathrm{true}})\rangle$, i.e., $\{|m_{k}\rangle\langle m_{k},|m_{0}\rangle=|\psi_{\vec{x}_{\mathrm{true}}}\rangle\}$ is possible to be an optimal measurement to attain the quantum Cramér-Rao bound, as shown in reference [151, 136]. The calculation in this appendix basically coincides with the appendix in reference [136].

Since $\{|m_{k}\rangle\langle m_{k}|\}$ contains the information of the true value, in practice one need to use the estimated value of $\vec{x}$ (denoted as $\hat{\vec{x}}$) to construct $|m_{0}\rangle=|\psi(\hat{\vec{x}})\rangle$ to perform this measurement, then improve the accuracy of $\hat{\vec{x}}$ adaptively. Thus, it is reasonable to assume $|\psi_{\vec{x}_{\mathrm{true}}}\rangle=|m_{0}\rangle+\sum_{x_{a}}\delta x_{a}|\partial_{x_{a}}\psi\rangle|_{\vec{x}=\hat{\vec{x}}}$. The probability for $|m_{k}\rangle\langle m_{k}|$ is $p_{k}=|\langle\psi|m_{k}\rangle|^{2}$, which gives the CFIM as below

[TABLE]

At the limit $\hat{\vec{x}}\rightarrow\vec{x}_{\mathrm{true}}$, it is

[TABLE]

For the $k=0$ term,

[TABLE]

and $\lim_{\vec{x}\rightarrow\vec{x}_{\mathrm{true}}}\langle\psi|m_{0}\rangle=1$. Thus, the CFIM is

[TABLE]

Due to the fact

[TABLE]

one can have

[TABLE]

where

[TABLE]

To let $\mathcal{I}(\vec{x}_{\mathrm{true}})=\mathcal{F}(\vec{x}_{\mathrm{true}})$, $Q$ has to be a zero matrix. Since the diagonal entry

[TABLE]

in which all the terms within the summation are non-negative. Therefore, its value is zero if and only if

[TABLE]

Furthermore, this condition simultaneously makes the non-diagonal entries of $Q$ vanish. Thus, it is the necessary and sufficient condition for $Q=0$, which means it is also the necessary and sufficient condition for $\mathcal{I}=\mathcal{F}$ at the point of true value. One may notice that the limitation in equation (28) is a $0/0$ type. Thus, we use the formula $|\psi_{\vec{x}_{\mathrm{true}}}\rangle=|m_{0}\rangle+\sum_{x_{j}}\delta x_{j}|\partial_{x_{j}}\psi\rangle|_{\vec{x}=\hat{\vec{x}}}$ to further calculate above equation. With this formula, one has

[TABLE]

$\langle\partial_{x_{j}}\psi|m_{k}\rangle$ cannot generally be zero for all $x_{j}$ since all $\partial_{x_{j}}\psi\rangle$ are not orthogonal in general. Thus, equation (28) is equivalent to

[TABLE]

Appendix J Gradient in GRAPE for Hamiltonian estimation

J.1 Gradient of CFIM

The core of GRAPE algorithm is to obtain the expression of gradient. The dynamics of the system is described by

[TABLE]

For the Hamiltonian estimation under control, the dynamics is

[TABLE]

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. To perform the algorithm, 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.

For a set of probability distribution $p(y|\vec{x})=\mathrm{Tr}(\rho\Pi_{y})$ with $\{\Pi_{y}\}$ a set of POVM. The gradient of $\mathcal{I}_{ab}$ at target time $T$ reads [172]

[TABLE]

where

[TABLE]

and $\mathcal{M}^{(1)}_{j}$, $\mathcal{M}^{(2)}_{j,a(b)}$ and $\mathcal{M}^{(3)}_{j,a(b)}$ are Hermitian operators and can be expressed by

[TABLE]

The notation $A^{\times}(\cdot):=[A,\cdot]$ is a superoperator. $\delta_{jm}$ is the Kronecker delta function. $D^{j^{\prime}}_{j}$ is the propagating superoperator from the $j$th time point to the $j^{\prime}$th time with the definition $D^{j^{\prime}}_{j}:=\prod^{j^{\prime}}_{i=j}\exp(\Delta t\mathcal{E}_{i})$ for $j\leq j^{\prime}$. We define $D^{j^{\prime}}_{j}=\leavevmode\hbox{\small 1\kern-4.3pt\normalsize 1}$ for $j>j^{\prime}$. $\rho_{j}=D^{j}_{1}\rho(0)$ is the quantum state at $j$th time.

For a two-parameter case $\vec{x}=(x_{0},x_{1})$, the objective function can be chosen as $F_{\mathrm{eff}}(T)$ according to corollary 3.1.3, and the corresponding gradient is

[TABLE]

For the objective function

[TABLE]

the gradient reads

[TABLE]

J.2 Gradient of QFIM

Now we calculate the gradient of the QFIM. Based on the equation

[TABLE]

we can obtain

[TABLE]

Similarly, we have

[TABLE]

Next, from the definition of QFIM, the gradient for $\mathcal{F}_{ab}$ at target time $T$ reads [172]

[TABLE]

Combing equations (13), (14) and (15), one can obtain

[TABLE]

Substituing the specific expressions of $\frac{\delta\rho(T)}{\delta V_{k}(j)}$ given in reference [172], one can obtain the gradient of $\mathcal{F}_{ab}$ as below

[TABLE]

The gradient for the diagonal entry $\mathcal{F}_{aa}$ reduces to the form in reference [172], i.e.,

[TABLE]