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

Quantum estimation involving multiple parameters remains an important problem of both theoretical and practical interest. In this work, we study the problem of simultaneous estimation of two parameters that are respectively associate with the length and direction of the Bloch vector for identically prepared qubit states that is confined to a plane, where in order to obtain the optimal estimation precision for both parameters, collective measurements on multiple qubits are necessary. Upon treating $N$ qubits as an ensemble of spin-1/2 systems, we show that simultaneous optimal estimation for both parameters can be attained asymptotically with a simple collective measurement strategy -- first, we estimate the length parameter by measuring the populations in spaces corresponding to different total angular momentum values $j$, then we estimate the direction parameter by performing a spin projection onto an optimal basis. Furthermore, we show that when the state is nearly pure, for sufficiently but not arbitrarily large $N$, most information will be captured in the largest three $j$-subspaces. Then, we study how the total angular-momentum measurement can be realized by observing output signatures from a Bell multiport setup, either exactly for $N=2,3$, or approximately when the qubits are nearly pure for other $N$ values. We also obtain numerical results that suggest that using a Bell multiport setup, one can distinguish between projection onto the $j=N/2$ and $j=N/2-1$ subspaces from their respective interference signatures at the output.


Introduction
An important aspect in physical and information sciences is the acquisition of optimal knowledge about some quantities of interest, through the acts of measurement on the information carriers or probes. Quantum parameter estimation and quantum metrology study, in particular, the ultimate attainable precision in estimating the parameters of interest, under the constraints set by quantum theory. Estimation of a relative phase, such as one in an optical interferometer for gravitational wave sensing [1][2][3] or atomic states for frequency measurements and magnetic field sensing [4][5][6][7], temperature in quantum systems [8][9][10], as well as the spatial extent of composite light sources [11][12][13][14] arXiv:2109.07430v2 [quant-ph] 24 Mar 2022 in imaging tasks, are examples in which the ultimate precise estimators are highly desired.
While it is quite typical to have just a single parameter of interest, more generally, one might encounter quantum estimation problems that concern genuinely multiple parameters, and for which simultaneous or joint-parameter estimation is called for [15][16][17]. They include, for instance, estimation of phases and decoherence strength [18,19], numerous parameters for unitary operations [20][21][22], and all or some parameters characterizing the spatial extent of light sources, such as the centroid, brightness, separation and orientation [23][24][25]. Understandably, finding and implementing the optimal measurement for joint estimation of multiple parameters is generally more challenging than in the case of single-parameter estimation. In particular, joint estimation need not be "compatible", i.e., unless certain conditions are met, using no matter what measurement and estimation strategy, one can never simultaneously estimate all the parameters as precisely as one would optimally achieve when estimating just one parameter at a time, while assuming the other parameters are known and fixed. This is not surprising, as the measurement bases needed to attain the optimal precision for different parameters need not correspond to complementary observables [26,27].
Many important aspects and interesting results on multi-parameter estimation have been established. They include: the general conditions for "compatible" joint estimation [28], the mathematical conditions for the measurement operators being optimal in estimating all the parameters [29,30], and bounds on estimation precision and their evaluations [31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47]. However, given a specific multi-parameter quantum estimation problem, the measurement that attains the optimal joint estimation precision is still not explicitly deducible and recognizable, not to say its implementation scheme. Finding them remains thus a problem of both theoretical and practical interest.
In this work, we perform a case study on joint estimation of parameters that respectively define the length and the direction of the Bloch vector for an ensemble of identically prepared qubit states. In particular, the Bloch vector is assumed to be confined to a known plane, such that we are here essentially studying a two-parameter estimation problem: one for specifying the length, another one for specifying the angle within the plane. Within the larger context for multiparameter estimation problem, here are looking at a specific model which is known as asymptotically classical [48], and moreover, that the two parameters are globally orthogonal, such that their joint estimation is compatible in the large ensemble limit. Despite its simplicity and being specific, however, this elementary model has significance for various physical problems, such as in quantum sensing and quantum imaging. For example, joint estimation of a phase shift and the amplitude of phase diffusion [18,19], states parameters in twolevel quantum mechanical systems such as the purity and the phase for the polarization degrees of freedom of light [49], and the separation and centroid of two binary light sources which serves as a model to discuss resolution limits in imaging [50], can be fittingly mapped into our model under appropriate circumstances. Hereinafter, we will not make explicit references to the physical meaning of the parameters, keeping the estimation problem at the abstract level. We also take notice of two related and detailed work on studying optimal full estimation of qubit states [36,51], but with a somewhat different analysis and focus than what we will have here.
The organization of this work is as follows. In Sec. 2, we review the basic ingredients of quantum parameter estimation that we will use in this work. Readers that are familiar with the relevant concepts and results in quantum parameter estimation can then proceed directly to Secs. 3 and 4, where we formally introduce and study our twoparameter estimation problem. As a result, upon treating the qubits as a collection of spin-1/2 elementary systems, we show how the compatible joint estimation can be achieved in the large ensemble limit by the global or collective measurement of the total angular momentum squared operator, followed by a commuting, local projective measurement. Furthermore, we look deeper into the special case of nearly-pure qubit estimation, and highlight the link with the problem of sub-Rayleigh resolution imaging. We then discuss how such estimation scheme can be carried out in principle, either exactly or approximately, and by using either passive unitary setup or more general quantum circuits. From the obtained numerical results, we propose that the projection onto the two largest total angular momentum values could be achieved using a Bell multiport setup. Finally, we discuss briefly a few points in Sec. 5 and conclude in Sec. 6.

Preliminaries
For the sake of completeness, in this section we review the basic ingredients of quantum parameter estimation theory that will be relevant for this work. We start with the singleparameter estimation scenario, then followed by multi-parameter estimation scenario.

Single-parameter estimation
In this scenario, our task is to estimate the parameter, labeled 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 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 1 , k 2 , · · · , k ν } as our data D, from which we construct an estimate ‡θ N (D), using for example the maximum-likelihood principle [52,53]. ‡ The estimatorθ N depends of course on ν as well. Here we write out only the dependence on N but not ν explicitly to avoid overloading the notation.
In this work, we will consider three specifications. Firstly, our choice of estimation precision quantifier is the mean squared error (MSE), ∆ 2θ N , which is the expected squared difference between the true value of the parameter and its estimate. That is, where L(D|θ) is the likelihood of observing the data D, given the true parameter value θ. That is, for any D = {k 1 , k 2 , · · · , k ν } which contains ν occurrences of the outcome , Secondly, we consider locally unbiased estimators, for which the MSE will be lower bounded by the Cramér-Rao bound (CRB) [54,55]: where are respectively the Fisher information (FI) and the quantum FI (QFI)-the largest FI upon optimizing the choice of measurement, withṗ (θ) = ∂ θ p (θ). 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, Eq. (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] where the Hermitian symmetric-logarithmic derivative (SLD) L θ is defined implicitly by where {a, b} = 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 F N ;θ = N F 1;θ , and the SLD will have a singleparticle operator or an independent-sum structure, i.e., θ . 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 total ∆ 2θ N ≥ F −1 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 total ∆ 2θ N ≥ λ(F 1;θ ) −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 total ∆ 2θ N ≥ c N (F 1;θ ) −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.

Joint-parameter estimation
We move on to review the estimation of multiple parameters θ = {θ 1 , θ 2 , · · · , θ 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θ N = {θ N ;1 ,θ N ;2 , · · · ,θ N ;K } 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 C N , whose (i, j)-th matrix element is given by where p (θ) = Tr{ρ N (θ)M } and L(D|θ) = p (θ) ν with ν = ν, whereas the FI is generalized to the FI matrix F N ;θ , whose (i, j)-th matrix element reads It then follows that we have the multi-parameter CRB, a matrix inequality, which reads § Equivalently, Eq. (9) can be written as a numerical bound νTr{C N W } ≥ Tr{F −1 N ;θ W } 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 (RLD) 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] In the Holevo bound, · 1 is the trace norm, Z is the matrix with elements is a set of Hermitian matrices satisfying Tr{X i ∂ρ N (θ) ∂θ j } = δ i,j and Tr{ρ N (θ)X i } = 0 ∀i. Meanwhile, in the multi-parameter QFI CRB, we have the QFI matrix F N ;θ , with matrix elements F N ;i,j : where the SLDs are defined analogously as in Eq. (6), i.e., ∂ θ i ρ N (θ) : Note that the diagonal elements of the QFI matrix are exactly the QFI for the individual parameters as in single-parameter estimation, i.e., F N ;i,i = F N ;θ i , while the diagonal elements of the covariance matrix, are, of course, the MSE for the individual parameters, i.e., νC N ;i,i = ν∆ 2θ i;N . Moreover, both the Holevo bound and QFI matrix is additive, i.e., the Holevo bound for ρ N (θ) = ρ(θ) ⊗N is N −1 times of that for ρ(θ) [36,39], and By the same reasoning as in the single-parameter cases, the first inequality in Eq. (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][37][38][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 Tr{ρ N (θ)[L θ i , L θ j ]} = 0 ∀i = j is satisfied, where [a, b] = ab − ba is the commutator [19,28,61,62]. Finally, should a more stringent condition Tr{ρ N (θ)L θ i L θ j } = 0 ∀i = j be met, we have F N ;i,j = 0 ∀i = 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 νC N ;i,i = ν∆ 2θ i;N = 1/F N ;i,i = 1/F N ;θ 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 Tr{ρ N (θ)L θ i L θ j } = 0 ∀i = j is then a sufficient condition their compatibility. In this sense, our definition of compatibility here follows that of Ref. [28], which is different from those for example in Refs [44,[65][66][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.

Two-parameter estimation for qubit states
In this section, we apply the general formalism described above to the estimation of the length and direction of Bloch vector for a qubit state. Then, we show how one can asymptotically achieve compatible joint estimation for the two parameters, by a measuring an observable that can be thought of as the total angular momentum squared operator, followed by a commuting, local projective measurement. Furthermore, we look deeper into the case of nearly-pure qubit estimation, and reveal a connection with the problem of superresolution imaging in the sub-Rayleigh limit.

The qubit model
Consider a qubit described by the state ρ(θ), which in the familiar Bloch-sphere representation as visualized in Fig. 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 where σ is the usual Pauli operator vector. The two parameters of interest of estimation are θ = {ε, ϕ}. Note that in Refs. [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. It is helpful to think of the qubit as a spin-1/2 system, and so equivalently ρ(θ) can be expressed as where tanh β = s(ε) and σ , = x, y, z, is the Pauli-operator. With Eq. (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 with the generalization to N -spin angular momentum, i.e., J = 1 is the Pauli operator vector for the i-th qubit. Note that, as can be seen immediately from Eqs. (12,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 ρ 0 polarized in the +z direction, using In this work, we will work on estimation of θ in the state ρ N (θ) as given in Eq. (13) regardless of its interpretation. Nevertheless, we will return for a short discussion in Sec. 5 on comparison of estimation precision in the corresponding quantum metrology scenario, where instead of using just product input states for ρ 0 , one may consider the use of entangled states as well.

Single-parameter estimation
First, let us work out the estimation precision limit, should we just estimate one of the parameter in the state ρ N (ε, ϕ) in Eq. (13), with another parameter being actually known.
3.2.1. Estimation of ε With ϕ being known and fixed, it is straightforward to find that and Then, as e −iσyϕ/2 σ z e iσyϕ/2 =ê s · σ, from Eq. (15), the optimal precision F N ;ε can be attained by performing the projective measurement specified by the directions ±ê s , i.e., } for each of the qubit, which has the respective probability p ± (ε) = 1±s(ε) 2 of obtaining them, and so ε = s −1 (p + (ε) − p − (ε)). Note that in this case the optimal measurement does not depend on the parameter ε itself. An unbiased estimator (in the asymptotic ν → ∞ limit)ε N that attains the QFI F N ;ε for any N is the maximum-likelihood estimator (MLE) [52,53]. More precisely, denote ν (i) ± as the relative frequencies of obtaining the outcome of 1±ês·σ (i) 2 for the ith qubit over the ν repetitions, such that ν 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 √ νν ± converges to a Gaussian peaked at √ νp ± (ε)), the MLE can be conveniently written aŝ which coincides with the linear inversion estimator. Otherwise, the maximization in Eq. (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.

Estimation of ϕ
With ε being known and fixed instead, we find that where the optimal precision F N ;ϕ can be obtained with the projective measurement onto the direction specified by ±ê s for each qubit, whereê s = cos ϕê x −sin ϕê 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 ±ê s withê s = cos αê x −sin αê z where eventually α → ϕ, such that the respective probablity is given by p ± (ϕ) = 1±s(ε) sin(ϕ−α) 2 , and then ). In order to attain the QFI F N ;ϕ in the asymptotic limit ν → ∞ for any N we may consider the MLE as defined bŷ where now N is the sum of the relative frequencies for the projection onto . Similar to the discussion earlier, the MLE will coincide with the linear inversion estimator withφ N (D) := ϕ + sin −1 ( ν + −ν − s(ε) ) after taking α → ϕ, when ν + −ν − s(ε) is bounded between zero and one, which will be mostly the case as ν → ∞. As a remark, F N ;ϕ can only be attained by the said projective measurement onto the direction ±ê s with α = ϕ for any s(ε) that is strictly not equal to zero or unity; for s(ε) = 1, arbitrary values of α would also attain F N ;ϕ .

Joint-parameter estimation
First, from a quick observation where the state ρ N (θ) (13) and the two SLDs (15,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 multiparameter 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 Eq. (16) and Eq. (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, aŝ e s ·ê s = 0, they correspond to mutually unbiased bases [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 2 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 Eq. (13), it is inviting for us to consider measuring the angular momentum squared operator, J 2 , which is invariant under rotation, for estimation of ε. Indeed, as the eigenvalues or spectrum of the state ρ(θ), which is 1±s(ε) 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 2 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 Refs. [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 2 measurement, by applying the well-known error-propagation formula [70] for ν → ∞, which reads (see Appendix A for details of the evaluation) Eq. (21) applies to most estimatorε N that takes in the mean of J 2 (as estimated from the data by J 2 (D) = N/2 j=N/2− N/2 j(j +1)ν j , where ν j is the relative frequency for the projection onto the j-subspace of J 2 over ν repetitions) as input, and precisely estimates ε when the estimated mean is exactly the true mean value, i.e., J 2 (D) = Tr{ρ N (θ)J 2 }. Then, we check that in the large N limit, ν∆ 2ε N has the asymptotic behaviour where O(1/N 2 ) stands for terms that are of order 1/N 2 and beyond. We thus have a clear demonstration that as N gets ever larger, ν∆ 2ε i.e., indeed the J 2 measurement permits construction of estimator which asymptotically estimates ε as precise as using the optimal projective measurement when assuming the parameter ϕ is known (c.f. Eqs. (15,16)). Importantly, although ν∆ 2ε N → F −1 N ;ε only asymptotically, the estimation of ε is global here, in the sense that the observable J 2 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 2 commutes with L ϕ -one can explicitly confirm with Eq. (18), we may still effectively measure in the L ϕ basis, and obtain, in the ν → ∞ limit for any N , using for example the MLEφ N introduced in Eq. (20) in Sec. 3.2.2. In Fig. 2, we illustrate the convergence of the weighted sum of two the estimation precision, 1 2 ν∆ 2ε N + 1 2 ν∆ 2φ N , with the J 2 and L ϕ measurement strategy, to the ultimate bound 1 2 F −1 N ;ε + 1 2 F −1 N ;ϕ , for the chosen parametrization of s(ε) = 1 − ε 2 /8.

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 ε 0 with s(ε 0 ) = 1, such that we can perform the Taylor expansion in powers of δ 1 around ε 0 for some ε = ε 0 − δ. 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 − ε 2 /8 with ε 1, a choice that is motivated by a case study in quantum imaging superresolution problem: As reported in Ref. [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 2 measurement works in this ε 1 regime. From the theory of angular momentum, see textbooks like Refs. [71,72], we know that the full N -spin Hilbert space, H total , can be decomposed into direct sum of orthogonal m e a s u r e L ε a n d L ϕ s e p a r a t e ly m e a s u r e J 2 a n d L ϕ subspaces specified by two numbers, j and g(j): 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 Refs. [73] and [74] , Denote P j as the projector onto the j-angular momentum space, i.e., the space µ j g(j)=1 H j,g(j) . Then, with the quantum state ρ N (ε, ϕ) of Eq. (13), evidently ρ N (ε, ϕ) = N/2 j=N/2− N/2 P j ρ N (ε, ϕ)P j , and the probability that There is a typo in the equation (A4) of Ref. [74].
the N -qubit ensemble is found to have the total angular momentum value j is given by Notice that according to expectations p j does not depend on ϕ. As ε 1, we expand Eq. (25) in powers of ε, and the first two leading terms are where a j,N = 1 16 (1 − j − N 2 − δ j,0 ). More precisely, as the leading order of the series increases in powers of N ε 2 , we see that the truncated series expansion in Eq. (26) should approximate Eq. (25) sufficiently well if N ε 2 1. Then, from Eq. (26) we see that the leading contribution for the FI, up to N ε 2 , is from the three subspaces with j = N/2, N/2 − 1, and N/2 − 2. That is, and with f j = 1 such that where η(·) is the unit step function, and here we define η(0) ≡ 1. For N ≥ 4, the FI up to N ε 2 is then equal to Upon comparing with the QFI in Eq. (16) that now reads 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 2 operator carry most of the information about ε when ε 1, and N ε 2 1. In particular, from Eq. (27) and Eq. (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 Eq. (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 ε 2 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 ε 2 > 1 and the asymptotic series approximation Eq. (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 2 measurement still provides the asymptotic optimal precision as given by Eq. (21) and Eq. (22). All these can be summarized in Fig. 3, where we present the FIs associated the the three largest angular momentum modes, both for the exact values and the approximate ones as given by Eqs. (30)(31)(32), as well as their joint contribution.

Realization of the J 2 and L ϕ measurement
Let us now discuss the implementation of the J 2 and L ϕ measurement. While measuring the local operator L ϕ should be straightforward, measuring the global operator J 2 is more challenging. For the simplest scenario with N = 2, exact implementation of such a scheme is known for bosonic system: measurement of the J 2 basis is realized using the Hong-Ou-Mandel (HOM) interference effect [76,77], whereas the L ϕ measurement is performed in its local basis, after the interference. In fact, a recent experiment demonstrating superresolution estimation of the centroid and separation between two mutually incoherent sources in the sub-Rayleigh limit [78]-a task which can be well approximated as the estimation of ε and ϕ in our model [50]-utilizes essentially the same idea.
For the sake of clarity as well as introducing the notation, below, we first review the known results of how the setup of HOM interferometer followed by projective measurements at the end works for the N = 2 case. Note that here, our study  Their (exact) joint contribution is depicted by the black solid line, which we compare against with the whole information obtained from the J 2 measurement (black dashed line)-essentially the inverse of ∆ 2ε N in Eq. (21). Despite up to N = 300 in this case where the three modes still account for most of the information obtained with J 2 measurement, we see that it is decreasing as N ε 2 starts to get larger than one, which is a manifestation of the fact that information carried in the lower angular momentum modes are now becoming significant.
include both bosons and fermions. Then, new in this work we show how this can be generalized to the N = 3 case, before moving on to briefly discussing the more complicated scenarios with larger N , where the J 2 measurement is approached from a quantum circuit perspective, with ancilla qubits involved. Lastly, we look closer into the case of nearly-pure qubit estimation, where in conjunction with the discussion in Sec. 3.3.2, non-exact but approximate implementations that focus on distinguishing between just the j = N/2 and j = N/2 − 1 subspaces will be considered. In what follows, for convenience of the notation, we will use {h, v} to specify a basis for the qubit system, such that, e.g., σ z = |h h| − |v v|.

N = 2
A HOM interferometer, schematically depicted in Fig. 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 U linking the input and output creation operators, i.e., The HOM in (a) can be seen as a special case of a Bell multiport with N = 2, where the Hadamard transformation is now generalized to the discrete Fourier transform. For N = 3, the measurement statistics (including coincidences) at (between) the various output ports allow us to implement the J 2 and L ϕ measurement exactly; see Eqs. (41) and (42). For N ≥ 4, however, Bell multiport cannot implement J 2 perfectly, as not all interference signatures for different j values are now distinct. Nevertheless, we suggest that the interference signatures for j = N/2 and j = N/2 − 1 will be distinct, which we can then make use of in the limit of nearly-pure qubit estimation, as most detection events in this case are contained in this two subspaces; see Sec. 3.3.2. where the operator a † α,i (b † α,i ) specifies the creation of a qubit with the state α = {h, v} at the input (output) port i = {1, 2}. As usual, we have β,j } = 0 for fermions, which denies the possibility of having two fermions with same degrees of freedom at the same output. In this notation, where |vac is the vacuum state, a α,i = (a † α,i ) † . 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: j = 1 : all two qubits exit at the same ports; j = 0 : all two qubits exit at distinct ports, [bosons] (38) or j = 1 : all two qubits exit at distinct ports; j = 0 : all two qubits exit at the same ports.
[fermions] (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 { 1+ê s ·σ 2 , 1−ê s ·σ 2 } 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.

N = 3
A natural question arises: Can we still implement the J 2 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 Fig. 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 Eq. (36), is the discrete Fourier transform, i.e., 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 2 has the following unique signatures at the output ports: } at the outputs, all the outcome statistics indeed allow us to implement exact J 2 and L ϕ measurements with a Bell tritter for N = 3.

N ≥ 4
For N ≥ 4, exact realization of the measurement of J 2 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 Eq. (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 U on the creation operators set {a † α, } with fixed α is already hard [85,86]-a central identity in boson sampling [87,88], finding a transformation 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 2 and L ϕ . Indeed then, upon identifying the J 2 (and L ϕ common) eigenbasis as the Schur basis [89]-the union of bases for both irreducible representations of the globallysymmetric 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][90][91][92]. In particular, the whole circuit can be implemented with O(N 3 ) 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 3 ) error, the full Schur transformation circuit can be realized with an error γ using overall O(N 4 log(N/γ)) universal gates, as well as O(log N ) ancilla [91,92].

Nearly-pure qubit estimation: Approximate realization
Consider now the situation in Sec. 3.3.2, i.e., nearly-pure qubit estimation. In particular, we stick with the parameterization of s(ε) = 1 − ε 2 /8, ε 1. As we have discussed in Sec. 3.3.2, in the limit of ε 1 and N ε 2 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 Eq. (13) differently, namely, up to terms of O(N ε 2 ), where w 0 (ε) = 1 − N ε 2 16 , w 1 (ε) = N ε 2 16 , τ 0 (ϕ) = e −iJyϕ |h · · · h h · · · h|e iJyϕ , and τ 1 (ϕ) = e −iJyϕ N i=1 τ 1,i e iJyϕ /N with τ 1,i = |v v| (i) N j( =i) |h h| (j) . Ignoring differences beyond O(N ε 2 ) with the exact ρ N (ε, ϕ), we may then consider the state as given by the last line of Eq. (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 S 1 . Furthermore, our numerical results show that S 1 and the signature set for τ 0 (ϕ), 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 S 0 is then while the rest of signatures have probability Over ν repetitions then, should we obtain the relative frequencies ν 0 and ν 1 = 1 − ν 0 for observing signatures in S 0 and S 1 \ S 0 respectively, we may estimate ε by linear inversion estimator as defined bŷ which in this case is always physical, and is unbiased for the approximate state at last line of Eq. (43). As a reminder to the reader, by obtaining ν 0 and ν 1 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 Fig. 5, we present a numerical simulation for the approximate estimation of different true values of ε with the state ρ N (θ) as given by Eq. (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 S 0 and S 0 \ S 1 increases, making the linear estimator less accurate. Evidently, this effect increases as ε increases as well. It turns out that the probability of getting the signature sets S 0 and S 1 \ S 0 are respectively equal to p N/2 and p N/2−1 , see Eqs. (44) and (45) versus Eqs. (27) and (28). Is this correspondence accidental? Evidently, |h · · · h = j = N 2 , m z = N 2 , and so τ 0 (ϕ) = PN , and one can show easily that Tr{PN 2 τ 1 (ϕ)} = 1 N . Hence, just like Pr(S 0 ) and Pr(S 1 \ S 0 ), we have p N/2 ≈ w 0 (ε) + w 1 (ε)(1/N ) = 1 − (N − 1)ε 2 /16 and p N/2−1 ≈ w 1 (ε)(1 − 1/N ) = (N − 1)ε 2 /16, indeed. It is, therefore, inviting for us to associate the 1/N overlap between S 1 and S 0 as originated from the 1/N fraction of τ 1 (ϕ) in the PN 2 subspace, and therefore identify S 0 as the signature set for j = N/2, while S 1 \ S 0 as the signature set for j = N/2 − 1.
We take note that 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 τ 1,i = |v v| (i) N j( =i) |h h| (j) , i = 1, · · · , N . Similarly, S 0 is obtained from considering just τ 0 , 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 S 0 and S 1 \ S 0 do unambiguously distinguish between the projection onto the two j-subspaces respectively, we check explicitly the signatures for all the eigenstates, |j, m z , g(j) , 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 S 0 , while all the kets with j = N/2 − 1 have signatures in S 1 \ S 0 . We then end this section by suggesting the following from our numerical observations: Given N spin-1/2s with one spin entering simultaneously from each of the input ports of an N -Bell multiport, the interference signatures for the j = N/2 and j = N/2 − 1 subspaces are distinguishable.
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 2 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.

Discussion
Ideally, we would like to have a passive implementation of the J 2 measurement for any valid range of ε and N . Unfortunately however, we only manage to do so for N = 2, 3 with Bell multiport, and at the moment be content with approximate realizations in the nearly-pure qubit regime, where projection onto the j = N/2 and N/2 − 1 subspaces is sufficient. To proceed further, instead of Bell multiport with the corresponding discrete Fourier transform then, one may perhaps consider more general complex Hadamard transform [94]. Of course, more generally, measuring J 2 and L ϕ need not be the only scheme that achieves the multi-parameter QFI CRB; there might be some other schemes as well, with different implementation method altogether. We demonstrated one such scheme here and will leave these other possibilities open for now.
Regarding the interference signatures in S 0 and S 1 with N -Bell multiport in the above section, we also observe an interesting feature for the case of bosonic qubits.
With N bosons and N outputs in the Bell multiport, there are in total |S| := 2N − 1 N different signatures possible. Then, we find that apart from some specific N , the total number of signatures in S 1 , |S 1 |, is exactly |S|, which again confirms that there are no room for other j = N/2 values to have distinguishable signature set. Up to N = 12 as we have verified, the few exceptional N values for which |S 1 | = |S| are found to be N = 6, 10, 12, with |S 1 |/|S| = 449 462 , 90358 92378 , 1321717 1352078 respectively. Incidentally, N = 6, 10, 12 are the first three elements in the sequence of numbers which are not integer power of a prime number, and it is a well known open problem in mathematical physics and quantum information whether there exists N + 1 mutually unbiased bases (MUB) for Hilbert space with dimension of such numbers [68,95]. Given that one of the important tools in studying MUB is the complex Hadamard matrices for which discrete Fourier transform is one of such [68,94], it is then perhaps not surprising to have such a correspondence between the non-unity ratio of |S 1 |/|S| and the MUB problem. However, exactly how or why they are related would require a separate analysis, and is beyond the scope of this work.
Finally, as we have put forward in Sec. 3.1, the qubit model we studied in this work can be seen as a particular case of quantum metrology or quantum channel estimation problem, in which we would like to estimate the phase parameter in an unitary encoding channel as well as the dephasing parameter in the dephasing channel, using the Nproduct input state (see Eq. (14)). More generally, one may consider entangled input states (such as one-or two axis-spin squeeze states) and general measurement strategies instead, and obtain better estimation precisions for the two parameters. Then, one may show that, when assuming one parameter is known and we estimate only the second parameter, their optimal attainable MSE in the large N ensemble limit are known to be respectively. Note that we have used the parametrization η = s(ε) = ε here. For their joint estimation, there are no rigorous analytical results proving that the estimation of these two parameters are compatible such that these asymptotic optimal MSE can be attained simultaneously, although there are some numerical measure J 2 and L ϕ with separable states optimal scheme using two-axis spin-squeezed states Figure 6. ξ, the equal-weight sum of joint estimation MSE for the unitary phase parameter ϕ and the dephasing strength parameter ε, normalized by their respective optimal asymptotic precision (attainable with entangled states in general) when assuming another parameter is known and fixed (Eqs. (48,49)), for different number of probes N . For the black solid curve, we consider estimation scheme using separable states with J 2 and L ϕ measurement introduced in Sec. 3, which has a limit of ξ → 1 2 2−ε 2 1−ε 2 as N → ∞; in the plot, we use ε = 0.9. Except in the limit of very strong noise, ε → 1, we see that we cannot attain asymptotic compatible joint estimation for both parameters, i.e., ξ → 1. For comparison, in the red solid curve (which is also the red curve of Fig. 4 in Ref. [28]), the corresponding ξ is plotted for an optimal scheme with entangled input states such as the two-axis spin-squeezed state, which shows the trend that might possibly lead to ξ → 1 as N → ∞. evidence suggested that this might be true [28]. Still, for comparison, let us now consider our joint estimation scheme using product input states with the collective J 2 and L ϕ measurement, in which Eq. (23) gives ∆ 2φ N = 1/(N ε 2 ) and Eq. (21) gives ∆ 2ε N ∼ (1 − ε 2 )/N as N → ∞. This shows that we are able to attain the ultimate estimation precision ∆ 2ε N,met for ε as N → ∞, but unfortunately however there is a factor of 1/(1 − ε 2 ) larger than ∆ 2φ N,met -we therefore cannot attain compatible optimal estimation of these two parameters using our scheme with product input states. For illustration, in Fig. 6 we depict (black solid line) the quantity ξ, which is the equal-weight sum of the MSEs for our joint estimation scheme using product input states with the collective J 2 and L ϕ measurement, Eqs. (21,23), normalized by the respective asymptotic optimal MSE Eqs. (48,49). Asymptotically, we have ξ → 1 2 1 1−ε 2 + 1 = 1 2 2−ε 2 1−ε 2 . For comparison, we also plot the corresponding known (numerical) results of optimal joint estimation for ϕ and ε using the two-axis spinsqueezed states (red solid line, which is also the red curve of Fig. 4 in Ref. [28]), which suggest that ξ → 1 might be attainable. Additional remark: We realize that there is a recent work by J. O. de Almeida et al. [96] on estimation of the separation and relative intensity of two incoherent point sources, which explores similar idea of collective "angular momentum" measurements.

Conclusion
In summary, in this work, we have studied the joint estimation of the length parameter ε and the direction parameter ϕ of the Bloch vector for qubit states, a model which has relevance to various physical problems, e.g., superresolution quantum imaging. We show that these two parameters can be simultaneously estimated with optimal precisions asymptotically, using the collective measurement of the angular-momentum squared operator, followed by a local projective measurement. We also discuss how such measurement scheme can be realized using Bell multiport, either exactly for N = 2, 3, or approximately when the qubits are nearly pure for other N values. Finally, we propose that with a Bell multiport setup, the interference signatures in the j = N/2 and j = N/2 − 1 angular-momentum spaces will be distinct.

Acknowledgments
The author would like to express his gratitude to Konrad Banaszek, Chandan Datta, Marcin Jarzyna, and Jan Ko lodyński for their insightful comments and discussions. The author would also like to extend his appreciation to Konrad Banaszek for his encouragement and support throughout this work. The author also thanks Marcin Jarzyna for providing the data for reference in plotting Fig. 6. This work is part of the project "Quantum Optical Technologies" carried out within the International Research Agendas Programme of the Foundation for Polish Science co-financed by the European Union under the European Regional Development Fund.
Appendix A. Precision of estimation for ε with J 2 In order to help the readers verify the expression in Eq. (21), we list down explicitly here the various constituent terms. With ρ N (ε, ϕ) as in Eq. (13), and the notation A ≡ Tr{ρ N (ε, ϕ)A} for the operator A, we have, by a simple rotation, or equivalently a re-definition of the (x, y, z)-directions, . After traversing the HOM interferometer, these states are respectively transformed as From Eq. (B.2), we thus see that the projection onto the respective j-eigenspace of J 2 can be identified with the unique signatures as stated in Eq. (38). Similar analysis can be done for fermions to obtain Eq. (39).
We thus see that the signatures for j = 3 2 after traversing the Bell tritter are either all qubits exit at the same output ports, or all exit at different output ports.

(C.3)
The m z = 1 2 state, upon introducing u = e 2πi/3 , is also 1 which has the signature of two qubits exit at the same output ports, and one at another. The same conclusion holds for the m z = − 1 2 state. Lastly, for the eigenkets in the second j = 1 2 -subspace, we have j = 1 2 , m z = 1 2 , g(j) = 2 = 1 √ 2 (a † h,1 a † h,2 − a † v,1 a † h,2 )a † h,3 |vac , The m z = 1 2 state is also i which again has the signature of two qubits exiting at the same output ports, and one at another, and the same conclusion holds for the m z = − 1 2 state. Therefore, overall, the projections onto j = 1 2 subspaces all have the same interference signature, which is perfectly distinguishable from the signatures associated with the j = 3 2 subspace as summarized in Eq. (41). Similar analysis can be done for fermions to obtain Eq. (42).