# On the existence of unbiased resilient estimators in discrete quantum
systems

Javier Navarro<EMAIL_ADDRESS>BCAM - Basque Center for Applied
Mathematics, Mazarredo 14, 48009 Bilbao, Spain  Ricard Ravell Rodríguez
<EMAIL_ADDRESS>BCAM - Basque Center for Applied Mathematics, Mazarredo
14, 48009 Bilbao, Spain  Mikel Sanz<EMAIL_ADDRESS>BCAM - Basque Center
for Applied Mathematics, Mazarredo 14, 48009 Bilbao, Spain  Department of
Physical Chemistry, University of the Basque Country UPV/EHU, Apartado 644,
48080 Bilbao, Spain  EHU Quantum Center, University of the Basque Country
UPV/EHU, Apartado 644, 48080 Bilbao, Spain  IKERBASQUE, Basque Foundation for
Science, Plaza Euskadi 5, 48009, Bilbao, Spain

###### Abstract

Cramér-Rao constitutes a crucial lower bound for the mean squared error of an
estimator in frequentist parameter estimation, albeit paradoxically demanding
highly accurate prior knowledge of the parameter to be estimated. Indeed, this
information is needed to construct the optimal unbiased estimator, which is
highly dependent on the parameter. Conversely, Bhattacharyya bounds result in
a more resilient estimation about prior accuracy by imposing additional
constraints on the estimator. Initially, we conduct a quantitative comparison
of the performance between Cramér-Rao and Bhattacharyya bounds when faced with
less-than-ideal prior knowledge of the parameter. Furthermore, we demonstrate
that the $n^{th}$order classical and quantum Bhattacharyya bounds cannot be
computed –given the absence of estimators satisfying the constraints– under
specific conditions tied to the dimension $m$ of the discrete system.
Intriguingly, for a system with the same dimension $m$, the maximum non-
trivial order $n$ is $m-1$ in the classical case, while in the quantum realm,
it extends to $m(m+1)/2-1$. Consequently, for a given system dimension, one
can construct estimators in quantum systems that exhibit increased robustness
to prior ignorance.

##  I Introduction

Quantum metrology studies the limits in precision by estimating physical
parameters encoded in quantum systems. Quantum properties such as entanglement
and squeezing are used to surpass classical estimation strategies [1, 2, 3].
Applications like quantum phase estimation [4, 5, 6, 7], quantum-enhanced
position and velocity estimation [8, 9, 10], quantum illumination [11, 12,
13], quantum thermometry [14, 15, 16, 17], and channel discrimination [18]
among others have boosted the field. Parameter estimation gives the
theoretical framework to establish the optimal performance of particular
estimation protocols and its main figure of merit is the Cramér-Rao bound
(CRB) [19, 20].

Paradoxically, to achieve the ultimate precision established by the CRB an
infinite number of measurements have to be performed or the prior value of the
parameter has to be known. Both scenarios are unlikely or experimentally
unfeasible. In the finite sample regime, CRB is frequently unattainable [21,
22, 23]. Several bounds as Hammersley–Chapman–Robbins bound [24], the family
of Bhattacharyya bounds (BhBB) [25] or Barankin [26, 27] bounds, and more
generally Abel bounds [28], which include the previous bounds as special
cases, were introduced. These bounds succeed in establishing a sensible limit
for the achievable precision [29], especially in noisy estimation problems
[27].

Recently, there has been interest in developing tighter bounds for quantum
parameter estimation tasks. Gessner et al. [30] derived lower bounds on the
variance of estimators in quantum metrology by increasing the unbiasedness
constraints of the estimator, i.e., in the same spirit as their classical
analogs were obtained. They presented a hierarchy of increasingly tight bounds
that include the quantum CRB (QCRB) [31] at the lowest order and quantum
Barankin, Bhattacharyya, and Abel bounds as other particular cases.

In this article, we study the classical and the quantum Bhattacharyya bounds
(QBhBB). We compare the performance of these bounds with the CRB when the
prior knowledge of the parameter is not sharp deriving a lower bound for the
mean squared error (MSE). We illustrate this result for the Mach-Zehnder
interferometer by comparing the actual MSE of the estimator with BhB and CRB
and we see a clear advantage when using the BhB. Finally, we address the
impossibility of computing classical BhB under certain conditions when dealing
with discrete probability distributions. We prove the necessary and sufficient
conditions for this bound to exist and we show how these conditions get
relaxed for the QBhB.

This article is organized as follows: in sec. II, we introduce estimation
theory, present the BhB and motivate them using the Mach-Zehnder
interferometer setting. In sec. III, we study the cases when the BhB cannot be
computed, i.e. the bound is divergent for both, classical probability
distributions and quantum systems. Finally, in sec. IV, the conclusions are
presented.

##  II Estimation theory

To extract information from a physical system measurements have to be
performed. The measurement outcome is modeled as a random variable following a
probability distribution that depends on the state of the system. The
information is encoded in the random variable and estimation theory is the
branch of statistics that deals with the optimal protocols and the ultimate
limit for the information extraction. The CRB establishes a limit on the
precision while estimating a magnitude or a parameter. Consider that the
outcome of the experiment, denoted by the random variable $X$, follows the
family of probability distribution $P_{\theta}(x)$, that is parameterized by
the parameter $\theta$ which is being estimated. An estimator
$\tilde{\Theta}(x)$ is a function of the random measurement outcome that
estimates the value of the parameter $\theta$. The performance of an estimator
is characterized by its bias $b(\theta)$, mean squared error (MSE) and
variance $(\Delta\tilde{\Theta})^{2}$,

$\displaystyle b(\theta)$
$\displaystyle=\int(\theta-\tilde{\Theta}(x))P_{\theta}(x)\>dx,$ (1)
$\displaystyle(\Delta\tilde{\Theta})^{2}$
$\displaystyle=\int(\langle\tilde{\Theta}\rangle_{\theta}-\tilde{\Theta}(x))^{2}P_{\theta}(x)\>dx,$
(2) MSE $\displaystyle=\int(\theta-\tilde{\Theta}(x))^{2}P_{\theta}(x)\>dx.$
(3)

The mean value is defined as $\langle f\rangle_{\theta}=\int
f(x)P_{\theta}(x)\>dx$. An estimator $\tilde{\Theta}_{\theta}(x)$ is locally
unbiased on the point $\theta$ when
$\langle\tilde{\Theta}_{\theta}\rangle_{\theta}=\theta$ and
$\frac{db(\theta)}{d\theta}|_{\theta=\theta}=0$. Note that for unbiased
estimators MSE and the variance coincide. Unbiased Cramér-Rao, usually named
as the Cramér-Rao bound (CRB), states that for any unbiased estimator

$(\Delta\tilde{\Theta}_{\theta})^{2}\geq\frac{1}{F_{C}},$ (4)

where $F_{C}$ is the classical Fisher information and is given by

$F_{C}=\int_{x\in X_{+}}\left(\frac{\partial
P_{\theta}}{\partial\theta}\right)^{2}\frac{1}{P_{\theta}(x)}\>dx,$ (5)

where $X_{+}$ is the space of events $x$ with nonzero probability
$P_{\theta}(x)>0$.

The CRB is always saturable in the asymptotic limit or when the true value of
$\theta$ is known [32]. This is not the case in realistic scenarios. In
quantum parameter estimation [1], a family of quantum states $\rho_{\theta}$
parametrized by a parameter $\theta$ plus a particular POVM measurement
$\\{\Pi_{x}\\}$ substitute the family of probability distributions,
$P_{\theta}(x)$. The lower bound on the MSE of the estimator is now determined
by the quantum Cramér-Rao bound (QCRB) [1, 31, 14], its version for locally
unbiased estimators reads

$(\Delta\tilde{\Theta}_{\theta})^{2}\geq\frac{1}{F_{Q}},$ (6)

where $F_{Q}=\operatorname{Tr}[\rho_{\theta}L^{2}]$ is the quantum Fisher
information. The symmetric logarithmic derivative (SLD) is the operator $L$
satisfying

$\frac{d\rho_{\theta}}{d\theta}=\frac{\rho_{\theta}L+L\rho_{\theta}}{2}.$ (7)

From Eq.(7) it is deduced that the SLD is a hermitian operator and its
eigenvectors can form an orthonormal basis. In the derivation of the QCRB, a
general POVM $\\{\Pi_{x}\\}$ was considered but it was shown that a projection
on the basis of eigenvectors of $L$ is the optimal measurement that saturates
the QCRB [14]. As in the classical estimation theory, the QCRB is always
saturable in the asymptotic regime. Furthermore, it has recently been shown to
be non-saturable within a Bayesian perspective [33, 23].

###  II.1 Bhattacharyya bounds and their motivation

Bhattacharyya bounds were introduced in the context of classical parameter
estimation problems to give a lower bound for the variance of the estimate
[25]. The BhB gives a tighter bound than the CRB by including higher-order
derivatives of the probability distribution in the calculation. It is shown
that BhB converge to the variance of the best unbiased estimator –the
estimator that is unbiased for any point in the region of the allowed values
of the parameter–, when the sampling distribution is a member of an
exponential family of distributions [34]. The $n^{th}$order BhB bound gives
the lowest value of the variance for an estimator that fulfills
$\frac{d^{i}b(\theta)}{d\theta^{i}}|_{\theta=\theta_{0}}=0$ for $1\leq i\leq
n$. This estimator is unbiased in a larger region than the one obtained for
the CRB, as will be shown later, and its variance is closer to the real
variance of a reasonable parameter estimation problem, namely where the real
value of the parameter is unknown, and is given by [30],

$(\Delta\tilde{\Theta}_{\theta_{0}})^{2}\geq\max_{\bm{a}}\frac{\left(\bm{a}^{\top}\bm{\lambda}\right)^{2}}{\bm{a}^{\top}\bm{C}\bm{a}}=\bm{\lambda}^{\top}\bm{C}^{-1}\bm{\lambda}.$
(8)

The last equality is valid only when $\bm{C}^{-1}$ exists; whereas if this is
not fulfilled, one needs to resort to the maximization problem. The vectors
$\bm{\lambda}$ and $\bm{a}$ have $n$ entries and $\bm{C}$ is a $n\times n$
matrix, where $n$ is the order of the bound we are computing. They are defined
as [30]

$\displaystyle\bm{\lambda}=\begin{pmatrix}1\\\ 0\\\ \vdots\\\
0\end{pmatrix},\>\>\>\>\>\>\>\bm{C}_{kl}=\int_{x\in
X_{+}}\frac{\partial^{k}P_{\theta_{0}}(x)\partial^{l}P_{\theta_{0}}(x)}{P_{\theta_{0}}(x)},$
(9)

where

$\partial^{k}P_{\theta_{0}}(x)=\frac{\partial^{k}P_{\theta}(x)}{\partial\theta^{k}}\bigg{|}_{\theta=\theta_{0}}.$

We would like to point out the similarities between the $n^{th}$order BhB
bound and the multiparameter Cramér-Rao bound (MCRB) [35] in a multiparameter
estimation problem with $n$ parameters. The computational effort to compute
both bounds is the same and the BhB can be codified as a particular case of
the CRB in a multiparameter estimation setting. The BhB can be extended to
quantum estimation theory [30] in a coherent way that includes them as a
specific case of a larger set of bounds. Doing so, the QBhB state that the
variance of any unbiased estimator with null $n$ first bias derivatives
satisfies

$(\Delta\hat{\Theta}_{\theta_{0}})^{2}\geq\max_{\bm{a}}\frac{\left(\bm{a}^{\top}\bm{\lambda}\right)^{2}}{\bm{a}^{\top}\bm{Q}\bm{a}}=\bm{\lambda}^{\top}\bm{Q}^{-1}\bm{\lambda}.$
(10)

Again, the last equality holds when $\bm{Q}^{-1}$ exists. The $n\times n$
matrix $\bm{Q}$ is defined as

$\displaystyle\bm{Q}_{kl}=\operatorname{Tr}\left(\frac{d^{k}\rho_{k}}{d\theta^{k}}L_{l}\right),$
(11)

where $L_{l}$ is the generalization of the symmetric logarithmic derivative
that satisfies

$\displaystyle\frac{d^{l}\rho_{\theta}}{d\theta^{l}}=\frac{\rho_{\theta}L_{l}+L_{l}\rho_{\theta}}{2}.$
(12)

###  II.2 An example: Bhattacharyya bound in Mach-Zehnder interferometer

Historically more attention has been paid to the Barakin bounds [26, 27]. The
Bhattacharyya bounds have not been computed in any physical problem before.
Here, we compute BhB and the estimators achieving the bound for the
paradigmatic Mach-Zehnder interferometer. We show that when the prior is not
sharp, and we move away from the point for which the unbiased estimators were
constructed, the BhB bound offers an improvement to the achievable precision.

$\varphi_{a}$$\varphi_{b}$$a$$b$$a^{\prime}$$b^{\prime}$$\ket{r}$$\ket{r}$$n_{a}$$n_{b}$
Figure 1: Mach-Zehnder interferometer setup

####  II.2.1 Mach-Zehnder

In the Mach-Zehnder conﬁguration, depicted in Fig. (1), two modes of the
electromagnetic field interfere on a balanced beamsplitter, the output beams
acquire a relative phase $\theta=\varphi_{a}-\varphi_{b}$, and interfere again
on the second beamsplitter. Finally, the photon numbers $n_{a}$ and $n_{b}$
are measured at the output ports. We consider the input state $\ket{r,r}$
where both input modes are in the Fock state with $r$ photons. It has been
shown that this state achieves the Heisenberg limit [36] while estimating the
relative phase $\theta$. The evolution induced by the interferometer is given
by the unitary $U=e^{\frac{\theta}{2}(a^{\dagger}b-b^{\dagger}a)}$ [6]. The
probability distribution $P_{\theta}(2q)$ is then

$P_{\theta}(2q)=|\langle
r-q,r+q|e^{\frac{\theta}{2}(a^{\dagger}b-b^{\dagger}a)}|r,r\rangle|^{2}$ (13)

where $2q=n_{a}-n_{b}$ is the photon counting difference. The width of this
probability distribution is $|q|<r\theta$. If $\theta\ll 1$ the probability
distribution can be well approximated by [36]

$P_{\theta}(2q)=J^{2}_{q}(r\theta),$ (14)

where $J_{q}(x)$ are the Bessel functions. We compute the CRB, the
$2^{nd}$order BhB, and the estimators achieving each bound in
$\theta_{0}=10^{-3}$. In Fig. (LABEL:PlotMSE), we can see the comparison
between the MSE performed by these two estimators in an interval of $\theta$
that contains $\theta_{0}$. One can conclude that when the prior knowledge is
not sharp, the BhB offers a more sensible bound i.e. MSEBh remains flat in an
interval of $\theta$ while MSECR increases rapidly when
$\theta\neq\theta_{0}$. In addition, the protocol that saturates the BhB has a
smaller MSE. So it might be advisable to use this bound as a figure of merit
instead of the CRB for some specific scenarios. Sll

###  II.3 Analytical motivation

In this section, we derive a lower bound for the MSE and show that it depends
on the bias. We first define $\hat{\Theta}_{\theta_{0}}(x)$ as an estimator
unbiased at $\theta_{0}$. We compute the MSE of this estimator in the point
$\theta^{\prime}$

$\displaystyle\begin{split}\text{MSE}(\theta^{\prime})&=\int\left[\theta^{\prime}-\tilde{\Theta}_{\theta_{0}}(x)\right]^{2}P_{\theta^{\prime}}(x)dx.\end{split}$
(15)

We now use the Cauchy-Schwarz inequality to find a lower bound for Eq. (15).
Doing so, one gets

$\displaystyle\left(\int f(x)g(x)dx\right)^{2}\leq\int f^{2}(x)dx\int
g^{2}(x)dx.$ (16)

Making the choice
$f(x)=(\theta^{\prime}-\tilde{\Theta}_{\theta_{0}})\sqrt{P_{\theta^{\prime}}(x)}$
and $g(x)=\sqrt{P_{\theta^{\prime}}(x)}$, we obtain

$\displaystyle\left(\int(\theta^{\prime}-\tilde{\Theta}_{\theta_{0}})P_{\theta^{\prime}}(x)dx\right)^{2}\leq\int(\theta^{\prime}-\tilde{\Theta}_{\theta_{0}})^{2}P_{\theta^{\prime}}(x)dx\int
p(x|\theta)dx$ (17)

The first integral is the square of the bias, the second is the MSE, and the
last integral is just the unity. So we can conclude that

$\displaystyle b^{2}(\theta^{\prime})\leq\text{MSE}(\theta^{\prime}).$ (18)

The squared bias gives a lower bound for the MSE. Since the derivatives of the
bias for the BhB estimator are zero, this estimator will have a smaller lower
bound for the MSE for an interval
$\theta_{0}-\epsilon<\theta^{\prime}<\theta_{0}+\epsilon$. Unfortunately, the
saturability of (18) highly depends on the particular probability distribution
and it is not possible in general. Nevertheless, in Appendix (A) we prove that
if the probability distribution satisfies certain conditions there always
exists an interval where the $n^{th}$order BhB performs at least as well as
the CRB.

##  III Conditions of the existence of estimators

Up to this point, we presented the usefulness of the BhB for experimentally
friendly scenarios. However, one cannot always construct estimators that
fulfill the conditions needed to lower bound their MSE by the BhB. This is
true especially if one deals with discrete probability distributions. In this
section, we give the conditions for the existence of the estimators for finite
support probability distributions. We conclude that in some cases, higher-
order BhB might not be computable or irrelevant. The largest BhB computable –
and that gives us new information– is related to the dimension of the support
of the probability distribution. We also show that the BhB diverges if and
only if no estimator is satisfying the unbiasedness conditions.

We present our results for both classical and quantum parameter estimation
settings. The problem of quantum parameter estimation encompasses two key
optimization challenges: (i) the optimization across all conceivable
observables or, more broadly, the optimization across all feasible POVMs that
can be executed; and (ii) the optimization across all potential estimators
$\tilde{\Theta}$ that can be derived from the outcomes of measurements. While
optimization (ii) is addressed by the classical Cramér-Rao bound, optimization
(i) introduces a distinctly quantum dimension to the problem. We will
demonstrate the implications of this extra freedom to select a POVM.

###  III.1 Classical

We consider a family of N-point support probability functions $P_{\theta}(x)$
where $P_{\theta}(x_{i})\geq 0$ for $x_{i}\in\\{x_{1},..,x_{N}\\}$. The family
of probability distributions is parametrized by the parameter $\theta$ that we
would like to infer.
The most general way to write an estimator $\tilde{\Theta}(x)$ is,

$\tilde{\Theta}(x)=\begin{cases}\tilde{\theta}_{1}&\quad\text{if }x=x_{1},\\\
\tilde{\theta}_{2}&\quad\text{if }x=x_{2},\\\ \vdots&\quad\vdots\\\
\tilde{\theta}_{N}&\quad\text{if }x=x_{N}.\\\ \end{cases}.$

We want to compute an unbiased estimator in the point $\theta_{0}$ and impose
the first $n$ derivatives of the bias to be zero in the point $\theta_{0}$.
Given that $b(\theta)=\langle\tilde{\Theta}\rangle_{\theta}-\theta$ the
conditions are

$\displaystyle\begin{split}\langle\tilde{\Theta}\rangle_{\theta_{0}}&=\theta_{0},\\\
\frac{d\langle\tilde{\Theta}\rangle_{\theta_{0}}}{d\theta}\bigg{|}_{\theta=\theta_{0}}&=1,\\\
\frac{d^{l}\langle\tilde{\Theta}\rangle_{\theta_{0}}}{d\theta^{l}}\bigg{|}_{\theta=\theta_{0}}&=0.\>\>\>\>\>\>\>\>\>l=2,3,...,n.\\\
\end{split}$ (19)

Since we are considering discrete probability distributions the mean value can
be written as finite sum
$\langle\tilde{\Theta}\rangle_{\theta_{0}}=\sum_{i}P_{\theta_{0}}(x_{i})\tilde{\Theta}(x_{i})$.
Defining $\tilde{\Theta}(x_{i})=\tilde{\theta_{i}}$, Eq. (LABEL:condiciones)
can be written as a system of equations $A\vec{x}=\vec{b}$ where

$A=\begin{pmatrix}P_{\theta_{0}}(x_{1})&P_{\theta_{0}}(x_{2})&\dots&1-\sum_{i}^{N-1}P_{\theta_{0}}(x_{i})\\\
\partial P_{\theta_{0}}(x_{1})&\partial
P_{\theta_{0}}(x_{2})&\dots&-\sum_{i}^{N-1}\partial P_{\theta_{0}}(x_{i})\\\
\partial^{2}P_{\theta_{0}}(x_{1})&\partial^{2}P_{\theta_{0}}(x_{2})&\dots&-\sum_{i}^{N-1}\partial^{2}P_{\theta_{0}}(x_{i})\\\
\vdots&\vdots&\ddots&\vdots\\\
\partial^{n}P_{\theta_{0}}(x_{1})&\partial^{n}P_{\theta_{0}}(x_{2})&\dots&-\sum_{i}^{N-1}\partial^{n}P_{\theta_{0}}(x_{i})\\\
\\\ \end{pmatrix},\quad\quad\vec{x}=\begin{pmatrix}\tilde{\theta}_{1}\\\
\tilde{\theta}_{2}\\\ \vdots\\\
\tilde{\theta}_{N}\end{pmatrix},\quad\quad\vec{b}=\begin{pmatrix}\theta_{0}\\\
1\\\ 0\\\ \vdots\\\ 0\end{pmatrix},$ (20)

where we defined
$\partial^{l}P_{\theta_{0}}(x)\equiv\frac{\partial^{l}P_{\theta}(x)}{\partial\theta^{l}}|_{\theta=\theta_{0}}$.
The system has $n+1$ equations and $N$ variables. We consider the case where
the first $N-1$ derivatives are linearly independent. This implies that if
$n\leq N-1$, there is always a solution. If $n>N-1$, then the higher
derivatives are no longer independent, and they can be written as a
combination of the first $N-1$ linearly independent derivatives. In this case,
the system has no solution or the solution is trivial, that is to say, the
$n^{th}$order BhB is equal to $N-1$ order BhB for all $n>N-1$. It is then
concluded that given a $N$ points probability distribution is advisable to
compute BhB bound up to order $N-1$.

In the appendix (B), the conditions for the existence of the estimator are
given. It is also shown that the bound diverges if and only if the estimator
does not exist.

###  III.2 Quantum

We extend the previous results to determine when computing QBhB gives more
information about the estimation problem. In quantum estimation theory the
measurement is a POVM $\\{\Pi_{i}\\}$, where $i$ are the possible outcomes,
each one with probability
$P_{\theta}(i)=\operatorname{Tr}(\Pi_{i}\rho_{\theta})$. The classical
processing of the measurement outcomes is accounted by the estimator
$\tilde{\Theta}(x)$ that associates an estimate $\tilde{\Theta}(i)$ to the
outcome $i$. Every estimate $\tilde{\Theta}(i)$ has a probability given by
$P_{\theta}(i)$. Both tasks, measurement, and classical processing, can be
treated simultaneously as a hermitian operator
$\hat{\Theta}=\sum_{i}\tilde{\Theta}(i)\Pi_{i}$ [22]. The conditions of Eq.
(LABEL:condiciones) are now written in terms of the operator $\hat{\Theta}$ as

$\displaystyle\begin{pmatrix}\operatorname{Tr}[\rho_{\theta_{0}}\hat{\Theta}]\\\
\operatorname{Tr}[d\rho_{\theta_{0}}\hat{\Theta}]\\\
\operatorname{Tr}[d^{2}\rho_{\theta_{0}}\hat{\Theta}]\\\ \vdots\\\
\operatorname{Tr}[d^{n}\rho_{\theta_{0}}\hat{\Theta}]\\\
\end{pmatrix}=\begin{pmatrix}\theta_{0}\\\ 1\\\ 0\\\ \vdots\\\ 0\\\
\end{pmatrix},$ (21)

where again

$d^{k}\rho_{\theta_{0}}=\frac{d^{k}\rho_{\theta}}{d\theta^{k}}\bigg{|}_{\theta=\theta_{0}}.$

We rewrite these conditions as

$\displaystyle\begin{pmatrix}(\rho_{\theta_{0}})_{ij}\hat{\Theta}_{ji}\\\
(d\rho_{\theta_{0}})_{ij}\hat{\Theta}_{ji}\\\
(d^{2}\rho_{\theta_{0}})_{ij}\hat{\Theta}_{ji}\\\ \vdots\\\
(d^{n}\rho_{\theta_{0}})_{ij}\hat{\Theta}_{ji}\\\
\end{pmatrix}=\begin{pmatrix}\theta_{0}\\\ 1\\\ 0\\\ \vdots\\\ 0\\\
\end{pmatrix},$ (22)

where we are using the Einstein summation convention and $i,j=1,2,...,N$. Now
we have a system of equations with complex coefficients for the complex
variables $\hat{\Theta}_{ij}$. Since the estimator operator must be hermitian,
not all the variables $\hat{\Theta}_{ij}$ are independent. Instead of $N^{2}$
independent variables, there are now $\frac{N(N+1)}{2}$ independent variables.
In that way, we can vectorize the system of equations as follows

$\begin{pmatrix}(\rho_{\theta_{0}})_{1}&(\rho_{\theta_{0}})_{2}&\dots&(\rho_{\theta_{0}})_{\frac{N(N+1)}{2}}\\\
(d\rho_{\theta_{0}})_{1}&(d\rho_{\theta_{0}})_{2}&\dots&(d\rho_{\theta_{0}})_{\frac{N(N+1)}{2}}\\\
(d^{2}\rho_{\theta})_{1}&(d^{2}\rho_{\theta})_{2}&\dots&(d^{2}\rho_{\theta})_{\frac{N(N+1)}{2}}\\\
\vdots&\vdots&\ddots&\vdots\\\
(d^{n}\rho_{\theta})_{1}&(d^{n}\rho_{\theta})_{2}&\dots&(d^{n}\rho_{\theta})_{\frac{N(N+1)}{2}}\\\
\end{pmatrix}\begin{pmatrix}\hat{\Theta}_{1}\\\ \hat{\Theta}_{2}\\\ \vdots\\\
\hat{\Theta}_{\frac{N(N+1))}{2}}\end{pmatrix}=\begin{pmatrix}\theta_{0}\\\
1\\\ 0\\\ \vdots\\\ 0\end{pmatrix}.$ (23)

So we have a system of equations equivalent to the one for a classical
probability distribution with $\frac{N(N+1)}{2}$ support points. As before if
this system of equations cannot be solved the QBhB variance will go to
infinity. On the other hand, if the first $\frac{N(N+1)}{2}-1$ derivatives are
linearly independent, for all $n>\frac{N(N+1)}{2}-1$ one has that
$n^{th}$order QBhB is equal to $\frac{N(N+1)}{2}-1$ order QBhB. The difference
between the quantum and classical cases is that in the former we have the
freedom to choose a measurement. This extra freedom implies that the highest
relevant QBhB bound –which gives a different bound than lower orders– is
lifted in the quantum case because $N-1<\frac{N(N+1)}{2}-1$.

Let us elaborate more on this because at first glance the result may seem to
be at odds with the classical case. After projective measurements, one deals
with $N$ point probability distributions. Therefore, the results derived in
the classical setting must also apply, so it might seem that there will be a
contradiction to what we stated here. To solve this apparent contradiction, we
explain why $N$ order QBh bound need not be equal to the $N-1$ one. The key
point is that to compute the $N$ and the $N+1$ order QBhB, different
measurements have to be performed, giving different probability distributions
$P_{N}$ and $P_{N+1}$. This extra freedom in the measurements is what gives
the quantum systems more non-trivial information, or in other words, the
maximal order for which QBhB gives different information is higher than for
classical scenarios. Thus, since $P_{N}$ and $P_{N+1}$ are different
probability distributions, there is no discrepancy with what we presented in
section (III.1).

Figure 3: QCR and QBhB as a function of $\lambda$. One sees that they differ
and as expected QBhB bound will be higher than QCR. In addition, we plot these
values of $\lambda$ because the plot is symmetric around $\lambda=0.5$ and
also divergent for this same value –since the state will be
$\rho\propto\mathbb{I}$ and won’t contain any information on the parameter–.

###  III.3 Example of the quantum case

In this section, we show a case in which we see the lifting in the constraints
produced by the quantum nature of the problem, meaning that we see that the
quantum bounds are not trivial –giving the same value as bounds from lower
orders– for $n>N-1$, being $n$ the order of the BhB employed and $N$ the
number of outcomes.

Consider that our state $\rho(\theta)$ is obtained through the following
unitary operation $U=e^{-i\theta^{2}H}$ where $H=\sigma_{x}$ is the $x$-Pauli
matrix, and we would like to estimate $\theta$. The number of outcomes here is
$N=2$. Our initial state is given by

$\rho(0)=\frac{1}{2}\mathbb{I}+\frac{(2\lambda-1)}{2}\sigma_{z},$ (24)

where $\lambda\in(0,1)$ and

$\rho(\theta)=U\rho(0)U^{\dagger}=\left(\begin{array}[]{cc}\left(\lambda-\frac{1}{2}\right)\cos\left(2\theta^{2}\right)+\frac{1}{2}&\frac{1}{2}i(2\lambda-1)\sin\left(2\theta^{2}\right)\\\
i(1-2\lambda)\sin\left(\theta^{2}\right)\cos\left(\theta^{2}\right)&\frac{1}{2}\left((1-2\lambda)\cos\left(2\theta^{2}\right)+1\right)\\\
\end{array}\right).$ (25)

We can proceed and calculate $L_{1}$, and $L_{2}$ following the definitions
presented in previous sections. Doing some algebra one ends up with the QBhB
matrix which reads

$\bm{Q}_{BhB^{2}}=\left(\begin{array}[]{cc}16\theta^{2}(1-2\lambda)^{2}&16\theta(1-2\lambda)^{2}\\\
16\theta(1-2\lambda)^{2}&\frac{16(1-2\lambda)^{2}\left((\lambda-1)\lambda-4\theta^{4}\right)}{(\lambda-1)\lambda}\\\
\end{array}\right),$ (26)

where the first element corresponds to the QFI. To finally get the BhB, we
take the inverse of the matrix and look at the first element of it. In Fig.
(3) we can see the lower bound on the variance of the estimator when using the
QCR bound and QBhB. Thus, computing the quantum Bhattacharyya bound gives new
non-trivial information even though $n>N-1$.

##  IV Conclusions

Since long ago, it has been clear in estimation theory that while the CRB is
extremely useful and illustrative is not always meaningful. For that reason,
different bounds were constructed and introduced in the classical regime. On
the other hand, until very recently, there were no works extending these
bounds to the quantum regime [30].

In this article, we first motivated the use of Bhattacharyya bounds using a
simple example. Essentially, when one does not know the prior value of the
parameter exactly, but knows that is within a small range, the BhB are more
desirable. We then gave necessary and sufficient conditions for an estimator
that satisfies the BhB to exist. The existence conditions are related to the
support of the probability distribution which determines the freedom we have
to construct an estimator. When we computed the $n$-th order BhB bound we
found a limit for the variance of an estimator that satisfies $n+1$
unbiasedness conditions. When the support of the probability distribution is
small, we might not have enough freedom to construct an estimator satisfying
the desired constraints. Namely, the estimator exists if and only if the
maximization presented in Eq. (8) can be calculated (non-divergent). Moreover,
following a similar reasoning, we could anticipate when computing higher
orders in the BhB gives no extra information. Specifically, we saw that for
$N$ outcomes, it is “sufficient” to calculate the $n^{th}$order BhB bound
where $n\leq N-1$. However, in the quantum case, where we have more freedom in
the measurements, the order is lifted to $n\leq N(N+1)/2-1$ implying that in
the quantum regime, the estimators are more resilient. Finally, it would be of
theoretical and practical interest to extend the results presented here to
continuous probability distributions as well as to biased estimators. For the
latter case, we think it will be conceptually difficult to compare the
different bounds such as the BhB and CRB.

###### Acknowledgements.

We acknowledge fruitful discussions with Manuel Gessner and Maximilian
Reichert. The authors acknowledge financial support from OpenSuperQ+100 (Grant
No. 101113946) of the EU Flagship on Quantum Technologies, as well as from the
EU FET-Open project EPIQUS (Grant No. 899368), also from Project Grant No.
PID2021-125823NA-I00 595 and Spanish Ramón y Cajal Grant No. RYC-2020-030503-I
funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”
and “ERDF Invest in your Future,” this project has also received support from
the Spanish Ministry of Economic Affairs and Digital Transformation through
the QUANTUM ENIA project call –Quantum Spain, and by the EU through the
Recovery, Transformation and Resilience Plan– NextGenerationEU within the
framework of the Digital Spain 2026 Agenda, we acknowledge funding from Basque
Government through Grant No. IT1470-22 and the IKUR Strategy under the
collaboration agreement between Ikerbasque Foundation and BCAM on behalf of
the Department of Education of the Basque Government.

## References

  * [1] Matteo GA Paris.  Quantum estimation for quantum technology.  International Journal of Quantum Information, 7(supp01):125–137, 2009.
  * [2] Géza Tóth and Iagoba Apellaniz.  Quantum metrology from a quantum information science perspective.  Journal of Physics A: Mathematical and Theoretical, 47(42):424006, 2014.
  * [3] Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone.  Advances in quantum metrology.  Nature photonics, 5(4):222–229, 2011.
  * [4] A gravitational wave observatory operating beyond the quantum shot-noise limit.  Nature Physics, 7(12):962–965, 2011.
  * [5] Junaid Aasi, Joan Abadie, BP Abbott, Richard Abbott, TD Abbott, MR Abernathy, Carl Adams, Thomas Adams, Paolo Addesso, RX Adhikari, et al.  Enhanced sensitivity of the ligo gravitational wave detector by using squeezed states of light.  Nature Photonics, 7(8):613–619, 2013.
  * [6] Rafal Demkowicz-Dobrzański, Marcin Jarzyna, and Jan Kołodyński.  Quantum limits in optical interferometry.  Progress in Optics, 60:345–435, 2015.
  * [7] Manuel Gessner, Luca Pezzè, and Augusto Smerzi.  Sensitivity bounds for multiparameter quantum metrology.  Physical review letters, 121(13):130503, 2018.
  * [8] Quntao Zhuang and Jeffrey H Shapiro.  Ultimate accuracy limit of quantum pulse-compression ranging.  Physical review letters, 128(1):010501, 2022.
  * [9] Maximilian Reichert, Roberto Di Candia, Moe Z Win, and Mikel Sanz.  Quantum-enhanced doppler lidar.  npj Quantum Information, 8(1):147, 2022.
  * [10] Maximilian Reichert, Quntao Zhuang, and Mikel Sanz.  Heisenberg-limited quantum lidar for joint range and velocity estimation.  arXiv preprint arXiv:2311.14546, 2023.
  * [11] Si-Hui Tan, Baris I Erkmen, Vittorio Giovannetti, Saikat Guha, Seth Lloyd, Lorenzo Maccone, Stefano Pirandola, and Jeffrey H Shapiro.  Quantum illumination with gaussian states.  Physical review letters, 101(25):253601, 2008.
  * [12] Mikel Sanz, U Las Heras, Juan José García-Ripoll, Enrique Solano, and R Di Candia.  Quantum estimation methods for quantum illumination.  Physical review letters, 118(7):070803, 2017.
  * [13] Maximilian Reichert, Quntao Zhuang, Jeffrey H. Shapiro, and Roberto Di Candia.  Quantum illumination with a hetero-homodyne receiver and sequential detection.  Phys. Rev. Appl., 20:014030, Jul 2023.
  * [14] Mohammad Mehboudi, Anna Sanpera, and Luis A Correa.  Thermometry in the quantum regime: recent theoretical progress.  Journal of Physics A: Mathematical and Theoretical, 52(30):303001, 2019.
  * [15] Ricard Ravell Rodríguez, Mohammad Mehboudi, Michał Horodecki, and Martí Perarnau-Llobet.  Strongly coupled fermionic probe for nonequilibrium thermometry.  New Journal of Physics, 26(1):013046, 2024.
  * [16] Mohammad Mehboudi, Aniello Lampo, Christos Charalambous, Luis A Correa, Miguel Ángel García-March, and Maciej Lewenstein.  Using polarons for sub-nk quantum nondemolition thermometry in a bose-einstein condensate.  Physical review letters, 122(3):030403, 2019.
  * [17] Felix Binder, Luis A Correa, Christian Gogolin, Janet Anders, and Gerardo Adesso.  Thermodynamics in the quantum regime.  Fundamental Theories of Physics, 195:1–2, 2018.
  * [18] L Davidovich and RL de Matos Filho.  Quantum illumination strategy for parameter estimation.  Physical Review A, 108(4):042612, 2023.
  * [19] Harald Cramer.  Mathematical methods of statistics.  Princeton mathematical series. Princeton Univ. Press, Princeton, NJ, 1954.
  * [20] C. Radhakrishna Rao.  Information and the Accuracy Attainable in the Estimation of Statistical Parameters, pages 235–247.  Springer New York, New York, NY, 1992.
  * [21] Jesús Rubio, Paul Knott, and Jacob Dunningham.  Non-asymptotic analysis of quantum metrology protocols beyond the cramér–rao bound.  Journal of Physics Communications, 2(1):015027, 2018.
  * [22] Johannes Jakob Meyer, Sumeet Khatri, Daniel Stilck França, Jens Eisert, and Philippe Faist.  Quantum metrology in the finite-sample regime.  arXiv preprint arXiv:2307.06370, 2023.
  * [23] Marcin Jarzyna and Rafał Demkowicz-Dobrzański.  True precision limits in quantum metrology.  New Journal of Physics, 17(1):013010, 2015.
  * [24] Douglas G Chapman and Herbert Robbins.  Minimum variance estimation without regularity assumptions.  The Annals of Mathematical Statistics, pages 581–586, 1951.
  * [25] A Bhattacharyya.  On some analogues of the amount of information and their use in statistical estimation.  Sankhyā: The Indian Journal of Statistics, pages 1–14, 1946.
  * [26] R McAulay and E Hofstetter.  Barankin bounds on parameter estimation.  IEEE Transactions on Information Theory, 17(6):669–676, 1971.
  * [27] Luc Knockaert.  The barankin bound and threshold behavior in frequency estimation.  IEEE Transactions on Signal Processing, 45(9):2398–2401, 1997.
  * [28] Jonathan S Abel.  A bound on mean-square-estimate error.  IEEE Transactions on Information Theory, 39(5):1675–1680, 1993.
  * [29] Athanasios Papoulis and S Unnikrishna Pillai.  Probability, random variables and stochastic processes. Sec 8.2.  2002\.
  * [30] Manuel Gessner and Augusto Smerzi.  Hierarchies of frequentist bounds for quantum metrology: From cramér-rao to barankin.  Phys. Rev. Lett., 130:260801, Jun 2023.
  * [31] Samuel L Braunstein and Carlton M Caves.  Statistical distance and the geometry of quantum states.  Physical Review Letters, 72(22):3439, 1994.
  * [32] Rafał Demkowicz-Dobrzański, Wojciech Górecki, and Mădălin Guţă.  Multi-parameter estimation beyond quantum fisher information.  Journal of Physics A: Mathematical and Theoretical, 53(36):363001, 2020.
  * [33] Wojciech Górecki, Rafał Demkowicz-Dobrzański, Howard M Wiseman, and Dominic W Berry.  $\pi$-corrected heisenberg limit.  Physical review letters, 124(3):030501, 2020.
  * [34] BJN Blight and PV Rao.  The convergence of bhattacharyya bounds.  Biometrika, 61(1):137–142, 1974.
  * [35] Jing Liu, Haidong Yuan, Xiao-Ming Lu, and Xiaoguang Wang.  Quantum fisher information matrix and multiparameter estimation.  Journal of Physics A: Mathematical and Theoretical, 53(2):023001, 2020.
  * [36] MJ Holland and K Burnett.  Interferometric detection of optical phase shifts at the heisenberg limit.  Physical review letters, 71(9):1355, 1993.
  * [37] Tom M Apostol.  Calculus, Volume 1, Sec. 7.6.  John Wiley & Sons, 1991.

##  Appendix A Analytical comparison MSE${}_{BhB^{n}}$ and MSECR

Given a family of probability distribution $P_{\theta}(x)$ parameterized by a
continuous parameter $\theta\in\Theta$ we show that if
$\exists\>\>|\partial^{n^{\prime}}P_{\theta}(x)|<M\>\>\forall\>(n^{\prime},\theta,x)$
where $M$ is a real constant and satisfies,

$\displaystyle\int
P_{\theta_{0}}(x)\bigg{|}\frac{\partial^{i}P_{\theta_{0}}(x)}{P_{\theta_{0}}(x)}\bigg{|}^{k}\>dx<\infty\quad\forall\>i,k$
(27)

then there is an interval where using the estimator saturating $n^{th}$order
BhB performs at least as well as the CR. Condition (27) ensures the existence
of the elements of the Bhattacharyya matrix and is satisfied by the most
common families of probability distributions, for instance, the Beta, Cauchy,
Exponential, Gamma, and Normal or Gaussian among others. Lets compute the MSE,

$\displaystyle\text{MSE}_{BhB^{n}}(\theta)=\int(\theta-\tilde{\Theta}(x))^{2}P_{\theta}(x)dx$
(28)

The optimal form for the estimator saturating the $n^{th}$order BhB bound in
$\theta_{0}$ is
$\tilde{\Theta}_{\theta_{0}}(x)=\theta_{0}+\bm{g}^{\top}(x)\>\bm{C}^{-1}\>\bm{\lambda}$
where
$\bm{g}^{\top}(x)=(\partial^{1}P_{\theta_{0}}(x),\partial^{2}P_{\theta_{0}}(x),\dots,\partial^{n}P_{\theta_{0}}(x))$.
Eq. (28) can be written,

$\displaystyle\text{MSE}_{BhB^{n}}(\theta)=\int(\theta-\theta_{0}-\frac{1}{P_{\theta_{0}}(x)}\>\>\bm{g}^{\top}(x)\>\bm{C}^{-1}\>\bm{\lambda})^{2}P_{\theta}(x)dx$
(29)
$\displaystyle=(\theta-\theta_{0})^{2}-2(\theta-\theta_{0})\int\frac{1}{P_{\theta_{0}}(x)}\>\>(\bm{g}^{\top}(x)\>\bm{C}^{-1}\>\bm{\lambda})P_{\theta}(x)dx+\int\left(\frac{1}{P_{\theta_{0}}(x)}\>\>\bm{g}^{\top}(x)\>\bm{C}^{-1}\>\bm{\lambda}\right)^{2}P_{\theta}(x)dx$
(30)

The probability distribution is Taylor expanded up to order $n^{\prime}$,
without losing generality we consider $n^{\prime}\geq n$,

$\displaystyle
P_{\theta}(x)=P_{\theta_{0}}(x)+\sum_{j=1}^{n^{\prime}}\frac{\partial^{j}P_{\theta}(x)}{j!}(\theta-\theta_{0})^{j}+R^{n^{\prime}}_{\theta}(x)$
(31)

where $R^{n^{\prime}}_{\theta}(x)$ is the remainder up to order $n^{\prime}$.
Substituting Eq. (31) and simplifying Eq. (30) can be written as

$\displaystyle\delta^{2}-2\delta\sum_{j=1}^{n^{\prime}}\frac{\delta^{j}}{j!}\int\frac{\bm{g}^{\top}(x)\partial^{j}P_{\theta_{0}}(x)}{P_{\theta_{0}}(x)}dx\>\>\bm{C}^{-1}\>\bm{\lambda}+\sum_{j=0}^{n^{\prime}}\frac{\delta^{j}}{j!}\int\frac{\partial^{j}P_{\theta_{0}}(x)}{P_{\theta_{0}}^{2}(x)}(\bm{g}^{\top}(x)\>\bm{C}^{-1}\>\bm{\lambda})^{2}dx+R_{BhB^{n}}^{\prime}$
(32)

where $\delta=\theta-\theta_{0}$, $\partial_{0}P_{\theta_{0}}=P_{\theta_{0}}$
and $R_{BhB^{n}}^{\prime}$ is,

$\displaystyle R_{BhB^{n}}^{\prime}$
$\displaystyle=\int\left(\delta-\frac{1}{P_{\theta_{0}}(x)}\>\>\bm{g}^{\top}(x)\>\bm{C}^{-1}\>\bm{\lambda}\right)^{2}R^{n^{\prime}}_{\theta}(x)dx$
(33)
$\displaystyle=-2\delta\int\frac{1}{P_{\theta_{0}}(x)}\>\>(\bm{g}^{\top}(x)\>\bm{C}^{-1}\>\bm{\lambda})R^{n^{\prime}}_{\theta}(x)dx+\int\left(\frac{1}{P_{\theta_{0}}(x)}\>\>\bm{g}^{\top}(x)\>\bm{C}^{-1}\>\bm{\lambda}\right)^{2}R^{n^{\prime}}_{\theta}(x)dx.$
(34)

Note that the first integral in Eq. (32) is the definition of the BhB matrix
(9) so it can be simplified,

$\displaystyle\delta^{2}-2\delta\sum_{j=1}^{n^{\prime}}\frac{\delta^{j}}{j!}\bm{C}_{j,i}\>\>\bm{C}^{-1}_{i,1}+\sum_{j=0}^{n^{\prime}}\frac{\delta^{j}}{j!}\bm{S}_{i,j,l}\bm{C}^{-1}_{i,1}\bm{C}^{-1}_{l,1}+R_{BhB^{n}}^{\prime}$
(35)
$\displaystyle=\delta^{2}-2\delta^{2}-2\delta\sum_{j=n}^{n^{\prime}}\frac{\delta^{j}}{j!}\bm{C}_{j,i}\>\>\bm{C}^{-1}_{i,1}+\sum_{j=0}^{n^{\prime}}\frac{\delta^{j}}{j!}\bm{S}_{i,j,l}\bm{C}^{-1}_{i,1}\bm{C}^{-1}_{l,1}+R_{BhB^{n}}^{\prime}$
(36)

where we are using the Einstein convention –there are implicit sums over $i,l$
indices– and have defined the tensor $\bm{S}$

$\displaystyle\bm{S}_{i,j,l}=\int\frac{\partial^{i}P_{\theta_{0}}(x)\>\partial^{k}P_{\theta_{0}}(x)\>\partial^{l}P_{\theta_{0}}(x)}{P^{2}_{\theta_{0}}(x)}dx.$
(37)

Note that the elements $\bm{S}_{i,j,l}$ exists because (27). Analogously for
the estimator saturating CRB we have

$\displaystyle\text{MSE}_{CR}$
$\displaystyle=\delta^{2}-2\delta\sum_{j=1}^{n^{\prime}}\frac{\delta^{j}}{j!}\frac{1}{F_{C}}\int\frac{\partial^{1}P_{\theta_{0}}(x)\partial^{j}P_{\theta_{0}}(x)}{P_{\theta_{0}}(x)}dx+\sum_{j=0}^{n^{\prime}}\frac{\delta^{j}}{j!}\frac{1}{F_{C}^{2}}\int\frac{\partial^{1}P_{\theta_{0}}\partial^{1}P_{\theta_{0}}\partial^{j}P_{\theta_{0}}}{P^{2}_{\theta_{0}}}dx+R_{CR}^{\prime}$
(38)
$\displaystyle=\delta^{2}-2\delta\sum_{j=1}^{n^{\prime}}\frac{\delta^{j}}{j!}\frac{1}{F_{C}}\bm{C}_{j,1}+\sum_{j=0}^{n^{\prime}}\frac{\delta^{j}}{j!}\frac{1}{F_{C}^{2}}\bm{S}_{1,1,j}+R_{CR}^{\prime},$
(39)

where we used that $\bm{C}_{1,1}$ is the classical Fisher information $F_{C}$.
The remainder term can be estimated [37],

$\displaystyle q(x)\>\frac{\delta^{{n^{\prime}}+1}}{({n^{\prime}}+1)!}\leq
R^{n^{\prime}}_{\theta}(x)\leq
Q(x)\>\frac{\delta^{{n^{\prime}}+1}}{({n^{\prime}}+1)!},$ (40)

where

$\displaystyle q(x)\leq\partial^{{n^{\prime}}+1}P_{\theta}(x)\leq
Q(x),\quad\quad\forall\theta.$ (41)

We use the functions $Q(x)$ and $q(x)$ to bound the $n^{th}$order BhB error
$R^{\prime}_{BhB^{n}}$ and the CR error $R^{\prime}_{CR}$ respectively

$\displaystyle R^{\prime}_{BhB^{n}}\leq\int
Q(x)\>\frac{\delta^{{n^{\prime}}+1}}{({n^{\prime}}+1)!}\left(\delta-\frac{1}{P_{\theta_{0}}(x)}\>\>\bm{g}(x)\>\bm{C}^{-1}\>\bm{\lambda}\right)^{2}dx=F(\delta),$
(42) $\displaystyle R^{\prime}_{CR}\geq\int
q(x)\>\frac{\delta^{{n^{\prime}}+1}}{({n^{\prime}}+1)!}\left(\delta-\frac{1}{F_{C}}\frac{\partial
P_{\theta_{0}}(x)}{P_{\theta_{0}}(x)}\right)^{2}dx=G(\delta),$ (43)

where both functions are polynomials of order $n^{\prime}+3$ in the variable
$\delta$, i.e. $F,G\in P_{n^{\prime}+3}(\delta)$. To compare the performance
of the CRB-inspired protocol and the BhB one we integrate the MSE difference
along an interval around $\theta_{0}$ i.e.
$\theta_{0}-\frac{\Delta}{2}\leq\theta\leq\theta_{0}+\frac{\Delta}{2}$,

$\displaystyle\int_{\theta_{0}-\frac{\Delta}{2}}^{\theta_{0}+\frac{\Delta}{2}}\left(\text{MSE}_{CR}(\theta)-\text{MSE}_{BhB^{n}}(\theta)\right)d\theta$
(44)
$\displaystyle\geq\int_{-\frac{\Delta}{2}}^{\frac{\Delta}{2}}\left(\delta^{2}-2\delta\sum_{j=1}^{n}\frac{\delta^{j}}{j!}\frac{1}{F_{C}}\bm{C}_{j,1}+\sum_{j=0}^{n}\frac{\delta^{j}}{j!}\frac{1}{F_{C}^{2}}\bm{S}_{1,1,j}+G(\delta)+\delta^{2}-\sum_{j=0}^{n}\frac{\delta^{j}}{j!}\bm{S}_{i,j,l}\bm{C}^{-1}_{i,1}\bm{C}^{-1}_{l,1}-F(\delta)\right)d\delta$
(45)
$\displaystyle=\int_{-\frac{\Delta}{2}}^{\frac{\Delta}{2}}\left(-2\delta\sum_{j=2}^{n}\frac{\delta^{j}}{j!}\frac{1}{F_{C}}\bm{C}_{j,1}+\sum_{j=0}^{n}\frac{\delta^{j}}{j!}\frac{1}{F_{C}^{2}}\bm{S}_{1,1,j}+G(\delta)-\sum_{j=0}^{n}\frac{\delta^{j}}{j!}\bm{S}_{i,j,l}\bm{C}^{-1}_{i,1}\bm{C}^{-1}_{l,1}-F(\delta)\right)d\delta=H(\Delta).$
(46)

The rhs of Eq. (45) is bounded by $H(\Delta)$ which is a polynomial of order
$n^{\prime}+4$ in $\Delta$. Choosing $n^{\prime}$ so that $n^{\prime}+4$ is
odd, $H(\Delta)$ has at least one real root, showing that there is an interval
$-\frac{\Delta}{2}\leq\delta\leq\frac{\Delta}{2}$ where employing the BhB
estimator is as good as employing the CR estimator.

##  Appendix B Existence of the estimator

We study the solution of the system $A\vec{x}=\vec{b}$ with

$A=\begin{pmatrix}P_{\theta_{0}}(x_{1})&P_{\theta_{0}}(x_{2})&\dots&1-\sum_{i}^{N-1}P_{\theta_{0}}(x_{i})\\\
\partial P_{\theta_{0}}(x_{1})&\partial
P_{\theta_{0}}(x_{2})&\dots&-\sum_{i}^{N-1}\partial P_{\theta_{0}}(x_{i})\\\
\partial^{2}P_{\theta_{0}}(x_{1})&\partial^{2}P_{\theta_{0}}(x_{2})&\dots&-\sum_{i}^{N-1}\partial^{2}P_{\theta_{0}}(x_{i})\\\
\vdots&\vdots&\ddots&\vdots\\\
\partial^{n}P_{\theta_{0}}(x_{1})&\partial^{n}P_{\theta_{0}}(x_{2})&\dots&-\sum_{i}^{N-1}\partial^{n}P_{\theta_{0}}(x_{i})\\\
\end{pmatrix},\quad\quad\vec{x}=\begin{pmatrix}\tilde{\theta}_{1}\\\
\tilde{\theta}_{2}\\\ \vdots\\\
\tilde{\theta}_{N}\end{pmatrix},\quad\quad\vec{b}=\begin{pmatrix}\theta_{0}\\\
1\\\ 0\\\ \vdots\\\ 0\end{pmatrix}.$ (47)

We point out that the matrix has dimensions $(n+1)\times N$ so if $n+1>N$
there are more equations than variables and a solution might not exist. Also,
it is important to stress that the derivatives inside the matrix are all
evaluated at point $\theta_{0}$ but we do not write it for the sake of
clarity.

Rouché-Frobenius theorem states that the system has a solution if and only if
$\operatorname{rank}(A|\vec{b})=\operatorname{rank}(A)$. If all the rows of
$A$ are linearly independent matrix $(A|\vec{b})$ wouldn’t have a larger rank
than $A$. Then, we are interested in the case where $A$ has linearly dependent
rows. Notice that the first row is linearly independent of the others so the
system $A\vec{x}=\vec{b}$ has no solution if and only if
$A^{\prime}\vec{x}=\vec{\lambda}$ has no solution. From now on we study the
latter system

$A^{\prime}=\begin{pmatrix}\partial P_{\theta_{0}}(x_{1})&\partial
P_{\theta_{0}}(x_{2})&\dots&-\sum_{i}^{N-1}\partial P_{\theta_{0}}(x_{i})\\\
\partial^{2}P_{\theta_{0}}(x_{1})&\partial^{2}P_{\theta_{0}}(x_{2})&\dots&-\sum_{i}^{N-1}\partial^{2}P_{\theta_{0}}(x_{i})\\\
\vdots&\vdots&\ddots&\vdots\\\
\partial^{n}P_{\theta_{0}}(x_{1})&\partial^{n}P_{\theta_{0}}(x_{2})&\dots&-\sum_{i}^{N-1}\partial^{n}P_{\theta_{0}}(x_{i})\\\
\end{pmatrix},\quad\quad\vec{x}=\begin{pmatrix}\tilde{\theta}_{1}\\\
\tilde{\theta}_{2}\\\ \vdots\\\
\tilde{\theta}_{N}\end{pmatrix},\quad\quad\vec{\lambda}=\begin{pmatrix}1\\\
0\\\ \vdots\\\ 0\end{pmatrix}.$ (48)

The matrix $A^{\prime}$ has dimensions $n\times N$ and since the last column
is linearly dependent (is the sum of the $N-1$ first columns) there are at
most $N-1$ linearly independent rows. As we mentioned before linearly
dependent rows are a necessary condition for the system not having a solution.
For that reason, we suppose that there is a row $m$ that can be written as,

$\begin{pmatrix}\partial^{m}P_{\theta_{0}}(x_{1})&\partial^{m}P_{\theta_{0}}(x_{2})&\dots&-\sum_{i}^{N-1}\partial^{m}P_{\theta_{0}}(x_{i})\\\
\end{pmatrix}=\sum_{l}\alpha_{l}\begin{pmatrix}\partial^{l}P_{\theta_{0}}(x_{1})&\partial^{l}P_{\theta_{0}}(x_{2})&\dots&-\sum_{i}^{N-1}\partial^{l}P_{\theta_{0}}(x_{i})\\\
\end{pmatrix}.$ (49)

Without losing generality we can consider that $m\neq 1$. Now we consider the
$m$ row of the augmented matrix $(A^{\prime}|\vec{\lambda})$,

$\begin{pmatrix}\partial^{m}P_{\theta_{0}}(x_{1})&\partial^{m}P_{\theta_{0}}(x_{2})&\dots&-\sum_{i}^{N-1}\partial^{m}P_{\theta_{0}}(x_{i})&0\\\
\end{pmatrix}.$ (50)

If the $m$ row of the augmented matrix Eq. (50) is linearly independent from
the rows of $(A^{\prime}|\vec{\lambda})$ then
$\operatorname{rank}(A^{\prime}|\vec{\lambda})>\operatorname{rank}(A^{\prime})$
and there is no solution.

Thus we investigate whether it is linearly dependent or not. According to our
assumption, we can write Eq. (50) as,

$\begin{pmatrix}\partial^{m}P_{\theta_{0}}(x_{1})&\dots&-\sum_{i}^{N-1}\partial^{m}P_{\theta_{0}}(x_{i})&0\\\
\end{pmatrix}=\sum_{l}\alpha_{l}\begin{pmatrix}\partial^{l}P_{\theta_{0}}(x_{1})&\dots&-\sum_{i}^{N-1}\partial^{l}P_{\theta_{0}}(x_{i})&0\\\
\end{pmatrix}$ (51)

but the first element of the sum,

$\alpha_{1}\begin{pmatrix}\partial
P_{\theta_{0}}(x_{1})&\dots&-\sum_{i}^{N-1}\partial P_{\theta_{0}}(x_{i})&0\\\
\end{pmatrix}$ (52)

is not a row of the matrix $(A^{\prime}|\vec{\lambda})$ so Eq. (50) is
linearly independent when $\alpha_{1}\neq 0$. Then, we conclude that the
system given in Eq. (47) has no solution if and only if,

$\begin{pmatrix}\partial^{m}P_{\theta_{0}}(x_{1})&\partial^{m}P_{\theta_{0}}(x_{2})&\dots&-\sum_{i}^{N-1}\partial^{m}P_{\theta_{0}}(x_{i})\\\
\end{pmatrix}=\sum_{l}\alpha_{l}\begin{pmatrix}\partial^{l}P_{\theta_{0}}(x_{1})&\partial^{l}P_{\theta_{0}}(x_{2})&\dots&-\sum_{i}^{N-1}\partial^{l}P_{\theta_{0}}(x_{i})\\\
\end{pmatrix},$ (53)

where $\alpha_{1}\neq 0$.

Now we look at the Bhattacharyya bound

$(\Delta\hat{\Theta}_{\theta_{0}})^{2}\geq\max_{\bm{a}}\frac{\left(\bm{a}^{\top}\bm{\lambda}\right)^{2}}{\bm{a}^{\top}\bm{C}\bm{a}}.$
(54)

We study the case when the BhB variance is unbounded. It is unbounded if and
only if there is a vector $\bm{a}^{\prime}$ that satisfies
$\bm{C}\cdot\bm{a}^{\prime}=0$ and $\bm{a}^{\prime\top}\bm{\lambda}\neq 0$.
This vector $\bm{a}^{\prime}$ exists if and only if the next system of
equations

$\begin{pmatrix}\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial
P_{\theta_{0}}(x_{i})\partial
P_{\theta_{0}}(x_{i})&\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial
P_{\theta_{0}}(x_{i})\partial^{2}P_{\theta_{0}}(x_{i})&\dots&\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial
P_{\theta_{0}}(x_{i})\partial^{n}P_{\theta_{0}}(x_{i})\\\
\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial^{2}P_{\theta_{0}}(x_{i})\partial
P_{\theta_{0}}(x_{i})&\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial^{2}P_{\theta_{0}}(x_{i})\partial^{2}P_{\theta_{0}}(x_{i})&\dots&\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial^{2}P_{\theta_{0}}(x_{i})\partial^{n}P_{\theta_{0}}(x_{i})\\\
\vdots&\vdots&\ddots&\vdots\\\ \sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial
P_{\theta_{0}}(x_{i})\partial^{n}P_{\theta_{0}}(x_{i})&\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial^{n}P_{\theta_{0}}(x_{i})\partial^{2}P_{\theta_{0}}(x_{i})&\dots&\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial^{n}P_{\theta_{0}}(x_{i})\partial^{n}P_{\theta_{0}}(x_{i})\\\
\end{pmatrix}\begin{pmatrix}k\\\ a_{2}\\\ \vdots\\\ a_{n}\end{pmatrix}=0,$
(55)

has a solution and also $k\neq 0$ because the condition that
$\bm{a}^{\prime\top}\bm{\lambda}\neq 0$ has to be fulfilled too.

The system of equations can be rewritten as,

$\displaystyle\begin{pmatrix}\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial
P_{\theta_{0}}(x_{i})\partial^{2}P_{\theta_{0}}(x_{i})&\dots&\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial
P_{\theta_{0}}(x_{i})\partial^{n}P_{\theta_{0}}(x_{i})\\\
\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial^{2}P_{\theta_{0}}(x_{i})\partial^{2}P_{\theta_{0}}(x_{i})&\dots&\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial^{2}P_{\theta_{0}}(x_{i})\partial^{n}P_{\theta_{0}}(x_{i})\\\
\vdots&\ddots&\vdots\\\
\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial^{n}P_{\theta_{0}}(x_{i})\partial^{2}P_{\theta_{0}}(x_{i})&\dots&\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial^{n}P_{\theta_{0}}(x_{i})\partial^{n}P_{\theta_{0}}(x_{i})\\\
\end{pmatrix}\begin{pmatrix}a_{2}\\\ \vdots\\\
a_{n}\end{pmatrix}=k\begin{pmatrix}\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial
P_{\theta_{0}}(x_{i})\partial P_{\theta_{0}}(x_{i})\\\
\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial
P_{\theta_{0}}(x_{i})\partial^{2}P_{\theta_{0}}(x_{i})\\\ \vdots\\\
\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial
P_{\theta_{0}}(x_{i})\partial^{n}P_{\theta_{0}}(x_{i})\end{pmatrix}.$

Rouché-Frobenius theorem states that the system has a solution if and only if
$\operatorname{rank}(A|\vec{b})=\operatorname{rank}(A)$. This is equivalent to
say that $\vec{b}$ is a linear combination of columns of the matrix $(A)$,
i.e.

$\displaystyle\begin{pmatrix}\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial
P_{\theta_{0}}(x_{i})\partial P_{\theta_{0}}(x_{i})\\\ \vdots\\\
\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial
P_{\theta_{0}}(x_{i})\partial^{n}P_{\theta_{0}}(x_{i})\end{pmatrix}=\sum_{l\neq
1}\alpha_{l}\begin{pmatrix}\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial^{l}P_{\theta_{0}}(x_{i})\partial
P_{\theta_{0}}(x_{i})\\\ \vdots\\\
\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial^{l}P_{\theta_{0}}(x_{i})\partial^{n}P_{\theta_{0}}(x_{i})\end{pmatrix}.$

Rearranging the previous expression,

$\displaystyle\begin{pmatrix}\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial^{m}P_{\theta_{0}}(x_{i})\partial
P_{\theta_{0}}(x_{i})\\\ \vdots\\\
\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial^{m}P_{\theta_{0}}(x_{i})\partial^{n}P_{\theta_{0}}(x_{i})\end{pmatrix}=\sum_{l=1,l\neq
m}\alpha^{\prime}_{l}\begin{pmatrix}\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial^{l}P_{\theta_{0}}(x_{i})\partial
P_{\theta_{0}}(x_{i})\\\ \vdots\\\
\sum_{i}\frac{1}{P_{\theta_{0}}(x_{i})}\partial^{l}P_{\theta_{0}}(x_{i})\partial^{n}P_{\theta_{0}}(x_{i})\end{pmatrix}.$
(56)

where $\alpha_{1}^{\prime}\neq 0$. That is, there is a column $m$ that is a
combination of the first column and other columns of $\bm{C}$ . Note that the
elements of the $l$ column of the rhs of Eq. (56) can be written as,

$\displaystyle\begin{pmatrix}\partial^{l}P_{\theta_{0}}(x_{1})&\partial^{l}P_{\theta_{0}}(x_{2})&\dots&\partial^{l}P_{\theta_{0}}(x_{N})\\\
\end{pmatrix}.\begin{pmatrix}\frac{1}{P_{\theta_{0}}(x_{1})}\partial^{j}P_{\theta_{0}}(x_{1})\\\
\frac{1}{P_{\theta_{0}}(x_{2})}\partial^{j}P_{\theta_{0}}(x_{2})\\\ \vdots\\\
\frac{1}{P_{\theta_{0}}(x_{N})}\partial^{j}P_{\theta_{0}}(x_{N})\\\
\end{pmatrix},$

where $j$ goes from 1 to $n$ for each element of the column. Defining

$\displaystyle A_{m}$
$\displaystyle=\begin{pmatrix}\partial^{m}P_{\theta_{0}}(x_{1})&\partial^{m}P_{\theta_{0}}(x_{2})&\dots&\partial^{m}P_{\theta_{0}}(x_{N})\\\
\end{pmatrix},$ (57) $\displaystyle B_{j}$
$\displaystyle=\begin{pmatrix}\frac{1}{P_{\theta_{0}}(x_{1})}\partial^{j}P_{\theta_{0}}(x_{1})\\\
\frac{1}{P_{\theta_{0}}(x_{2})}\partial^{j}P_{\theta_{0}}(x_{2})\\\ \vdots\\\
\frac{1}{P_{\theta_{0}}(x_{N})}\partial^{j}P_{\theta_{0}}(x_{N})\\\
\end{pmatrix}.$ (58)

Equation (56) can be written as,

$\displaystyle A_{m}\cdot B_{j}$
$\displaystyle=\sum_{l=1}\alpha_{l}\>A_{l}\cdot B_{j}\quad\forall j.$ (59)

This implies

$\displaystyle A_{m}$ $\displaystyle=\sum_{l=1}\alpha_{l}\>A_{l}.$ (60)

Substituting the definitions given in Eq. (57) in Eq. (60) we get

$\begin{pmatrix}\partial^{m}P_{\theta_{0}}(x_{1})&\partial^{m}P_{\theta_{0}}(x_{2})&\dots&-\sum_{i}^{N-1}\partial^{m}P_{\theta_{0}}(x_{i})\\\
\end{pmatrix}=\sum_{l}\alpha_{l}\begin{pmatrix}\partial^{l}P_{\theta_{0}}(x_{1})&\partial^{l}P_{\theta_{0}}(x_{2})&\dots&-\sum_{i}^{N-1}\partial^{l}P_{\theta_{0}}(x_{i})\\\
\end{pmatrix},$ (61)

We have proved that BhB is unbounded if and only if Eq. (61) is true. We note
that Eq. (61) is the same as Eq. (53). This is the same as saying that there
is no estimator satisfying the conditions. As stated before, the matrix of Eq.
(47) has at most $N$ linearly independent rows associated with $N$ independent
conditions. Consider the case where the first $N-1$ higher order derivatives
are linearly independent, then $n^{th}$order BhB bound will be equal to the
$(N-1)^{th}$ BhB bound if $n\geq N-1$.

###  B.1 Bhattacharyya bounds for pure states

In this section, we show that contrary to what was recently claimed in the
literature [30], Bhattacharyya bounds for pure states are computable and can
be different from the Cramer-Rao bound in general. For this aim, we first
explore why the conclusion exposed in [30], namely that for pure states
Bhattacharyya bounds are the same as Cramer-Rao bound is faulted. After that,
we present a simple counterexample of the conclusion. The proof in [30] is
based on the fact that for pure states, $L_{i}$ (which represent the analogous
of SLD but for higher derivatives) are bounded. From this, the authors
concluded that the Bhattacharyya matrix is then filled with unbounded terms
and when taking the optimization problem, one only looks at the first entry of
the matrix which corresponds to the quantum Fisher Information. As exposed in
[Yating], even when the operators $L_{i}$ are unbounded the product of $\rho
L_{i}L_{j}$ may not be so, and this is what appears in the BhBattacarya
matrix. The authors there argue that the usual Cramer-Rao bound (and for
extension all Bhattacharyya bounds) is only constructed for bounded operators.
For that reason, they “rederive” the Cramer-Rao bound for a newly defined
operator $Q=L\sqrt{\rho}$. In the same fashion, one could do the same with
Bhattacharyya bounds and write the Bhattacharyya matrix as a function of
$Q_{i}$ operators. [RR: This we have to explicitly do here. ]This is shown in
the next example. Consider a state
$|\psi\rangle=\sqrt{x}|0\rangle+\sqrt{1-x}|1\rangle$ where $x$ is the
parameter one wants to estimate. From this state one can construct $Q_{i}$
operators that fulfill

$\displaystyle\frac{1}{2}(\sqrt{\rho}Q^{\dagger}_{j}+Q_{j}\sqrt{\rho})=\frac{\partial^{j}\rho}{\partial
x^{j}},$ (62) $\displaystyle\sqrt{\rho}Q_{j}-Q^{\dagger}_{j}\sqrt{\rho}=0,$
(63)

where the second equation comes from the hermiticity of $L_{j}$ [Yating].
After simple algebra one gets that (up to second order in Bhattacharyya
bounds) in the basis $\\{|\psi\rangle,|\psi^{\perp}\rangle\\}$ the operators
read

$Q_{1}=\begin{pmatrix}\langle\psi|\frac{\partial\rho}{\partial
x}|\psi\rangle&0\\\ 2\langle\psi^{\perp}|\frac{\partial\rho}{\partial
x}|\psi\rangle&0\end{pmatrix},$ (64)

and

$Q_{2}=\begin{pmatrix}\langle\psi|\frac{\partial^{2}\rho}{\partial
x^{2}}|\psi\rangle&0\\\ 2\langle\psi^{\perp}|\frac{\partial^{2}\rho}{\partial
x^{2}}|\psi\rangle&0\end{pmatrix},$ (65)

where $\rho=|\psi\rangle\langle\psi|$. From this, one can rewrite the quantum
Bhattacharyya matrix (up to the second order) in terms of $Q_{j}$ as

$Q_{BhB}=\begin{pmatrix}\operatorname{Tr}{Q^{\dagger}_{1}Q_{1}}&\frac{1}{2}\left(\operatorname{Tr}{Q^{\dagger}_{1}Q_{2}}+\operatorname{Tr}{Q^{\dagger}_{2}Q_{1}}\right)\\\
\frac{1}{2}\left(\operatorname{Tr}{Q^{\dagger}_{1}Q_{2}}+\operatorname{Tr}{Q^{\dagger}_{2}Q_{1}}\right)&\operatorname{Tr}{Q^{\dagger}_{2}Q_{2}}\end{pmatrix}.$
(66)

In this way, one can compute the second orther Bhattacharyya bound as the
first element of the inverse of $Q_{BhB}$, i.e.
$\Delta_{BhB}=Q^{-1}_{BhB,11}$. In Fig. (4) one can clearly see that the
second order Bhattacharyya bound and CR bound differ.

Figure 4: Plot of the variances for different values of $x$. The upper one
corresponds to the QBhB bound and the other to CR bound.

L_1 and L_2 formulas: Doing so, one gets

$\displaystyle L_{1}$
$\displaystyle=\left(\begin{array}[]{cc}4\theta(1-2\lambda)\sin\left(2\theta^{2}\right)&4i\theta(2\lambda-1)\cos\left(2\theta^{2}\right)\\\
-4i\theta(2\lambda-1)\cos\left(2\theta^{2}\right)&4\theta(2\lambda-1)\sin\left(2\theta^{2}\right)\\\
\end{array}\right),$ (69) $\displaystyle L_{2}$
$\displaystyle=\left(\begin{array}[]{cc}\frac{4(2\lambda-1)\left(\theta^{2}\left(\cos\left(2\theta^{2}\right)-2\lambda+1\right)-(\lambda-1)\lambda\sin\left(2\theta^{2}\right)\right)}{(\lambda-1)\lambda}&\frac{4i(2\lambda-1)\left((\lambda-1)\lambda\cos\left(2\theta^{2}\right)+\theta^{2}\sin\left(2\theta^{2}\right)\right)}{(\lambda-1)\lambda}\\\
-\frac{4i(2\lambda-1)\left((\lambda-1)\lambda\cos\left(2\theta^{2}\right)+\theta^{2}\sin\left(2\theta^{2}\right)\right)}{(\lambda-1)\lambda}&\frac{4(2\lambda-1)\left((\lambda-1)\lambda\sin\left(2\theta^{2}\right)-\theta^{2}\left(\cos\left(2\theta^{2}\right)+2\lambda-1\right)\right)}{(\lambda-1)\lambda}\\\
\end{array}\right).$ (72)

From this, one can construct the BhBattacharyaa matrix which in this case
reads

##  Appendix C Remarks to Ref. [30]

In this section, we raise some remarks to some claims made in Ref. [30].
Specifically, we present some arguments against the implication that if the
matrix $\bm{Q}$ appearing in Eq. (10) is singular, the bound is divergent.
Moreover, we also argue against the statement that for pure states all
Bhattacharyya bounds give the same value as the CR bound. We make our
arguments stronger using simple counterexamples of both propositions. Even
though our results are only concerned with the Bhattacharyya bound, we think
one could extend them for the other bounds presented in Ref. [30].

###  C.1 Invertibility of the matrix

In Ref. [30], Gessner and Smerzi argue that if the BhB matrix $\bm{C}$ is
singular (non-invertible) the corresponding BhB bound will be divergent. While
this may seem true at first inspection, we prove that is not true on a closer
look. We now present our objections using the quantum version of the BhB
bound, even though the conclusions also hold in the classical case. Moreover,
we restrict ourselves to the second order. Let us first invoke the formula
which gives the Bhattacharyya bound for the quantum case:

$(\Delta\hat{\Theta}_{\theta_{0}})^{2}\geq\max_{\bm{a}}\frac{\left(\bm{a}^{\top}\bm{\lambda}\right)^{2}}{\bm{a}^{\top}\bm{Q}\bm{a}}=\bm{\lambda}^{\top}\bm{Q}^{-1}\bm{\lambda}.$
(73)

as explained in Ref. [30]. It is easy to see, that if the matrix is singular,
one could always obtain zero in the denominator of Eq. (73). However, it does
not follow that $\left(\Delta\theta_{\mathrm{est}}\right)_{BhB}^{2}$ will be
zero, since it may be the case that the eigenvector with null eigenvalue of
$Q_{\mathrm{BhB}}$ may be orthogonal to $\bm{\lambda}$. This would imply that
$\left(\Delta\theta_{\mathrm{est}}\right)_{BhB}^{2}$ is zero over zero which
is undefined so one cannot say that the variance is unbounded. In addition, we
find illuminating to show the discrepancy with the aforementioned claim using
a simple example. Consider we have a state given by
$\rho=\theta|0\rangle\langle 0|+(1-\theta)|1\rangle\langle 1|$ where $\theta$
is the parameter we would like to estimate. It is clear that this state
fulfills $\frac{\partial\rho}{\partial\theta}\neq 0$ and
$\frac{\partial^{2}\rho}{\partial\theta^{2}}=0$. Then, $L_{2}=0$ and
$L_{1}=1/\theta|0\rangle\langle 0|+1/(\theta-1)|1\rangle\langle 1|$ and the
QBhB matrix will read

$\displaystyle\bm{Q}_{BhB^{2}}=\begin{pmatrix}\frac{1}{\theta(1-\theta)}&0\\\
0&0\\\ \end{pmatrix},$ (74)

and

$\left(\Delta\theta_{\mathrm{est}}\right)_{BhB}^{2}=\sup_{\bm{a}}\frac{\left(\bm{a}^{\top}\bm{\lambda}\right)^{2}}{\bm{a}^{\top}Q_{\mathrm{BhB}}\bm{a}}=\sup_{\bm{a}}\frac{a^{2}_{1}}{\frac{a^{2}_{1}}{\theta(1-\theta)}}=\frac{1}{\frac{1}{\theta(1-\theta)}},$
(75)

which is essentially the QCRB, since the QFI is equal to
$\frac{1}{\theta(1-\theta)}$. So, in this case, even if the QBhB matrix is
singular, there exists an estimator that saturates the QBhB even though its
value corresponds to the one given by the QCRB.

###  C.2 Bhattacharyya bounds for pure states

The literature claims that for pure states, Abel quantum bounds and in
particular $n^{th}$order QBhB bounds are always equal to QCRB [30]. From our
results, this would imply that there is always an estimator satisfying all the
unbiasedness conditions. But in fact, it is not true. In the following, we
give a simple counter-example and we point out why wrong conclusions are drawn
in [30]. Consider the diagonal state $\rho=\lambda|0\rangle\langle
0|+(1-\lambda)|1\rangle\langle 1|$, where in the limit $\lambda\rightarrow 0$
it corresponds to a pure state. The state evolves under the unitary
$U_{\theta}=e^{-i\theta\sigma_{x}}$. For this model the first and the third
derivatives are equal,
$\partial_{\theta}\rho_{\theta}=\partial^{3}_{\theta}\rho_{\theta}$ and

$\displaystyle\bm{Q}_{BhB^{3}}=\left(\begin{array}[]{ccc}4(1-2\lambda)^{2}&0&-16(1-2\lambda)^{2}\\\
0&-\frac{4(1-2\lambda)^{2}}{(\lambda-1)\lambda}&0\\\
-16(1-2\lambda)^{2}&0&64(1-2\lambda)^{2}\\\ \end{array}\right).$ (79)

Note that the first and the third rows are proportional. The bound in
$\lambda\rightarrow 0$ (pure state), is given by

$\left(\Delta\theta_{\mathrm{est}}\right)_{BhB}^{2}=\lim_{\lambda\rightarrow
0}\sup_{\bm{a}}\frac{\left(\bm{a}^{\top}\bm{\lambda}\right)^{2}}{\bm{a}^{\top}\bm{Q}_{\mathrm{BhB}}\bm{a}}.$
(80)

The supreme is achieved for $\bm{a}=(a_{1},0,\frac{a_{1}}{4})$ since it is an
eigenvector of $\bm{Q}_{\mathrm{BhB}}$ with eigenvalue $0$. In that case one
has

$\left(\Delta\theta_{\mathrm{est}}\right)_{BhB}^{2}=\lim_{\lambda\rightarrow
0}\sup_{a_{1}}\frac{a_{1}}{0}=\lim_{\lambda\rightarrow 0}\infty=\infty$ (81)

In this example the BhB bound diverges while the quantum Fisher is finite –and
hence the QCRB–, showing that for pure states the quantum Abel bounds and in
particular the QBhB bounds can be different than the QCRB and one may be
interested in computing them. The key is to consider the limit
$\lambda\rightarrow 0$ after the supremum in Eq. (73). In the literature, the
limit was taken before the supremum and it was observed that the matrix
$\bm{Q}$ had some divergent elements that would not contribute to the bound.


physics over the next decade.

# Future Neutrino Experiments

B. Jamieson
on behalf of the Hyper-Kamiokande Collaboration  University of Winnipeg,
Winnipeg, MB, Canada, R3B 2E9

###### Abstract

The discovery of neutrino oscillations and the large mixing parameters in the
Pontecorvo Maki Nakagawa Sakata matrix has opened a window to search for CP-
violation in neutrinos. Two long-baseline neutrino experiments, DUNE and
Hyper-Kamiokande, are being prepared in the next decade to search for this CP-
violation with intense beams of muon and anti-muon neutrinos. The JUNO
experiment in China will study neutrino oscillations at medium baselines with
several goals, including determining the mass ordering of the neutrinos by
using the interference of the solar and atmospheric neutrinos. Short-baseline
neutrino oscillation experiments have seen anomalies in their data that could
either be explained by sterile neutrinos or new nuclear effects. Several
experiments, including the Fermilab short baseline program, Prospect-II and
JUNO-TAO are being planned to understand the source of these anomalies. This
paper provides a short review of these exciting experiments that will lead to
new discoveries in neutrino physics over the next decade.

##  I Introduction

The title of this paper is too broad to be covered in a five-page review, and
therefore will not be a complete review. Instead, the paper will look at the
main physics sensitivities of four of the many future neutrino experiments.
First, the paper will focus on a review of the physics sensitivities of the
two large long-baseline neutrino oscillation experiments Deep Underground
Neutrino Experiment (DUNE) dune  and Hyper-Kamiokande (Hyper-K) hyperk .
Second, the JUNO juno  physics program will be reviewed, and finally, the
Fermilab Short Baseline (SBL) program will be discussed fermilabsbl1 ;
fermilabsbl2 .

First, we consider the Pontecorvo Maki Nakagawa Sakata (PMNS) matrix $U$
describing the mixing of neutrinos. The PMNS matrix can be written in terms of
mixing angles $\theta_{12}$, $\theta_{13}$, and $\theta_{23}$ describing the
strength of the oscillations, and CP-violating phases ($\delta$, $\alpha_{1}$,
and $\alpha_{2}$). Neutrino oscillation experiments are not sensitive to the
so-called Majorana phases ($\alpha_{1}$ and $\alpha_{2}$) and will be ignored
here. Neutrinos produced in a particular flavor ($\nu_{e}$, $\nu_{\mu}$,
$\nu_{\tau}$) propagate as a combination of mass states ($\nu_{1}$, $\nu_{2}$,
$\nu_{3}$), and the PMNS matrix is:

$\displaystyle\begin{pmatrix}\nu_{e}\\\ \nu_{\mu}\\\ \nu_{\tau}\\\
\end{pmatrix}=$ $\displaystyle\begin{pmatrix}1&0&0\\\ 0&c_{23}&s_{23}\\\
0&-s_{23}&c_{23}\end{pmatrix}\begin{pmatrix}c_{13}&0&s_{13}{\rm
e}^{-i\delta}\\\ 0&1&0\\\ -s_{13}{\rm
e}^{i\delta}&0&c_{13}\end{pmatrix}\times$
$\displaystyle\begin{pmatrix}c_{12}&s_{12}&0\\\ -s_{12}&c_{12}&0\\\
0&-0&1\end{pmatrix}\begin{pmatrix}\nu_{1}\\\ \nu_{2}\\\ \nu_{3}\\\
\end{pmatrix}$

where $c_{ij}=\cos{\theta_{ij}}$, and $s_{ij}=\sin{\theta_{ij}}$.

There are three main open questions in neutrino oscillation physics. The first
is whether the oscillation of $\nu_{\mu}\rightarrow\nu_{e}$ is the same as
$\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ indicating a CP-violating phase
$\delta\neq 0$? The second is whether $\nu_{1}$ is the lightest neutrino
giving us a Normal Hierarchy (NH), or $\nu_{3}$ is the lightest neutrino
giving us the Inverted Hierarchy (IH)? Finally, there is interest in knowing
whether the mixing of the atmospheric neutrinos is maximal giving
$\theta_{23}=45^{\circ}$ or if it falls below or above this angle
($\theta_{23}$ octant)? The next generation of long and medium baseline
experiments (Hyper-K, DUNE, and JUNO) have a good chance of answering these
open questions.

##  II DUNE

DUNE is a 1300 km baseline neutrino experiment under construction that
consists of a 1.2 MW $\nu_{\mu}$ or $\bar{\nu}_{\mu}$ neutrino beam at
Fermilab (the PIP-II beamline), a near detector complex, and four 17 kT liquid
Ar Time Projection Chambers (LArTPCs) located 1.5 km underground at the 1300
km baseline at Sanford Underground Research Facility. The unoscillated beam
flux prediction is shown in FIG. 1 along with the flux from other current
generation neutrino beams at Fermilab.

Figure 1: DUNE beam flux compared to current generation neutrino beams at
Fermilab (Figure from C. Wilkinson, DUNE collaboration). Figure 2: DUNE
disappearance sensitivity (left), and appearance sensitivity (right) for
different potential values of the neutrino oscillation parameters. Figures
prepared by the DUNE collaboration pickering .

The DUNE near detector complex will consist of a LArTPC, muon-spectrometer and
a beam monitor. The LArTPC and muon-spectrometer will be movable along 28.5 m
long rails, allowing the neutrino spectrum to be measured for off-axis angles
between zero and $\sim 3.3^{\circ}$ with peak neutrino energies from 0.5 GeV
to 2.0 GeV. Linear combinations of these neutrino spectra can be used to
obtain the lepton kinematics of the oscillated flux times cross-section at the
far detector. By using the same nuclear target (Ar) as the far detector, and
using this prism technique, the systematic uncertainties in neutrino
oscillation parameter determination can be reduced.

Recent talks from DUNE collaborators show plots of the appearance sensitivity
and disappearance sensitivity for different numbers of years running. The
appearance and disappearance sensitivities for DUNE are shown in FIG. 2.

##  III Hyper-Kamiokande

The Hyper-K experiment is a world-leading neutrino experiment, building on the
success of Super-Kamiokande and T2K. The experiment has a broad and ambitious
physics program covering many neutrino sources, as well as proton decay
measurements. The experiment will consist of upgraded 1.2 MW $\nu_{\mu}$ and
$\bar{\nu}_{\mu}$ beam at the Japan Proton Accelerator Research Complex
(J-PARC) on the East coast of Japan, a near detector complex including a new
Intermediate Water Cherenkov Detector (IWCD), and the Hyper-Kamiokande water
Cherenkov far detector, shown in Fig. 3. A key improvement to Hyper-K over
Super-K is the improved PMTs with nearly doubled quantum efficiency and the
reduction in the timing resolution (making the timing resolution twice as
good).

Figure 3: Conceptual drawing of the Hyper-Kamiokande water Cherenkov far
detector (Figure from Hyper-K Collaboration).

The Hyper-K detector will have an 8$\times$ larger fiducial mass over Super-
Kamiokande. Hyper-K will: be 71 m tall and 68 m in diameter; have a 1 m wide
veto region on the sides; and have a 2 m veto region on the top and bottom.
The 20,000 box-and-line 50 cm PMTs will provide 20% photo coverage with 1.5 ns
timing resolution and double the quantum efficiency of the Super-K PMTs.
Additional photo coverage will come from multi-PMT modules consisting of 19
eight-inch PMTs. The mPMT modules will provide improved position, timing, and
direction resolution to act as in-situ calibration of the 50 cm PMTs.

Figure 4: Hyper-K beam flux (bottom) compared to the oscillation probabilities
for $\nu_{\mu}$ disappearance (top), $\nu_{e}$ appearance (middle) (Figure
from Hyper-K collaboration).

The Hyper-K beam flux and oscillation probabilities as a function of neutrino
energy are shown in FIG. 4. By using a $2.5^{\circ}$ off-axis beam, the beam
flux peaks at the oscillation maximum, and it subtends a smaller range of
energies.

Upgrades to the T2K off-axis near detector which includes a new Super Fine
Grained Detector (SFGD) with over two million scintillator cubes between two
horizontally stacked TPC modules to allow measurements of neutrino interaction
products at larger angles from the beam direction in a magnetized detector.

A new Intermediate Water Cherenkov Detector (IWCD) at $\sim 1$ km from the
beam target will be built, which uses the same nuclear target as the far
detector, and uses the prism technique to measure different beam fluxes for
beams between $1^{\circ}$ and $3^{\circ}$ by raising and lowering an
instrumented region in a pit of water. IWCD will use finer-grained
photosensors in multi-PMT modules to improve timing and fiducial volume
resolution. By using linear combinations of the measurements at different off-
axis angles, either a monochromatic beam or a beam matching the oscillated
flux can be mimicked in the oscillation analysis to reduce the flux $\times$
cross-section uncertainties by directly using the measured lepton kinematics.

One of the key goals of Hyper-K is the discovery of a non-zero CP-violating
phase in neutrino oscillation. How quickly this discovery can be made depends
on improvements to systematic uncertainties as can be seen in the plot of the
fraction of $\delta_{CP}$ values excluding $\sin{\delta_{CP}}=0$ shown in FIG.
5. Without improving over the systematic uncertainties of T2K, there could be
the discovery of a non-zero $\delta_{CP}$ within 5 years of running for 50% of
the phase-space of values.

Figure 5: Hyper-K percentage of $\delta_{CP}$ values excluding
$\sin{\delta_{CP}}=0$, where the shaded bands show the effect of improving
systematic uncertainties (Figure from Hyper-K collaboration).

In addition to the CP-violation search, Hyper-K will improve measurements of
solar neutrinos, atmospheric neutrinos, and searches for nucleon decay.
Hyper-K will also act as a detector for supernova neutrinos, where of order
100,000 neutrinos from supernova within the distance to the galactic center
would be detected, allowing for discrimination between different models of
supernovas. Finally, Hyper-K will observe supernova relic neutrinos within 10
years of starting operations.

Hyper-Kamiokande began construction in 2021 with the access tunnel excavation
completed on schedule in June 2022. The main cavern excavation is underway.
The 50 cm PMT production is on schedule, and inspection and testing are
ongoing. Development of the PMT covers is in progress, as is the development
of multi-PMT modules with LED sources of different wavelengths to do detailed
detector calibration. Hyper-K is currently projected to begin operations in
2027.

##  IV JUNO

Jiangmen Underground Neutrino Observatory (JUNO) is a multi-purpose liquid
scintillator experiment in China that will be 700 m underground near Jiagnmen
city. JUNO is a 20 kton LAB-based liquid scintillator-based detector with 78%
PMT coverage allowing it to reach a better than $3\%/\sqrt{E({\rm MeV})}$
energy resolution, and calibrations to allow an energy scale uncertainty of
less than one percent. The reactor $\bar{\nu}_{e}$ will come from two reactor
complexes at a baseline of about 52.5 km totaling $\sim 27$ GWth of thermal
power. The detector will detect $\sim$60 reactor anti-neutrinos per day.

For the JUNO baseline, there is an interesting separation of the Normal
Ordering (NO) and Inverted Ordering as a phase shift in the oscillations as
seen in FIG. 6. The $\bar{\nu}_{e}$ disappearance oscillation probability is
approximated as

$\displaystyle P_{\bar{\nu}_{e}\rightarrow\bar{\nu}_{e}}=$ $\displaystyle
1-\cos^{4}{2\theta_{13}}\sin^{2}{2\theta_{12}}\sin^{2}{\frac{\Delta
m^{2}_{12}L}{4E}}$
$\displaystyle-\sin^{2}{2\theta_{13}}\left(\cos^{2}{\theta_{12}}\sin^{2}{\frac{\Delta
m^{2}_{31}L}{4E}}\right.$
$\displaystyle\left.+\sin^{2}{\theta_{12}}\sin^{2}{\frac{\Delta
m^{2}_{32}L}{4E}}\right).$

The first two terms give the oscillation from the solar only parameters, while
the last term adds in the effect from the atmospheric oscillations.

Figure 6: JUNO events per 2 days of data taking, showing no oscillation flux
(black) the effect of solar oscillations parameters (grey-dashed), and the
atmospheric oscillation parameters in the normal ordering (blue) and inverted
ordering (red).

JUNO will be able to measure the solar ($\sin^{2}{\theta_{12}}$, $\Delta
m^{2}_{21}$) and atmospheric ($|\Delta m^{2}_{32}|$) neutrino oscillation
parameters with a precision of 0.6% or better with six years of operation.
JUNO is likely to determine the neutrino mass ordering as it was designed to
do before the long-baseline experiments begin their measurements.

##  V Fermilab short baseline neutrinos

At short baselines, several neutrino experiments are showing unexpected
results at $L/E\sim 1$ km/GeV. The $\bar{\nu}_{e}$ disappearance has been
observed in reactor experiments abazajian  with a ratio
$R=N_{exp}/N_{calc}=0.934\pm 0.024$ at distances of less than 1 km across
several reactor experiments. In $\nu_{e}$ source calibrations of solar
neutrino experiments there is an observed disappearance with
$R=N_{exp}/N_{calc}=0.84\pm 0.05$ arxiv1906.01739 . The most famous examples
are the LSND and MiniBoone results which have an appearance of $\nu_{e}$ and
$\bar{\nu}_{e}$ neutrino beam experimentsminiboone . The observations could be
explained by sterile neutrinos with $\Delta m^{2}\sim 1$ eV2, an unaccounted-
for nuclear effect, or other interpretations beyond the standard model. The
Fermilab short baseline program will probe these questions using a set of
three liquid-argon TPCs which have excellent shower reconstruction and
electron / gamma-ray separation. The near Short Baseline Neutrino Detector
(SBND) is 110 m from the Booster Neutrino Beam’s target, MicroBooNE is 470 m
from the target, and ICARUS is 600 m from the target.

With three years of data from these three detectors, the sterile neutrino 3+1
parameter space can be covered at the 3-5 $\sigma$ level. With the detailed
interaction data of the LArTPCs nuclear effects can be studied with greater
certainty, and by using the same nuclear target for all three experiments,
systematic uncertainties can be kept to a minimum.

##  VI Conclusion

There are many exciting neutrino experiments under construction around the
world, many of which could not be covered in this short conference
proceedings. With these experiments, the next decade of measurements could see
us determine if there is CP violation in neutrinos, determine the mass
hierarchy of neutrinos, determine the octant of $\theta_{23}$, get closer to
measuring the absolute neutrino mass, and make progress on understanding short
baseline neutrino anomalies. With the higher precision experiments being
prepared we will likely discover new puzzles related to the properties of
neutrinos.

###### Acknowledgements.

This document is adapted from the “Instruction for producing FPCP2003
proceedings” by P. Perret and from eConf templates templates-ref . I would
like to thank Dr. Anne Schukraft for use of her slides on the Fermilab short
baseline program, Dr. Luke Pickering for his slides on DUNE, and Dr. Jinnan
Zhang for his slides on JUNO. I would also like to thank my collaborators on
Hyper-Kamikande for selecting me to present this talk, and for their tireless
efforts in preparing a world class experiment.

## References

  * (1) (DUNE Collaboration), “Long-baseline neutrino oscillation physics potential of the DUNE experiment,” Eur. Phys. J. C 80, 978 (2020).
  * (2) (Hyper-Kamiokande Collaboration), “Hyper-Kamiokande Design Report,” arXiv:1805.04163 (2018).
  * (3) (JUNO Collaboration), “JUNO physics and detector,” Progress in Particle and Nuclear Physics 123 (2022) 103927.
  * (4) (Fermilab SBL Collaboration), “A Proposal for a Three Detector Short-Baseline Neutrino Oscillation Program in the Fermilab Booster Neutrino Beam,” arXiV:1503.01520 (2015).
  * (5) Machado, P. A.N., Palmara, O., and Schmitz, D. W., “The Short-Baseline Neutrino Program at Fermilab,” Annu. Rev. Nucl. Part. Sci. 2019.69:363-387.
  * (6) K.N. Abazajian, and others, “Light sterile neutrinos: a white paper,” arXiv:1204.5379 (2012).
  * (7) L. Pickering, and DUNE collaboration, “DUNE long-baseline oscillations,” Lake Louise Winter Institution 2022 presentation.
  * (8) Sebastian Böser, Christian Buck, Carlo Giunti, Julien Lesgourgues, Livia Ludhova, Susanne Mertens, Anne Schukraft, Michael Wurm, “Status of Light Sterile Neutrino Searches,” arXiv:1906.01739 (2019).
  * (9) (MiniBooNE Collaboration), “Significant excess of electronlike events in the MiniBooNE short-baseline neutrino experiment,” Phys. Rev. Lett. 121 (2018), 221801.
  * (10) http://www.slac.stanford.edu/econf/editors/eprint-template/instructions.html