Multiparameter estimation for qubit states with collective measurements: a case study

In this scenario, our task is to estimate the parameter, labelled as θ, that is encoded in the qudit quantum state ρ(θ). Generally, we obtain information about, and hence an estimate for θ, from the measurement statistics of N_total independent and identical copies of ρ(θ). In particular, we consider N_total = νN, whereby we group N copies together as an N-qudit ensemble described by the state ρ_N (θ) = ρ(θ)^⊗N , and measure them with a probability-operator measurement (POM), specified by the set of N-qudit operators {M_ℓ }, with ${M}_{\ell }\geqslant 0\,\forall \ell ,\,{\sum }_{\ell }{M}_{\ell }={1}_{{d}^{N}}$. From a single measurement, a single outcome k will be recorded, with the probability of getting outcome k = ℓ given by the Born's rule, i.e., p_ℓ (θ) = Tr{ρ_N (θ)M_ℓ }. Then, upon repeating over ν rounds of measurement, we obtain a sequence of ν measurement outcomes {k₁, k₂, ⋯, k_ν } as our data D, from which we construct an estimate ¹ ${\hat{\theta }}_{N}(D)$, using for example the maximum-likelihood principle [52, 53].

In this work, we will consider three specifications. Firstly, our choice of estimation precision quantifier is the mean squared error (MSE), ${{\Delta}}^{2}{\hat{\theta }}_{N}$, which is the expected squared difference between the true value of the parameter and its estimate. That is,

$$\begin{equation}{{\Delta}}^{2}{\hat{\theta }}_{N}{:=}\mathbb{E}[{(\theta -{\hat{\theta }}_{N}(D))}^{2}]=\sum\limits _{D}L(D\vert \theta ){(\theta -{\hat{\theta }}_{N}(D))}^{2},\end{equation} \tag{ 1 }$$

where L(D|θ) is the likelihood of observing the data D, given the true parameter value θ. That is, for any D = {k₁, k₂, ⋯, k_ν } which contains ν_ℓ occurrences of the outcome ℓ,

$$\begin{equation}L(D\vert \theta )=\prod\limits _{\ell }{p}_{\ell }{(\theta )}^{{\nu }_{\ell }},\quad \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\enspace \sum\limits _{\ell }\,{\nu }_{\ell }=\nu .\end{equation} \tag{ 2 }$$

Secondly, we consider locally unbiased estimators, for which the MSE will be lower bounded by the Cramér–Rao bound (CRB) [54, 55]:

$$\begin{equation}\nu {{\Delta}}^{2}{\hat{\theta }}_{N}\,\geqslant \,\frac{1}{{F}_{N;\theta }}\,\geqslant \,\frac{1}{{\mathcal{F}}_{N;\theta }},\end{equation} \tag{ 3 }$$

where

$$\begin{equation}{F}_{N;\theta }={F}_{N;\theta }[{\rho }_{N}(\theta ),\left\{{M}_{\ell }\right\}]{:=}\sum\limits _{\ell }\frac{{\dot{p}}_{\ell }{(\theta )}^{2}}{{p}_{\ell }(\theta )},\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad {\mathcal{F}}_{N;\theta }={\mathcal{F}}_{N;\theta }[{\rho }_{N}(\theta )]{:=}\underset{\left\{{M}_{\ell }\right\}}{\mathrm{max}}{F}_{N;\theta },\end{equation} \tag{ 4 }$$

are respectively the Fisher information (FI) and the quantum FI (QFI)—the largest FI upon optimizing the choice of measurement, with ${\dot{p}}_{\ell }(\theta )={\partial }_{\theta }{p}_{\ell }(\theta )$. Lastly, we here focus on asymptotic quantum inference [35], where we study the ultimate precision obtainable in the limit of unlimited number of repetitions, ν → ∞. Realistically, we do not have unlimited repetitions of course, but the CRB will get sufficiently tight with large but finite ν. Then, the inequalities in the CRB, equation (3), are saturated asymptotically [35, 54], and we may turn our attention to the QFI as the figure of merit for the optimal attainable estimation precision.

Given ρ_N (θ), the QFI and a choice of measurement that achieves it is well known [56, 57]:

$$\begin{equation}{\mathcal{F}}_{N;\theta }[{\rho }_{N}(\theta )]=\mathrm{T}\mathrm{r}\left\{{\rho }_{N}(\theta ){L}_{\theta }^{2}\right\},\end{equation} \tag{ 5 }$$

where the Hermitian symmetric-logarithmic derivative (SLD) L_θ is defined implicitly by

$$\begin{equation}{\partial }_{\theta }{\rho }_{N}(\theta ){:=}\frac{1}{2}\left\{{L}_{\theta },{\rho }_{N}(\theta )\right\},\end{equation} \tag{ 6 }$$

where $\left\{a,b\right\}=ab+ba$ is the anti-commutator. Moreover, the QFI can be attained by choosing the projective measurement into the eigenstates of the SLD, though in general cases there could be other possible choices as well.

As a remark, for ρ_N (θ) = ρ(θ)^⊗N having a product state structure, it follows from additivity of QFI [58] that ${\mathcal{F}}_{N;\theta }=N{\mathcal{F}}_{1;\theta }$, and the SLD will have a single-particle operator or an independent-sum structure, i.e., ${L}_{\theta }={\sum }_{i=1}^{N}{\,L}_{\theta }^{(i)}$. Hence, local measurements are sufficient to attain the QFI. This is expected and can be understood intuitively: with uncorrelated information carriers, we will never obtain more than the sum of optimal information from each individual. In this case then, the division into ν repetitions and N-ensemble can be seen as superfluous, as whatever the ν and N separately are for fixed N_total = νN → ∞, they all correspond to the same situation of having N_total independent local measurements on N_total independent qudits i.e., ${N}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}{{\Delta}}^{2}{\hat{\theta }}_{N}\geqslant {\mathcal{F}}_{1;\theta }^{-1}$. More generally though, it does matter how we divide the N_total copies into ν repetitions of N-qudit ensemble, especially when we consider POM beyond local measurements on each qudits. Such a consideration is necessary, for example, when the local measurements are noisy due to some imperfect implementation, e.g., crosstalks mixing up the measurement outcomes locally, resulting in the effective POM elements no longer being rank-1 projectors. On one hand, insistence on using local measurements would thus have the achievable precision of estimating θ in ρ_N (θ) lower limited by ${N}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}{{\Delta}}^{2}{\hat{\theta }}_{N}\geqslant \lambda {({\mathcal{F}}_{1;\theta })}^{-1}$ for some constant λ > 1 that is characteristic of the noise. On the other hand, by allowing collective or joint measurements, interestingly, one can show that the measurement noise can be effectively negated, and now ${N}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}{{\Delta}}^{2}{\hat{\theta }}_{N}\geqslant {c}_{N}{({\mathcal{F}}_{1;\theta })}^{-1}$, with an N-dependent coefficient c_N that →1 in the large N limit [59]. In this situation then, it pays off to have larger N, so long as ν is still sufficiently large to keep the CRB tight. In this work we shall not consider complications due to imperfect measurements, though as we will see next, we still nevertheless need to consider collective measurement for joint estimation of multiple parameters, and keep the distinction between ν repetitions and grouping into N-qudit ensemble for good.

We move on to review the estimation of multiple parameters θ = {θ₁, θ₂, ⋯, θ_K } encoded in the state ρ( θ ). As in previous section, we consider N identical copies of the qudit grouped together and described by the state ρ_N ( θ ) = ρ( θ )^⊗N , and construct estimates of the parameters using asymptotic locally unbiased estimators ${\hat{\boldsymbol{\theta }}}_{N}=\left\{{\hat{\theta }}_{N;1},{\hat{\theta }}_{N;2},\cdots ,{\hat{\theta }}_{N;K}\right\}$ with the data D collected from ν repetitions of a measurement of ρ_N ( θ ) with some N-qudit POM {M_ℓ }. Then, the MSE is generalized to the covariance matrix ${\mathcal{C}}_{N}$, whose (i, j)th matrix element is given by

$$\begin{equation}{\mathcal{C}}_{N;i,j}=\mathbb{E}[({\theta }_{i}-{\hat{\theta }}_{N;i}(D))({\theta }_{j}-{\hat{\theta }}_{N;j}(D))]=\sum\limits _{D}\,L(D\vert \boldsymbol{\theta })({\theta }_{i}-{\hat{\theta }}_{N;i}(D))({\theta }_{j}-{\hat{\theta }}_{N;j}(D)),\end{equation} \tag{ 7 }$$

where p_ℓ ( θ ) = Tr{ρ_N ( θ )M_ℓ } and $L(D\vert \boldsymbol{\theta })={\prod }_{\ell }{p}_{\ell }{(\boldsymbol{\theta })}^{{\nu }_{\ell }}$ with ∑_ℓ ν_ℓ = ν, whereas the FI is generalized to the FI matrix F_{N; θ} , whose (i, j)th matrix element reads

$$\begin{equation}{F}_{N;i,j}=\sum\limits _{\ell }\frac{1}{{p}_{\ell }(\boldsymbol{\theta })}\frac{\partial {p}_{\ell }(\boldsymbol{\theta })}{\partial {\theta }_{i}}\frac{\partial {p}_{\ell }(\boldsymbol{\theta })}{\partial {\theta }_{j}}.\end{equation} \tag{ 8 }$$

It then follows that we have the multi-parameter CRB, a matrix inequality, which reads ²

$$\begin{equation}\nu {\mathcal{C}}_{N}\geqslant {F}_{N;\boldsymbol{\theta }}^{-1}.\end{equation} \tag{ 9 }$$

Equivalently, equation (9) can be written as a numerical bound $\nu \,\mathrm{T}\mathrm{r}\left\{{\mathcal{C}}_{N}W\right\}\geqslant \mathrm{T}\mathrm{r}\left\{{F}_{N;\boldsymbol{\theta }}^{-1}W\right\}$ for any non-negative K × K matrix W, known as the cost matrix or weight matrix. Then, while there are a rich choices of possible bounds one may consider, such as the right logarithmic derivative CRB [32], for simplicity of presentation, here we mention just two: first, the covariance matrix can be further lower limited by the Holevo bound, then by the multi-parameter QFI CRB [31, 34]:

$$\begin{equation}\nu \,\mathrm{T}\mathrm{r}\left\{{\mathcal{C}}_{N}W\right\}\geqslant \mathrm{T}\mathrm{r}\left\{{F}_{N;\boldsymbol{\theta }}^{-1}W\right\}\geqslant \underset{\left\{{X}_{i}\right\}}{\mathrm{min}}\left\{\mathrm{T}\mathrm{r}\left\{W\mathrm{R}\mathrm{e}Z\right\}+{\Vert}W\,\mathrm{I}\mathrm{m}\,Z{{\Vert}}_{1}\right\}\geqslant \mathrm{T}\mathrm{r}\left\{{\mathcal{F}}_{N;\boldsymbol{\theta }}^{-1}W\right\}.\end{equation} \tag{ 10 }$$

In the Holevo bound, ||⋅||₁ is the trace norm, Z is the matrix with elements Z_i,j = Tr{X_i X_j ρ_N ( θ )}, where ${\left\{{X}_{i}\right\}}_{i=1}^{K}$ is a set of Hermitian matrices satisfying $\mathrm{T}\mathrm{r}\left\{{X}_{i}\frac{\partial {\rho }_{N}(\boldsymbol{\theta })}{\partial {\theta }_{j}}\right\}={\delta }_{i,j}$ and Tr{ρ_N ( θ )X_i } = 0 ∀i. Meanwhile, in the multi-parameter QFI CRB, we have the QFI matrix ${\mathcal{F}}_{N;\boldsymbol{\theta }}$, with matrix elements ${\mathcal{F}}_{N;i,j}{:=}\frac{1}{2}\,\mathrm{T}\mathrm{r}\left\{{\rho }_{N}(\boldsymbol{\theta })\left\{{L}_{{\theta }_{i}},{L}_{{\theta }_{j}}\right\}\right\}$, where the SLDs are defined analogously as in equation (6), i.e., ${\partial }_{{\theta }_{i}}{\rho }_{N}(\boldsymbol{\theta }){:=}\frac{1}{2}\left\{{L}_{{\theta }_{i}},{\rho }_{N}(\boldsymbol{\theta })\right\}$. Note that the diagonal elements of the QFI matrix are exactly the QFI for the individual parameters as in single-parameter estimation, i.e., ${\mathcal{F}}_{N;i,i}={\mathcal{F}}_{N;{\theta }_{i}}$, while the diagonal elements of the covariance matrix, are, of course, the MSE for the individual parameters, i.e., $\nu {\mathcal{C}}_{N;i,i}=\nu {{\Delta}}^{2}{\hat{\theta }}_{i;N}$. Moreover, both the Holevo bound and QFI matrix is additive, i.e., the Holevo bound for ρ_N ( θ ) = ρ( θ )^⊗N is N⁻¹ times of that for ρ( θ ) [36, 39], and ${\mathcal{F}}_{N;\boldsymbol{\theta }}=N{\mathcal{F}}_{1;\boldsymbol{\theta }}$.

By the same reasoning as in the single-parameter cases, the first inequality in equation (10) can always be saturated [54] in the large repetition limit, ν → ∞. Remarkably then, in the large ensemble limit N → ∞, up to the first order convergence rate, i.e., the leading order of 1/N, the Holevo bound can always be attained as well (with some regularity conditions satisfied) [36–39, 60]. The last, weaker multi-parameter QFI CRB however, is only tight (asymptotically) and equal to the Holevo bound for asymptotically classical model [48], where the weak commutativity condition $\mathrm{T}\mathrm{r}\left\{{\rho }_{N}(\boldsymbol{\theta })[{L}_{{\theta }_{i}},{L}_{{\theta }_{j}}]\right\}=0\,\forall i\ne j$ is satisfied, where [a, b] = ab − ba is the commutator [19, 28, 61, 62]. Finally, should a more stringent condition $\mathrm{T}\mathrm{r}\left\{{\rho }_{N}(\boldsymbol{\theta }){L}_{{\theta }_{i}}{L}_{{\theta }_{j}}\right\}=0\,\forall i\ne j$ be met, we have ${\mathcal{F}}_{N;i,j}=0\,\forall i\ne j$, and the QFI matrix is diagonal. In this case then, the inverse of the QFI matrix can be done element-wise, and we have $\nu {\mathcal{C}}_{N;i,i}=\nu {{\Delta}}^{2}{\hat{\theta }}_{i;N}=1/{\mathcal{F}}_{N;i,i}=1/{\mathcal{F}}_{N;{\theta }_{i}}$ as N → ∞, i.e., for all the parameters, there exist a measurement scheme which allows us to estimate them simultaneously, each with an optimal asymptotic precision that is equal to that of the corresponding single-parameter estimation scenario. We shall refer to this as 'compatible' joint-parameter estimation, as in this case we can optimally estimate one particular parameter without affecting the other at all. As a remark, note that since the QFI matrix is non-negative, it can in fact be always diagonalized in a certain basis—equivalently, we can always define a new set of parameters θ '( θ ) from the initial parameters θ , where now their joint estimation is compatible. In the statistic literature, the subject of diagonalizing a FI matrix is also known as parameter orthogonality [63], and has also been discussed in the context of nuisance parameters in quantum estimation [64], where only a single parameter is of actual interest, and hence a reparametrization might be beneficial in constructing a good estimator for it. In general though, the new parameters θ ' might change their form for different true values of θ , as the diagonalization of the QFI matrix depends locally on θ . Here, we are concern with global parameter orthogonality, i.e., joint estimation of some given, fixed parameters θ for different possible true values, in which $\mathrm{T}\mathrm{r}\left\{{\rho }_{N}(\boldsymbol{\theta }){L}_{{\theta }_{i}}{L}_{{\theta }_{j}}\right\}=0\,\forall i\ne j$ is then a sufficient condition their compatibility. In this sense, our definition of compatibility here follows that of reference [28], which is different from those for example in references [44, 65–67], where compatibility is defined for saturation of the multi-parameter QFI CRB even if the QFI matrix is non-diagonal globally for all parameter values.

Consider a qubit described by the state ρ( θ ), which in the familiar Bloch-sphere representation as visualized in figure 1, has a Bloch vector s that is confined to the x–z plane, and deviates from the z-axis by some angle φ. We allow the length of the Bloch vector to be possibly further parametrized by the variable ɛ, such that

$$\begin{equation}\rho (\boldsymbol{\theta })=\rho (\varepsilon ,\varphi )=\frac{1}{2}\left(1+\boldsymbol{s}\cdot \boldsymbol{\sigma }\right),\quad \boldsymbol{s}=s(\varepsilon ){\hat{\mathrm{e}}}_{\boldsymbol{s}}(\varphi ),\quad {\hat{\mathrm{e}}}_{\boldsymbol{s}}(\varphi )=\mathrm{sin}\,\varphi {\hat{\mathrm{e}}}_{x}+\mathrm{cos}\,\varphi {\hat{\mathrm{e}}}_{z},\end{equation} \tag{ 11 }$$

where σ is the usual Pauli operator vector. The two parameters of interest of estimation are θ = {ɛ, φ}. Note that in references [36, 51] which study more general qubit model with up to three parameters, our qubit model here has also been studied as a sub-model with just two parameters, mainly from the prospect of analyzing the Holevo bound.

Zoom In Zoom Out Reset image size

It is helpful to think of the qubit as a spin-1/2 system, and so equivalently ρ( θ ) can be expressed as

$$\begin{equation}\rho (\varepsilon ,\varphi )={\mathrm{e}}^{-\mathrm{i}{\sigma }_{y}\varphi /2}\frac{{e}^{\beta {\sigma }_{z}}}{2\,\mathrm{cosh}\,\beta }{\mathrm{e}}^{\mathrm{i}{\sigma }_{y}\varphi /2},\end{equation} \tag{ 12 }$$

where tanh β = s(ɛ) and σ_ℓ , ℓ = x, y, z, is the Pauli-ℓ operator. With equation (12), N identical copies of the qubit in ρ( θ ) can then be treated as a collection of spin-1/2s, with the N-qubit state written concisely as

$$\begin{equation}{\rho }_{N}(\boldsymbol{\theta })={\rho }_{N}(\varepsilon ,\varphi )=\rho {(\varepsilon ,\varphi )}^{\otimes N}={\mathrm{e}}^{-\mathrm{i}{J}_{y}\varphi }\frac{{e}^{2\beta {J}_{z}}}{{(2\,\mathrm{cosh}\,\beta )}^{N}}{\mathrm{e}}^{\mathrm{i}{J}_{y}\varphi },\end{equation} \tag{ 13 }$$

with the generalization to N-spin angular momentum, i.e., $\boldsymbol{J}=\frac{1}{2}{\sum }_{i=1}^{N}{\boldsymbol{\sigma }}^{(i)}$, where σ ⁽ⁱ⁾ is the Pauli operator vector for the ith qubit. Note that, as can be seen immediately from equations (12) and (13), we can also treat our state parameter estimation problem as a quantum metrology or channel parameter estimation problem, where the state ρ_N ( θ ) is obtained by encoding the parameters θ into an input product state ρ₀ polarized in the +z direction, using N-independent unitary channel ${\mathcal{U}}_{\varphi }\sim \left\{{U}_{\varphi }={\mathrm{e}}^{\mathrm{i}{\sigma }_{y}\varphi /2}\right\}$ that each suffers a dephasing noise—${{\Lambda}}_{\eta }\sim \left\{\sqrt{(1+\eta )/2}1,\sqrt{(1-\eta )/2}{\sigma }_{y}\right\}$ with η = s(ɛ). That is, we can equivalently think of ρ_N ( θ ) in equation (13) as

$$\begin{equation}{\rho }_{N}(\boldsymbol{\theta })={{\Lambda}}_{\eta }^{\otimes N}{\circ}{\mathcal{U}}_{\varphi }^{\otimes N}[{\rho }_{0}],\quad {\rho }_{0}={\left(\frac{1}{2}(1+{\sigma }_{z})\right)}^{\otimes N}.\end{equation} \tag{ 14 }$$

In this work, we will work on estimation of θ in the state ρ_N ( θ ) as given in equation (13) regardless of its interpretation. Nevertheless, we will return for a short discussion in section 5 on comparison of estimation precision in the corresponding quantum metrology scenario, where instead of using just product input states for ρ₀, one may consider the use of entangled states as well.

First, let us work out the estimation precision limit, should we just estimate one of the parameter in the state ρ_N (ɛ, φ) in equation (13), with another parameter being actually known.

3.2.1. Estimation of ɛ

With φ being known and fixed, it is straightforward to find that

$$\begin{equation}{L}_{\varepsilon }=\frac{1}{1-s{(\varepsilon )}^{2}}\left[2{\mathrm{e}}^{-\mathrm{i}{J}_{y}\varphi }{J}_{z}{\mathrm{e}}^{\mathrm{i}{J}_{y}\varphi }-Ns(\varepsilon )1\right]\frac{\partial s(\varepsilon )}{\partial \varepsilon },\end{equation} \tag{ 15 }$$

and

$$\begin{equation}{\mathcal{F}}_{N;\varepsilon }=N\frac{1}{1-s{(\varepsilon )}^{2}}{\left(\frac{\partial s(\varepsilon )}{\partial \varepsilon }\right)}^{2}.\end{equation} \tag{ 16 }$$

Then, as ${\mathrm{e}}^{-\mathrm{i}{\sigma }_{y}\varphi /2}{\sigma }_{z}{\mathrm{e}}^{\mathrm{i}{\sigma }_{y}\varphi /2}={\hat{\mathrm{e}}}_{\boldsymbol{s}}\cdot \boldsymbol{\sigma }$, from equation (15), the optimal precision ${\mathcal{F}}_{N;\varepsilon }$ can be attained by performing the projective measurement specified by the directions $\pm {\hat{\mathrm{e}}}_{\boldsymbol{s}}$, i.e., $\left\{\frac{1+{\hat{\mathrm{e}}}_{\boldsymbol{s}}\cdot \boldsymbol{\sigma }}{2},\frac{1-{\hat{\mathrm{e}}}_{\boldsymbol{s}}\cdot \boldsymbol{\sigma }}{2}\right\}$ for each of the qubit, which has the respective probability ${p}_{\pm }(\varepsilon )=\frac{1\pm s(\varepsilon )}{2}$ of obtaining them, and so ɛ = s⁻¹(p₊(ɛ) − p₋(ɛ)). Note that in this case the optimal measurement does not depend on the parameter ɛ itself.

An unbiased estimator (in the asymptotic ν → ∞ limit) ${\hat{\varepsilon }}_{N}$ that attains the QFI ${\mathcal{F}}_{N;\varepsilon }$ for any N is the maximum-likelihood estimator (MLE) [52, 53]. More precisely, denote ${\nu }_{\pm }^{(i)}$ as the relative frequencies of obtaining the outcome of $\frac{1\pm {\hat{\mathrm{e}}}_{\boldsymbol{s}}\cdot {\boldsymbol{\sigma }}^{(i)}}{2}$ for the ith qubit over the ν repetitions, such that ${\nu }_{+}^{(i)}+{\nu }_{-}^{(i)}=1$ for all i. Then, our data is essentially D = {ν₊, ν₋} with ${\nu }_{\pm }{:=}{\sum }_{i=1}^{N}\frac{{\nu }_{\pm }^{(i)}}{N}$, and the MLE is defined by

$$\begin{equation}{\hat{\varepsilon }}_{N}(D){:=}\mathrm{arg}\underset{\varepsilon }{\mathrm{max}}\,L(D\vert \varepsilon )=\mathrm{arg}\underset{\varepsilon }{\mathrm{max}}\,{p}_{+}{(\varepsilon )}^{N\nu {\nu }_{+}}{p}_{-}{(\varepsilon )}^{N\nu {\nu }_{-}}.\end{equation} \tag{ 17 }$$

Note that, for data where ν₊ − ν₋ is non-negative (which will have overall likelihood approaching one as ν → ∞ by virtue of the central limit theorem such that $\sqrt{\nu }{\nu }_{\pm }$ converges to a Gaussian peaked at $\sqrt{\nu }{p}_{\pm }(\varepsilon )$), the MLE can be conveniently written as ${\hat{\varepsilon }}_{N}(D){:=}{s}^{-1}({\nu }_{+}-{\nu }_{-})$, which coincides with the linear inversion estimator. Otherwise, the maximization in equation (17) needs to be carried out explicitly, and the MLE will in general have a small bias (which vanishes as ν → ∞) as compared to the linear inversion estimator.

3.2.2. Estimation of φ

With ɛ being known and fixed instead, we find that

$$\begin{equation}{L}_{\varphi }=2s(\varepsilon ){\mathrm{e}}^{-\mathrm{i}{J}_{y}\varphi }{J}_{x}{\mathrm{e}}^{\mathrm{i}{J}_{y}\varphi },\end{equation} \tag{ 18 }$$

$$\begin{equation}{\mathcal{F}}_{N;\varphi }=Ns{(\varepsilon )}^{2},\end{equation} \tag{ 19 }$$

where the optimal precision ${\mathcal{F}}_{N;\varphi }$ can be obtained with the projective measurement onto the direction specified by $\pm {\hat{\mathrm{e}}}_{{\boldsymbol{s}}^{\prime }}$ for each qubit, where ${\hat{\mathrm{e}}}_{{\boldsymbol{s}}^{\prime }}=\mathrm{cos}\,\varphi {\hat{\mathrm{e}}}_{x}-\mathrm{sin}\,\varphi {\hat{\mathrm{e}}}_{z}$. Note that in this case the optimal measurement does depend on the parameter φ itself—a common feature in local estimation problem. In particular then, suppose we are performing measurement in the basis specified by $\pm {\hat{\mathrm{e}}}_{{\boldsymbol{s}}^{\prime }}$ with ${\hat{\mathrm{e}}}_{{\boldsymbol{s}}^{\prime }}=\mathrm{cos}\,\alpha {\hat{\mathrm{e}}}_{x}-\mathrm{sin}\,\alpha {\hat{\mathrm{e}}}_{z}$ where eventually α → φ, such that the respective probability is given by ${p}_{\pm }(\varphi )=\frac{1\pm s(\varepsilon )\mathrm{sin}(\varphi -\alpha )}{2}$, and then $\varphi =\alpha +{\mathrm{sin}}^{-1}(\frac{{p}_{+}(\varphi )-{p}_{-}(\varphi )}{s(\varepsilon )})$. In order to attain the QFI ${\mathcal{F}}_{N;\varphi }$ in the asymptotic limit ν → ∞ for any N we may consider the MLE as defined by

$$\begin{equation}{\hat{\varphi }}_{N}(D){:=}\mathrm{arg}\underset{\varphi }{\mathrm{max}}\,L(D\vert \varphi )=\mathrm{arg}\underset{\varphi }{\mathrm{max}}\,{p}_{+}{(\varphi )}^{N\nu {\nu }_{+}}{p}_{-}{(\varphi )}^{N\nu {\nu }_{-}},\end{equation} \tag{ 20 }$$

where now ${\nu }_{\pm }{:=}{\sum }_{i=1}^{N}\frac{{\nu }_{\pm }^{(i)}}{N}$ is the sum of the relative frequencies for the projection onto $\frac{1\pm {\hat{\mathrm{e}}}_{{\boldsymbol{s}}^{\prime }}\cdot {\boldsymbol{\sigma }}^{(i)}}{2}$. Similar to the discussion earlier, the MLE will coincide with the linear inversion estimator with ${\hat{\varphi }}_{N}(D){:=}\varphi +{\mathrm{sin}}^{-1}(\frac{{\nu }_{+}-{\nu }_{-}}{s(\varepsilon )})$ after taking α → φ, when $\frac{{\nu }_{+}-{\nu }_{-}}{s(\varepsilon )}$ is bounded between zero and one, which will be mostly the case as ν → ∞. As a remark, ${\mathcal{F}}_{N;\varphi }$ can only be attained by the said projective measurement onto the direction $\pm {\hat{\mathrm{e}}}_{{\boldsymbol{s}}^{\prime }}$ with α = φ for any s(ɛ) that is strictly not equal to zero or unity; for s(ɛ) = 1, arbitrary values of α would also attain ${\mathcal{F}}_{N;\varphi }$.

First, from a quick observation where the state ρ_N ( θ ) (13) and the two SLDs (15) and (18) are essentially all real in the same basis, we can immediately conclude that the weak commutativity condition holds. It then follows that the Holevo bound and the multi-parameter QFI CRB are equal, which dictate the first order convergence rate of the weighted trace of the covariance matrix as N → ∞. Moreover, one can readily verify that we have Tr{ρ_N ( θ )L_ɛ L_φ } = Tr{ρ_N ( θ )L_φ L_ɛ } = 0, and therefore the QFI matrix is diagonal, such that the two parameters are globally orthogonal and it is in principle possible to realize compatible joint estimation for them. That is, the optimal precisions in simultaneously estimating the two parameters in ρ_N (ɛ, φ) are given by their respective QFI as in the single-parameter scenarios; see equations (16) and (19).

Notice that, however, L_ɛ L_φ ≠ L_φ L_ɛ , which means that the optimal local measurement bases for estimation of ɛ and φ actually do not commute. In fact, as ${\hat{\mathrm{e}}}_{\boldsymbol{s}}\cdot {\hat{\mathrm{e}}}_{{\boldsymbol{s}}^{\prime }}=0$, they correspond to mutually unbiased bases (MUB) [27, 68]—acquiring optimal information about one parameter as such completely erases information about the other, and therefore we could never obtain simultaneous estimation of both parameters with optimal precision by measuring the N qubits separately. Therefore in order to attain still the multi-parameter QFI CRB, we must now consider collective or global measurements on the whole N-qubit ensemble instead.

3.3.1. Measurement with J ² and L_φ

We look for collective measurements that preserve information about the second parameter, say φ, when accessing information about the first parameter, say ɛ. Then, as φ serves as the angle of rotation generated by angular momentum J_y , see equation (13), it is inviting for us to consider measuring the angular momentum squared operator, J ², which is invariant under rotation, for estimation of ɛ. Indeed, as the eigenvalues or spectrum of the state ρ( θ ), which is $\frac{1\pm s(\varepsilon )}{2}$, is related exclusively to the parameter ɛ but not the parameter φ, the estimation of ɛ can also be viewed as estimation of the spectrum of ρ( θ ), where Keyl and Werner have shown that the J ² measurement on ρ_N ( θ ) produces estimates that have vanishing error as N → ∞ [69]. This fact has also been discussed in [17], and similar conclusions can also be observed in references [36, 51].

In spite of the insights and conclusions from the above references, let us here provide a slightly different but perhaps the most direct and explicit account on the MSE for estimation of ɛ with the J ² measurement, by applying the well-known error-propagation formula [70] for ν → ∞, which reads (see appendix A for details of the evaluation)

$$\begin{align}\hfill \nu {{\Delta}}^{2}{\hat{\varepsilon }}_{N}& =\frac{\mathrm{T}\mathrm{r}\left\{{({\boldsymbol{J}}^{2})}^{2}{\rho }_{N}(\boldsymbol{\theta })\right\}-\mathrm{T}\mathrm{r}{\left\{{\boldsymbol{J}}^{2}{\rho }_{N}(\boldsymbol{\theta })\right\}}^{2}}{\vert \frac{\partial }{\partial \varepsilon }\mathrm{T}\mathrm{r}\left\{{\boldsymbol{J}}^{2}{\rho }_{N}(\boldsymbol{\theta })\right\}{\vert }^{2}}\hfill \\ \hfill & =\frac{1}{N}\left[\frac{(3-2N)s{(\varepsilon )}^{4}+2(N-3)s{(\varepsilon )}^{2}+3}{2(N-1)s{(\varepsilon )}^{2}{\left(\frac{\partial s(\varepsilon )}{\partial \varepsilon }\right)}^{2}}\right].\hfill \end{align} \tag{ 21 }$$

Equation (21) applies to most estimator ${\hat{\varepsilon }}_{N}$ that takes in the mean of J ² (as estimated from the data by $\langle {\boldsymbol{J}}^{2}\rangle (D)={\sum }_{j=N/2-\lfloor N/2\rfloor }^{N/2}j(j+1){\nu }_{j}$, where ν_j is the relative frequency for the projection onto the j-subspace of J ² over ν repetitions) as input, and precisely estimates ɛ when the estimated mean is exactly the true mean value, i.e., ⟨ J ²⟩(D) = Tr{ρ_N ( θ ) J ²}. Then, we check that in the large N limit, $\nu {{\Delta}}^{2}{\hat{\varepsilon }}_{N}$ has the asymptotic behaviour

$$\begin{align}\hfill \nu {{\Delta}}^{2}{\hat{\varepsilon }}_{N}& =\frac{1}{N}\left[\frac{N}{N-1}(1-s{(\varepsilon )}^{2})+\frac{3s{(\varepsilon )}^{2}-6+3s{(\varepsilon )}^{-2}}{2(N-1)}\right]{\left(\frac{\partial s(\varepsilon )}{\partial \varepsilon }\right)}^{-2}\hfill \\ \hfill & \sim \frac{1}{N}(1-s{(\varepsilon )}^{2}){\left(\frac{\partial s(\varepsilon )}{\partial \varepsilon }\right)}^{-2}+O\left(\frac{1}{{N}^{2}}\right)\hfill \\ \hfill & ={\mathcal{F}}_{N;\varepsilon }^{-1}+O\left(\frac{1}{{N}^{2}}\right),\hfill \end{align} \tag{ 22 }$$

where O(1/N²) stands for terms that are of order 1/N² and beyond. We thus have a clear demonstration that as N gets ever larger, $\nu {{\Delta}}^{2}{\hat{\varepsilon }}_{N}\sim \frac{N}{N-1}{\mathcal{F}}_{N;\varepsilon }^{-1}\to {\mathcal{F}}_{N;\varepsilon }^{-1}$, i.e., indeed the J ² measurement permits construction of estimator which asymptotically estimates ɛ as precise as using the optimal projective measurement when assuming the parameter φ is known (cf equations (15) and (16)). Importantly, although $\nu {{\Delta}}^{2}{\hat{\varepsilon }}_{N}\to {\mathcal{F}}_{N;\varepsilon }^{-1}$ only asymptotically, the estimation of ɛ is global here, in the sense that the observable J ² is independent of the two parameters ɛ and ɛ, and so allows us to continue estimating φ without any potential lost of information as desired. That is, as J ² commutes with L_φ —one can explicitly confirm with equation (18), we may still effectively measure in the L_φ basis, and obtain, in the ν → ∞ limit for any N,

$$\begin{equation}\nu {{\Delta}}^{2}{\hat{\varphi }}_{N}={\mathcal{F}}_{N;\varphi }^{-1},\end{equation} \tag{ 23 }$$

using for example the MLE ${\hat{\varphi }}_{N}$ introduced in equation (20) in section 3.2.2.

In figure 2, we illustrate the convergence of the weighted sum of two the estimation precision, $\frac{1}{2}\nu {{\Delta}}^{2}{\hat{\varepsilon }}_{N}+\frac{1}{2}\nu {{\Delta}}^{2}{\hat{\varphi }}_{N}$, with the J ² and L_φ measurement strategy, to the ultimate bound $\frac{1}{2}{\mathcal{F}}_{N;\varepsilon }^{-1}+\frac{1}{2}{\mathcal{F}}_{N;\varphi }^{-1}$, for the chosen parametrization of s(ɛ) = 1 − ɛ²/8.

Zoom In Zoom Out Reset image size

3.3.2. Nearly-pure qubit estimation

Having established the general results for compatible joint estimation for ɛ and φ, let us look deeper into the scenario where the qubits are nearly pure, i.e., s(ɛ) ≈ 1. In particular, we focus on 'regular' parametrization where s(ɛ) is smooth around some value ɛ₀ with s(ɛ₀) = 1, such that we can perform the Taylor expansion in powers of δ ≪ 1 around ɛ₀ for some ɛ = ɛ₀ − δ. Then, while the conclusion below is generally true for any regular parametrization of s(ɛ) as well, for explicit illustration, let us consider specifically s(ɛ) = 1 − ɛ²/8 with ɛ ≪ 1, a choice that is motivated by a case study in quantum imaging superresolution problem: as reported in reference [50], estimating φ and ɛ in this parameterization could be seen as estimating the centroid and the separation between two close incoherent light sources (up to some multiplicative constant), respectively.

We shed some light as to how the J ² measurement works in this ɛ ≪ 1 regime. From the theory of angular momentum, see textbooks like references [71, 72], we know that the full N-spin Hilbert space, ${\mathcal{H}}_{\text{total}}$, can be decomposed into direct sum of orthogonal subspaces specified by two numbers, j and g(j):

$$\begin{equation}{\mathcal{H}}_{\text{total}}=\underset{j=N/2-\lfloor N/2\rfloor {}}{\overset{N/2}{\oplus }}\left(\underset{g(j)=1{}}{\overset{{\mu }_{j}}{\oplus }}{\mathcal{H}}_{j,g(j)}\right).\end{equation} \tag{ 24 }$$

Here, j = {N/2, N/2 − 1, ⋯, N/2 − ⌊N/2⌋} is nothing but the usual quantum number associated with the total angular momentum squared, while g(j) = {1, ⋯, μ_j − 1, μ_j } specifies the μ_j different subspaces for the same j value. Generally, the multiplicity factor μ_j depends on how the total angular momentum is composed; for our case of N spin-1/2s, it is known that, see for example references [73, 74], ³ ${\mu }_{j}=\left(\begin{matrix}\hfill N\hfill \\ \hfill N/2-j\hfill \end{matrix}\right)-\left(\begin{matrix}\hfill N\hfill \\ \hfill N/2-j-1\hfill \end{matrix}\right)=\frac{2j+1}{N+1}\left(\begin{matrix}\hfill N+1\hfill \\ \hfill N/2-j\hfill \end{matrix}\right)$. Denote ${\mathcal{P}}_{j}$ as the projector onto the j-angular momentum space, i.e., the space ${\oplus }_{g(j)=1}^{{\mu }_{j}}{\mathcal{H}}_{j,g(j)}$. Then, with the quantum state ρ_N (ɛ, φ) of equation (13), evidently ${\rho }_{N}(\varepsilon ,\varphi )={\sum }_{j=N/2-\lfloor N/2\rfloor }^{N/2}{\mathcal{P}}_{j}{\rho }_{N}(\varepsilon ,\varphi ){\mathcal{P}}_{j}$, and the probability that the N-qubit ensemble is found to have the total angular momentum value j is given by

$$\begin{align}\hfill {p}_{j}& =\mathrm{T}\mathrm{r}\left\{{\mathcal{P}}_{j}{\mathrm{e}}^{-\mathrm{i}{J}_{y}\varphi }{e}^{2\beta {J}_{z}}{\mathrm{e}}^{\mathrm{i}{J}_{y}\varphi }\right\}/{(2\,\mathrm{cosh}\,\beta )}^{N}={\mu }_{j}\sum\limits _{m=-j}^{j}{e}^{2\beta m}/{(2\,\mathrm{cosh}\,\beta )}^{N}\hfill \\ \hfill & ={\mu }_{j}\frac{\mathrm{sinh}\left(\beta (1+2j)\right)}{\mathrm{sinh}\beta }\frac{1}{{(2\,\mathrm{cosh}\,\beta )}^{N}}\hfill \\ \hfill & ={\mu }_{j}\frac{{\left[\varepsilon \sqrt{16-{\varepsilon }^{2}}\right]}^{-2j}\left[{(16-{\varepsilon }^{2})}^{1+2j}-{({\varepsilon }^{2})}^{1+2j}\right]}{16(1-{\varepsilon }^{2}/8)}\frac{{\varepsilon }^{N}}{{4}^{N}}{\left(1-\frac{{\varepsilon }^{2}}{16}\right)}^{N/2}.\hfill \end{align} \tag{ 25 }$$

Notice that according to expectations p_j does not depend on φ. As ɛ ≪ 1, we expand equation (25) in powers of ɛ, and the first two leading terms are

$$\begin{equation}{p}_{j}\approx {\mu }_{j}{\left(\frac{\varepsilon }{4}\right)}^{N-2j}\left[1+{a}_{j,N}{\varepsilon }^{2}\right],\end{equation} \tag{ 26 }$$

where ${a}_{j,N}=\frac{1}{16}(1-j-\frac{N}{2}-{\delta }_{j,0})$. More precisely, as the leading order of the series increases in powers of Nɛ², we see that the truncated series expansion in equation (26) should approximate equation (25) sufficiently well if Nɛ² ≪ 1. Then, from equation (26) we see that the leading contribution for the FI, up to Nɛ², is from the three subspaces with j = N/2, N/2 − 1, and N/2 − 2. That is,

$$\begin{equation}{p}_{\frac{N}{2}}\approx 1-\frac{N-1}{16}{\varepsilon }^{2},\end{equation} \tag{ 27 }$$

$$\begin{equation}{p}_{\frac{N}{2}-1}\approx (N-1)\frac{{\varepsilon }^{2}}{16}\left[1+\frac{(2-N)}{16}{\varepsilon }^{2}\right],\quad [\mathrm{d}\mathrm{e}\mathrm{fi}\mathrm{n}\mathrm{e}\mathrm{d}\,\mathrm{o}\mathrm{n}\mathrm{l}\mathrm{y}\,\mathrm{f}\mathrm{o}\mathrm{r}\,\,\,N\geqslant 2],\end{equation} \tag{ 28 }$$

$$\begin{equation}{p}_{\frac{N}{2}-2}\approx \frac{N(N-3){\varepsilon }^{4}}{512},\quad [\mathrm{d}\mathrm{e}\mathrm{fi}\mathrm{n}\mathrm{e}\mathrm{d}\,\mathrm{o}\mathrm{n}\mathrm{l}\mathrm{y}\,\mathrm{f}\mathrm{o}\mathrm{r}\,\,\,N\geqslant 4],\end{equation} \tag{ 29 }$$

and with ${f}_{j}=\frac{1}{{p}_{j}}{\left(\frac{\partial {p}_{j}}{\partial \varepsilon }\right)}^{2}$, we have

$$\begin{equation}{f}_{\frac{N}{2}}\approx \frac{{(N-1)}^{2}}{64}{\varepsilon }^{2},\end{equation} \tag{ 30 }$$

$$\begin{equation}{f}_{\frac{N}{2}-1}\approx \frac{N-1}{4}+\frac{3(N-1)}{64}\left(2-N-{\delta }_{N,2}\right){\varepsilon }^{2},\end{equation} \tag{ 31 }$$

$$\begin{equation}{f}_{\frac{N}{2}-2}\approx \frac{N(N-3)}{32}{\varepsilon }^{2},\end{equation} \tag{ 32 }$$

such that

$$\begin{align}\hfill {F}_{N;\varepsilon }& =\sum\limits _{j=N/2-\lfloor N/2\rfloor }^{N/2}{f}_{j}={f}_{\frac{N}{2}}+{f}_{\frac{N}{2}-1}+{f}_{\frac{N}{2}-2}+\left[\mathsf{\text{terms of}}\,\,\,O({\varepsilon }^{4})\mathsf{\text{and higher}}\right]\hfill \\ \hfill & \approx \frac{N-1}{4}+\left[\frac{{(N-1)}^{2}}{64}+\frac{3(N-1)}{64}(2-N-{\delta }_{N,2})+\frac{N(N-3)}{32}\eta (N-4)\right]{\varepsilon }^{2},\hfill \end{align} \tag{ 33 }$$

where η(⋅) is the unit step function, and here we define η(0) ≡ 1. For N ⩾ 4, the FI up to Nɛ² is then equal to

$$\begin{equation}{F}_{N;\varepsilon }=\frac{N}{4}\left(1+\frac{{\varepsilon }^{2}}{16}\right)-\left(\frac{1}{4}+\frac{5{\varepsilon }^{2}}{64}\right).\end{equation} \tag{ 34 }$$

Upon comparing with the QFI in equation (16) that now reads

$$\begin{equation}\enspace {\mathcal{F}}_{N;\varepsilon }=N\frac{4}{16-{\varepsilon }^{2}}=\frac{N}{4}\left(1+\frac{{\varepsilon }^{2}}{16}+\frac{{\varepsilon }^{4}}{256}+\cdots \,\right),\end{equation} \tag{ 35 }$$

we then confirm that in the large N limit, in decreasing order, the three j = N/2 − 1, N/2 − 2, N/2 subspaces of the J ² operator carry most of the information about ɛ when ɛ ≪ 1, and Nɛ² ≪ 1. In particular, from equations (27) and (28) we see that most detection events will be contributed by the j = N/2 subspace, followed by the j = N/2 − 1 subspace. It is however these rare detection events from the j = N/2 − 1 subspace that contain most information about ɛ, see equation (31); the abundant detection events from the j = N/2 subspace and the even rarer detection events from the j = N/2 − 2 subspace contribute only information with relative weight Nɛ² ≪ 1. The identification of a few specific angular momentum 'modes' that carry most information about the parameter in this regime is a reminiscent of the application of the spatial-mode demultiplexing technique for the superresolution imaging problem [23, 75]: at the close separation regime (corresponding to ɛ ≪ 1), most photons will be detected by the fundamental spatial mode—the transfer function of the imaging system (corresponding to j = N/2), followed by detection by the first 'excited' mode, which is the derivative of the transfer function with respect to the spatial dimension (corresponding to j = N/2 − 1). Similarly there, it is these rare detection events in the first 'excited' mode that carries most information and is mainly responsible for the superresolution capability of the spatial-mode demultiplexing technique. Note that, when N gets ever larger such that Nɛ² > 1 and the asymptotic series approximation equation (26) is no longer accurate, we expect now the information carried by the lower angular momentum modes will become increasingly important and needs to be taken into account, such that the J ² measurement still provides the asymptotic optimal precision as given by equations (21) and (22). All these can be summarized in figure 3, where we present the FIs associated the three largest angular momentum modes, both for the exact values and the approximate ones as given by equations (30)–(32), as well as their joint contribution.

Zoom In Zoom Out Reset image size

A HOM interferometer, schematically depicted in figure 4(a), consists of a 50:50 beam splitter (BS), with one qubit entering simultaneously from each of the input ports. The action of the BS can be summarized as a Hadamard transformation $\mathcal{U}$ linking the input and output creation operators, i.e.,

$$\begin{equation}\left(\begin{matrix}\hfill {b}_{\alpha ,1}^{{\dagger}}\hfill \\ \hfill {b}_{\alpha ,2}^{{\dagger}}\hfill \end{matrix}\right)=\mathcal{U}\left(\begin{matrix}\hfill {a}_{\alpha ,1}^{{\dagger}}\hfill \\ \hfill {a}_{\alpha ,2}^{{\dagger}}\hfill \end{matrix}\right)=\frac{1}{\sqrt{2}}\left(\begin{matrix}\hfill 1\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill -1\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {a}_{\alpha ,1}^{{\dagger}}\hfill \\ \hfill {a}_{\alpha ,2}^{{\dagger}}\hfill \end{matrix}\right),\end{equation} \tag{ 36 }$$

where the operator ${a}_{\alpha ,i}^{{\dagger}}({b}_{\alpha ,i}^{{\dagger}})$ specifies the creation of a qubit with the state α = {H, V} at the input (output) port i = {1, 2}. As usual, we have $\left[{a}_{\alpha ,i}^{{\dagger}},{a}_{\beta ,j}^{{\dagger}}\right]=\left[{b}_{\alpha ,i}^{{\dagger}},{b}_{\beta ,j}^{{\dagger}}\right]=0$ for bosons, and $\left\{{a}_{\alpha ,i}^{{\dagger}},{a}_{\beta ,j}^{{\dagger}}\right\}=\left\{{b}_{\alpha ,i}^{{\dagger}},{b}_{\beta ,j}^{{\dagger}}\right\}=0$ for fermions, which denies the possibility of having two fermions with same degrees of freedom at the same output. In this notation, we thus have ${\rho }_{N}(\varepsilon ,\varphi )={\otimes }_{i=1}^{N}{\rho }^{(i)}(\varepsilon ,\varphi )$, with

$$\begin{equation}{\rho }^{(i)}(\varepsilon ,\varphi )={\mathrm{e}}^{-\mathrm{i}{\sigma }_{y}\varphi /2}\left(\frac{1+s(\varepsilon )}{2}{a}_{\text{H},i}^{{\dagger}}\left\vert \text{vac}\right\rangle \!\left\langle \text{vac}\right\vert {a}_{\text{H},i}+\frac{1-s(\varepsilon )}{2}{a}_{\text{V},i}^{{\dagger}}\left\vert \text{vac}\right\rangle \!\left\langle \text{vac}\right\vert {a}_{\text{V},i}\right){\mathrm{e}}^{\mathrm{i}{\sigma }_{y}\varphi /2},\end{equation} \tag{ 37 }$$

where $\left\vert \text{vac}\right\rangle$ is the vacuum state, ${a}_{\alpha ,i}={({a}_{\alpha ,i}^{{\dagger}})}^{{\dagger}}$. One can then show that (see appendix B for explicit demonstration with the bosonic case) we have the following unique signatures at the output ports:

$$\begin{align}\hfill j=1& :\quad \text{all}\;\text{two}\;\text{qubits}\;\text{exit}\;\text{at}\;\text{the}\;\text{same}\;\text{ports;}\hfill \\ \hfill j=0& :\quad \text{all}\;\text{two}\;\text{qubits}\;\text{exit}\;\text{at}\;\text{distinct}\;\text{ports,}\,[\text{bosons}],\hfill \end{align} \tag{ 38 }$$

or

$$\begin{align}\hfill j=1& :\quad \text{all}\;\text{two}\;\text{qubits}\;\text{exit}\;\text{at}\;\text{distinct}\;\text{ports;}\hfill \\ \hfill j=0& :\quad \text{all}\;\text{two}\;\text{qubits}\;\text{exit}\;\text{at}\;\text{the}\;\text{same}\;\text{ports.}\,[\text{fermions}].\hfill \end{align} \tag{ 39 }$$

Note that these signatures are independent of the parameter φ, and hence information about it is preserved as it should. To then extract information about and estimate φ, we simply measure the projectors $\left\{\frac{1+{\hat{\mathrm{e}}}_{{\boldsymbol{s}}^{\prime }}\cdot \boldsymbol{\sigma }}{2},\frac{1-{\hat{\mathrm{e}}}_{{\boldsymbol{s}}^{\prime }}\cdot \boldsymbol{\sigma }}{2}\right\}$ at the end of each of the output ports, which of course corresponds to measuring the eigenspaces of L_φ . Note that, number-resolving detection technique is needed here to capture all the signatures.

Zoom In Zoom Out Reset image size

A natural question arises: can we still implement the J ² and L_φ measurement for the case of N ⩾ 3, using similar idea of observing different measurement statistics, including coincidence counts, at the various output ports? Here, we attempt to answer this question by considering a generalization of HOM interferometer to the Bell multiport setup [79, 80], as schematically depicted in figure 4(b), and can generally be realized with combinations of beam splitters, phase shifters, and mirrors [81]. In this case, one qubit simultaneously enters from the each of the N input of the Bell multiport, whose action, in generalization to equation (36), is the discrete Fourier transform, i.e.,

$$\begin{equation}{b}_{\alpha ,k}^{{\dagger}}=\sum\limits _{\ell =1}^{N}{\mathcal{U}}_{k\ell }{a}_{\alpha ,\ell }^{{\dagger}},\quad {\mathcal{U}}_{k\ell }=\frac{1}{\sqrt{N}}{\mathrm{e}}^{\frac{2\pi \mathrm{i}}{N}(k-1)(\ell -1)}.\end{equation} \tag{ 40 }$$

We work out explicitly the case for N = 3 here, where the Bell multiport is also known as the Bell tritter [79, 82] and has been realized experimentally for example in photonic systems [83, 84]. Leaving the details to appendix C, we find that the projection onto the j-eigenspace of J ² has the following unique signatures at the output ports:

$$\begin{align}\hfill j=\frac{3}{2}& :\quad \text{all}\;\text{three}\;\text{qubits}\;\text{exit}\;\text{at}\;\text{the}\;\text{same}\;\text{ports}\,\,\,\mathit{\text{or}}\,\,\,\text{distinct}\;\text{ports;}\hfill \\ \hfill j=\frac{1}{2}& :\quad \text{two}\;\text{qubits}\;\text{exit}\;\text{at}\;\text{the}\;\text{same}\;\text{ports,}\hfill \\ \hfill & \hspace{21.33955pt}\text{another}\;\text{at}\;\text{a}\;\text{different}\;\text{port,}\quad \quad [\text{bosons}],\hfill \end{align} \tag{ 41 }$$

or

$$\begin{align}\hfill j=\frac{3}{2}& :\quad \text{all}\;\text{three}\;\text{qubits}\;\text{exit}\;\text{at}\;\text{distinct}\;\text{ports;}\hfill \\ \hfill j=\frac{1}{2}& :\quad \text{two}\;\text{qubits}\;\text{exit}\;\text{at}\;\text{the}\;\text{same}\;\text{ports,}\hfill \\ \hfill & \hspace{21.33955pt}\text{another}\;\text{at}\;\text{a}\;\text{different}\;\text{port.}\quad \,\,\quad [\text{fermions}].\hfill \end{align} \tag{ 42 }$$

Thus, upon measuring the projectors $\left\{\frac{1+{\hat{\mathrm{e}}}_{{\boldsymbol{s}}^{\prime }}\cdot \boldsymbol{\sigma }}{2},\frac{1-{\hat{\mathrm{e}}}_{{\boldsymbol{s}}^{\prime }}\cdot \boldsymbol{\sigma }}{2}\right\}$ at the outputs, all the outcome statistics indeed allow us to implement exact J ² and L_φ measurements with a Bell tritter for N = 3.

For N ⩾ 4, exact realization of the measurement of J ² basis using interference signatures of the 'paths' in a passive, linear setup such as a HOM or Bell tritter, to the best knowledge of the author, is not available. For example, with the Bell multiport with N = 4 in equation (40), we do find that between the j = 2 and j = 1 eigenspaces the interference signatures are unique, but then there are no more room for another distinct signature for the j = 0 eigenspace. In fact, as computing the output probabilities for an arbitrarily given transformation $\mathcal{U}$ on the creation operators set $\left\{{a}_{\alpha ,\ell }^{{\dagger}}\right\}$ with fixed α is already hard [85, 86]—a central identity in boson sampling [87, 88], finding a transformation $\mathcal{U}$ that has output signatures that are different for each j-eigenspace, which now includes also inputs of different α, is even harder as N gets larger.

Let us hence briefly switch to discuss the problem from a quantum computation perspective, where instead of just a passive setup, we allow the usage of ancillary qubits and control gates to realize the measurement of the common eigenbasis of J ² and L_φ . Indeed then, upon identifying the J ² (and L_φ common) eigenbasis as the Schur basis [89]—the union of bases for both irreducible representations of the globally-symmetric unitary group as well as the permutation group—the task of projection onto this common eigenbasis can be mapped onto the outcomes of a general Schur transformation circuit, for which efficient constructions have been developed [89–92]. In particular, the whole circuit can be implemented with O(N³) two-level gates (unitaries that act non-trivially on two dimensions), which can then be approximately executed with the universal and fault-tolerant set of Clifford and T gates [93]. Then, should each of these two-level gates can be approximately executed with up to at most O(γ/N³) error, the full Schur transformation circuit can be realized with an error γ using overall O(N⁴ log(N/γ)) universal gates, as well as O(log N) ancilla [91, 92].

Consider now the situation in section 3.3.2, i.e., nearly-pure qubit estimation. In particular, we stick with the parameterization of s(ɛ) = 1 − ɛ²/8, ɛ ≪ 1.

As we have discussed in section 3.3.2, in the limit of ɛ ≪ 1 and Nɛ² ≪ 1, the rare j = N/2 − 1 detection events are the ones that contain most information about ɛ. We can see this from another angle, by writing equation (13) differently, namely, up to terms of O(Nɛ²),

$$\begin{align}\hfill {\rho }_{N}(\varepsilon ,\varphi )& ={\mathrm{e}}^{-\mathrm{i}{J}_{y}\varphi }{\left(\frac{1+s(\varepsilon )}{2}\left\vert \text{H}\right\rangle \!\left\langle \text{H}\right\vert +\frac{1-s(\varepsilon )}{2}\left\vert \text{V}\right\rangle \!\left\langle \text{V}\right\vert \right)}^{\otimes N}{\mathrm{e}}^{\mathrm{i}{J}_{y}\varphi }\hfill \\ \hfill & \approx \left(1-\frac{N{\varepsilon }^{2}}{16}\right){\mathrm{e}}^{-\mathrm{i}{J}_{y}\varphi }\left\vert \text{H}\cdots \text{H}\right\rangle \!\left\langle \text{H}\cdots \text{H}\right\vert {\mathrm{e}}^{\mathrm{i}{J}_{y}\varphi }\hfill \\ \hfill & \quad +\frac{{\varepsilon }^{2}}{16}{\mathrm{e}}^{-\mathrm{i}{J}_{y}\varphi }\left(\left\vert \text{H}\cdots \text{H}\text{V}\right\rangle \!\left\langle \text{H}\cdots \text{H}\text{V}\right\vert +\mathsf{p}\mathsf{e}\mathsf{r}\mathsf{m}\mathsf{u}\mathsf{t}\mathsf{a}\mathsf{t}\mathsf{i}\mathsf{o}\mathsf{n}\mathsf{s}\right){\mathrm{e}}^{\mathrm{i}{J}_{y}\varphi }\hfill \\ \hfill & \equiv {w}_{0}(\varepsilon ){\tau }_{0}(\varphi )+{w}_{1}(\varepsilon ){\tau }_{1}(\varphi ),\hfill \end{align} \tag{ 43 }$$

where ${w}_{0}(\varepsilon )=1-\frac{N{\varepsilon }^{2}}{16},{w}_{1}(\varepsilon )=\frac{N{\varepsilon }^{2}}{16}$, ${\tau }_{0}(\varphi )={\mathrm{e}}^{-\mathrm{i}{J}_{y}\varphi }\left\vert \text{H}\cdots \text{H}\right\rangle \!\left\langle \text{H}\cdots \text{H}\right\vert {\mathrm{e}}^{\mathrm{i}{J}_{y}\varphi }$, and ${\tau }_{1}(\varphi )={\mathrm{e}}^{-\mathrm{i}{J}_{y}\varphi }{\sum }_{i=1}^{N}{\tau }_{1,i}{\mathrm{e}}^{\mathrm{i}{J}_{y}\varphi }/N$ with ${\tau }_{1,i}=\left\vert \text{V}\right\rangle \!{\left\langle \text{V}\right\vert }^{(i)}{\otimes }_{j(\ne i)}^{N}\left\vert \text{H}\right\rangle \!{\left\langle \text{H}\right\vert }^{(j)}$.

Ignoring differences beyond O(Nɛ²) with the exact ρ_N (ɛ, φ), we may then consider the state as given by the last line of equation (43) instead, and suppose we send this state into the N-Bell multiport. Then, we verified explicitly up to N = 12 that the interference signatures for all the τ_1,i are the same, and we denote this set of signatures by ${\mathsf{S}}_{1}$. Furthermore, our numerical results show that ${\mathsf{S}}_{1}$ and the signature set for τ₀(φ), ${\mathsf{S}}_{0}$, are not exactly distinct, but they have an overlap, and the probability of having an interference signature from τ_1,i inside this overlapping subset is given by 1/N. The overall probability of observing ${\mathsf{S}}_{0}$ is then

$$\begin{equation}\mathsf{P}\mathsf{r}({\mathsf{S}}_{0})={w}_{0}(\varepsilon )+{w}_{1}(\varepsilon )\frac{1}{N}=1-\frac{(N-1){\varepsilon }^{2}}{16},\end{equation} \tag{ 44 }$$

while the rest of signatures have probability

$$\begin{equation}\mathsf{P}\mathsf{r}({\mathsf{S}}_{1}{\backslash}{\mathsf{S}}_{0})={w}_{1}(\varepsilon )\left(1-\frac{1}{N}\right)=\frac{(N-1){\varepsilon }^{2}}{16}.\end{equation} \tag{ 45 }$$

Over ν repetitions then, should we obtain the relative frequencies ν₀ and ν₁ = 1 − ν₀ for observing signatures in ${\mathsf{S}}_{0}$ and ${\mathsf{S}}_{1}{\backslash}{\mathsf{S}}_{0}$ respectively, we may estimate ɛ by linear inversion estimator as defined by

$$\begin{equation}\hat{\varepsilon }(D)=4\sqrt{\frac{{\nu }_{1}}{N-1}}=\sqrt{\frac{8(1+{\nu }_{1}-{\nu }_{0})}{N-1}},\end{equation} \tag{ 46 }$$

which in this case is always physical, and is unbiased for the approximate state at last line of equation (43). As a reminder to the reader, by obtaining ν₀ and ν₁ and then estimating ɛ, we have not disturbed any information about φ which we shall extract from L_φ measurement, as the interference signatures are invariant under an overall and equal transformation on all the qubits. In figure 5, we present a numerical simulation for the approximate estimation of different true values of ɛ with the state ρ_N ( θ ) as given by equation (13) (not the approximate state) with N = 4, 5, and 6, via the linear estimator (46). As N increases, we see that the biases in the estimator grow, which is unavoidable in this case, as now the number of angular momentum modes increases, and their contributions to both the set ${\mathsf{S}}_{0}$ and ${\mathsf{S}}_{0}{\backslash}{\mathsf{S}}_{1}$ increases, making the linear estimator less accurate. Evidently, this effect increases as ɛ increases as well.

Zoom In Zoom Out Reset image size

It turns out that the probability of getting the signature sets ${\mathsf{S}}_{0}$ and ${\mathsf{S}}_{1}{\backslash}{\mathsf{S}}_{0}$ are respectively equal to p_N/2 and p_N/2−1, see equations (44) and (45) versus equations (27) and (28). Is this correspondence accidental? Evidently, $\left\vert \text{H}\cdots \text{H}\right\rangle =\left\vert \enspace j=\frac{N}{2},{m}_{z}=\frac{N}{2}\right\rangle$, and so ${\tau }_{0}(\varphi )={\mathcal{P}}_{\frac{N}{2}}{\tau }_{0}(\varphi ){\mathcal{P}}_{\frac{N}{2}}$, such that $\mathrm{T}\mathrm{r}\left\{{\mathcal{P}}_{\frac{N}{2}}{\tau }_{0}(\varphi )\right\}=1$. Meanwhile, we have ${\tau }_{1}(\varphi )=({\mathcal{P}}_{\frac{N}{2}}+{\mathcal{P}}_{\frac{N}{2}-1}){\tau }_{1}(\varphi )({\mathcal{P}}_{\frac{N}{2}}+{\mathcal{P}}_{\frac{N}{2}-1})$, and one can show easily that $\mathrm{T}\mathrm{r}\left\{{\mathcal{P}}_{\frac{N}{2}}{\tau }_{1}(\varphi )\right\}=\frac{1}{N}$. Hence, just like $\mathsf{P}\mathsf{r}({\mathsf{S}}_{0})$ and $\mathsf{P}\mathsf{r}({\mathsf{S}}_{1}{\backslash}{\mathsf{S}}_{0})$, we have p_N/2 ≈ w₀(ɛ) + w₁(ɛ)(1/N) = 1 − (N − 1)ɛ²/16 and p_N/2−1 ≈ w₁(ɛ)(1 − 1/N) = (N − 1)ɛ²/16, indeed. It is, therefore, inviting for us to associate the 1/N overlap between ${\mathsf{S}}_{1}$ and ${\mathsf{S}}_{0}$ as originated from the 1/N fraction of τ₁(φ) in the ${\mathcal{P}}_{\frac{N}{2}}$ subspace, and therefore identify ${\mathsf{S}}_{0}$ as the signature set for j = N/2, while ${\mathsf{S}}_{1}{\backslash}{\mathsf{S}}_{0}$ as the signature set for j = N/2 − 1.

We take note that ${\mathsf{S}}_{1}$ has been so far not obtained from any arbitrary states in the j = N/2 − 1 subspace, but only those with m_z = N/2 − 1 via the states ${\tau }_{1,i}=\left\vert \text{V}\right\rangle \!{\left\langle \text{V}\right\vert }^{(i)}{\otimes }_{j(\ne i)}^{N}\left\vert \text{H}\right\rangle \!{\left\langle \text{H}\right\vert }^{(j)}$, i = 1, ⋯, N. Similarly, ${\mathsf{S}}_{0}$ is obtained from considering just τ₀, which is a particular state in the j = N/2 subspace with the specific value of m_z = N/2. To further confirm our assertion that ${\mathsf{S}}_{0}$ and ${\mathsf{S}}_{1}{\backslash}{\mathsf{S}}_{0}$ do unambiguously distinguish between the projection onto the two j-subspaces respectively, we check explicitly the signatures for all the eigenstates, $\left\vert \enspace j,{m}_{z},g(j)\right\rangle$, in each (j, g(j))-eigenspaces for j = N/2 and j = N/2 − 1 up to N = 6, and find that indeed all the kets with j = N/2 have signatures in ${\mathsf{S}}_{0}$, while all the kets with j = N/2 − 1 have signatures in ${\mathsf{S}}_{1}{\backslash}{\mathsf{S}}_{0}$. We then end this section by suggesting the following from our numerical observations:

$$\begin{align}\hfill & \text{Given}\,\,\,N\,\,\text{spin}-1\text{/}2\text{s}\;\text{with}\;\text{one}\;\text{spin}\;\text{entering}\;\text{simultaneously}\;\text{from}\;\text{each}\;\text{of}\;\hfill \\ \hfill & \text{the}\;\text{input}\;\text{ports}\;\text{of}\;\text{an}\,\,\,\mathit{\text{N}}-\text{Bell}\;\text{multiport,}\;\text{the}\;\text{interference}\;\text{signatures}\;\text{for}\;\hfill \\ \hfill & \text{the}\,\,\,j=N/2\,\,\,\text{and}\,\,\,j=N/2-1\text{subspaces}\;\text{are}\;\text{distinguishable.}\hfill \end{align} \tag{ 47 }$$

Put differently, we propose that an N-Bell multiport allows us to distinguish between projection onto the j = N/2 and j = N/2 − 1 subspaces of J ² from the outcome signatures. Regrettably, as the numerical analysis are basically just dry computation of probabilities and sorting of outcomes, we are not able to offer any graphical explanations or any easy-to-follow arguments for the proposed statement (47), other than to state we obtained the said results, and we welcome interested readers to approach the author for a copy of the numerical code.